;+
; Project     :	SOHO - CDS
;
; Name        :	XCOR_CDS
;
; Purpose     :	Coalign CDS images with generic FITS images
;
; Category    :	Class1, CDS, Analysis, Coordinates
;
; Explanation :	Does a two dimensional cross-correlation between a CDS raster
;		image and an image stored in a FITS file.  From this
;		cross-correlation, the relative pointing is determined.
;
;		The FITS file to correlate the CDS data to must meet the
;		following conditions:
;
;		* It must contain a simple 2D image in the main part of
;		  the FITS file--not in an extension.
;
;		* It must contain the time, either in the single keyword
;		  DATE_OBS in ISO/CCSDS format (e.g. EIT), or in the two
;		  keywords DATE-OBS and TIME-OBS (e.g. BBSO).
;
;		* It must contain the absolute pointing of the image, using the
;		  keywords CRPIX/CRVAL/CDELT.
;
;		* It must encompass the entire field of view of the CDS raster
;		  images, including a margin for pointing error.
;
;		The primary application for this routine is to coalign CDS
;		images with full disk images such as from EIT, MDI, Yohkoh, or
;		ground-based observatories.  Not all such images will have
;		valid CRPIX/CRVAL/CDELT keywords.  One way to add these
;		keywords is to use the routine IMAGE_TOOL.  The procedure would
;		be as follows:
;
;		1.  Copy over the FITS image file to your local directory, e.g.
;
;			cp $SYNOP_DATA/bbso/bbso_halph_fd_19960623_1622.fts .
;
;		2.  Run IMAGE_TOOL.
;		3.  Establish the coordinate system by fitting the limb.
;		4.  Select "Update FITS header" from the "File" pull-down menu.
;
;
; Syntax      :	XCOR_CDS, CDSFILE, IMAGEFILE, WINDOW  [, SUBIMAGE ]
;
; Examples    :	XCOR_CDS, 's3258r00', 'bbso_halph_fd_19960623_1622.fts', 6
;
; Inputs      :	CDSFILE	  = The name of the CDS FITS binary table file.
;			    Alternatively, it can be a quicklook data structure
;			    from READCDSFITS.
;
;		IMAGEFILE = The name of the FITS image file containing the data
;			    to coalign to.
;
;		WINDOW	  = The index of the CDS window to extract.  Valid
;			    values are 0 to N-1.
;
; Opt. Inputs :	None.
;
; Outputs     :	There are no required output parameters.  However graphics are
;		displayed during the processing, and information about the
;		alignment is printed to the screen.
;
; Opt. Outputs:	SUBIMAGE  = If passed, then the part of the FITS image file
;			    corresponding to the CDS raster image is extracted.
;			    The image will have the same dimensions as the CDS
;			    images.
;
; Keywords    :	OUTPUT	  = If passed, then contains the name of a file that
;			    information about the coalignment is appended to.
;			    The following information is written:
;
;				* The OPS_L and _R value from the FITS header.
;				* The X,Y pointing determined from the
;				  cross-correlation.
;				* The average OPS_L and _R values from the FITS
;				  binary table.
;				* The average sun sensor pitch and yaw values
;				  from the FITS binary table.
;
;			    The primary purpose of the OUTPUT keyword is to
;			    assist in the determination of the calibration of
;			    the OPS and sun sensor.
;
;		NOASK	  = If set, then the user is not prompted as to whether
;			    the data should be written to the OUTPUT file.
;
;		NOROTATE  = If set, then don't attempt to take solar rotation
;			    into account.  This feature is useful for aligning
;			    data at the limb.  However, the data should be
;			    taken very close in time if the fit is to be
;			    relevant.
;
;		NOLIMB	  = If set, and NOASK is also set, then images which go
;			    over the limb are not processed.
;
;		NOVIEW	  = If set, then the data windows are not displayed.
;
;		XCOR	  = Returns the calculated cross-correlation function.
;
;		PARAM	  = Returns the parameters of the fit to the
;			    cross-correlation function.
;
;		NCORR	  = Search size to use in calculating the cross
;			    correlation function.  Can be 1 or 2 numbers.
;			    Default is 20.
;
;		SC_ROLL	  = Roll angle during the observation.  This should be
;			    used when aligning data taken when the spacecraft
;			    is at a significant roll.  There are a couple of
;			    factors to keep in mind when using this keyword:
;
;				1.  This should only be used when coaligning
;				    with another SOHO dataset, such as EIT, so
;				    that both datasets are rotated by the same
;				    amount.
;
;				2.  The alignment is done in the rotated
;				    coordinate frame, rather than in the true
;				    frame.  The spacecraft roll angle is only
;				    used for determining the effects of solar
;				    rotation.
;
;		NOROLL	  = Keyword passed to READCDSFITS.  The default value
;			    is 1--i.e. ignore the spacecraft roll angle.  It is
;			    also set to 1 when the SC_ROLL keyword is used.
;
;		PIXEL_RANGE = Set to range of pixels to use within the
;			    subwindow to extract out the desired line.
;
; Calls       :	XCOR_CDS_GAUSS (included)
;
;		DATATYPE, READCDSFITS, FXREAD, UTC2TAI, FXPAR, UTC2STR,
;		STR2UTC, PB0R, ROT_XY, AVERAGE, CDS_CLEAN_IMAGE, PUT, SIGRANGE,
;		CROSS_CORR2, LSTSQR, CONTV, TRIM, ASK, SETFLAG, UNSETFLAG,
;		SETWINDOW, TVDEVICE
;
; Common      :	None.
;
; Restrictions:	None.
;
; Side effects:	Unless the /NOVIEW keyword is set, the two windows 0 and 2 are
;		created.  The TVDEVICE routine is used to channel images to
;		window 2.
;
; Prev. Hist. :	Based on an earlier unreleased version called XCOR_W_EIT.
;
; History     :	Version 1, 12-Dec-1996, William Thompson, GSFC
;		Version 2, 28-Jan-1997, William Thompson, GSFC
;			Corrected bug where contours were not plotted correctly
;			Put up both contours.  Smooth contours.
;			Added keywords XCOR and PARAM
;		Version 3, 04-Mar-1997, William Thompson, GSFC
;			Adjust correlation size to image size
;		Version 4, 05-May-1997, William Thompson, GSFC
;			Add print statement about pixel position.
;		Version 5, 07-May-1997, William Thompson, GSFC
;			Handle data that doesn't completely cover CDS field of
;			view.  Add keyword NCORR, allow different search sizes
;			in X and Y.
;		Version 6, 08-May-1997, William Thompson, GSFC
;			Make double pass through LSTSQR to improve fits.
;		Version 7, 13-Aug-1997, William Thompson, GSFC
;			Call TVDEVICE, /DISABLE before exiting
;		Version 8, 05-Sep-1997, William Thompson, GSFC
;			Allow either CDS filename or QLDS as input.
;		Version 9, 22-Oct-1997, William Thompson, GSFC
;			Fix display bug on wide windows.
;		Version 10, 28-Oct-1997, William Thompson, GSFC
;			Added keyword NOROTATE
;		Version 11, 17-Dec-1997, William Thompson, GSFC
;			Don't call UPD_CDS_POINT, as it is called automatically
;			by READCDSFITS.
;		Version 12, 29-Oct-1998, William Thompson, GSFC
;			Fixed bug where average OPS values weren't properly
;			calculated when some were missing.
;		Version 13, 23-Apr-1999, William Thompson, GSFC
;			Added keywords SC_ROLL and NOROLL.
;		Version 14, 27-Jul-1999, William Thompson, GSFC
;			Added keyword PIXEL_RANGE
;               Version 15, 05-Aug-2003, William Thompson, GSFC
;                       Take into account when SOHO is upside-down.
;
; Contact     :	WTHOMPSON
;-
;
;==============================================================================
;
	FUNCTION XCOR_CDS_GAUSS, POS, PARAM, PERR
