function sxt_composite,index,data,comp_index,dc_data=dc_data,sx=sx,     $
	sd=sd,satpix=satpix,plot=plot,                                  $
	smooth_factor=smooth_factor,sfd=sfd,		        	$
	dc_interpolate=dc_interpolate,qtest=qtest,register=register,    $
	despike=despike
;+
; NAME:
;   sxt_composite
; PURPOSE:
;   Prepare a composite SXT image from 2 or 3 Half or Quarter Res images.
; CALLING SEQUENCE EXAMPLES:
;   cimg = sxt_composite(index,data,sx=[16,17],comp_index)
;   cimg = sxt_composite(index,data,dc_data=dc_data,sx=[1,2,3],sd=[0,2,1])
;   cimg = sxt_composite(index,data,satpix=satpix,despike=7)
;   cimg = sxt_composite(index,data,satpix=satpix,/sfd,/dc_interp,/reg)
; INPUTS:
;   index = SXT index structure
;   data  = 3-d data cube (at least 2 images)
;
;  If data is byte and no dc_data then		leak_sub
;  If data is byte and dc_data    then		decompress and data-dc_data
;  If data is not byte type       then		assume already decompressed, 
;						-- and background subtracted
; OUTPUT:
;  Function result is a 2-d long array (linear scaling) unless data is float
;  type, then output is a 2-d float array.
; OPTIONAL OUTPUT:
;    comp_index = index structure for composite image, made from index
;	of longest exposure with history records appended.
; OPTIONAL INPUT KEYWORDS:
;   dc_data = Dark frame data cube, use these rather than leak_sub.
;   sx	   = Indices of the data array. If not supplied, compute values
;	     that result in ==> sx = [short, medium, long exposure]
;   sd     = Indices of the background array.  If not supplied, 
;	     then must supply   sd = [short, medium, long exposure]
;   satpix = Saturated pixel array (required if data is not byte type).
;	     Must be 3-d and correspond to data(*,*,*).
;   plot   = If set, plot the results as we go.  The results are displayed
;	     on a logarithmic scale.
;   smooth_factor = Boxcar smooth the bleed pixels.  If not defined, the
;	     default is to smooth over 5 pixels.
;   sfd    = /sfd will cause the image to come back scaled as an sfd image
;	     (per second/HR pixel, and scaled by sfd_comp).
;   dc_interpolate = If set = 0 leak_sub will not use DC interpolation.
;		[Default is to use DC interpolation.]
;   register = If set, do a whole pixel registration of short exposures to 
;	     the long exposure
;   despike = Calls de_spiker, despike=10 is a reasonable level in most cases.
; RESTRICTIONS:
;   Must be called separately for each composite output image.
; MODIFICATION HISTORY:
;    8-Nov-91, L. W. Acton (Wrote composite)
;   19-May-92, LWA  (Several more mods in betwee)
;   13-may-93, J. R. Lemen, Re-wrote.  Calls sxt_satpix
;    1-Dec-93, MDM - Corrected error with the use of "s2" (when it is -1)
;		- Modification on selection of replacement pixels.
;		  It used to check (a) that it was high signal (ge 90)
;				   (b) that it was not saturated
;		  Added check so that it replaced pixels in the longer 
;		  exposure only when the longer exposure was saturated
;    15-Dec-93, SLF - Fixed typo (stapix->satpix)
;    22-feb-94, JRL - Corrected the algorithm; 
;		    - added smooth_factor, sfd, dc_interpolate, register keywords
;    17-Apr-94, LWA - Fixed /sfd scaling option (1000.*  was needed)
;		    - Changed /plot to plot raw input images.
;		    - Added RESTRICTIONS statement to header.
;		    - Added despike=despike keyword.
;    23-jun-94, JRL - Output is floating if input is a fltarr;
;		    - Call sfd_comp with index=index to normalize properly
;    08-May-95, LWA - Fixed to return negative values.
;    02-Nov-96, LWA - Added output of index with history record.
;		    - Made dc_data a keyword intput.
;		    - Made dc_interpolate the default is dc_data not input.
;    04-Nov-96, LWA - Removed forgotten diagnostic.
;    10-Jan-97, LWA - Corrected calling sequence examples.
;-

progverno=3.00*1000
nparams = n_params()
if nparams lt 2 then begin
   doc_library,'sxt_composite'
   return,'**  Must have at least 2 input parameters'
endif

;-------------------------
;  Check input parameters
;-------------------------

sz0 = size(data)			; Data must be 3-d
if sz0(0) ne 3 then begin
   print,'***  Error in sxt_composite:  data must be 3-d'
   help,data
   return,-1
endif else nimg = sz0(3)

if n_elements(index) lt 2 then begin	; Must be at least 2 images
   print,'***  Error in sxt_composite:  index must have minimum length of 2'
   help,index
   return,-1
endif else nimg = nimg < n_elements(index)

if n_elements(sx) eq 0 then sx = indgen(nimg < 3)
nimg = n_elements(sx)
sx = sx(sort(gt_expdur(index(sx))))	; short, medium, long

