;+
; NAME:
;
; sumer_flatfield
;
; PURPOSE:
;
; To remove non-uniform instrument response from SOHO/SUMER
; observations.
;
; The environment variable SUMER_FLATFIELD must be set or the flat
; field file must be specified explicitly.
;
; CATEGORY:
;
; Data handling.
;
; CALLING SEQUENCE:
;
; sumer_flatfield,index,data
; sumer_flatfield,index,data,flatfile_file=flatfield_file
; sumer_flatfield,index,data,/reverse_ff
; 
; INPUTS:
;
; The index and data structures from rd_sumer.
;
; OPTIONAL INPUTS:
;
; flatfield_file: The full pathname of a file containing a flat field
;                 image. The default is to use the flat field image
;                 taken closest to, but before, the observation. 
;       
; KEYWORD PARAMETERS:
;
; reverse_ff: Reverses the flat fielding.
;
; search: Search for flatfield images. The search is usally done only
;         during initialization.
;
; OUTPUTS:
;
; The input data is overwritten with the flat fielded data.
;
; OPTIONAL OUTPUTS:
;
; None.
;
; COMMON BLOCKS:
;
; flatfield_com: stores the list of available flatfield images. The
;                list can be updated by using the search keyword.
;
; SIDE EFFECTS:
;
; Warns the user if no flat field images where taken previous to the
; observation being flat fielded.
;
; RESTRICTIONS:
;
; Can not be used on data subjected to compression schemes 6 and
; above, data that has been binned, or data that has already been
; flatfielded.
;
; PROCEDURE:
;
; Extracts the appropriate subset of the flat field image and
; convolves it with the input data. If a flat field image file is not
; specified, the flat field image taken closest to the observation is
; used.
;
; EXAMPLE:
;
; IDL> filename = '$SUMER_DATA/sum_960513_085055.fits'
; IDL> rd_sumer,index,data
; IDL> sumer_flatfield,index,data
;
; MODIFICATION HISTORY:
;
;   HPW 24-JUL-1996: Based on Terry Kucera's flatfield routine written
;                    25-JUN-1996.
;
;   HPW 03-SEP-1996: Modified the routine to update the history field
;                    in the index structure when flat fielding is
;                    applied.
;
;   HPW 27-NOV-1996: Added the flatfield common block which stores the
;                    list of available flat field image files.
;
;   HPW 24-JUN-1997: Modified to use the structure returned by the
;                    routine mk_ffdb. Added verbose keyword. Modified
;                    the history variable to include the name of the
;                    flat field data used.
;
;   HPW 18-AUG-1997: Fixed bug related to specifying the flat field
;                    file name.
;
;   HPW 18-DEC-1997: Modified the routine to handle files created from
;                    both real time data and cd-rom data.
;
;   HPW 29-DEC-1997: Renamed keyword "reverse" "reverse_ff."
;
;   DMZ 28-Dec-1998: Renamed to SUMER_FLATFIELD, and improved error checking
;
;   TAK 12-Aug-1999: Data changed to floating point to prevent round 
;                    off errors.
;
;   HPW 16-Aug-1999: Changed time format used to determine closest
;                    flat field image.

pro sumer_flatfield,index,data,reverse_ff=reverse_ff,$
                  flatfield_file=flatfield_file,$
                  search=search,verbose=verbose,err=err

;; return to main level on error
on_error,1

TRUE  = 1b
FALSE = 0b
err=''

common flatfield_com,ffdb

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

;; check compression scheme
comp = sgt_compress(index)
if (comp ge 6) then begin
    err='Inappropriate data compression scheme for flat fielding.
    message,err,/cont
    return
endif

;; check to see if data has been binned
bin = sgt_bin(index)
mm  = where(bin ne 1 AND bin ne 0,count)
if (count gt 0) then begin
    message,'Binned data cannot be flat fielded.',/cont
    return
endif
;; check to see if data has been flat fielded already
if keyword_set(reverse_ff) then reverse_ff = TRUE else reverse_ff = FALSE
n_ff   = 0  		;; number of exposures to flat field
n_cols = n_tags(data)	;; number of columns of data
for i=0,n_cols-1 do begin
  match = where(index.(i).spectrum.flatfield eq 0,count)
  n_ff  = n_ff + count
endfor

if ((n_ff eq 0) and not(reverse_ff)) then begin
  message,'Data has already been flat fielded.',/CONTINUE
  return
end
  
;; ----------------------------------------------------------------------
;; Read flatfield image
;; ----------------------------------------------------------------------