;
	X = POS.X - PARAM(2)
	Y = POS.Y - PARAM(3)
	Q = (X/PARAM(4))^2 + (Y/PARAM(5))^2 + X*Y*PARAM(6)
	PERR = 0
	RETURN, PARAM(0) + PARAM(1)*EXP(-Q)
	END
;
;==============================================================================
;
	PRO XCOR_CDS, CDSFILE, IMAGEFILE, WINDOW, SUBIMAGE, OUTPUT=OUTFILE, $
		NOASK=NOASK, NOLIMB=NOLIMB, NOVIEW=NOVIEW, XCOR=C,	$
		PARAM=PARAM, NCORR=NCORR, NOROTATE=NOROTATE, NOROLL=K_NOROLL, $
		SC_ROLL=SC_ROLL, PIXEL_RANGE=IPIX
;
	ON_ERROR, 2
;
;  Check the input parameters.
;
	IF N_PARAMS() LT 3 THEN MESSAGE,	$
		'Syntax:  XCOR_CDS, CDSFILE, IMAGEFILE, WINDOW  [, SUBIMAGE ]'
;
	INPUT_TYPE = DATATYPE(CDSFILE,1)
	IF (N_ELEMENTS(CDSFILE) NE 1) OR ((INPUT_TYPE NE 'String') AND	$
		(INPUT_TYPE NE 'Structure')) THEN MESSAGE,		$
		'CDSFILE must be a scalar string or structure.'
	IF (N_ELEMENTS(IMAGEFILE) NE 1) OR (DATATYPE(IMAGEFILE,1) NE 'String')$
		THEN MESSAGE, 'IMAGEFILE must be a scalar string.'
