;+
; NAME:
;       IVM_CALCAL
;
; PURPOSE:
;       Compute the instrument matrix for the IVM, from a standard
;       ICAL... calibration data file.
;
; CATEGORY:
;       Data analysis
;
; CALLING SEQUENCE:
;       ivm_calcal,calfile,d0,flatfh,bigmat,nn,npts,calm
;
; INPUTS:
;       calfile - Name of the calibration data file
;            d0 - name of the data camera dark file
;        flatfh - Name of the data camera flat file
;          calm - name of file to which output matrix is written
;            nn - dimension of images, assumed square
;          npts - number of wavelength pts in spectrum
;
; OPTIONAL INPUTS:
;       
; KEYWORD PARAMETERS:
;          scat - scattered-light parameters for ivm_scatter
;      keep_int - if set, don't delete intermediate files
;         fourf - if set, assume four-frame modulation sequence
;         csmth - width of smoothing function to apply to cal matrix
;
; OUTPUTS:
;       bigmat - inverted calibration matrix, nn x nn x 4 x 4.
;
; RESTRICTIONS:
;
; PROCEDURE:
;
; MODIFICATION HISTORY:
;       Mon Apr 15 10:01:26 1996, Don Mickey <mickey@ahinahina>
;               Added doc header. This routine mostly breaks out the
;               cal-matrix code from rivm, which is too long.
;       Thu Aug 24 13:53:58 2000, Don Mickey
;               Read oblk from raw data file; pass it to calmat.
;	October 30, 2000  Barry LaBonte
;	Oct 3, 2002  BJL  Replace SPAWN with REMOVE_FILE.
;	Dec 30, 2009 ELW  Uses modified demod + calmat -- passes array between
;-

PRO ivm_calcal,calfile,d0,flatfh,bigmat,nn,npts,calm,$
               scat=scat,keep_int=keep_int,fourf=fourf,csmth=csmth

hm=2176                         ;This will break, only on a few *old* files.
                                ;It's the size of magic header + oblk.

setenv,'SHELL=/bin/sh'		; ELW-20091230 keep remove_file sane

calf=calfile+'_df'
cals=calfile+'_df_st'
;calm=calfile+'_instr'

oblk = ivm_rd_head(calfile)
; nw   = oblk.n_image_sets
nw = 10

get_lun,lun1
get_lun,lun2

a=assoc(lun1,intarr(nn,nn),hm)	;raw data/cal files

aaa=assoc(lun1,fltarr(nn,nn))	;filtered flat file

;a) flat the cal's
openr,lun1,flatfh
c=aaa(npts-1)                   ;get last (continuum) flat, for cal data
close,lun1
openr,lun1,calfile
openw,lun2,calf
print,'Flat cal: ',calfile,' ',systime()
ncalfrm = 59
if (keyword_set(fourf)) then ncalfrm=39
;; Scattered light correction
b = a(3) - d0                   ;pick some frame, calc scat light
IF keyword_set(scat) then begin
    print,'ivm_calcal:  Using scattered-light correction.'
    IF n_elements(scat) GE 2 then $
      bb = ivm_scatter(b,cscat,lev=scat(0),pow=scat(1)) $
    ELSE bb = ivm_scatter(b,cscat,lev = scat)
    cscat = temporary(cscat)/total(b) ;scale to total
endif

for j=0,ncalfrm do begin
    b=a(j) - d0
    IF keyword_set(scat) then begin
        btot = total(b)         ; For scale.
        bs = b - cscat*btot     ; Subtract scaled offset, same shape for all.
        b = bs * btot/total(bs) ; Rescale to preserve photon count.
    endif 
    writeu,lun2,c*b
end
close,lun1&close,lun2
print,'End Flat cal: ',systime()


;b) demodulate the cal's

;ELW 20091223 - need to tidy it up and test it
FileIn=1
FileOut=1 ;for debugging purposes
; 'cals' can be eliminated [by setting fileoutput=0] if we use 'caldat' instead
;                               in    out        in          out
demod,nn,nn,nw,FileIn,FileOut,idatra,caldat,rawfile=calf,stkfile=cals,fourf=fourf
;				/no_flt, /head

idatra=0	; array from calf not needed?... for now 20091223

;if (keyword_set(fourf)) then begin
;    demod_2,nn,nn,calf,cals     ; 4 frame input -> 4 fram uncal Stokes out
;endif else begin
;    demod,nn,nn,calf,cals       ; creates 6frames -> 4 Stokes uncalibrated
;endelse
;
IF NOT KEYWORD_SET(keep_int) THEN BEGIN
    print,'Removing ',calf
;  STOP
;    REMOVE_FILE, calf
    file_delete, calf
ENDIF
;IF NOT KEYWORD_SET(keep_int) THEN REMOVE_FILE,calf

;c) create the calibration matrix:
; ELW-20091223 - calmat modified to take caldat (array input) in place of cals
;CALMAT, bigmat,cals,oblk,smooth=csmth ;creates big cal matrix

CALMAT, oblk,caldat,bigmat, SMOOTH=csmth ;creates big cal matrix
print,'  Averaged cal matrix:'
print, calmatsum(bigmat,nn,nn)	;prints averaged matrix

; Save the instrument matrix
OPENW, lun1,calm
WRITEu, lun1,bigmat
CLOSE, lun1
PRINT, 'Cal written to disk: ',systime()

; Invert matrix and print average for monitoring
tmp = REBIN( bigmat, nn/4, nn/4, 4, 4)
CINV, tmp,nn/4,nn/4
PRINT,CALMATSUM(tmp,nn/4,nn/4)

IF NOT KEYWORD_SET(keep_int) THEN BEGIN
    print,'Removing ',cals
;  STOP
;    REMOVE_FILE,cals
    file_delete, cals
ENDIF

FREE_LUN,lun1
FREE_LUN,lun2

RETURN
END