;+
; PRO IVM_SPUR_MOD, oblk,fileout, vb, geomr=geomr,geomc=geomc,datr=datr,datc=datc,$
;		fgeomr=fgeomr, fgeomc=fgeomc, fdatr=fdatr, fdatc=fdatc
;
; DESCRIPTION:  Computes spurious modulation vector of geometry images and
;	corrects geometry images.  Corrects both geometry and data images for
;	extinction.
;
; INPUT PARAMETERS:
;	oblk = IVM data file header block structure
;	filein  = flag for whether or not to read from input file
;	fileout = flag for whether or not to write output file
;	vb = verbose flag (0 is errors-only, 1 is minimal)
;    Input arrays can be either 3D or 4D:
;	geomr = Input warped geometry images array  FLTARR(ncol,nrow, nimages)
;	datr  = Input ___data images array   	    FLTARR(ncol,nrow, nimages)
;    If none of the array keywords above are set, file input is used:
;	fgeomr = Filename of warped geometry images for input.
;	fdatr = filename of ___data images for input.
;
; OUTPUT PARAMETERS: 	[4D Output arrays: fltarr(ncol,nrow,nmod,nwave*nrepeat)]
;	datc   = Array of extinction corrected data images for output.
;	geomc  = Array of intensity corrected geometry images for output.
;	fdatc  = filename of extinction corrected data images for output.
;	fgeomc = Filename of intensity corrected geometry images for output.
;
; EXAMPLES: 	                         IN             OUT        IN        OUT
;   Array input:
;	IVM_SPUR_MOD, oblk,0,keepint, geomr=georb, geomc=geomp, datr=datr, datc=datp
;   File input:
;	IVM_SPUR_MOD, oblk,0,keepint, fgeomr=georb, fgeomc=geomp, fdatr=datr, fdatc=datp
;   Rivm call (file I/O if keepint [fileout] is set):
;	IVM_SPUR_MOD, oblk,0,keepint, geomr=igeomr, geomc=igeomc, datr=datreg,
;				datc=idatc, fgeomc=geomp, fdatc=datp
;
; METHOD:
;	There is a small leakage of the polarization modulation and 
;	spectral line shape into the Geometry images.
;	The Geometry passband includes the etalon order which
;	is on the line along with about 30 others.  The parabolic
;	wavelength shift is seen, but the "line" is not at the
;	same wavelength as the Data line, so this is not simply
;	parasitic light.
;
;	This procedure assembles a smooth function for 
;	removing the spurious modulation from the geometry 
;	images.  The remaining geometry image intensity variation
;	is caused by variable extinction, and is removed from
;	both the geometry and data images.
;
;	See TechNotes/IVM_note_20010614.html for more information.
;
; HISTORY:
;	Written  April 3, 1997  Barry LaBonte
;	Added extinction correction  May 29, 1997   BJL
;       Added array I/O, made file I/O optional + reduced - November 23, 2009 ELW
;       Logical operators used (~ || && replace NOT OR AND) - January 6, 2010 ELW
;       File input - May 27, 2010 ELW
;       Got rid of parens in arrays - Oct 15, 2010 ELW
;-

PRO IVM_SPUR_MOD, oblk,filein, fileout,vb,geomr=geomr,geomc=geomc,datr=datr,datc=datc,$
		fgeomr=fgeomr, fgeomc=fgeomc, fdatr=fdatr, fdatc=fdatc

compile_opt idl2
if (vb ge 1) then print,'IVM_SPUR_MOD start: ',systime()
if (vb ge 2) then print,' ** Removing Spurious Extinction and Modulation **'

; Get the array dimensions
ncol    = oblk.spix
nrow    = oblk.ppix
nmod    = oblk.nframes
nwave   = oblk.n_image_sets
nrepeat = oblk.ntimes
nwr     = nwave*nrepeat
nimages = nmod * nwr

; The model is:
;	nmod transmissions, spatially and spectrally invariant 
;	nrepeat-1 order polynomial time variation of extinction,
;		spatially and spectrally invariant
;	nwave spectral intensities, spatially invariant.
;
; NOTE: This should be improved to account for the spatial variation
;	of the spectral variation.

; First pass  --------------------------------------------------------
; Use medians to avoid bad images from clouds, etc.

; Geom input ('fgeomr' or 'geomr')
if (filein) then begin
    if (~keyword_set(fgeomr)) then begin
	print, 'ivm_spur_mod: ERROR! Must define fgeomr input file. Stopping'
	stop
    endif
    geomr = FLTARR(ncol,nrow, nimages)
    OPENR, lun1, fgeomr, /GET_LUN
    readu, lun1, geomr
    CLOSE, lun1 & FREE_LUN, lun1