;
	IF N_ELEMENTS(OUTFILE) NE 0 THEN BEGIN
	    IF (N_ELEMENTS(OUTFILE) NE 1) OR (DATATYPE(OUTFILE,1) NE $
		'String') THEN MESSAGE, 'OUTFILE must be a scalar string.'
	ENDIF
;
;  Check the NCORR parameter, if passed.
;
	CASE N_ELEMENTS(NCORR) OF
	    1: M_CROSS = NCORR
	    2: M_CROSS = NCORR
	    ELSE: M_CROSS = 20
        ENDCASE
;
;  Determine the default value for NOROLL.  If SC_ROLL is passed, or the NOROLL
;  keyword was not passed in, then set NOROLL=1.  Otherwise, if the user didn't
;  pass in a value, then leave it undefined, and let READCDSFITS determine the
;  default value for roll.
;
	IF N_ELEMENTS(SC_ROLL) EQ 1 THEN BEGIN
	    NOROLL = 1
	END ELSE IF N_ELEMENTS(K_NOROLL) EQ 1 THEN BEGIN
	    NOROLL = K_NOROLL
	END ELSE NOROLL = 1
;
;  Create the windows used for this program.  Use window 2 for the images, and
;  window 0 for the plots.
;
	IF NOT KEYWORD_SET(NOVIEW) THEN BEGIN
	    WINDOW, 0
	    WINDOW, 2
	    SETWINDOW, 0
	    TVDEVICE, 2
	ENDIF
;
;  Set the missing pixel flag value for the image processing and display
;  routines.  Keep the old setting, so that one can go back.
;
	OLDFLAG = !IMAGE.MISSING
	SETFLAG, MISSING=-1
;
;  Read in the data files.
;
	IF INPUT_TYPE EQ 'Structure' THEN A = CDSFILE ELSE BEGIN
	    A = READCDSFITS(CDSFILE,PRESELECT=WINDOW,NOROLL=NOROLL)
	    IF DATATYPE(A,1) NE 'Structure' THEN MESSAGE,	$
		    'Unable to read file ' + CDSFILE
	ENDELSE
;
	IF NOT FILE_EXIST(IMAGEFILE) THEN MESSAGE,	$
		'File ' + IMAGEFILE + ' does not exist'
	MESSAGE = ''
	FXREAD,IMAGEFILE,TEMP,HEADER,ERRMSG=MESSAGE
	IF MESSAGE NE '' THEN MESSAGE, MESSAGE