; Make sure all the images have the same resolution:
if gt_res(index(sx(0))) ne total(gt_res(index(sx)))/ nimg then begin
    message,'The images must all have the same resolution'
    print,'Resolutions present = ',gt_res(index(sx),/str)
    return,-1
endif

if keyword_set(dc_data) then begin		; Assume we have dc images
   if n_elements(sd) eq 0 then begin
	sd = indgen(n_elements(sx))
	print,'*** Warning:  I will supply dark frame indices: sd = ',sd
   endif   
endif

; ----
; Are data byte type?
; ----

d_type = sz0(sz0(0)+1)
if d_type eq 1 then begin	; Yes - must decompress, background subtract
if NOT keyword_set(dc_interpolate) then dc_interpolate=1         ;LWA 2/11/96
;  if nparams eq 3 then begin	; Background was supplied
  if keyword_set(dc_data) then begin   ; Background was supplied  LWA 2/11/96
; Removed >0 in following 2 statements.   LWA 5/8/95.
      img = (sxt_decomp(data(*,*,sx)) - sxt_decomp(dc_data(*,*,sd),/silent))
  endif else img = leak_sub(index(sx),data(*,*,sx),dc_interpolate=dc_interpolate)
  t0 = gt_expdur(index(sx))
endif else t0 = gt_expdur(index(sx))	; t0 = [short, med, long]

; ----
;  Set up the plotting window and display the input images.
; ----

xsize = n_elements(data(*,0,0))
ysize = n_elements(data(0,*,0))
if keyword_set(plot) then begin
  set_plot,'x'			; Make sure
  wind_x = !d.x_size
  wind_y = !d.y_size
  if (wind_x lt 2*xsize) or (wind_y lt 2*ysize) then $
	wdef,lun,2*xsize,2*ysize,/free,/uright  else $
  	wshow & erase
  if (d_type eq 1) then tvscl,data(*,*,sx(0)),0,ysize $
        else tvscl,alog(data(*,*,sx(0))>1),0,ysize
  if (d_type eq 1) then tvscl,data(*,*,sx(1)),xsize,ysize $
        else tvscl,alog(data(*,*,sx(1))>1),xsize,ysize
  if (nimg eq 3) then begin
    if (d_type eq 1) then tvscl,data(*,*,sx(2)),xsize,0 $
          else tvscl,alog(data(*,*,sx(2))>1),xsize,0
  endif
endif

; -----
; Set up the saturated pixel array
; -----
if n_elements(satpix) eq 0 then begin
   satpix = sxt_satpix(index,data,sx) 
   s_sat  = indgen(n_elements(sx))
endif else s_sat = sx

; ----
; Work out the shifts for the short and medium exposures
; ----

if keyword_set(register) then begin
  cen = sxt_cen(index(sx))
  dx   = (cen(0,0:nimg-2) - cen(0,nimg-1)) / (2^gt_res(index(sx(0))))	; They are the same res
  dy   = (cen(1,0:nimg-2) - cen(1,nimg-1)) / (2^gt_res(index(sx(0))))	; They are the same res
  dx   = fix(dx + 0.5*sgn(dx))	; Number of binned pixels to shift the short/med image by
  dy   = fix(dy + 0.5*sgn(dy))	; positive dx ==> long exposure (nimg-2) is to the right
				; ==> when extracting from short and med exposures, will subtract
				;     from the subscripts
  if abs(dx(0))+abs(dy(0)) gt 0 then 							$
	message,string('Shifting Short  exposure by ',dx(0),',',dy(0),' pixels',	$
	format='(2(a,i4),a)'),/info
  if nimg eq 3 then 									$
     if abs(dx(1))+abs(dy(1)) gt 0 then							$
	message,string('Shifting Medium exposure by ',dx(1),',',dy(1),' pixels',	$
	format='(2(a,i4),a)'),/info
endif

; ---------------------------------------------------------------------------
; Composite Algorithm
; ---------------------------------------------------------------------------
; Algorithm for 3 image case:
;
;1.  Replace the saturated pixels in the medium exposure with the
;    corresponding pixels in the short exposure.
;2.  Replace saturated pixels in long exposure with corresponding pixels
;    from the combined medium/short exposure == cimg (combined image).
;3.  Make a temporary copy of the combined image and smooth with a
;    boxcar width of 5  == timg (temporary image).
;4.  cimg(s2) = timg(s2) where s2 = saturated pixels in the long exposure.
;5.  Let cimg(m1) = short_exposure(m1) where m1 is given by:
;    m1 = where((short_exposure ge 90) and (short_expousre is not saturated) and
;	         (long_exposure is saturated))
;
; Algorithm for 2 image case:
;   Same as above except Step 1 is omitted.
; ---------------------------------------------------------------------------

; Set up the output array (cimg2) and working array (cimg1):
typx = data_chk(data,/type) > 3			; Long or float if float on input
cimg1 = make_array(xsize,ysize,type=typx)	; short or combined short/med
cimg2 = make_array(xsize,ysize,type=typx)	; Final image