endif else begin
    r=data_4Dto3D(geomr, ncol, nrow, nmod, nwr) 	; Assure geomr is 3D
    if (r ne 1) then begin
        print, 'ivm_spur_mod:  unable to reform input data.'
        return
    endif
endelse

intenr = FLTARR( nimages )

FOR i=0,nimages-1 DO BEGIN
	intenr[i] =  MEDIAN( geomr[*,*,i] )
ENDFOR
intenr = REFORM( intenr, nmod, nwave, nrepeat, /overwrite)

; Average over modulation and wavelength
tvect = FLTARR(nrepeat)
FOR i=0,nrepeat-1 DO tvect[i] = MEDIAN( intenr[*,*,i] )

; Fit for variation with "time" (actually sequence index)
; Must remove here so it will remain in the Geometry for
; correcting the Data.
x = ( FINDGEN(nrepeat) + 0.5 ) * nmod * nwave - 0.5
res = POLY_FIT( x, tvect, nrepeat-1 )
xt = FINDGEN(nimages)
timevar = REPLICATE( res[0,0], nimages)
FOR i=1,nrepeat-1 DO timevar = timevar + res[0,i] * xt^i
tnorm = REBIN( timevar, 1 )
timevar = FLOAT( timevar / tnorm[0] )
inten = intenr / timevar

; Average over wavelength and repeats
mvect = FLTARR(nmod)
FOR i=0,nmod-1 DO mvect[i] = MEDIAN( inten[i,*,*] )

; Average over modulation and repeats
; This should be spatially variant
wvect = FLTARR(nwave)
FOR i=0,nwave-1 DO wvect[i] = MEDIAN( inten[*,i,*] )

; Invert
cvect = REFORM( (1./mvect) # (1./wvect), nmod*nwave )

; Normalize
fnorm = REBIN( cvect, 1 )
cvect = cvect / fnorm[0]

; Now determine the extinction variation
cv = REFORM( REBIN( cvect, nmod*nwave, nrepeat ), nimages )
intenc = REFORM( intenr, nimages ) * cv
inorm = REBIN( intenc, 1 )
intenc = intenc / inorm[0]

; Second pass  --------------------------------------------------------
; ELW 20091123 - get rid of extra file read - reuse geomr from earlier
; OPENR, lun1, fgeomr, /GET_LUN
; aa = ASSOC( lun1, FLTARR(ncol,nrow) )

; Data input ('fdatr' or 'datr')
if (filein) then begin
    if (~keyword_set(fdatr)) then begin
	print, 'ivm_spur_mod: ERROR! Must define fdatr input file. Stopping'
	stop
    endif
    datr = FLTARR(ncol,nrow, nimages)
    OPENR, lun3, fdatr, /GET_LUN
    READU, lun3, datr
    CLOSE, lun3 & FREE_LUN, lun3
endif else begin
    r=data_4Dto3D(datr,  ncol, nrow, nmod, nwr)    ; Make sure datr is 3D
    if (r ne 1) then begin
        print, 'ivm_spur_mod:  unable to reform input data.'
        return
    endif
endelse

;==============================================
; Compute output arrays, make all I/O arrays 4D
;==============================================
geomc = FLTARR(ncol,nrow, nimages)
datc  = FLTARR(ncol,nrow, nimages)

FOR i=0,nimages-1 DO BEGIN
    geomc[*,*,i] = geomr[*,*,i] * cv[i] / intenc[i]
    datc[*,*,i]  = datr[*,*,i] / intenc[i]
ENDFOR

if keyword_set(datc)  then r=data_3Dto4D(datc,  ncol, nrow, nmod, nwr)
if keyword_set(geomc) then r=data_3Dto4D(geomc, ncol, nrow, nmod, nwr)
if keyword_set(datr)  then r=data_3Dto4D(datr,  ncol, nrow, nmod, nwr)
if keyword_set(geomr) then r=data_3Dto4D(geomr, ncol, nrow, nmod, nwr)

;==============================================
; Write output files, if 'fileout' is set
;==============================================
if (fileout) then begin
    OPENW,  lun2, fgeomc, /GET_LUN
    WRITEU, lun2, geomc
    CLOSE, lun2 & FREE_LUN, lun2

    OPENW,  lun4, fdatc, /GET_LUN
    WRITEU, lun4, datc
    CLOSE, lun4 & FREE_LUN, lun4
endif

if (vb ge 2) then begin
    ; Save answers
    fnorm = REBIN( mvect, 1 )
    PRINT, ' Geometry spurious modulation '
    PRINT, ' Modulation transmissions: ', mvect/fnorm[0]
    fnorm = REBIN( wvect, 1 )
    PRINT, ' Wavelength transmissions: ', wvect/fnorm[0]
    PRINT, ' Extinction transmissions: ', intenc
endif
if (vb ge 2) then print,'IVM_SPUR_MOD end: ',systime()

END