;
;  Store the image in a bigger array with zeroes around the edge.
;
	SZB = SIZE(TEMP)
	IF SZB(0) NE 2 THEN MESSAGE, 'Reference image must be two dimensional'
	B = FLTARR(SZB(1)+2, SZB(2)+2)
	B(1,1) = TEMP
;
;  Get the times of each observation.  Use the middle of each observation.
;
	CDS_TIME = (UTC2TAI(A.HEADER.DATE_OBS) + UTC2TAI(A.HEADER.DATE_END))/2.
;
;  The image time could be stored in the CCSDS-format keyword DATE_OBS (e.g
;  EIT), or it could be in the standard FITS keyword DATE-OBS, and in the
;  not-so-standard keyword TIME-OBS (e.g. BBSO).
;
	IMAGE_TIME = FXPAR(HEADER,'DATE_OBS')
	IF DATATYPE(IMAGE_TIME,1) NE 'String' THEN BEGIN
	    IMAGE_DATE = FXPAR(HEADER,'DATE-OBS')
	    IMAGE_TIME = FXPAR(HEADER,'TIME-OBS')
	    IF (DATATYPE(IMAGE_DATE,1) NE 'String') OR	$
		    (DATATYPE(IMAGE_TIME,1) NE 'String') THEN MESSAGE,	$
		    'FITS file ' + IMAGEFILE +	$
		    ' does not contain a proper date/time'
;
;  Sometimes the dates and times have embedded blanks, which must be removed.
;
	    IMAGE_DATE = STRCOMPRESS(IMAGE_DATE, /REMOVE_ALL)
	    IMAGE_TIME = STRCOMPRESS(IMAGE_TIME, /REMOVE_ALL)
	    IMAGE_DATE = UTC2STR(STR2UTC(IMAGE_DATE,/DMY), /DATE_ONLY)
	    IMAGE_TIME = IMAGE_DATE + ' ' + IMAGE_TIME
	ENDIF
	IMAGE_TIME = UTC2TAI(IMAGE_TIME) + FXPAR(HEADER,'EXPTIME') / 2.
;
;  Extract the XCEN and YCEN values from the CDS header, and rotate to the time
;  of the other observation.
;
	ERROR = ''
	R = SQRT(A.HEADER.XCEN^2 + A.HEADER.YCEN^2)
	R0 = PB0R(CDS_TIME,/SOHO,/ARCSEC)
	XCEN = A.HEADER.XCEN
	YCEN = A.HEADER.YCEN
	IF (R LT R0(2)) AND NOT KEYWORD_SET(NOROTATE) THEN BEGIN
	    IF N_ELEMENTS(SC_ROLL) EQ 1 THEN BEGIN
		XX = XCEN
		YY = YCEN
		COS_R = COS(SC_ROLL/!RADEG)
		SIN_R = SIN(SC_ROLL/!RADEG)
		XCEN = XX * COS_R - YY * SIN_R
		YCEN = YY * COS_R + XX * SIN_R
	    ENDIF
	    RESULT = ROT_XY(XCEN, YCEN, TSTART=CDS_TIME, TEND=IMAGE_TIME, $
		    ERROR=ERROR, /SOHO)
	    IF ERROR EQ '' THEN BEGIN
		IF N_ELEMENTS(SC_ROLL) EQ 1 THEN BEGIN
		    XCEN = RESULT(0) * COS_R + RESULT(1) * SIN_R
		    YCEN = RESULT(1) * COS_R - RESULT(0) * SIN_R
	        END ELSE BEGIN
		    XCEN = RESULT(0)
		    YCEN = RESULT(1)
		ENDELSE
	    END ELSE BEGIN
		XCEN = A.HEADER.XCEN
		YCEN = A.HEADER.YCEN
	    ENDELSE
	ENDIF
;
;  If the NOASK and NOLIMB keywords were both set, and the observation goes
;  over the limb, then don't process the data.
;
	IF KEYWORD_SET(NOASK) AND KEYWORD_SET(NOLIMB) AND	$
		(SQRT((ABS(XCEN)+A.HEADER.IXWIDTH/2.)^2 +	$
		(ABS(YCEN)+A.HEADER.IYWIDTH/2.)^2) GT R0(2)) THEN GOTO, FINISH