; ---------------------------------------------------------------------------
; This sets up m1 (see step 5).
; m1 = where((short_exposure ge 90) and (short_exposure is not saturated) and
;            (long_exposure is saturated))
; ---------------------------------------------------------------------------
if n_elements(img) ne 0 then begin
   cimg1(0,0) = img(*,*,nimg-2) 	; Short (3- or 2-image)
   cimg2(0,0) = img(*,*,nimg-1)		; Long (3-image) or Med (2-image)
   m1 = where((img(*,*,0) ge 90) and 			$; High signal in short, 
		(satpix(*,*,s_sat(0)) eq 0) and 	$; but short not saturated, &
		(satpix(*,*,s_sat(nimg-1)) eq 1),ngood)	 ; saturated in long exposure
   if ngood gt 0 then  m1_data = (img(*,*,0))(m1)	 ; Use the short exposure
endif else begin
;Removed >0 in following 2 statements.  LWA 5/08/95.
   cimg1(0,0) = data(*,*,sx(nimg-2))
   cimg2(0,0) = data(*,*,sx(nimg-1))
   m1 = where((data(*,*,sx(0)) ge 90) and		$; High signal in short,
		(satpix(*,*,s_sat(0)) eq 0) and 	$; but short not saturated, &
		(satpix(*,*,s_sat(nimg-1)) eq 1),ngood)	 ; saturated in long exposure
   if ngood gt 0 then  m1_data = (data(*,*,sx(0)))(m1)	 ; Use the short exposure
endelse

; ---------------------------------------------------------------------------
; This is Step 1 for 3-image case:  Insert shortest into the medium exposure
; ---------------------------------------------------------------------------
if nimg eq 3 then begin
  s1 = where(satpix(*,*,s_sat(nimg-2)),npix)	; s1 = saturated pixels in med. exposure
  if npix gt 0 then begin
    if n_elements(img) ne 0 then 	$
	cimg1(s1) = ((img(*,*,nimg-3))(s1)    * t0(nimg-2)/t0(nimg-3)+.5) else $
;Removed >0 in following statement.  LWA 5/08/95.
	cimg1(s1) = (((data(*,*,sx(nimg-3)))(s1))*t0(nimg-2)/t0(nimg-3)+.5)
  endif
endif
delvarx,img					; Finished with img array

; ---------------------------------------------------------------------------
; This is Step 2: Replace saturated pixels in long exposure with 
;                 short (2-image) or combined med/short (3-image)
; ---------------------------------------------------------------------------
s2 = where(satpix(*,*,s_sat(nimg-1)),npix)	; Saturated pixels in long exposure
if npix gt 0 then 		$		; Replace saturated pixels in long exposure
  cimg2(s2) = cimg1(s2) * t0(nimg-1) / t0(nimg-2) + .5
delvarx,cimg1					; Finished with short/med array

; ---------------------------------------------------------------------------
; This section represents Steps 3 and 4:
; ---------------------------------------------------------------------------
; Make smoothed array and use it to replace bleed pixel
if n_elements(smooth_factor) eq 0 then smooth_factor = 5
if smooth_factor gt 1 then begin
  if fix(smooth_factor/2)*2 eq smooth_factor then message,'Warning: Smooth_factor is not an odd number = '+strtrim(smooth_factor,2),/info
  if (npix gt 0) then begin			;MDM put conditional 1-Dec-93
    timg = smooth(cimg2,smooth_factor)		; Boxcar smooth over smooth_factor pixels
    cimg2(s2) = timg(s2)
  endif						; npix gt 0
endif

; ---------------------------------------------------------------------------
; Implement Step 5:
; ---------------------------------------------------------------------------
if ngood gt 0 then cimg2(m1) = (m1_data * t0(nimg-1)/t0(0) + .5)

if keyword_set(despike) then cimg2=de_spiker(cimg2,despike)
if keyword_set(plot) then tvscl,alog(cimg2>1) 
if keyword_set(sfd) then cimg2 = sfd_comp(cimg2,index=index(sx(nimg-1))) ; Scale as an sfd image 
if keyword_set(qtest) then stop 

;---------------------------------------------------------------------------
; This section builds the comp_index structure.
; ** Still need to include index.his.Q_COMPRESSION once I find out how
;	is it used.   LWA 2-Nov-96 **
;---------------------------------------------------------------------------
comp_index=index(sx(nimg-1))			;Longest exposure.
his_index, /enable
if nimg eq 2 then begin
   time_compos=[gt_time(index(sx(0))),gt_time(index(sx(1))),0]
   day_compos=[gt_day(index(sx(0))),gt_day(index(sx(1))),0]
endif else begin
   time_compos=[gt_time(index(sx(0))),gt_time(index(sx(1))),$
	gt_time(index(sx(2)))]
   day_compos=[gt_day(index(sx(0))),gt_day(index(sx(1))),gt_day(index(sx(2)))]
endelse

his_index, comp_index, 0, 'q_composite', progverno*1000
his_index, comp_index, 0, 'time_compos', time_compos
his_index, comp_index, 0, 'day_compos', day_compos

return,cimg2 
end