FUNCTION destr_bilin2,data,Slit,RefPix,Det,Quiet=Quiet

;
;+
;Project:
;       SOHO - SUMER
;Name:
;                 destr_bilin
;Purpose:
;      Apply distortion correction to sumer data
;   ++ derived from destretchn.pro of T. Moran modif. by L.A. Smaldone
;               on Feb,6 1997
;                                         
;Use:
;       data_out = destr_bilin(data,slit,RefPix,detector)
;Inputs
;      Data - data to be destretched. 1 2-d array
;      Slit - Number of Slit used
;      RefPix - Refernce Pixel for the wavelength. Should be either
;	the actual reference pixel or the reference co-pixel +1
;	(the latter is in the binary fits header)   
;      Det - Detector type, 'a' or 'b'
;Keyword Input:
;      Quiet - supress the general information printed out at start of program
;Keyword Output:
;      delta - the delta arrray
;Output
;       Data_Out - destretched data
;Restrictions:
;      This program should NOT be used on data which has been binned or
;	subject to compression schemes 6 and up.
;
; The delta x and y  files must be in a directory called by the logical
; name 'sumer_calib'.
; For example, to set the logical name from within IDL, 
; IDL> setenv,'SUMER_CALIB = files_directory ' , where the delta files are
; contained in files_directory. 
;
;Written:
;       Tom Moran,    July 15, 1996
;	Changes Oct, 3 1996:
;       Corrected in case where full 1024 pixels used
;       Comon block added & some of file reading logic changed
;       Quiet keyword introduced
;	Changes Dec 3, 1996:
;	Program can now destretch data from either detector A or B
;	Files called by Destretch must now be in directory defined
;	    by logical SUMER_CALIB. These files are:
;	    del_x_at.dat, del_x_bt.dat, del_y_at.dat, del_y_bt.dat, 
;	    size_fac.dat
;	Deals correctly with cases in which 120" slit is
;	    used with 360 Y-direction formats
;	Changes Feb 6, 1997 by L.A. Smaldone :
;            substitute spline interpolations by bilinear interpolation to
;                       speed-up the routine
;-
;      

common deltafile_com, Del_X_file_old,Del_Y_file_old,i_del_x,i_del_y, i_pix_size


;determin size
dsize=size(data)
dY =  dsize(2)
dW = dsize(1)

up = 299
center = 179
down = 59

case slit of
     1: SlitLoc = center
     2: SlitLoc = center
     3: SlitLoc = up
     4: SlitLoc = center
     5: SlitLoc = down
     6: SlitLoc = up
     7: SlitLoc = center
     8: SlitLoc = down
     9: SlitLoc = center
 else: message,'Problem with slit'
endcase

case slit of
     1: f_t = 0
     2: f_t = 0
     3: f_t = 1
     4: f_t = 0
     5: f_t = 1
     6: f_t = 1
     7: f_t = 0
     8: f_t = 1
     9: f_t = 0
endcase

if dY eq 360 then YRange = [0,359] else YRange=SlitLoc+[-dY/2+1,dY/2]

;wavelength direction range

if dW eq 1024 then Wrange=[0,1023] else begin
   if dW mod 2 eq 0 then WRange = [refpix-1,refpix-1]+[-dW/2+1,dW/2] $
                 else WRange = [refpix-1,refpix-1]+[-dW/2,dW/2]
endelse

;;;;;;;;;;;;;;;;;;; check if WRange is inside [0,1023] ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

 if WRange(0) lt 0 then data1=data(dW-1-WRange(1):dW-1,YRange(0):YRange(1)) else $
 if WRange(1) gt 1023 then data1=data(0:1023-WRange(0),YRange(0):YRange(1)) else data1=data
 WRangeff=WRange>0
 WRangeff=WRangeff<1023
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
      
; read in delta x,y arrays
; Test to determine if SUMER_CALIB is defined
dir = getenv('SUMER_CALIB')
;;dir = 'contrib:[moran]'
  if dir eq '' then message,$
	'Environment variable SUMER_CALIB must be defined'
			;Get names of Destretch data files
  case STRUPCASE(det) of
     'A': begin
         Del_x_file = concat_dir(dir,'del_x_at.dat') 
         Del_y_file = concat_dir(dir,'del_y_at.dat')
          end
     'B': begin
         Del_x_file = concat_dir(dir,'del_x_bt.dat')
         Del_y_file = concat_dir(dir,'del_y_bt.dat')
          end
     else:  message,"Invalid detector: must be 'A' or 'B'" 
  endcase
  size_fac_file = concat_dir(dir,'size_fac.dat')
				;Test to see if files exist
  test = loc_file(Del_x_file,count=c)
  if c eq 0 then message,"Can't find file "+Del_x_file
  test = loc_file(Del_y_file,count=c)
  if c eq 0 then message,"Can't find file "+Del_y_file
  test = loc_file(size_fac_file,count=c)
  if c eq 0 then message,"Can't find file "+Size_fac_file
 				;Define the following if they don't exist
  if n_elements(Del_X_file_old) eq 0 then Del_X_file_old = ''
  if n_elements(Del_Y_file_old) eq 0 then Del_Y_file_old = ''
 

