;+
; NAME: IVM_PARASITIC
;
; PURPOSE: remove the scattered ``parasitic'' light from IVM images,
;	i.e. that light leaked in from adjacent Fabry-Perot orders.
;	Removes continuum from data.
;	Also removes ghost images from data.
;
; CATEGORY: data reduction
;
; CALLING SEQUENCE:
;	ivm_parasitic, oblk, fileoutput, vb, geomin=geomin, datain=datain,
;		dataout=dataout, geomf=geomf, dataf=dataf, datoutf=datoutf, $
;		     flat=flat, header=header, gdark=gdark, ddark=ddark
;
; INPUTS:
;	oblk = data file header structure
;	fileoutput = 1/0 flag according to whether file output should be
;			written/not to file specified by 'datoutf' keyword.
;	vb = 1/0 flag according to whether or not to be verbose
;
; KEYWORD PARAMETERS:
;
;    Array I/O:
;	geomin = input geometry array (ncol,nrow,nmod), processed through warp
;	datain = input data array, processed through register
;	dataout = output array for corrected data
;	geomout = output array for corrected geometry
;
;    File I/O:
;	geomf = input geometry file name, processed through warp
;	dataf = input data file name, processed through register
;	datoutf = output filename for corrected data (used if fileoutput=1)
;	geooutf = output filename for corrected geometry (used if fileoutput=1)
;
;    Optional:
;	mask = If set, use mask
;	flat = If set, these are flats,  have no modulation.
;	header = If set,  skip this many bytes when reading
;		for raw flat data.
;	ddark,gdark = Dark images for dealing with raw
;		 (un-biascorrected) images, namely flats.
;
;
; COMMON BLOCKS: none.
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
; EXAMPLE:
;    File input: 	['datoutf' file is written only if 'keepint' is set]
; 	IVM_PARASITIC, oblk, keepint, geomf=geomp, dataf=datp, datoutf=datrp
;    Array input:
;	IVM_PARASITIC, oblk, keepint, geomin=igeomc, datain=idatc, $
;				dataout=idatrp, datoutf=datrp
; MODIFICATION HISTORY:
;	DLM ??? written
;	KDL 10/96 added more commentary, made it general for use with
;		flat-field data as well.
;	BJL 96.12.11 Added ghost image correction
;	BJL April 6, 2000 Update call to IVM_GHOST
;	KDL 17July2006 added check to return fraction =0.0 if not Fe data.
;	  05/2010 NB this is not the case any more, non-Fe (e.g. Na) data use fraction of 0.04.
;	ELW 20091202 Added array I/O; less-excessive file I/O is optional.
;	KDL 05/2010: realized that because geom-camera's parasitic light was never getting 
;		removed, it still had a spectral variation to it, and
; 		later that just got folded back into things like deblur.  This is, I think,
;		the underlying cause of the funny dots from the deblur, believe it or not.
;	ELW 09/2012: Added mask
;-

PRO ivm_parasitic, oblk, fileinput, fileoutput, vb, $
		datain=datain, dataout=dataout, dataf=dataf, datoutf=datoutf, $
		geomin=geomin, geomout=geomout, geomf=geomf, geooutf=geooutf, $
                flat=flat, header=header, gdark=gdark, ddark=ddark, mask=mask

if (vb ge 1) then print,'IVM_PARASITIC start:  ',systime()
if (vb ge 2) then print, '  ** Correcting Spectral Parasitic Light **'

; Parasitic light fractions
fraction = IVM_PARA_FRAC(oblk)
if (vb ge 1) then print,'    Parasitic light fraction:  ', fraction

; Data dimensions from header
ncol = long(oblk.spix)
nrow = long(oblk.ppix)
nmod = long(oblk.nframes)
nwavel = long(oblk.n_image_sets)
nrepeat = long(oblk.ntimes)
nimage = nwavel * nrepeat

