function sxt_composite,index,data,comp_index,unc_out=unc_out,           $
	sx=sx,satpix=satpix,unc_in=unc_in,plot=plot,nofill=nofill,      $
        sfd=sfd, dc_interpolate=dc_interpolate,qtest=qtest,             $
        register=register,clean=clean,fillpix=fillpix,                  $
        threshold=threshold,debug=debug, fillmap=fillmap
;+
;
;*************************************************************************
;*******This version returns an uncertainty array, plus photon noise******
;*******DMcK  12-Jan-2001*************************************************
;*************************************************************************
;
; 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,comp_index)
;   cimg = sxt_composite(index,data,[comp_index,unc=unc,sx=sx,           $
;        satpix=satpix,plot=plot,nofill=nofill,sfd=sfd,                  $
;        dc_interpolate=dc_interpolate,qtest=qtest,register=register,    $
;        clean=clean,fillpix=fillpix,threshold=threshold,debug=debug])
; 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, 
;						background subtracted, and
;						FLOATING-POINT.
; 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 KEYWORD OUTPUT:
;   unc_out = array of 1-sigma uncertainties in each pixel.  Uses
;            sxt_dn_uncert for a temperature of 2.3 MK.  If data are
;	     byte type, then this array will come from an internal call
;	     call to SXT_PREP.  If data are not byte type, and no UNC_IN
;	     is passed in, then this will _only_ be the photon noise; if
;	     data are not byte type and UNC_IN is passed in, then they
;	     will be added in quadrature.
;   fillpix = Boolean array showing interpolated pixels.o
;   fillmap = fill map where 1: filled from short
;                            2: filled from medium
;
; 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 are not byte type).
;	     Must be 3-d and correspond to data(*,*,*).
;   unc_in = Array of uncertainties of the input images, e.g., from
;	     SXT_PREP.  It is assumed that this UNC_IN does _not_ come
;            from SXT_DN_UNC.  (Only used if input data are not byte type).
;   plot   = If set, plot the results as we go.  The results are displayed
;	     on a logarithmic scale.
;   nofill = set this keyword to PREVENT filling bleeds w/ interpolated data.
;   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
;   clean  = calls sxt_clean with defaults to despike output image.
;   threshold = if data are byte type this is the signal level in the
;	     shorter exposure below which the saturated regions of the
;            longer exposures are filled in with interpolated data for
;            statistical appearance sake.  (Default is 15 DN.)
;   debug  = stops program at end for debugging purposes.
; RESTRICTIONS:
;   Must be called separately for each composite output image.
;   It is assumed that any input uncertainty comes from SXT_PREP, but
;   _not_ from SXT_DN_UNC.  Photon noise will be added!
; METHOD:
; ---------------------------------------------------------------------------
;1.  Replace the saturated pixels in the medium exposure with the
;    corresponding above-threshold time-normalized pixels from the 
;    short exposure.
;2.  Replace the saturated pixels in the medium exposure with interpolated
;    corresponding below-threshold time-normalized pixels from the
;    short exposure.
;
;    This is completes the process for the 2 image case.
;
;3.  For the 3 image case replace the saturated pixels in the long exposure
;    with the corresponding above-threshold time-normalized pixels from
;    step 2.
;4.  Replace saturated pixels in long exposure with interpolated
;    corresponding below-threshold time-normalized pixels from step 2.
; ---------------------------------------------------------------------------
; 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.
;    20-Apr-99, LWA - Added /float and /orbit_correct to leak_sub call.
;                   - Version 3.01
;    14-Jul-99, BNH - history record wasn't being updated from LEAK_SUB.
;                   - Bumped version to 3.02.
;    11-Apr-00, LWA - Complete rewrite and improvement of logic.
;		    - Added unc, smoothpix and threshold keywords.
;		    - Changed version to 4.00.
;    11-Apr-00, LWA - Corrected typo and bumped version to 4.01.
;                   - Added keyword debug and bumped version to 4.02.
;    10-Oct-00, LWA - Replaced smoothing with sxt_suture.pro.
;		    - Changed version to 5.00.
;		    - This version does NOT return revised uncertainty.
;    12-Jan-01, DMcK & LWA - Carried the uncertainties through the 
;		      compositing process, and added option of input
;		      UNC_IN.  Also cut out deadwood and fixed some typos,
;		      and changed default threshold to 15.
;		      Changed version to 6.00.
;    15-Jan-01, SLF - add FILLMAP (flags for SSC uncertainty info storage)
;-

progverno=6.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
;-------------------------

if NOT keyword_set(threshold) then threshold = 15 $
	else threshold = threshold

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 

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 

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

; Turn on history for this index record.
comp_index=index(sx(nimg-1))
his_index, /enable

; 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

; ----
; Run SXT_PREP to process the input data as required.
; ----
d_typ = sz0(sz0(0)+1)
if d_typ eq 1 then begin
   sxt_prep,index(sx),data(*,*,sx),iout,img,prep_unc,satp,$
;	/second_order_leak, $
	/dc_interpolate,/dc_orbit_correct,$
	register=register,ref_image=index(sx(nimg-1)),$
        sxt_cleanx=clean,$   
	/float
   t0 = gt_expdur(index(sx))