if not(keyword_set(flatfield_file)) then begin

  ;; construct flat field image database

  if (datatype(ffdb) ne 'STC') or keyword_set(SEARCH) then begin
   ffdb = mk_sumer_dbase_ff(verbose=verbose,err=err)
   if err ne '' then begin
    message,err,/cont
    return
   endif
  endif      

  ;; find closest flat field image
  this_detector = sgt_detector(index,/TEXT)
  this_time     = anytim(sgt_time(index)) 
  match = where(this_detector eq ffdb.detector,count)
  if count eq 0 then begin
     err="No flat field files taken with detector "+this_detector
     message,err,/cont
     return
  endif
  diff = min(abs(this_time-ffdb(match).utime),pos_min)
  flatfield_file = ffdb(match(pos_min)).filename

  fn  = str_replace(flatfield_file,"\","/")
  fn  = str_sep(fn,"/")
  fn  = fn(n_elements(fn)-1)
  fmt = '(2a25,a15)'
  if keyword_set(VERBOSE) then begin
    print,"Using the following flat field data:"
    print,format=fmt,fn,ffdb(match(pos_min)).time,ffdb(match(pos_min)).detector
  endif

  ;; read file
  rd_sumer,flatfield_file,index_ff,data_ff,data_cols=[1]

endif else begin

  flatfield_files = findfile(flatfield_file)
  exist           = file_exist(flatfield_file)

  if (exist) then begin
    rd_sumer,flatfield_files(0),index_ff,data_ff,data_cols=[1]
  endif else begin
    err='Specified flat field image not found.'
    message,err,/cont
    return
  endelse

  fn  = str_sep(flatfield_file,"/")
  fn  = fn(n_elements(fn)-1)
  fmt = '(2a25,a15)'
  tim = fmt_tim(filename2time(flatfield_files(0)))
  det = sgt_detector(index_ff,/text)
  if keyword_set(VERBOSE) then begin
    print,"Using the following flat field data:"
    print,format=fmt,fn,tim,det
  endif

endelse

;; ----------------------------------------------------------------------
;; CHECK ORIENTATION OF THE EXPOSURE
;;
;; Real time data (FITS) has north on bottom while cd-rom data (FTS)
;; has north on top. If the data file and the flat field file are not
;; from the same source, the flat field file must be reversed for the
;; orientations to agree. Note that the data source is determined from
;; the file name written in the index.
;; ----------------------------------------------------------------------

len      = strlen(index_ff.gen.filename)
ff_ext   = strupcase(strmid(index_ff.gen.filename,len-3,3))
len      = strlen(index.gen.filename)
data_ext = strupcase(strmid(index.gen.filename,len-3,3))

if ff_ext ne data_ext then data_ff.(0) = reverse(data_ff.(0),2)

;; ----------------------------------------------------------------------
;; Convolve flatfield image with data
;; ----------------------------------------------------------------------

for i=0,n_cols-1 do begin

  dsize = size(data.(i))
  if dsize(0) eq 2 then n_exp = 1 else n_exp = dsize(3)
  dy = dsize(2) & dw = dsize(1)

  ;; get reference pixels
  refpix   = sgt_refpix(index.(i))
  w_refpix = refpix(0)			
  y_refpix = refpix(1)

  ;; wavelength direction range
  if (dw eq 1024) then begin
    w_range = [0,1023]
  endif else begin
    if (dw mod 2 eq 0) then begin
      w_range = w_refpix - 1 + [-dw/2+1,dw/2]	
    endif else begin
      w_range = w_refpix - 1 + [-dw/2,dw/2]		
    endelse  
  endelse

  ;; spatial direction range
  if (dy eq 360) then begin
    y_range = [0,359]
  endif else begin
    y_range = y_refpix + [-dy/2+1,dy/2]		
  endelse

  ;; extract subset of flatfield image

  ff = float(data_ff.(0)(w_range(0):w_range(1),y_range(0):y_range(1)))

  for j=0,n_exp-1 do begin
    flatfield = index.(i).spectrum(j).flatfield
    history   = index.(i).spectrum(j).history
    if history eq '' then sep = '' else sep = ','
    if ((flatfield) and (reverse_ff)) then begin
      ;; reverse flat fielding of data
      data.(i)(*,*,j) = float(data.(i)(*,*,j)*ff)
      ;; set flatfield "off" in index
      index.(i).spectrum(j).flatfield = FALSE
      index.(i).spectrum(j).history   = history+sep+'flatfield-reverse:'+fn
    endif
    if (not(flatfield) and not(reverse_ff)) then begin
      ;; flat field data
      data.(i)(*,*,j) = float(data.(i)(*,*,j)/ff)
      ;; set flatfield "on" in index
      index.(i).spectrum(j).flatfield = TRUE
      index.(i).spectrum(j).history   = history+sep+'flatfield:'+fn
    endif
  endfor

endfor

return
end