;
;  Determine the coordinate system for the IMAGE data.
;
	XPIX = FXPAR(HEADER,'CRPIX1')  &  YPIX = FXPAR(HEADER,'CRPIX2')
	XVAL = FXPAR(HEADER,'CRVAL1')  &  YVAL = FXPAR(HEADER,'CRVAL2')
	RSUN0 = PB0R(IMAGE_TIME,/SOHO,/ARCSEC)
	SOLAR_R = FXPAR(HEADER,'SOLAR_R')
	IF SOLAR_R GT 0 THEN BEGIN
		XDEL = RSUN0(2) / SOLAR_R
		YDEL = XDEL
	END ELSE BEGIN
		XDEL = FXPAR(HEADER,'CDELT1')
		YDEL = FXPAR(HEADER,'CDELT2')
	ENDELSE
;
;  If the images are upside-down, then change the sign of XDEL,YDEL.
;
        IF GET_SOHO_ROLL(IMAGE_TIME) EQ 180 THEN BEGIN
            XDEL = -XDEL
            YDEL = -YDEL
        ENDIF
;
;  Extract the window data, and average over wavelength.  Remove cosmic rays.
;
	AA = GT_WINDATA(A,WINDOW)
	IF N_ELEMENTS(IPIX) EQ 2 THEN AA = AA(IPIX(0):IPIX(1),*,*)
	AA = AVERAGE(AA,1,MISSING=-1)
	CDS_CLEAN_IMAGE, AA, /NOMISSING
	SZ = SIZE(AA)
;
;  From the CDS header, determine the relative scales.
;
	XSIZE = A.HEADER.IXWIDTH / SZ(1)
	YSIZE = A.HEADER.IYWIDTH / SZ(2)
	XSCALE = XSIZE / ABS(XDEL)
	YSCALE = YSIZE / ABS(YDEL)
;
;  If feature tracking is off, then adjust the scales accordingly.
;
	IF (A.HEADER.TRACKING EQ 0) AND NOT KEYWORD_SET(NOROTATE) THEN BEGIN
	    CDS_END = UTC2TAI(A.HEADER.DATE_END)
	    RESULT = ROT_XY(A.HEADER.XCEN, A.HEADER.YCEN, TSTART=CDS_TIME, $
		    TEND=CDS_END, ERROR=ERROR, /SOHO)
	    XCOR = RESULT(0)-A.HEADER.XCEN
	    XSIZE = (A.HEADER.IXWIDTH + 2*XCOR) / SZ(1)
	    XSCALE = XSIZE / ABS(XDEL)
	ENDIF
;
;  Calculate the pixel numbers containing the CDS data range, and extract the
;  subimage.  Add 1 to the pixel numbers to accommodate the fact that the image
;  was increased in size.
;
	IC = XPIX + (XCEN-XVAL)/XDEL
	JC = YPIX + (YCEN-YVAL)/YDEL
	I0 = IC - XSCALE * (SZ(1)-1) / 2.
	J0 = JC - YSCALE * (SZ(2)-1) / 2.
	BB = POLY_2D(B,[I0+1,0,XSCALE,0],[J0+1,YSCALE,0,0],0,SZ(1),SZ(2))
;
;  Use the same missing pixel flag value.
;
	W = WHERE(BB EQ 0, COUNT)
	IF COUNT GT 0 THEN BB(W) = -1
;
;  Display the images
;
	IF NOT KEYWORD_SET(NOVIEW) THEN BEGIN
	    PUT, SIGRANGE(AA,FRAC=0.99), 1, 2, SCALE=[XSIZE,YSIZE], /ADJUST, $
		    GEOMETRY=GEOM
	    PUT, SIGRANGE(BB,FRAC=0.99), 2, 2, SCALE=[XSIZE,YSIZE], /ADJUST, $
		    /DATA
	ENDIF