;if arrays have not already been read or the data file names are different
;	from last time, read in correction arrays.

  if (n_elements(i_del_x) eq 0)   or  $
          (Del_X_file_old ne Del_X_file) or (Del_Y_file_old ne Del_Y_file) $
    then begin
    i_del_x = intarr(1024,360)
    i_del_y = intarr(1024,360)
    i_pix_size = intarr(1024,360)

    del_x = fltarr(1024,360) 
    del_y = fltarr(1024,360)
    pix_size = fltarr(1024,360)

    print,'Reading in del_x, del_y files - takes about 1. minute'
    print,' '

     close,1
     openr,1,Del_X_file
     readf,1, i_del_x
     close,1

     openr,1,Del_Y_file 
     readf,1, i_del_y
     close,1

     openr,1,size_fac_file
     readf,1, i_pix_size
     close,1
    
     Del_X_file_old = Del_X_file
     Del_Y_file_old = Del_Y_file
  endif

   del_x = i_del_x/1000.
   del_y = i_del_y/1000.
   pix_size = i_pix_size/1000.

    i_frame = fltarr(1024,360)
    o_frame = fltarr(1024,360)
    x_axis = indgen(1024)
    y_axis = indgen(360)

    i_frame(WRangeff(0):WRangeff(1),YRange(0):YRange(1)) = data1

    i1 = (wrange(0) - 10) > 0
    i2 = (wrange(1) + 10) < 1023
    j1 = (yrange(0) - 10) > 0
    j2 = (yrange(1) + 10) < 359

    sub_frame = i_frame(i1:i2,j1:j2)

    sub_x = x_axis(i1:i2)
    sub_y = y_axis(j1:j2)

    sub_dely = del_y(i1:i2,j1:j2)
    sub_del = del_x(i1:i2,j1:j2)
    sub_size = pix_size(i1:i2,j1:j2)

    sub_t_frame = sub_frame

   c = ''

; contruct the new frame


      for j=0,j2-j1-1 do begin
         sub_del(*,j)= sub_x + sub_del(*,j)-i1
      endfor 

     for i=0,i2-i1-1 do begin

       sub_dely(i,*)= sub_y + sub_dely(i,*)-j1

     endfor
       sub_t_frame = bilin_interp(sub_frame,sub_del,sub_dely,0)

; correcting for pixel size
     sub_t_frame = sub_t_frame * sub_size
; creating full frame with subframe enclosed' 
      o_frame(i1:i2,j1:j2) = sub_t_frame

      o_frame(0:4,7) = 0.
      o_frame = o_frame > 0.
     
      o_frame(*,0:5) = 0.
      o_frame(*,356:359) = 0.

;return,o_frame   

;;;;;;;;;; cut out the original image format   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

data_out=fltarr(dW,dY)
if WRange(1) gt 1023 then $
data_out(0:1023-WRange(0),YRange(0):YRange(1))=o_frame(WRangeff(0):WRangeff(1),YRange(0):YRange(1)) else $
if WRange(0) lt 0 then $
data_out(dW-1-WRange(1):dW-1,YRange(0):YRange(1))=o_frame(WRangeff(0):WRangeff(1),YRange(0):YRange(1)) else $
      data_out=o_frame(WRange(0):WRange(1),YRange(0):YRange(1))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;      
return,data_out
end
      
;    ss = float(s)


; optimized version
;    ss00 = ss(x0,y0)
;    ss01 = ss(x0,y1)
;    ssfx = (ss(x1,y0)-ss00)*fx
;    ans  = ss00 + ssfx + (ss01-ss00 + (ss(x1,y1)-ss01)*fx - ssfx)*fy
;    endelse

;return, ans
;end

