PRO astrom, raref, decref, idgs, ximg, yimg, astr, err, fl, $
	NOQUESTIONS=NoQues, DUMPTYPE=DumpType
;+
; NAME:
;	ASTROM
; PURPOSE:
;	Computes parameters describing a linear fit to a plate    
;	solution using the algorithm supplied by R. S. HARRINGTON, USNO
;	Use the procedure ASTROMIT for interactive use of this program
;
; CALLING SEQUENCE:
;	astrom, raref, decref, idgs, ximg, yimg, astr, err, [ fl ,
;				/NOQUESTIONS, DUMPTYPE= ]
;
; INPUTS:
;	raref - Array of RA'S of astrometric reference stars (DEGREES)
;	decref- Array of DEC'S of astrometric reference stars (DEGREES)
;	idgs  - Array of indices of selected reference stars.  IDGS must
;		contain the same or fewer elements than RAREF and DECREF.
;		If the total # of reference stars is N, then to select
;		all the reference stars, set IDGS = INDGEN(N)
;	ximg  - Array of image source X COORDS.                  
;	yimg  - Array of image source Y COORDS.
;	astr  - A structure, of type ASTROMETRY, containing the initial guess
;		at the plate solution.
;
; OPTIONAL INPUTS:
;	fl    - Camera Focal Length (in pixel units).  If omitted, a focal
;		length of 1.7E5 (UIT camera) will be assumed
; KEYWORDS:
;	NOQUESTIONS - If set, user will not be asked 'Ren/Del/Cont'
;	DUMPTYPE    - If 0, printed information will NOT be dumped.  If 1,
;		printed information will be dumped to ASTROM.DMP in the
;		same format as the screen.  If 2, a higher precision and
;		sexigesimal format will be dumped to file ASTROM.DMP.
;		DumpType=2 is noticeably slower...
; OUTPUTS:
;	astr  - The input ASTROMETRY structure is filled with the just-computed
;		plate solution.
;	err   - (output) error in star location (measured - calc) in arcsec   
;		Same number of elements as IDGS
;
; REVISION HISTORY:
;	Written  by J. K. HILL, STX CORP.                        7/11/86  
;	Converted to IDL, made interactive by R. H. CORNETT, STX 5/29/87
;	03-SEP-90 Fixed Error calculation and made it dump in to file. 
;	   (Deutsch)
;	Added option to delete or rename output file, J. Isensee, 
;	   ST Systems Corp.
;	Added Calling sequence message, N. Collins, STX, 11/20/90
;	20-AUG-91 Fixed Error calculation by removing extra SQRT inserted in
;	   line 161 probably when the ERR variable was changed to ERROUT by 
;	   someone.  Therefore, all error values were wrong, probably since 
;	   11/20/90.  Also changed variable ERROUT back to ERR because only 
;	   SOME of the ERRs had been changed, which gave more errors. 
;	   (E. Deutsch)
;	20-AUG-91 Added /NoQues parameter to avoid Del/Ren Prompt. (E. Deutsch)
;	24-AUG-91 Added DumpType= keyword and corresponding code.  (E. Deutsch)
;	24-AUG-91 Changed Dumpfile Unit=8 to /get_lun              (E. Deutsch)
; 	Fixed the returned ERR array.  The square of the errors was being
;	   returned; a square root is now taken prior to returning.  The 
;	   displayed error analysis is not affected.  M. R. Greason, 
;	   Hughes STX, 10 December 1992.
;	Converted to the new, image-distortion-related, astrometry scheme.
;	Reformatted the procedure and beefed up the internal documentation.
;	   M. R. Greason, Hughes STX, 6 January 1993.
;	Fixed variable-name error.  RSH, HSTX, 25-Mar-93
;	New astrometry software WBL HSTX  Feb 94
;       Changed formatting for dumptype = 2, NRC HSTX Jul 95
;	Changed star index format from I3 to I4.  Modified to access Astro2
;	   distortion correction.  MRG, HSTX, Aug 95
;- 
on_error,2
;
;			Check parameters.
;
npar = n_params()
IF (npar LT 7) THEN begin
 	print, 'Syntax - astrom, raref, decref, idgs, ximg, yimg, astr, err, 
	print,'					[ fl, DUMPTYPE =, NOQUES = ]
	RETURN