;
;  Calculate the cross correlation function, and find the maximum.
;
	M_CROSS = M_CROSS < [(SZ(1)-2), (SZ(2)-2)]
	C = CROSS_CORR2(AA,BB,M_CROSS,XC,YC,MISSING=-1)
	N_CROSS_X = N_ELEMENTS(XC)
	N_CROSS_Y = N_ELEMENTS(YC)
	W = WHERE(C EQ MAX(C))
	XMAX = XC(W(0) MOD N_CROSS_X)
	YMAX = YC(W(0) / N_CROSS_X)
	IF NOT KEYWORD_SET(NOVIEW) THEN		$
		PLOT_IMAGE, C, ORIGIN=[XC(0),YC(0)], SCALE=1
;
;  Fit a bilateral gaussian to the cross correlation function and determine the
;  best fit.
;
	POS = REPLICATE({X: 0.0, Y: 0.0}, N_CROSS_X, N_CROSS_Y)
	POS.X = XC # REPLICATE(1., N_CROSS_Y)
	POS.Y = REPLICATE(1., N_CROSS_X) # YC
	PARAM0 = [1.D-2, MAX(C), XMAX+0.01, YMAX+0.01, 6.D0, 6.D0, 1.D-4]
	PARAM = PARAM0
;
;  Make two iterations.  The first time, determine the weights from the
;  correlation values.  The second time, use the previous fit to modify the
;  weights.
;
	WGT = C(*)
	FOR ITER = 0,1 DO BEGIN
	    LSTSQR, POS(*), C(*), 'XCOR_CDS_GAUSS', PARAM, PERR,	$
		    MAX_ITER=100, WEIGHT=WGT^4, N_ITER=N_ITER, CHISQR=CHISQR
	    IF ITER EQ 0 THEN BEGIN
		IF (ABS(PARAM(2)-XMAX) GT 3) OR (ABS(PARAM(3)-YMAX) GT 3) $
			THEN PARAM = PARAM0
		WGT = WGT * XCOR_CDS_GAUSS(POS(*), PARAM)
	    ENDIF
	ENDFOR
;
;  If the parameter errors were not calculated, then set them all to 999.
;
	IF N_ELEMENTS(PERR) EQ 0 THEN BEGIN
		N_ITER = 100
		PERR = REPLICATE(999,N_ELEMENTS(PARAM))
	ENDIF
;
;  Correct the parameter errors for the CHISQR value.
;
	PERR = PERR * SQRT(CHISQR)
;
;  Display the fit.
;
	IF NOT KEYWORD_SET(NOVIEW) THEN BEGIN
		OPLOT,PARAM(2:2),PARAM(3:3),PSYM=1,SYMSIZE=3,COLOR=0
		SCALE=[XSIZE,YSIZE]
		CONTV, SIGRANGE(SMOOTH(AA,3)), ORIGIN=-PARAM(2:3)*SCALE, $
			SCALE=SCALE
		SETIMAGE, GEOM(0), GEOM(1), GEOM(2), GEOM(3)
		CONTV, SIGRANGE(SMOOTH(BB,3)), ORIGIN=PARAM(2:3)*SCALE, $
			SCALE=SCALE
		EMPTY
	ENDIF
;
;  Calculate the pointing inferred from the other image.
;
	XCEN = (IC - XPIX)*XDEL + XVAL - XSIZE*PARAM(2)
	YCEN = (JC - YPIX)*YDEL + YVAL - YSIZE*PARAM(3)
;
;  Rotate the pointing to the correct time.
;
	XROT = XCEN
	YROT = YCEN
	R = PB0R(CDS_TIME,/SOHO,/ARCSEC)
	R = R(2)
	IF (R^2 GT (XCEN^2 + YCEN^2)) AND NOT KEYWORD_SET(NOROTATE) THEN BEGIN
	    XX = XCEN
	    YY = YCEN
	    IF N_ELEMENTS(SC_ROLL) EQ 1 THEN BEGIN
		XX = XCEN * COS_R - YCEN * SIN_R
		YY = YCEN * COS_R + XCEN * SIN_R
	    ENDIF
	    ERROR = ''
	    RESULT = ROT_XY(XX, YY, TSTART=IMAGE_TIME, TEND=CDS_TIME, $
		    ERROR=ERROR, /SOHO)
	    IF ERROR EQ '' THEN BEGIN
		IF N_ELEMENTS(SC_ROLL) EQ 1 THEN BEGIN
		    XROT = RESULT(0) * COS_R + RESULT(1) * SIN_R
		    YROT = RESULT(1) * COS_R - RESULT(0) * SIN_R
	        END ELSE BEGIN
		    XROT = RESULT(0)
		    YROT = RESULT(1)
		ENDELSE
	    ENDIF
	ENDIF
