;+
;Project:
;       SOHO - SUMER
;Name:
;       SUM_Flatfield
;Purpose:
;      Apply flat field correction to sumer data
;Use:
;       data_out = sum_flatfield(data,slit,RefPix,FFdat)
;Inputs
;      Data - data to be flat fielded. Can be an array of images with
;	the images number as the 3rd dimension if all the images are
;	from the same spot on the detector
;      Slit - Nuber of Slit used
;      RefPix = 1023 - pixpos(header_data) + 1
;      FFdat = variable containing flat field matrix.
;Keyword input:
;      Reverse - undo Flat Field
;Keyword Output:
;      FF - the flat field
;Output
;       Data_Out - Flat Fielded data
;Calls:
;       None
;Common:
;       None
;Restrictions:
;      This program should NOT be used on data which has been binned or
;       subjected to compression schemes 6 and up.
;Written:
;       Terry Kucera,    June 25, 1996
;Modified:
;	Udo Schuehle,  May 11, 1998
;	Use of MPAE-restore files
;	Udo Schuehle, May 10, 2001
;	Accept large image formats with short slits
;       Dietmar Germerott, January 22, 2013
;         Correct handing of out of detector areas
;-
;:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

Function sum_flatfield,data,Slit,RefPix,FFdat,Reverse=Reverse,Flat=FF

;determine size of data;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

dsize=size(data)
if dsize(0) eq 2 then nexp = 1 else nexp = dsize(3)
dY =  dsize(2)
dW = dsize(1)

      data=float(data)
      data_out=fltarr(dW,dY,nexp)

;determine subframe of image;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

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

if dY GE 360 THEN SlitLoc = center	; accept large format with short slit
					;spatial image range:
YRange=[SlitLoc,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

fftmp=ffdat

;---------------------------------------------------------------------------
;  handle areas outside detector e.g. refpix 20 of 50x360 image
;---------------------------------------------------------------------------

IF wrange(0) LT 0 THEN BEGIN
  woff=abs(wrange(0))
  wrange=wrange+woff
  fftmp=fltarr(1024+woff,360)   ;increase array
  fftmp(*)=1; set ff mask to 1
  fftmp(woff:woff+1023,*)=ffdat ;and shift
ENDIF

IF wrange(1) GT 1023 THEN BEGIN
  woff=wrange(1)-1023

  fftmp=fltarr(1024+woff,360)   ;increase array
  fftmp(*)=1; set ff mask to 1
  fftmp(0:1023,*)=ffdat
ENDIF
						;select subframe of flat field
FF=FFtmp(WRange(0):WRange(1),YRange(0):YRange(1))

;perform flat field correction or undo;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

for  i=0,nexp-1 do begin
  if keyword_set(Reverse) then data_out(*,*,i) = data(*,*,i) * FF $
  else data_out(*,*,i) = data(*,*,i) / FF
  
endfor

return,data_out

end