# QREGISTER: Quick register of images (replaces individual macros)
#  Written by A. C. Phillips at Lick Observatory, Jun 94

procedure qregister (input, output, reference, coordin, inversion, dispim)

file	input {prompt = "input image"}
file	output {prompt = "output image"}
file	reference {prompt = "reference image"}
file	coordin {prompt = "coord.list in reference stars"}
bool	inversion {prompt = "Is input inverted wrt. reference?"}
bool	dispim {prompt = "Display image?"}

begin
{

# define local variables
int	nx, ny, nxr, nyr
real	xstar1, ystar1, xstar2, ystar2
file	inp, ref, out
file	qtmp0, qtmp1, qtmp2, qtmp3

# Construct names
inp = input
ref = reference
out = output

if (access (inp//".imh"))
	inp = input//".imh"
else if (access (inp//".fit"))	
	inp = input//".fit"
else if (access (inp//".hhh"))	
	inp = input//".hhh"
if (access (ref//".imh"))
	ref = reference//".imh"
else if (access (ref//".fit"))	
	ref = reference//".fit"
else if (access (ref//".hhh"))	
	ref = reference//".hhh"

if (! access (ref)) {
	beep
        print ""
        print ("Error: "//ref//" does not exist ! -- exiting")
        bye
}
if (! access (inp)) {
	beep
        print ""
        print ("Error: "//inp//" does not exist ! -- exiting")
        bye
}
if (access (out) || access (out//".imh") || access (out//".hhh")) {
	beep
        print ""
        print ("Error: "//out//" already exists ! -- exiting")
        bye
}

qtmp0 = mktemp ("tmpq")
qtmp1 = mktemp ("tmpq")
qtmp2 = mktemp ("tmpq")
qtmp3 = mktemp ("tmpq")

# Get the dimnesions of the images:
hselect (inp, "naxis[12]", yes) | scan (nx, ny)
hselect (ref, "naxis[12]", yes) | scan (nxr, nyr)
#
# Get the initial centerings for REFERENCE image:
#
print ("... working on coordinate lists\n")
fitpsf (ref, 7., func="elgauss", coord=coordin, inter=no, veri=no, out=qtmp1)
txdump (qtmp1, "xcen,ycen", yes, >qtmp2)
delete (qtmp1, veri-, def+)
fitpsf (ref, 7., func="elgauss", coord=qtmp2, inter=no, veri=no, out=qtmp1)
delete (qtmp2, veri-, def+)
print ("\n\tREFER (check that all stars appear real):")
print ("ID   XCEN   YCEN   XSIG   YSIG   AMPL")
txdump (qtmp1, "lid,xcen,ycen,xsigma,ysigma,amplitude", yes, >"STDOUT")
txdump (qtmp1, "xcen,ycen", yes, >qtmp3)
delete (qtmp1, veri-, def+)
#
# Now get input image coords:
#
if (dispim)
	display (inp, zs+)

print ("Mark 1st reference star:\n")
print (imcur) | scan (xstar1, ystar1)
print ("Mark 2st reference star:\n")
print (imcur) | scan (xstar2, ystar2)
xytran (qtmp3, qtmp2, xstar1, ystar1, xstar2, ystar2, inversion=inversion, ck+, nxt=nx, nyt=ny, tol=2., common_out=qtmp0)
delete (qtmp3, veri-, def+)
#
# Get coordinates for FINAL REFER and INPUT lists:
#
fitpsf (ref, 7., func="elgauss", coord=qtmp0, inter=no, veri=no, out=qtmp1)
delete (qtmp0, veri-, def+)
print ("\nID  XCEN   YCEN   XSIG   YSIG   AMPL", >qtmp3)
print ("\tREFER("//ref//") PSF Params:", >>qtmp3)
txdump (qtmp1, "lid,xcen,ycen,xsigma,ysigma,amplitude", yes,>>qtmp3)
txdump (qtmp1, "xcen,ycen", yes, >qtmp0)
delete (qtmp1, veri-, def+)
fitpsf (inp, 7., func="elgauss", coord=qtmp2, inter=no, veri=no, out=qtmp1)
delete (qtmp2, veri-, def+)
txdump (qtmp1, "xcen,ycen", yes, >qtmp2)
delete (qtmp1, veri-, def+)
fitpsf (inp, 7., func="elgauss", coord=qtmp2, inter=no, veri=no, out=qtmp1)
print ("\n\tINPUT (check that all stars appear real):")
print ("ID  XCEN   YCEN   XSIG   YSIG   AMPL")
txdump (qtmp1, "lid,xcen,ycen,xsigma,ysigma,amplitude", yes, > "STDOUT")
txdump (qtmp1, "xcen,ycen", yes, >>qtmp0)
delete (qtmp1, veri-, def+)
delete (qtmp2, veri-, def+)
table (qtmp0, nc=2, >qtmp2)
geoscal (qtmp2, nx=nxr, ny=nyr, >qtmp1)
print ("Writing transformed image...\n")
geotran (inp, out, qtmp1, qtmp2, interp="spline3")
delete (qtmp0, veri-, def+)
delete (qtmp1, veri-, def+)

fitpsf (out, 7., func="elgauss", coord=qtmp2, inter=no, veri=no, out=qtmp1)
print ("\tOUPUT("//out//") PSF Params:", >>qtmp3)
txdump (qtmp1, "lid,xcen,ycen,xsigma,ysigma,amplitude", yes,>>qtmp3)
delete (qtmp1, veri-, def+)
delete (qtmp2, veri-, def+)

type (qtmp3)
delete (qtmp3, veri-, def+)

}
end