;
;  Print out the information about the result.
;
	PRINT, ''
	PRINT, 'OPS_L/R = ', A.HEADER.OPS_L, A.HEADER.OPS_R
	PRINT, 'X/Y in FITS header = ', A.HEADER.XCEN, A.HEADER.YCEN
	PRINT, 'X Calculated = ', XROT, ' +/- ', PERR(2)*XSIZE
	PRINT, 'Y Calculated = ', YROT, ' +/- ', PERR(3)*YSIZE
	PRINT, 'Cross correlation coefficient = ' +	$
		TRIM(INTERPOLATE(C,PARAM(2)+M_CROSS(0),PARAM(3)+M_CROSS(1)))
;
;  Recalculate the pixel numbers containing the CDS data range,
;
	IC = XPIX + (XCEN-XVAL)/XDEL
	JC = YPIX + (YCEN-YVAL)/YDEL
	I0 = IC - XSCALE * (SZ(1)-1) / 2.
	J0 = JC - YSCALE * (SZ(2)-1) / 2.
	PRINT, 'CDS image starts at pixel (' + TRIM(I0) + ',' + TRIM(J0) + ')'
	PRINT, ''
;
;  If requested, extract again the subimage.
;
	IF N_PARAMS() EQ 4 THEN SUBIMAGE =	$
		POLY_2D(B,[I0+1,0,XSCALE,0],[J0+1,YSCALE,0,0],0,SZ(1),SZ(2))
;
;  If the output filename was passed, then ask if the user wishes to write this
;  entry to the output file.
;
	IF N_ELEMENTS(OUTFILE) EQ 1 THEN BEGIN
	    IF KEYWORD_SET(NOASK) THEN BEGIN
		IF N_ITER EQ 100 THEN ANSWER = 'N' ELSE ANSWER = 'Y'
	    END ELSE BEGIN
	        ASK, 'Do you want to write this to the output file? ', ANSWER
	    ENDELSE
	    IF ANSWER EQ 'Y' THEN BEGIN
;
;  Get the OPS and sun sensor values.
;
		W = WHERE(A.PITCHDATA NE A.PITCHDESC.MISSING)
		PITCH = AVERAGE(SENSOR_CALIB(A.PITCHDATA(W),	$
			A.PTCHBITSDATA(W), 0))
		W = WHERE(A.YAWDATA NE A.YAWDESC.MISSING)
		YAW   = AVERAGE(SENSOR_CALIB(A.YAWDATA(W),	$
			A.YAW_BITSDATA(W), 1))
		W = WHERE(A.OPS_LDATA NE A.OPS_LDESC.MISSING)
		OPS_L = AVERAGE(A.OPS_LDATA(W))
		W = WHERE(A.OPS_RDATA NE A.OPS_RDESC.MISSING)
		OPS_R = AVERAGE(A.OPS_RDATA(W))
;
	        OPENU, UNIT, OUTFILE, /APPEND, /GET_LUN
	        PRINTF, UNIT, A.HEADER.OPS_L, A.HEADER.OPS_R, XROT, YROT, $
		        OPS_L, OPS_R, PITCH, YAW
	        FREE_LUN, UNIT
	    ENDIF
	ENDIF
;
;  Reinstate the old missing pixel flag value, disable the TV device, and
;  return.
;
FINISH:
	IF OLDFLAG.SET THEN SETFLAG,MISSING=OLDFLAG.VALUE ELSE	$
		UNSETFLAG, /MISSING
	TVDEVICE, /DISABLE
	RETURN
	END