;------------------------------------------------------
function destretch,data,Slit,RefPix,Quiet=Quiet
;
;+
;Project:
;       SOHO - SUMER
;Name:
;                 destretch
;Purpose:
;      Apply distortion correction to sumer data
;Use:
;       data_out = destretch(data,slit,RefPix)
;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)   
;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.
;Written:
;       Tom Moran,    July 15, 1996
;	Changes Oct, 3 1996:
;       corrested in case where full 1024 pixels used
;       comon block added & some of file reading logic changed
;       Quiet keyword introduced
;-


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


if not keyword_set(Quiet) then begin
   print,' '
   print,' '
   print,'                        Destretch'
   print,' '
   print,' '
   print,'Written by Tom Moran, GSFC, July, 1996'
   print,' '
   print,'The procedure corrects for wavelength and spatial distortions '
   print,'in the SUMER A detector'
   print,' '
   print,'Wavelength sigma = 0.1 pixels - Maximum deviations from the mean '
   print,'are about 0.25 pixels '
   print,' '
   print,'Spatial sigma = 0.3 pixels - Maximum deviations from the mean '
   print,'are about 0.6 pixels  
   print,' '
   print,'The correction is valid for spatial pixels 20 to 340
   print,' '
   print,'More data on uncertainties will be provided as it becomes available'
   print,' '
   print,' '
   print,'To use pixel shift frames without destretching the data, readf'
   print,'in del_x_os.dat, del_x_e.dat, del_y_os.dat as integer arrays - '
   print,'intarr(1024,360).  Then convert them to float arrays and '
   print,'divide by 1000.  One can then take line center pixel positions'
   print,'found by analyzing an uncorrected frame and correct the line'
   print,'center position with the shift found in the del frames.  del_x_os'
   print,'contains shifts for the 300 arcsecond slit, and del_x_e contains'
   print,'shifts for the 120 arcsecond slit on the top and bottom of the array' 
   print,' '
endif

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

up = 300
center = 180
down = 60

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

YRange=SlitLoc+[-dY/2,dY/2-1]

;wavelength direction range

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

;YRange(0) = 0
;YRange(1) = 359
;WRange(0) = 0
;WRange(1) = 1023

; read in delta x,y arrays

  dir = 'science:[software.moran]'
  if (f_t eq 0) then Del_X_file =dir+'del_x_os.dat' $
     else  Del_X_file =dir+'del_x_e.dat' 
  Del_Y_file = dir+'del_y_os.dat'
  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,dir+'size_fac.dat'
     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(WRange(0):WRange(1),YRange(0):YRange(1)) = data

    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_o_frame = o_frame(i1:i2,j1:j2)
 
   sub_t_frame = sub_o_frame
    sub_t_frame(*,*) = 0.    

   c = ''

; contruct the new frame

      print,' Performing wavelength correction - takes about 1.5 minutes'
      print,' '

      for j=0,j2-j1-1 do begin
         old_pix = sub_x + sub_del(*,j)
         sub_o_frame(*,j) = spline(sub_x,sub_frame(*,j), old_pix,2.01)
      endfor 

     print,' Performing spatial correction - takes about 1.5 minutes'      
     print,' '
     for i=0,i2-i1-1 do begin

;      for i=500,512 do begin
;        print,'i = ', i    
;        plot,sub_dely(i,*)
;        read,' sub_dely, hit return ', c

       old_pix = sub_y + sub_dely(i,*)

;       for jj=0,359 do print,jj,' ',sub_y(jj),' ',old_pix(jj)

;       read,'hit return ', c

       sub_t_frame(i,*) = spline(sub_y,sub_o_frame(i,*), old_pix,2.01)

;       for jj=0,359 do print,jj,' ',sub_o_frame(i,*),' ',sub_t_frame(i,*)

;       read,'hit return ', c

;       window,0,xsize=500,ysize=400
;       set_xy
;       print,i
;       plot,sub_o_frame(i,*)
;       oplot,sub_t_frame(i,*)
;       read,'hit return ', c
       
     endfor

     print,'correcting for pixel size'
     print,' '

     sub_t_frame = sub_t_frame * sub_size

     print,'creating full frame with subframe enclosed' 
     print,' '
      o_frame(i1:i2,j1:j2) = sub_t_frame

      o_frame(0:4,7) = 0.
      o_frame = o_frame > 0.

      window,3,xsize=1024,ysize=360
 
      data_out = o_frame
     
      data_out(*,0:5) = 0.
      data_out(*,356:359) = 0.

      tvscl,data_out
 
return,data_out   
  
end