ENDIF
;
;				Camera focal length is 2.0D4 for the
;				schwarzschild camera, 1.7D5 for UIT.
;
IF (npar LT 10) THEN fl = 1.7d5
;
NoQues = keyword_set(NoQues)
IF (n_elements(DumpType) EQ 0) THEN DumpType = 1
;
;			Define numerical constants.
;
one = 1.0d0
r90 = !Dpi / 2.d
r270 = 3.d * r90
r360 = !Dpi * 2.d
r2d = 180.0d0/!Dpi
;
;			Other initialization.
;
xycen = astr.crpix - 1.
crvalr = astr.crval/r2d
sdo = sin(crvalr(1))
cdo = cos(crvalr(1))
nmatch = n_elements(idgs)       ; Number of selected reference stars.
;
;			Extract specified target ID's .
;
ra_ref  = double(raref(idgs)) / r2d
dec_ref = double(decref(idgs)) / r2d
x_img   = double(ximg(idgs))
y_img   = double(yimg(idgs))
;
;			Adjust reference star measures & print out all data on
;			star.
;
x = x_img / fl
y = y_img / fl
;
xymean = [total(x),total(y)] / nmatch	; Mean position of reference stars.
;
;			Store coordinate differences from plate center.
;
IF (crvalr(0) GT r270) THEN BEGIN
	select = where (ra_ref LT r90, nselect)
	IF (nselect GT 0) THEN ra_ref(select) = ra_ref(select) + r360
ENDIF ELSE IF (crvalr(0) LT r90) THEN BEGIN
	select = where (ra_ref GT r270, nselect)
        IF (nselect GT 0) THEN ra_ref(select) = ra_ref(select) - r360
ENDIF
;
a = ra_ref - double(crvalr(0))
d = dec_ref
sd = sin(d) & cdx = cos(d)
sa = sin(a) &  ca = cos(a)
denom = (sd * sdo) + (cdx * cdo * ca)
;
;			New A, D are in fact XSI, ETA (radians).
;
a = cdx * sa / denom
d = (sd * cdo - cdx * sdo * ca) / denom  
admean = [total(a),total(d)] / nmatch
;
;			Compute linear plate constants.
;			All measures and standard coordinates are first
;			referred to their respective centroids, eliminating the
;			constant plate constant.  no assumptions are made
;			regarding rigid rotations or perpendicular axes.
;    
a = a - admean(0)
d = d - admean(1)
x = (xymean(0) - x) + a
y = (xymean(1) - y) + d 
;               
sxx    = total( a * a )
syy    = total( d * d )
sxy    = total( a * d )
sxmx   = total( x * a )
sxmy   = total( x * d )
symx   = total( y * a )
symy   = total( y * d )
denom1 = (sxx * syy) - (sxy * sxy)
;
;			Scale constants.
;
ax = (sxmx * syy - sxmy * sxy) / denom1
by = (symy * sxx - symx * sxy) / denom1
;
;			Rotation constants.
;
bx = (sxmy * sxx - sxmx * sxy) / denom1
ay = (symx * syy - symy * sxy) / denom1
denom2 = one - ax - by + (ax * by) - (ay * bx)
;
;			Compute fit sigmas (radians).
;			We're not happy with the values this part gets!
;
IF (nmatch EQ 3) THEN BEGIN
	rmsa = 0.0d0
	rmsd = 0.0d0
ENDIF ELSE BEGIN
	rmsa = total((x - ax * a - bx * d) ^ 2 )
	rmsd = total((y - by * d - ay * a) ^ 2 )
	rmsa = sqrt(rmsa / (nmatch-3.d0))
	rmsd = sqrt(rmsd / (nmatch-3.d0)) 
