function leak_sub, index, data, index_out, $
        udata_in=udata_in, udata_out=udata_out, $ ;jmm, 12-feb-95
	update_index=update_index, $
	leak_only=leak_only, leak_index=leak_index0, save=save, qdebug=qdebug, $
	dc_interpolate=dc_interpolate, dc_scalar=dc_scalar, float=float, $
	force_darksub=force_darksub, orbit_correct=orbit_correct, $
	noleak=noleak, second_order=second_order, synthetic_sfc=synthetic_sfc, $
        grill=grill, nogrill=nogrill, force_standard_sfc=force_standard_sfc, $
        loud=loud
;
;+
;NAME:
;	leak_sub
;PURPOSE:
;	Subtract the Scattered light image from input FFI images.  This is
;	needed for Thin-Al and AlMg images taken after 13-Nov-92. 
;CALLING SEQUENCE:
;	data_out = leak_sub(index, data, index_out)
;	data_out = leak_sub(index, data, /update_index)
;INPUT:
;	index	- The index structure for each image
;	data	- The data array (NX x NY x N).  It should have had the
;		  dark current already subtracted.  If the input is
;		  byte type, then DARK_SUB is called.
;OUTPUT:
;	data_out- The leak (scattered light) subtracted image
;	index_out- The modified index header
;OPTIONAL KEYWORD INPUT:
;       udata_in - the uncertainty in the data array, passed through to DARK_SUB
;       udata_out(INPUT/OUTPUT) - The uncertainty in the output array, passed out of DARK_SUB,
;                   This must be set to a variable before the call, i.e.
;                   YOU (the user) have a statement like "udata_out=1"
;                   before calling LEAK_SUB with the udata_out keyword, the calculation
;                   will not be done if a zero or undefined value is passed in.
;                   Note that leak_only data "officially" has no uncertainty,
;                   So don't use this with the leak_only keyword, you will get an array of zeros
;	leak_only- If set, the output is only the scattered light
;		  images that go with each index.
;	save	- If set, store the data in a common block to avoid 
;		  having to read it from memory many times.
;	dc_interpolate - If set, do integration time interpolation on
;		  the dark current subtraction images.
;	dc_scalar - If set, use a scalar value for the dark current
;		  subtraction
;	force_darksub - If set, then call DARK_SUB even if the data type 
;		  is non-byte.  This is needed with restore low-8 data.
;       noleak - if set, don't subtract leak(!) [default:  noleak=0]
;		 (seems silly in a routine called
;                leak_sub, but needed for leak bybass from callers
;                like sxt_prep...)
;       synthetic_sfc - if set, generate and apply synthetic SFC
;                       [ hook for P.G.Shirts routines ]
;       force_standard_sfc - if set, override synthetic_sfc keyword 
;                            and environmental (used for recursive calls
;                            when synthetic SFC not available or invalid
;       loud - if set, print some ~useless albeit historical info
;
;OPTIONAL KEYWORD OUTPUT:
;	leak_index- Return the index of the reference images
;METHOD:
;	1. Call DARK_SUB if "data" is byte type
;	2. For PFI, extract the relevant portion of the images
;	3. Scale the scattered light image for exposure duration.
;	4. Subtract off the scattered light
;HISTORY:
;	Written 7-Jan-93 by M.Morrison
;	14-Jan-93 (MDM) - Modification to perform leak subtraction 
;			  properly for quarter resolution images.
;	 3-Feb-93 (MDM) - Modified to correct AlMg leak images also.
;			- Added /SAVE option to store the leak image
;			  in a common block - avoids reading it multiple
;			  times, but uses memory up.
;	 9-mar-93 (JRL) - Fixed a bug in case no leak images are available.
;	 7-Apr-93 (MDM) - Modified to work with PFI images
;			  Adjusted code to do REBIN only once per uniq resolution
;			  (faster/more memory efficient)
;			- Added check to see if leak calibration images exist.
;			- Corrected an error that leak subtraction would not
;			  have been performed if the neutral density filter was
;			  in place!!
;V2.0	 6-Jun-93 (MDM) - Added version number, and started calls to HIS_INDEX
;V2.1	 9-Jul-93 (MDM) - Corrected an error caused when a smaller image is 
;			  imbedded in a larger array (256x256 inside a 512x512xN)
;V2.11	 3-Aug-93 (MDM) - Modification for handling unassembled ORs
;V3.00	 4-Aug-93 (MDM) - Modified to use EXT_SUBIMG2 which is much faster than
;			  rebinning the whole image.
;			- Added INDEX_OUT parameter
;			- Do not do a "fix" of data when it is non-byte type
;			  when creating the DATA_OUT variable
;V3.01	26-Aug-93 (MDM) - Modified to use ALIGN1IMG instead of EXT_SUBIMG2
;			- Added /UPDATE_INDEX option
;	26-Aug-93 (MDM) - Replaced Al.1 image sfc921127.0947 with sfc930819.1607
;				** V3.xx was never put on-line **
;V4.00	27-Aug-93 (MDM) - Replaced code with the essence of DARK_SUB
;V4.01  11-Oct-93 (MDM) - Made sure that the output index was passed to 
;			  HIS_INDEX so that the .HIS structure was appended
;V4.02	13-Oct-93 (MDM) - Added /DC_INTERPOLATE and /DC_SCALAR options
;			- Adjusted code since the DARK_SUB history was not being
;			  preserved
;V4.03	 8-Feb-95 (MDM) - Added /FLOAT keyword to be passed through to DARK_SUB
;	12-Feb-95 (jmm) - Added the UDATA keywords
;V4.04	27-Feb-95 (MDM) - Put McTiernan modifications on-line
;V4.05	23-Mar-95 (MDM) - Added /force_darksub
;V4.06	 9-May-95 (MDM) - Modified how the /FLOAT option worked (it was always
;			  doing a FIX before on the leak image).
;			- Use ROUND instead of FIX when not going to floating
;			  output
;V4.07   11-Feb-98 (SLF) - add ORBIT_CORRECT (pass to dark_sub)
;        26-Oct-98 (SLF) - add /NOLEAK keyword and function
;V4.08   10-Jan-2000 - S.L.Freeland add /SECOND_ORDER keyword and function
;                      Add history bit
;V4.09   23-Aug-2000 - Add /SYNTHETIC_SFC keyword and function
;                      Add history bit
;V4.10    8-Jan-2001 - S.L.Freeland - add Al.1 Grill correction option
;                      (no longer applied to SFCs)
;V4.20   13-Nov-2001 - S.L.Freeland - remove a call to 'fix' - why it
;                      was there at all is TBD (ask mons...)
;V4.21   14-Jan-2002 - S.L.Freeland - add /FORCE_STANDARD_SFC 
;V4.22    7-Mar-2003 - S.L.Freeland - add /LOUD and made default quieter
;V4.23   30-May-2003 - S.L.Freeland - merge LWA Long Word Overflow fix
;                                     (val=...) per Oct 16 2002 email
;V5.0    24-Jul-2003 - S.L.Freeland - do grill correction for both synthetic
;                                     and standard in same place  
;V5.1    21-Jan-2005 - S.L.Freeland - avoid global side affect with 
;                                     SYNTHETIC_SFC (possible cube processing
;                                     problem with toggle during recursion)
;V5.2    28-Jan-2005 - S.L.Freeland - continue Syn Vs Term mods - fix
;                                     logic flaw introduced in 28-jan-mod
;V5.3     2-Feb-2005 - S.L.Freeland - fixed typo...
;V5.4    21-Jan-2005 - L. W. Acton  - Added code to handle early and late
;                      straylight drifts.
; 		       Boxmessage, 'Synthetic requested but ...' only for /loud. 
;	 28-Mar-2005 - LWA - Commented out early straylight drift commands.
;	 18-May-2006 - LWA - Removed . in 2^15. to make compatible with
;			IDL 6+ on linux machines.
;	 20-Jul-2009 - LWA - Special leak processing for 24-Oct-95 eclipse.
;	 27-Jul-2009 - LWA - Revised special 24-Oct-95 leak processing.
;			Added special leak processing for 17-Apr-96 eclipse.
;	 30-Jul-2009 - LWA - Special leak processing for 26-Feb-98 eclipse.
;V5.5	 31-Jul-2009 - LWA - Updated progverno to 5.5.
;			Commented out all special eclipse processing, moved
;			it to sxt_deleak.pro.
;	29-Jun-2012 - LWA - Corrected noleak in header.
;-

common get_leak_image, lindex, ldata
common get_leak_image2, called
second_order=keyword_set(second_order)
;
;progverno = 5.4*1000
progverno = 5.5*1000
;
n = n_elements(index)
index_out = index
his_index, index_out
;
quncert = keyword_set(udata_out) ;jmm, 12-feb-95
noleak=keyword_set(noleak)
yesleak=1-noleak
force_standard_sfc=keyword_set(force_standard_sfc)
use_syn=keyword_set(synthetic_sfc) or get_logenv('sxt_syn_sfc') ne ''
use_syn=use_syn * (1-force_standard_sfc) ; override synthetic?
force_standard_sfc=0   ; !!only used by recurive call - reset on exit
orbit_correct=keyword_set(orbit_correct)
dc_interpolate=keyword_set(dc_interpolate) or orbit_correct ; FORCE ON
;                                                           ; for orbit corr
nogrill=keyword_set(nogrill)
grill=1-keyword_set(nogrill)
old_grill=get_logenv('leak_sub_oldgrill') ne ''   ; default:use new Grills
;
; Needed later on for straylight drift (LWA, 1/21/2005)
; Fitted slopes in filter number order.  Second 5 are with ND filter.
lateslope=[1.97e-7,4.35e-8,2.64e-7,5.19e-8,1.71e-7,$
        3.00e-6,2.08e-6,2.08e-6,2.27e-6,2.42e-7]
;early_dt=ssw_deltaT('12-JAN-93  07:51:32', '24-OCT-93  01:33:10',/sec)
late_dt=ssw_deltaT(' 6-FEB-99  06:15:46', '14-DEC-01  01:55:09',/sec)
;
if (keyword_set(leak_only)) then begin
    nx = max(gt_shape(index, /x))
    ny = max(gt_shape(index, /y))
    if (keyword_set(float)) then data_out = fltarr(nx,ny,n) $
			else data_out = intarr(nx, ny, n)
    if(quncert) then udata_out = bytarr(nx, ny, n) ;jmm, 12-feb-95
end else begin
    siz = size(data)
    nx = siz(1)
    ny = siz(2)
    typ = siz( siz(0)+1 )

    if (typ eq 1) or (keyword_set(force_darksub)) then begin
       udata_out = quncert      ;jmm, 12-feb-95, udata_out has to be set to be passed into DARK_SUB
       data_out = dark_sub(index, data, index_out, interpolate = dc_interpolate, $
                           force_scalar = dc_scalar, float = float, $
                           orbit_correct=orbit_correct, $ ; slf 11-feb-98
                           udata_in = udata_in, udata_out = udata_out) ;jmm, 12-feb-95
    end else begin
       data_out = data
       if(quncert) then begin   ;jmm, 12-feb-95
          if(keyword_set(udata_in)) then udata_out = udata_in $
            else udata_out = byte(make_array(size = size(data_out)) + 0) ;pass out an array of zeros if nothing is passed in and no dark subtraction.
       end
    end
end
;
if (keyword_set(qdebug)) then print, 'Finished making output array'
;
;---------------------------------------- Deterimine the optimal order of the images to do the subtraction
;
sel_leak_image, index, leak_files, dset
;
xcorn = fix(gt_corner(index, /lower, /x)/4)     ;0-255
ycorn = fix(gt_corner(index, /lower, /y)/4)
res = fix(gt_res(index))

val = xcorn +  ycorn*256UL + res*256UL*256UL  + ulong64(dset)*256UL*256UL*64       ;dset specifies different leak image

;      low 8-bits 0:7      second 8-bits  8:15          third 8-bits     fourth 8-bits
;---- do all FR in order of DPE, then all HR, then all QR
;
ss = sort([val])
;
last_val = -999
;
qdebug=keyword_set(qdebug)
loud=keyword_set(loud) or qdebug
if yesleak then begin 

for iimg=0,n-1 do begin
    i = ss(iimg)		;do the subtraction in order of DPE so that
    ;
    synstatus=0
    if (dset(i) ne -1) then begin
	case 1 of 
           use_syn: begin
		   delvarx,leak_index0,leak_data
              mk_syn_sfc, index(i), leak_index0, leak_data , synstatus=synstatus
              if not synstatus then begin 
                 if loud then box_message,'Synthetic requested but bad or no coefficients.'
		 get_leak_image,index(i),leak_index0,leak_data,qdebug=qdebug,/save
		 endif else get_leak_image, index(i), leak_index=leak_index0, leak_data=leak_data
            endcase
            else: if (val(i) ne last_val) then get_leak_image, index(i), leak_index0, leak_data, qdebug=qdebug, /save
        endcase	    

        if (keyword_set(qdebug)) then print, 'Using leak image ', fmt_tim(leak_index0), ' for val=', val(i)
	

	factor = (gt_expdur(index(i)) /  ([gt_expdur(leak_index0),2691.09])(synstatus) )
        if get_logenv('synsfc_debug') ne '' then stop,'syndebug'
	if (keyword_set(qdebug)) then print, 'LEAK_SUB: Exposure normalization factor: ', factor
        algrill=sxt_get_grill(index(i),status,loud=loud, geny=1-(old_grill)) * grill
; Note: algrill will be all zeros if non Al.1 or problem with
; grill file access (so null effect in those cases)
;
; Additions here by LWA to handle straylight drift.  (2/11/2005)
        if get_logenv('synsfc_debug') ne '' then goto, SKIPDRIFT
        totleak=leak_data * factor 
;        imgtime=int2secarr(index(i),'12-JAN-93')
;        leaktime=int2secarr(leak_index0,'12-JAN-93')
        imgtime=int2secarr(index(i),'6-Feb-99 06:15:46')
        leaktime=int2secarr(leak_index0,'6-Feb-99 06:15:46')
        case 1 of
           gt_res(index(i)) eq 0 : slpfact=0.25
           gt_res(index(i)) eq 1 : slpfact=1.0
           gt_res(index(i)) eq 2 : slpfact=4.0
        endcase
        slpfact=slpfact*gt_expdur(index(i))/1000.       ;Exposure correction
;        if imgtime lt 0 then begin
;           if gt_filta(index(i)) eq 6 then begin        ;ND filter case
;              slp=slpfact*earlyslope(gt_filtb(index(i))-2+5)
;           endif else begin
;              slp=slpfact*earlyslope(gt_filtb(index(i))-2)
;           endelse
;           totleak=totleak+((imgtime-leaktime)*slp)
;        endif
;        if imgtime gt 1.9148495e+08 then begin  ;Start of late drift period, 6-FEB-99 06:15:46
        if imgtime gt 0 then begin  ;Start of late drift period, 6-FEB-99 06:15:46
           if gt_filta(index(i)) eq 6 then begin        ;ND filter case
              slp=slpfact*lateslope(gt_filtb(index(i))-2+5)
           endif else begin                             ;FiltA open case
              slp=slpfact*lateslope(gt_filtb(index(i))-2)
           endelse
           totleak=totleak+((imgtime-leaktime)*slp)
        endif
        SKIPDRIFT :
        if (keyword_set(leak_only)) then begin
           if get_logenv('synsfc_debug') ne '' then data_out(0,0,i) = leak_data * factor $
              else data_out(0,0,i) = totleak
        endif else begin
            siz = size(leak_data)
            nx2 = siz(1)-1
            ny2 = siz(2)-1
	    ;;data_out(0,0,i) = data_out(0:nx2,0:ny2,i) - fix(leak_data * factor+algrill)
            if get_logenv('synsfc_debug') ne '' then begin
	       if (data_type(data_out) ge 4) then data_out(0,0,i) = data_out(0:nx2,0:ny2,i) - (leak_data * factor+algrill) $
	    			else data_out(0,0,i) = data_out(0:nx2,0:ny2,i) - round(leak_data * factor + algrill)
            endif else begin
	       if (data_type(data_out) ge 4) then data_out(0,0,i) = data_out(0:nx2,0:ny2,i) - (totleak+algrill) $
	    			else data_out(0,0,i) = data_out(0:nx2,0:ny2,i) - round(totleak + algrill)
            endelse
; End of straylight drift additions and changes.
	    his_index, index_out, i, 'time_leak', leak_index0.his.time_leak
	    his_index, index_out, i, 'day_leak',  leak_index0.his.day_leak
            pver=progverno
            if second_order then begin
	       if gt_res(index_out(i)) ge 1 and $
		  gt_pfi_ffi(index_out(i),/ffi) then begin
	            if n_elements(called) eq 0 or keyword_set(qdebug) then $
		       box_message,'Applying SECOND ORDER leak correction'
	               so_data=sxt_deleak(index_out(i), data_out(*,*,i), leak_index0, leak_data)
	               data_out(0,0,i)=so_data
;                       pver = progverno + 2^15.          ; set history flag
                       pver = progverno + 2^15          ; Removed dot. LWA 5/18/2006 
	       endif 
	    endif  
            pver=pver+([0.,2^14.])(use_syn)        ; set history flag
	    his_index, index_out, i, 'q_leak_sub', pver
       endelse
   end

;   ; Special leak handling for 24-Oct-95 eclipse.  LWA  20-Jul-09/27-Jul-09
;   ; First transit.
;   if anytim(index_out(i)) gt anytim('24-OCT-95  03:22:00') and anytim(index_out(i)) $
;      lt anytim('24-OCT-95  03:48:00') and gt_pfi_ffi(index_out(i)) eq 1 then begin
;      nsfc=total(leak_data(240:253,240:253))
;      ndat=total(data_out(240:253,240:253,i))
;      data_out(0:255,0:255,i)=data_out(0:255,0:255,i)-((ndat/nsfc)*leak_data)
;   endif
;   ; Second transit.
;   if anytim(index_out(i)) gt anytim('24-OCT-95  05:34:00') and anytim(index_out(i)) $
;      lt anytim('24-OCT-95  05:53:00') and gt_pfi_ffi(index_out(i)) eq 1 then begin
;      nsfc=total(leak_data(10:23,10:23))
;      ndat=total(data_out(10:23,10:23,i))
;      data_out(0:255,0:255,i)=data_out(0:255,0:255,i)-((ndat/nsfc)*leak_data)
;   endif
;   ; Special leak handling for 17-Apr-96 eclipse.  LWA  27-Jul-09
;   if anytim(index_out(i)) gt anytim('17-Apr-96 22:04') and anytim(index_out(i)) $
;      lt anytim('18-Apr-96 00:08') and gt_pfi_ffi(index_out(i)) eq 1 then begin
;      nsfc=total(leak_data(240:253,240:253))
;      ndat=total(data_out(240:253,240:253,i))
;      data_out(0:255,0:255,i)=data_out(0:255,0:255,i)-((ndat/nsfc)*leak_data)
;   endif
;   ; Special leak handling for 26-Feb-98 eclipse.  LWA  30-Jul-09
;   if anytim(index_out(i)) gt anytim('26-Feb-98 18:50:34') and anytim(index_out(i)) $
;      lt anytim('26-Feb-98 18:59:14') and gt_pfi_ffi(index_out(i)) eq 1 then begin
;      nsfc=total(leak_data(10:23,10:23))
;      ndat=total(data_out(10:23,10:23,i))
;      data_out(0:255,0:255,i)=data_out(0:255,0:255,i)-((ndat/nsfc)*leak_data)
;   endif
;   ; Special leak handling for 16-Feb-99 eclipse.  LWA  30-Jul-09
;   if (anytim(index_out(i)) gt anytim('16-Feb-99 08:57:00') and $
;      anytim(index_out(i)) lt anytim('16-Feb-99 09:06:00'))  or $
;      (anytim(index_out(i)) gt anytim('16-Feb-99 09:07:00') and $
;      anytim(index_out(i)) lt anytim('16-Feb-99 09:09:00')) and $
;      gt_pfi_ffi(index_out(i)) eq 1 then begin
;      nsfc=total(leak_data(480:506,20:46))
;      ndat=total(data_out(480:506,20:46,i))
;      data_out(0:511,0:511,i)=data_out(0:511,0:511,i)-((ndat/nsfc)*leak_data)
;   endif

    last_val = val(i)
    called=1
endfor
endif      ; end of Leak inhibit

;
if (not keyword_set(save)) then begin
    lindex = 0	;clear the GET_LEAK_SUB common block variables
    ldata = 0
end
;
if (keyword_set(qstop)) then stop
if (keyword_set(update_index)) then index = index_out
;stop
return, data_out
end