;+
; NAME: ivm_calstk -- Handle wavelength variation of IVM calibration
;
; INPUT PARAMETERS:
;	oblk = raw oblk data header structure from IMGR* file
;	filein  = Use file input if set, else array
;	fileout = Use file output if set, else array
;	vb   = verbosity [0/unset = none, 1 = start/end times]
;    Array:
;	iss = Instrumental Stokes spectra - input array (nx,ny,4,nw) 
;	caldbarr = name of cal variation array (used if caldbfile not set)
;    File:
;	sfile = name of instrumental Stokes spectra, *_st file
;	cfile = name in instrument matrix file, *_instr file
;	caldbfile = name of cal variation database file, ivm_calvar_db.* file
;
; OUTPUT PARAMETERS:
;	css = Calibrated Stokes spectra - output array (nx,ny,4,nw)
;
; OPTIONAL KEYWORDS:
;	bigmat = If set, use this as the continuum instrument matrix.
;		Default reads it from cfile.
;	outf = If set, write 'css' array to the specified file.
;	signflip = If set, flip V sign
;
; HISTORY:
;	Written  October 24, 2000  Barry LaBonte
;	Bug fix for zero division  July 11, 2003  BJL
;	Added array I/O, RAM-centric operation  December 23, 2009  ELW
;	Booleanized (NOT=~ AND=&&)  January 8, 2010 ELW
;	File input parameter  June 1, 2010 ELW
;	caldbarr input parameter added  August 31, 2010 ELW
;-

PRO ivm_calstk, oblk, filein,fileout, vb, bigmat, css, $
		iss=iss, outf=outf, sfile=sfile, cfile=cfile,  $
		caldbfile=caldbfile, caldbarr=caldbarr, signflip=signflip

    if (vb ge 1) then PRINT, 'IVM_CALSTK start: ' + SYSTIME()

    ; Get data dimensions from oblk header
    nx = oblk.spix
    ny = oblk.ppix
    nw = oblk.n_image_sets	; Number of wavelengths

    ; Get the wavelength scale, milliAngstroms
    w = oblk.fp_set[0:nw-1]
    w = w - MIN(w)
    fpma = IVM_FP_STEP( oblk.f_wavelength * 10000. )
    w = temporary(w) * fpma
    wix = SORT(w)
    dw = DERIV( temporary(w[wix]) )
    zero = WHERE( dw EQ 0., nzero )
    IF( nzero GT 0 ) THEN dw(zero) = MEDIAN(dw)
    bix = SORT( wix )

    ; Read in the continuum instrument matrix if 'cfile' set
    IF( KEYWORD_SET(cfile) ) THEN BEGIN
	bigmat = FLTARR( nx, ny, 4, 4 )
	OPENR, lun, cfile, /GET_LUN
	READU, lun, bigmat
	CLOSE, lun & free_lun, lun
    ENDIF

    ; If signflip flag set, postmultiply by fixmat to correct
    if (keyword_set(signflip)) then begin
	fixmat = [[1.,0.,0.,0.], $
		  [0.,1.,0.,0.], $
		  [0.,0.,1.,0.], $
		  [0.,0.,0., -1.]] 
	FOR i = 0,3 DO BEGIN		; Taken from calmat.pro
	    mcol = bigmat[*,*,*,i] 	; This is a row, despite the name
	    FOR j=0,3 DO BEGIN
		temp = FLTARR(nx,ny)
		FOR k=0,3 DO BEGIN
		    temp = temp + mcol[*,*,k] * fixmat[j,k]
		END
		bigmat[0,0,j,i] = temp
	    END
	END
	if (vb ge 1) then PRINT, 'ivm_calstk:  v-sign flip accounted for.'
    endif

    ; Read instrumental Stokes spectra from file 'sfile' if 'iss' is not set
    if (filein) then begin
	if (~keyword_set(sfile)) then begin
	    print,'ivm_calstk: ERROR! Input filename "sfile" not set. Stopping.'
	    stop
	endif
	iss = FLTARR( nx, ny, 4, nw )
	OPENR, lun, sfile, /GET_LUN
	READU, lun, iss
	CLOSE, lun & free_lun, lun
    endif

    ; Get the cal variation coefficients
    if (keyword_set(caldbfile)) then begin 		    ; Read from file
	cstr = RD_IVM_DB( caldbfile )
	coef = FLTARR(2, 4, 4 )
	coef[*,1:3,1:3] = cstr.calvar
    endif else if (n_elements(caldbarr) eq 18) then begin   ; Set from parameter
	coef = FLTARR(2, 4, 4 )
	coef[*,1:3,1:3] = caldbarr
    endif else begin
	print, 'ivm_calstk:  CALDBFILE or CALDBARR must be defined.  Stopping.'
	stop
    endelse

    ; Form the first and second derivatives of the normalized intensity
    ; w.r.t. wavelength
    d1 = FLTARR( nx, ny, nw )
    d2 = FLTARR( nx, ny, nw )
    icon = iss[*,*,0,nw-1]		; Continuum intensity
    FOR j=0,ny-1 DO BEGIN
	FOR i=0,nx-1 DO BEGIN
	    inten = REFORM( iss[i,j,0,wix] ) / ( icon[i,j]  > 1. )
	    inten = MEDIAN( inten, 3 )
	    tmp = DERIV( inten ) / dw
	    d1[i,j,bix] = tmp
	    d2[i,j,bix] = DERIV( tmp ) / dw
	ENDFOR
    ENDFOR

    ; Walk through wavelengths, generate correction to instrument matrix 
    css = FLTARR( nx, ny, 4, nw)

    FOR k=0,nw-1 DO BEGIN
	wim = bigmat
	; Only correct polarized elements
	FOR j=1,3 DO BEGIN
	    FOR i=1,3 DO BEGIN
		wim[*,*,i,j] = wim[*,*,i,j] + coef[0,i,j]*d1[*,*,k] $
					+ coef[1,i,j]*d2[*,*,k]
	    ENDFOR
	ENDFOR

	; Walk through the FOV
	FOR j=0,ny-1 DO BEGIN
	    FOR i=0,nx-1 DO BEGIN
		; Invert the instrument matrix
		wim[i,j,*,*] = INVERT( REFORM( wim[i,j,*,*] ) )
		; Calibrate the Stokes data
		stkinstr =  REFORM( iss[i,j,*,k] )
		css[i,j,0,k] = stkinstr[0]
		stkinstr[0] = 1.0
		stkcal = stkinstr # REFORM( wim[i,j,*,*] )
		css[i,j,1:3,k] = stkcal[1:3]
	    ENDFOR
	ENDFOR
    ENDFOR

    ; Write the calibrated data to 'outf' if both it, and 'fileout' are set
    if (fileout && keyword_set(outf)) then begin
	OPENW, lun, outf, /get_lun
	WRITEU, lun, css
	CLOSE, lun & FREE_LUN, lun
    endif

    if (vb ge 3) then PRINT, '  ivm_calstk end: ' + SYSTIME()
END