;+
; PRO ivm_flatdata, infile, darkfile, flatfile, outfile, outim, hstr, fileout
;
; Dark subtracts and flat field divides images in an IVM file,
; interpolating in wavelength and time as needed.
;
; INPUT PARAMETERS:
;	  infile = Filename of image file to flat, string.  May be Data
;			or Geometry camera, magnetogram or calibration.
;	darkfile = Filename of appropriate dark current file, string.
;	flatfile = Filename of daily flat interpolation file, string,
;			output from IVM_DAY_FLAT, format IFIn_ccyymmdd.
;	 outfile = Name of file to write flatted data, if "fileout" is set.
;	    hstr = OBLK header structure for datafile.
;	 fileout = 1/0 according to whether or not to store in 'outfile' file
;
; OUTPUT PARAMETERS:
;	   outim = Output data, of the form  FLTARR(nx,ny,nmod*gnw*gnrepeat)
;	 maskout = Output mask, so as to steer clear of flat's puka
;
; HISTORY:
;	Replaces FLDAT.
;	Written  June 20, 2001  Barry LaBonte
;	Use /DOUBLE keyword in MIN_CURVE_SURF  May 12, 2003  BJL
;	Handle case of 1 flat  July 8, 2003 BJL
;	File output optional; added outim/fileout params - 20091120 ELW
;	Mask computed based on flat, added to parameters - 20120809 ELW
;-

PRO ivm_flatdata, infile, darkfile, flatfile, outfile, outim, hstr, fileout, maskout

; Set up
gnx = hstr.spix
gny = hstr.ppix
gnw = hstr.n_image_sets
gnmod = hstr.nframes
gnrepeat = hstr.ntimes
nimages=gnmod*gnw*gnrepeat

; Read in the dark value
IVM_GT_DARK, darkfile, dadark, gedark
; Which one to use?
IF( STRPOS( infile, 'IMG0' ) GE 0 ) THEN dark = dadark ELSE dark = gedark

; Read in flat information.
IVM_RD_DAYFLAT, flatfile, ngood, nw, nflats, nx, ny, coef, wavel, $
	days, secs, pca, fday

; Hours since 16 UT.
hrs = FLOAT(secs)/3600. + (days - days(0))*24. - 16.
inhrs = REBIN( TRANSPOSE(hrs), nw, nflats )

; Get the "location" of this gram
gw = hstr.fp_set( 0:gnw-1 )
tstr = ANYTIM2INTS( FILE2TIME( hstr.filename, /YOHKOH ) )
ghrs = FLOAT(tstr.time)/3600000. + (tstr.day - tstr(0).day)*24. - 16.
gx = gw
gy = REPLICATE( ghrs, gnw )

; If there is only 1 flat, just take the coefficients as is
IF( nflats EQ 1 ) THEN BEGIN
    newcoef = reform( coef )
ENDIF ELSE BEGIN
    ; Interpolate principal components to the (wavelength,time) of this gram
    newcoef = FLTARR(gnw, ngood)
    ;if (wavel(29,0) lt wavel(28,0)) then stop
    FOR i=0,ngood-1 DO BEGIN
        if nflats GT 1 then begin
	    newcoef(*,i) = MIN_CURVE_SURF( coef(*,*,i), wavel, inhrs, $
					   XPOUT=gx, YPOUT=gy,/DOUBLE )
         endif else begin   ; TRM 2004-03-11 ... no time interp
	    ss = sort(wavel)
	    ; newcoef(*,i) = interpol((reform(coef(*,0,i)))[ss], wavel[ss],gx)
	    newcoef(*,i) = interpol((coef(*,0,i))[ss], wavel[ss],gx)
         endelse
    ENDFOR
ENDELSE

; Define I/O, flat vars
inim  = INTARR(gnx, gny, nimages)
outim = FLTARR(gnx, gny, nimages)
flat  = FLTARR(nx, ny)
nhb = 2176		; Header bytes to skip over

; Read input file
OPENR, luni, infile, /GET_LUN
POINT_LUN, luni, nhb
READU, luni, inim
CLOSE, luni & FREE_LUN, luni

;; Walk through the spectral repeats
FOR k=0,gnrepeat-1 DO BEGIN

    ; Walk through wavelengths
    FOR n=0,gnw-1 DO BEGIN

	; Compute start/end index into array, for chunk
	istart = gnmod*(n + k*gnw)  &  iend = gnmod*(n + k*gnw + 1) - 1

	; Dark subtract
	outemp = float(inim(*,*,istart:iend)) - dark

	; Build flat
	replicate_inplace, flat, 0.	; Reinitialize 'flat' to zero
	FOR j=0,ngood-1 DO flat = flat + pca(*,*,j) * newcoef(n,j)
	flat = fday * (flat + 1.)

	; Apply flat
	FOR i=0,gnmod-1 DO outim(*,*,istart+i) = outemp(*,*,i) * flat

    ENDFOR
ENDFOR

; Compute mask based on fday
ruffmask=bytarr(nx,ny)				; Define rough mask:
ruffmask[nx/16:15*nx/16,ny/16:15*ny/16]=1	;   Mask edges only
ruffmask[3*nx/10:7*nx/10, 0:ny/6]=0		;   Mask lower puka area
ruffmask[2*nx/5:3*nx/5, 7*ny/8:15*ny/16]=0	;   Mask upper puka area
maskout = bytarr(nx,ny)				; Define mask
max_cntr = max(fday[where(ruffmask eq 1)])	; Get max of center portion
notmsk = bytscl(fday, max_cntr, max(fday))	; Scale to define pukas
maskout[where(notmsk eq 0)] = 1			; Set mask based on above
ss=intarr(3,3)+1
maskout = erode(maskout,ss)			; Erode, and we're done

; Output data if 'fileout' specified
if (fileout) then begin
    ; Write result
    OPENW, luno, outfile, /GET_LUN
    WRITEU, luno, outim
    CLOSE, luno & FREE_LUN, luno
endif

RETURN
END