endif else begin
   if n_elements(satpix) eq 0 then begin   ;Check for satpix array
      print,'***  Error in sxt_composite: MUST PROVIDE SATPIX FOR NON-BYTE DATA'
      return,-1
   endif
  unc_out=unc_out
   if keyword_set(unc_out) then begin        ; If UNC_OUT is desired, then
    if n_elements(unc_in) eq 0 then begin    ; are there input UNC_IN ?
      print,'***  SXT_COMPOSITE: There are no input uncertainties.  ***'
      print,'    The output uncertainties will ONLY include photon noise'
    endif else begin
      print,'    The input uncertainties will be added to photon noise,',$
	' in quadrature.'
    endelse
   endif
   sxt_prep,index(sx),data(*,*,sx),iout,img,register=register,$
	sxt_cleanx=clean,$    
        ref_image=index(sx(nimg-1)),/float
   t0 = gt_expdur(index(sx))	; t0 = [short, med, long]
   satp = satpix(*,*,sx)
endelse

; ----
;  Set up the plotting window and display the input images.
; ----
xsize = sz0(1)			;LWA  4/10/00
ysize = sz0(2)			;LWA  4/10/00
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_typ eq 1) then tvscl,data(*,*,sx(0)),0,ysize $
        else tvscl,alog(data(*,*,sx(0))>1),0,ysize
  if (d_typ 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_typ eq 1) then tvscl,data(*,*,sx(2)),xsize,0 $
          else tvscl,alog(data(*,*,sx(2))>1),xsize,0
  endif
endif

; -----
; Prepare the uncertainty array
; -----
if d_typ eq 1 then begin
   uncert = sxt_dn_unc(iout,img,te=alog10(2.3e6),/float,prep_unc)
endif else begin
   if n_elements(unc_in) eq 0 then begin
      uncert = sxt_dn_unc(iout,img,te=alog10(2.3e6),/float)
   endif else begin
      uncert = sxt_dn_unc(iout,img,te=alog10(2.3e6),/float,unc_in(*,*,sx))
   endelse   
endelse

; -----
; Prepare the composite images
; -----
   cimg0 = img(*,*,0)
   cimg1 = img(*,*,1)
   ii0 = where(cimg0 gt threshold and satp(*,*,1) eq 1,nii0)
   cimgA = img(*,*,1)
   cuncA = uncert(*,*,1)
   if nii0 gt 0 then begin 
	cimgA(ii0) = (img(*,*,0))(ii0)*t0(1)/t0(0)
	cuncA(ii0) = (uncert(*,*,0))(ii0)*t0(1)/t0(0)
   endif
   ii00 = where(cimg0 le threshold and satp(*,*,1) eq 1,nii00)
   if nii00 gt 0 then begin
      if NOT keyword_set(nofill) then begin
         fill_imgA = sxt_suture(img(*,*,1),satp(*,*,1),uncert(*,*,1),$
		unc_out=uncert_outA)
	 cimgA(ii00) = fill_imgA(ii00)
	 cuncA(ii00) = uncert_outA(ii00)
      endif else begin
         cimgA(ii00) = (img(*,*,0))(ii00)*t0(1)/t0(0)
	 cuncA(ii00) = (uncert(*,*,0))(ii00)*t0(1)/t0(0)
      endelse
   endif
   cimg = cimgA
   cunc = cuncA

   if nimg eq 3 then begin
      ii1 = where(cimgA gt threshold and satp(*,*,2) eq 1,nii1)
      cimgB = img(*,*,2)
      cuncB = uncert(*,*,2)
      if nii1 gt 0 then begin
	 cimgB(ii1) = cimgA(ii1)*t0(2)/t0(1)
	 cuncB(ii1) = cuncA(ii1)*t0(2)/t0(1)
      endif
      ii10 = where(cimgA le threshold and satp(*,*,2) eq 1,nii10)
      if nii10 gt 0 then begin
         if NOT keyword_set(nofill) then begin
            fill_imgB  = sxt_suture(img(*,*,2),satp(*,*,2),uncert(*,*,2),$
		unc_out=uncert_outB)
            cimgB(ii10) = fill_imgB(ii10)
	    cuncB(ii10) = uncert_outB(ii10)
         endif else begin
            cimgB(ii10) = cimgA(ii10)*t0(2)/t0(1)
	    cuncB(ii10) = cuncA(ii10)*t0(2)/t0(1)
         endelse
      endif
      cimg = cimgB
      cunc = cuncB
   endif

unc_out = cunc		; 12-Jan-01 DMcK

; -----
; Prepare the array showing which pixels have been filled.
; -----
fillpix = bytarr(sz0(1),sz0(2))
if nii00 gt 0 and NOT keyword_set(nofill) then fillpix(ii00) = 1
fillmap=fillpix
if nimg eq 3 then begin
   if nii10 gt 0 and NOT keyword_set(nofill) then begin
      fillpix(ii10) = 1
      fillmap(ii10) = 2
   endif
endif

if keyword_set(plot) then tvscl,safe_log10(cimg)
if keyword_set(sfd) then cimg = sfd_comp(cimg,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 **
;---------------------------------------------------------------------------
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

if keyword_set(debug) then stop

return,cimg 

end