ENDELSE
;
;			Finish computing the plate solution.
;
;				The CD matrix.
;
cd = replicate(0.d, 2, 2) + [one-by, ay, bx, one-ax] / (denom2 * fl)
astr.cd = r2d * cd
;
;				The CRPIX array.
;
xymean = xymean * fl
v = -admean + cd # xymean
det = cd(0,0) * cd(1,1) - cd(1,0) * cd(0,1)	; Determinant of CD.
astr.crpix(0) = (v(0) * cd(1,1) - v(1) * cd(0,1)) / det + 1.
astr.crpix(1) = (v(1) * cd(0,0) - v(0) * cd(1,0)) / det + 1.
;
;				The CRVAL array.  Use XY2AD to convert
;				the plate center to (RA,DEC).
;
astflg = strpos(astr.ctype(0), 'UIT')		; Identifier flag: >=0 -> UIT
IF (astflg GE 0) THEN uit_xy2ad, xycen(0), xycen(1), astr, ac, dc $
                 ELSE xy2ad, xycen(0), xycen(1), astr, ac, dc
astr.crval = [ac, dc]
astr.crpix = xycen + 1.		; Restore original center position
;
;			Compute and display the errors of the fit.
;
;				Start with the display titles.
;
IF (DumpType GT 0) THEN openw, dmp, 'astrom.dmp', /GET_LUN
IF (DumpType EQ 1) THEN BEGIN
	printf, dmp, '                         STAR POSITIONS and ERRORS'
	printf, dmp, '  I     X      Y       RA (calc) Dec      RA (cat)  Dec     d(RA) d(Dec) d(tot)'
ENDIF
IF (DumpType EQ 2) THEN BEGIN
	printf, dmp, '[',!stime,']                    Star Positions and Errors'
	printf, dmp, 'Idx        X       Y            RA   (calc)   Dec             RA   (cat)    Dec         d(RA)   d(Dec)    d(tot)'
	printf, dmp, '----    ------- -------    --------------------------    --------------------------    ------- -------   --------'
ENDIF
;
print, ' '
print, '                         STAR POSITIONS and ERRORS'
print, '  I     X      Y       RA (calc) Dec      RA (cat)  Dec     d(RA) d(Dec) d(tot)'
;
;				Use XY2AD to compute star positions.
;
IF (astflg GE 0) THEN uit_xy2ad, x_img, y_img, astr, asdeg, dsdeg $
                 ELSE xy2ad, x_img, y_img, astr, asdeg, dsdeg 
as = asdeg / r2d
ds = dsdeg / r2d
;
;				Compute the residuals, in arcseconds.
;
era = (as - ra_ref) * cos(dec_ref) * (r2d * 3600.)
erd = (ds - dec_ref) * (r2d * 3600.)
err = (era ^ 2) + (erd ^ 2)
;
;				Display the info relating to each star.
;
ra_ref = ra_ref*r2d    &   dec_ref = dec_ref*r2d
fmt = '(I4,1X,2F7.1,1X,2F9.4,1X,2F9.4,1X,3F7.2)'
fmth = '(I4,2X,2F8.1,1X,2a30,4x,f7.4,1x,f7.4,3x,f8.4)'
FOR k = 0, nmatch-1 DO BEGIN
	print, FORM=fmt, idgs(k), x_img(k), y_img(k), asdeg(k), dsdeg(k), $
		ra_ref(k), dec_ref(k), era(k), erd(k), sqrt(err(k))
	IF (DumpType EQ 1) THEN printf, dmp, FORM=fmt, idgs(k), x_img(k), $
		y_img(k), asdeg(k), dsdeg(k), ra_ref(k), dec_ref(k), $
		era(k), erd(k), sqrt(err(k))
 	IF (DumpType EQ 2) THEN BEGIN 
                printf, dmp, form=fmth, idgs(k), x_img(k), $
		y_img(k), adstring(asdeg(k), dsdeg(k), 2), $
		adstring(ra_ref(k), dec_ref(k), 2), era(k), erd(k), $
		sqrt(err(k))
        ENDIF