;+
; NAME:
;
; sumer_destretch
;
; PURPOSE:
;
; Calls Tom Moran's destretch routine which does sub-pixel corrections
; for optical distortion for SOHO/SUMER observations.
;
; CATEGORY:
;
; Data analysis.
;
; CALLING SEQUENCE:
;
; sumer_destretch,index,data
; 
; INPUTS:
;
; The index and data structures returned from rd_sumer.
;
; OPTIONAL INPUTS:
;
; None.
;       
; KEYWORD PARAMETERS:
;
; verbose: if set an explanation of destretch will be printed 
;       and  program will keep user informed of progess
;
; OUTPUTS:
;
; The input data is overwritten with the corrected data.
;
; OPTIONAL OUTPUTS:
;
; None.
;
; COMMON BLOCKS:
;
; None.
;
; SIDE EFFECTS:
;
; None.
;
; RESTRICTIONS:
;
; The corrected data may correspond to a slightly different region of
; the detector than was input.
;
; PROCEDURE:
;
; Calls destretch. If a subsection of the detector is passed as input,
; the array is embedded in a full detector array and that is processed
; with destretch. Information about the processing is noted in the
; index structures using the OPTDIST and HISTORY tag names.
;
; EXAMPLE:
;
; IDL> file = "sum_960513_085055.fits"
; IDL> rd_sumer,file,index,data,data_cols=[1]
; IDL> flatfield,index,data
; IDL> sumer_destretch,index,data
;
; MODIFICATION HISTORY:
;
;   HPW 03-SEP-1996: 
;
;   TAK 13-NOV-1996: Verbose keyword controls the DESTRETCH QUIET
;                    keyword.
;
;   HPW 19-NOV-1996: Added check on detector.
;
;   HPW 10-OCT-1997: Fixed bug related to sgt_range.
;
;   HPW 18-DEC-1997: Modified the routine to handle files created from
;                    both real time data and cd-rom data.
;
;   DMZ 29-DEC-1998: renamed to SUMER_DESTRETCH
;-

pro sumer_destretch,index,data,verbose=verbose,show_exps=show_exps

;; check for binning
bin = sgt_bin(index)
mm  = where(bin ne 1 and bin ne 0,count)
if (count gt 0) then begin
  message,'Binned data should not be destretched.',/INFORMATIONAL
endif

;; check detector
detector = sgt_detector(index,/text)
if (detector ne 'A' and detector ne 'B') then begin
  message,'Data is from '+detector,/informational
  message,'Only data from detector A or B may processed with DESTRETCH.',/info
  return
endif

;; parameters
dims   = sgt_dims(index)
slit   = sgt_slit(index)
refpix = sgt_refpix(index)
slit   = sgt_slit(index)

n_col = n_tags(data)
for ic = 0,n_col-1 do begin

  n_exp = n_elements(data.(ic)(0,0,*))
  for ie=0,n_exp-1 do begin

    ;; compute optical distortion correction - check orientation of data
    time_1   = systime(1)
    data_in  = data.(ic)(*,*,ie)
    len      = strlen(index.gen.filename)
    data_ext = strupcase(strmid(index.gen.filename,len-3,3))
    if data_ext eq "FTS" then data_in = reverse(data_in,2)
    data_out = destr_bilin2(data_in,slit,refpix(0,ic),detector,/QUIET)
    if data_ext eq "FTS" then data_out = reverse(data_out,2)
    method   = 'destr_bilin'

    if keyword_set(VERBOSE) then begin
      time = systime(1) - time_1
      mm   = fix(time/60)
      ss   = round(time mod 60)
      print,' Exposure ',string(format='(i4)',ie+1),$
          ' of ',string(format='(i4)',n_exp)+' processed in '+$
          trim(mm,'(i2.2)')+":"+trim(ss,'(i2.2)')
    endif
    
    ;; overwrite old uncorrected data
    wr = (sgt_range(index,/spectral))(*,ic)
    yr = (sgt_range(index,/spatial))(*,ic)
    data.(ic)(*,*,ie) = data_out ;; (wr(0):wr(1),yr(0):yr(1))

    ;; note processing in the index structure
    history = index.(ic).spectrum(ie).history
    if history eq '' then sep = '' else sep = ','
    index.(ic).spectrum(ie).history = history + sep + method
    index.(ic).spectrum(ie).optdist = method

    if keyword_set(SHOW_EXPS) then begin
      print,wr,yr
      nx = n_elements(data_in(*,0))
      ny = n_elements(data_in(0,*))
      if ic eq 0 and ie eq 0 then window,xs=nx,ys=ny*3,/free
      loadct,3
      img1 = data_in
      img2 = data_out ;; (wr(0):wr(1),yr(0):yr(1))
      diff = img1 - img2
      tvscl,alog10(img1 > 1),0
      tvscl,alog10(img2 > 1),1
      tvscl,diff,2
    endif

  endfor
  
endfor

return
end