; Left off here ELW20120921 - coding mask input
; if (keyword_set(mask) &&

IF( KEYWORD_SET(flat) ) THEN BEGIN
	nmod = 1 & nrepeat = 1
ENDIF

IF( KEYWORD_SET(header) ) THEN head=header ELSE head=0
IF( ~KEYWORD_SET(gdark) ) THEN gdark = 0.
IF( N_ELEMENTS(gdark) GT 1 ) THEN gdk = REBIN(gdark,1,1) ELSE gdk = gdark
IF( ~KEYWORD_SET(ddark) ) THEN ddark = 0.
IF( N_ELEMENTS(ddark) GT 1 ) THEN ddk = REBIN(ddark,1,1) ELSE ddk = ddark

; Open input for intensity scale factor - read from file if fileinput set
if (fileinput && keyword_set(geomf)) then begin
    ; aa = ASSOC( lun1, FLTARR(ncol,nrow,nmod), head )
    geomin=FLTARR(ncol,nrow,nmod, nimage)
    OPENR, lun1, geomf, /GET_LUN
    if (head gt 0) then POINT_LUN, lun1, head
    readu, lun1, geomin
    CLOSE, lun1 & FREE_LUN, lun1
endif else if ~keyword_set(geomin) then begin
    print, 'ivm_parasitic: must set geomf keyword.'
    return
endif
if (fileinput && keyword_set(dataf)) then begin
    datain = fltarr(ncol,nrow, nmod, nimage)
    openr, lun2, dataf, /GET_LUN
    if (head gt 0) then POINT_LUN, lun1, head
    readu, lun2, datain
    close, lun2 & free_lun, lun2
endif else if ~keyword_set(datain) then begin
    print, 'ivm_parasitic: must set either datain keyword.'
    return
endif

; Scale geom images so that signal level is the same as data images
;  at continuum wavelength.
gc0 = 0.
dc0 = 0.
FOR j=1,nrepeat DO BEGIN
    ; Was:   b = aa(j*nwavel - 1)
;   b = geomin(*,*,*, j*nwavel - 1)
    b = geomin(64:447,64:447,*,j*nwavel - 1) ; KDL 20120920 should use mask, trying this as proxy for testing
    gc0 = gc0 + TOTAL(b) - gdk
    ; Was:   b = aaa(j*nwavel - 1)
;   b = datain(*,*,*, j*nwavel - 1)
    b = datain(64:447,64:447,*, j*nwavel - 1)
    dc0 = dc0 + TOTAL(b) - ddk
ENDFOR
scalf = dc0 / gc0
scalf = scalf(0)

; Compute correction image for each wavelength step, from mean
;	geom-camera image at that step. Rescale, then mult by fraction.
;	Subtract the correction from the data-cam set at that step.

dataout = FLTARR(ncol,nrow,nmod,nimage)
geomout = FLTARR(ncol,nrow,nmod,nimage)

geomr = geomin - gdark				;Geometry, corrected for dark
datar = datain - ddark				;Data-camera, corrected for dark

; 05/2010 KDL:  compute variation of geom camera here w/r/t the mean level, remove it below.
; 09/2012 KDL: below is probably safe w/r/t computing geomavg, but for limb data maybe
; not.  Best would be to use the (geometry??) mask for this...but it's *probably* ok...
geomavg=reform(rebin(geomr[ncol/4:3*ncol/4,nrow/4:3*nrow/4,*,*],1,1,nmod,nimage))
avgeomavg=rebin(geomavg,1)

FOR j=0,nimage-1 DO BEGIN
    FOR i=0,nmod-1 DO BEGIN
	; Correct data-camera for parasitic
	datac = datar[*,*,i,j] - geomr[*,*,i,j] * fraction * scalf
	; Correct for ghosts at the same time
	dataout(*,*,i,j) = IVM_GHOST(datac, oblk)
	; Correct geom-camera for parasitic 
        nfact=geomavg[i,j]/avgeomavg(0) 
        geomout[*,*,i,j]=GEOMIN[*,*,i,j]/nfact
    ENDFOR
ENDFOR

; Open output file 'datoutf' and write if set, else set 'dataout' array only
if fileoutput && keyword_set(datoutf) then begin
    OPENW, lun4, datoutf, /GET_LUN
    writeu, lun4, dataout
    CLOSE, lun4 & FREE_LUN, lun4
endif

; Open output file 'geooutf' and write if set, else set 'geomoutf' array only
if fileoutput && keyword_set(geooutf) then begin
    OPENW, lun4, geooutf, /GET_LUN
    writeu, lun4, geomout
    CLOSE, lun4 & FREE_LUN, lun4
endif

if (vb ge 2) then print,'  ivm_parasitic end:  ',systime()
END