ENDFOR
;
;				Display the summary info.
;
print, 'Using n=' + strn(nmatch) + ' stars:'
print, FORM='(A29,T58,3F7.2)', $
	'RMS positional error(arcsec):', sqrt(total(era ^ 2) / nmatch), $
	sqrt(total(erd ^ 2) / nmatch), sqrt(total(err) / nmatch)
print, FORM='(A11,F6.3,3X,A19,F6.3,3X,A21,F6.3)', 'RMS Error: ', $
	sqrt(total(err)/nmatch), 'RMS Error/sqrt(n): ', $
	sqrt(total(err)/nmatch) / sqrt(nmatch), 'RMS Error/sqrt(n-3): ', $
	sqrt(total(err)/nmatch) / sqrt((nmatch - 3) > 1)
;
IF (DumpType EQ 1) THEN BEGIN
	printf, dmp, FORM='(A29,T58,3F7.2)', $
		'RMS positional error(arcsec):', sqrt(total(era^2)/nmatch), $
		sqrt(total(erd^2)/nmatch), sqrt(total(err)/nmatch)
	printf, dmp, FORM='(A11,F6.3,3X,A19,F6.3,3X,A21,F6.3)','RMS Error: ', $
		sqrt(total(err)/nmatch),'RMS Error/sqrt(n): ', $
		sqrt(total(err)/nmatch)/sqrt(nmatch),'RMS Error/sqrt(n-3): ', $
		sqrt(total(err)/nmatch)/sqrt((nmatch-3)>1)
ENDIF
;
IF (DumpType EQ 2) THEN BEGIN
	printf, dmp, FORM='(A20,T53,A29,T87,f7.4,1x,f7.4,1X,F10.4)', $
		'Using n='+strn(nmatch)+' Stars       ', $
		'RMS positional error(arcsec):',sqrt(total(era^2)/nmatch), $
		sqrt(total(erd^2)/nmatch),sqrt(total(err)/nmatch)
	printf, dmp, FORM='(A11,F8.5,3X,A19,F8.5,3X,A21,F8.5)','RMS Error: ', $
		sqrt(total(err)/nmatch),'RMS Error/sqrt(n): ', $
		sqrt(total(err)/nmatch)/sqrt(nmatch),'RMS Error/sqrt(n-3): ', $
		sqrt(total(err)/nmatch)/sqrt((nmatch-3)>1)
ENDIF
IF (DumpType GT 0) THEN free_lun, dmp	; Close the file.
;
;			The error should be returned in arcseconds, 
;			not arcseconds^2.
;
err = sqrt(err)
;
;			Give the user a chance to manipulate the output
;			summary file.
;
IF ((NoQues EQ 0) AND (DumpType GT 0)) THEN BEGIN
;
;				What does the user want to do with the file?
;
	print, ' '
	print, '____________________________________________________________'
	print, ' '
	print, 'Above results have been printed to file ASTROM.DMP'
	print, ' '
	answ = ' '
	read, 'Type D to DELETE, R to RENAME, Return to KEEP AS IS:  ', answ
	answ = strupcase(strmid(strtrim(answ, 2), 0, 1))
	IF (answ EQ 'D') THEN BEGIN
;
;				Delete it?
;
		IF (!version.os EQ "vms") THEN spawn, 'delete astrom.dmp;*' $
		                          ELSE spawn, 'rm astrom.dmp'
;
	ENDIF ELSE IF (answ EQ 'R') THEN BEGIN
;
;				Rename it?
;
		print, ' '
		newname = ''
		read, 'Enter new file name: ', newname
		print,' '
		IF !version.os EQ "vms" THEN cmnd = 'rename ' $
		                        ELSE cmnd = 'mv '
		spawn, cmnd + 'astrom.dmp ' + newname
;
	ENDIF
ENDIF
;
RETURN
END