;+
;NAME:
;	RIVM
;PURPOSE:
;	Reads in raw dark files, flat files, cal files, and data files
;	and reduces them -- everything short of inversion.
;	This is the driver program for IVM data reduction
;	Note: Rivm is an outgrowth of ivmreduce, enough changes to
;	give it a new name.
;CATEGORY:
;	Data analysis
;CALLING SEQUENCE:
;	rval=rivm(calfile,datfile)
;INPUTS:
;	calfile = string, name of calseq raw file
;	datfile = string, name of raw datafile
;OPTIONAL INPUTS:
;KEYWORDS:
;       save_binarr = If eq 1, then write output arrays to binary files
;       REDCAL = if set, then calibration has already been reduced,etc.
;       REDFLAT = if set, then flats have already been filtered, etc.
;	NO_REG - if set, does not attempt to coalign images using
;		geometry data.
;       REGPARM = [bx,by,hx,hy,no_lim].  Optional parameters for coalignment.
;                 Usually not necessary to explicitely set.
;                 bx,by: size of coalignment window (default: 128,128).
;                 hx,hy: centerpoint of coalignment window in image frame.
;                        (default: center of image).
;		  no_lim: default=0 (5-pix shift limit).  Useful to set to 
;		  1 if guider was not working well and large shifts
;		  are needed
;       NO_LIM - Useful to set to 1 if guider was not working well
;                and large shifts are needed.  REGPARM overrides this
;                keyword. 
;       SHREFIDX - Image index for reference image in shifter_2d.
;       DETREND - if set, take limb-darkening out of correlation subimages
;                 in shifter_2d and in geometry/data coalignment procedure
;       GIF - if set, write a Bl_qnd image.
;	OLD_DATA = if set, prompts for header sizes of d0,f0,calfile,datfile
;			default is to assume the currently stable number
;	D0FIL = string, name of dark file for camera 0
;	F0FIL = string, name of flat file for camera 0
;	FileIO - If set, use file input and output
;	KEEP_INT - If set, saves intermediate files by setting "fileoutput" to
;		1 for each procedure.  Array I/O will still be used, except
;		where "FileIO" is set.  Final output is saved regardless of
;		setting.  Needs of order 850Mb free.
;	NPOINTS = number of spectral point - default=oblk->n_image_sets - 1
;	NFUS = number of fabry=perot steps between npts - default=from header
;	SHUT - if file named then will use that file to subtract the
;              shutter light leakage as stored in that file
;	PARASITIC - if set, correct for parasitic light: leakage from
;		adjacent orders of the FP into the spectrum.  Only for 
;		630.25 as of 07/17/2006
;	SPURIOUS -- if set, remove spurious extinction (usually clouds) and
;		modulation.  Now unbundled from 'parasitic' correction.
;	SCAT - if set, amount of scattered-light correction to apply to all
;		(flat, cal, data) raw images, using ivm_scatter.pro.
;		Scat may be a scalar or a two-element array containing
;		[level,power] for ivm_scatter.
;	UNBLUR - if set, use geometry camera images to correct data
;                for seeing distortions. Set unblur to a value >3 to
;                select a tile size.
; 	DESTRETCHONLY - if set and /unblur is set, corrects large-scale
;		spatial distortions only, and does not correct defocus or blur
;       WARP - if set, warp geometry images to coalign with data images.
;		WILL be performed if UNBLUR keyword is set.
;       MATCHED - if set, the coalignment between the geometry and
;                 data cameras has already been determined. The keyword
;                 should be a string containing the name of the file
;                 which has the 4x2 transformation matrix.
;	FOURSTEP = if set, assume IVM V2 observations, four-step modulation
;		sequence. Implemented Feb 1995. Default if nframes=4 in
;		data file header.
;	CALSMOOTH = width of smoothing filter to be applied to cal matrices.
;		Default is 16. Zero turns it off.
;	SEEING = if set, calls 'ivm_sub.pro after doing the deblur.
;		Generally recommended for any gram where there is variable
;		seeing, i.e. most non-morning-survey data.
;		Recommended to set '/unbias' keyword, too.
;	UNBIAS	Calls 'ivm_unbias.pro' to eliminate any	leftover intensity 
;		biases, as manifest by offsets in histograms of each 
;		polarization image, up to 0.5% of I.  Also sets /unbias
;		keyword in xbgram7, which subtracts the appropriate 
;		continuum-polarization	image from each of the 
;		wavelength-polarization images.  NOT recommended if
;		there is very little quiet-sun in the field of view.
;	MASK = Output (nn x nn) byte array image of of 1's and 0's:
;		 1 are presumed okay to invert, 0 not (incl pukas, edges). 
;	NDERIV = number of derivative terms to use in destretch/deblur
;		code.  Default = 3.  Only other choice = 5.
;	DAYFLAT = Names of daily flat files from IVM_DAY_FLAT, string
;		array(2), name format IFIn_ccyymmdd
;	KEEP512 = if 1, then keep 512 output files, even if (binstokes && ~keep_int)
;	BINSTOKES = if 1, then flag, if > 1, then amount to bin stokes
;		spectra/images by prior to inverting.  Useful occasionally
;		for 512^2 data.  Will not remove unbinned file even if 
;		keep_int not set.
;	NO_GEOMREFLAT = if set, will not call 'ivm_geomreflat' which removes
;		residual fringing in the geometry camera.
;	SODIUM = if set, skip telluric and unbias steps
;	VERBOSE = Verbosity keyword -- values are as follows:
;		    0 [or unset] - print errors-only
;		    1 - print start/end times of blocks
;		    2 - print memory usage statistics
;		    3 - excessive verbosity, where available
;	JUMP = string, one of the following, to skip into this program.
;		'shift':	goto, DOSHIFT
;		'warp':		goto, DOWARP
;		'para':		goto, DOPARA
;		'spur':		goto, DOPARA
;		'scat':		goto, DOSCAT
;		'blur':		goto, DOBLUR
;		'seeing':	goto, DOSEE
;		'avg':		goto, DOAVG
;		'demod':	goto, DODEMOD
;		'calib':	goto, DOCAL
;		'shutsub':	goto, DOSHUTSUB [commented out - ELW-200912]
;		'unbias':	goto, DOUNBIAS
;OUTPUTS:
;	Returns status:  0 is normal, -1 indicates an error.
;
;OPTIONAL OUTPUTS:
;	Display format images of Intensity, B_total, Polar_angle, and
;		Azimuth, both full size and thumbnail. Separate files, named
;		for active region and obs time.
;
;COMMON BLOCKS:
;	None
;SIDE EFFECTS:
;	A journal file is written to disk for each call, corresponding to
;	the input file.
;RESTRICTIONS:
;	need lots of disk space! About 850Mb for the intermediate files
;	if keep_int is set, or of order 200Mb if not;
;PROCEDURE:
;MODIFICATION HISTORY:
;   V1.0
;	original algorithm by DLM, 1991
;	many modules by DLM, TRM  1991-1992
; 	modified 6/4/92 to ignore camera 1 (geom camera) images, due to
;		light leak...
; 	modified 6/23/92 to not include registration of images
; 		(fast guider ok!)
;	put into readable/understandable form: KDL, June 26,1992
;	change call to xbgram3 for: unbias, 6302.5, dl option. DLM 92.07.07
;	changed headersize offsets for dc & flat - needs to be updated
;		when stablized: KDL 8/92
;	modified so can exit cleanly w/o calling xbgram3 8/92 KDL
;	modified to include possible shutter light-leak correction 9/92 KDL
;	modified to accept old data of varying header sizes 3/3/93
;	Calls split_fil to separate data, geom images; parflat to reduce
;		flat data (instead of getflats); shifter and shifterd to
;		coregister data frames from geom sequence. 4/93 DLM
;	Allow for repeats of spectrum (ntimes) 6/93 DLM
;       calls shifter_2d (cross-corr. in 2 dim) instead of shifter  1/94 GC
;       moved the case for "jump" keyword. Now allows to read in flat
;               and/or cal if already done.    1/94 GC
;       TRM Oct 94: Modified to use xbgram5 instead of xbgram4
;	DLM Feb 95: Add fourstep keyword. If set call demod_2. Also
;		fitmap, to verify I-profile fit. And xbgram6.
;	DLM Apr 95: Add calsmooth keyword, for median smoothing of
;		calib matrix. Default is 16-px median smooth.
;	DLM Apr 95: Make deleting intermediate files the default.
;	DLM Feb 96: Implement better scattered-light correction based
;		on psf derived from May 1994 eclipse data. Changes to
;		this procedure, parfl2, fldat.
;  V2.0
;	DLM Mar 96: True breakthrough. Add capability to measure the
;	        instantaneous seeing, using the geometry images.
;	        New keywords 'unblur' and 'matched'.
;	DLM Apr 96: Add parasitic spectral light subtraction, reorder
;		procedures. Do scattered light correction to data
;		after parasitic, rather than in fldat.
;	JPW May 96: Add REGPARM keyword
;	DLM May 96: Add DETREND keyword
;	KDL Oct 96: Changed Matching algorithm to be automatic, not
;		relying on setpts & manual intervention
;  V3.0
; 	BJL ?? added in ability to do multiple platforms 
;	KDL Nov 97-Jan 98: new xbgram7/bjeff7 included, with correct calls 
;		to voigt profiles (instead of gaussians).  Found problem with
;		ABCOAL_BJL, didn't resolve (using original), and 
;		forcing warping to be done if deblur is going to be performed.
;		Bugs in DEBLUR algorithms found & fixed.   Added keywords
;		UNBIAS, MASK, SEEING, FITFUNCT.	Wavelength-vector is 
;		calculated in xbgram7 unless npts or nfus are specified.  
;		Separated out shutter-subtraction, unbiasing, seeing in 
;		jump keywords.  Added more documentation, tidied things up.
;  V3.4
;	DLM Aug 2000: Changed csmth default from 16 to nn/16.
;  V3.5
;	BJL  October 2000: Changed to use IVM_CALSTK to handle wavelength
;		dependence of calibration.
;	BJL  June 13, 2000: Added /NDERIV keyword.
;  V4.0
;	BJL  June 22, 2001: Added DAYFLAT keyword, change to use
;		IVM_FLATDATA to handle temporal variation of flat field.
;  	KDL  Oct 2001: revamped DEBLUR  to wrapper-program IVM_DESEE
;		and IVM_DESTRETCHDEBLUR
;		removed NDERIV (always 3 now)
;		re-wrote inversion area to use quick-look, jls, or triplet
;		according to new 'INVERT' keyword (triplet default)
;		YES, got all the other 'quicklook' options fixed, too.
;	BJL July 29 2002: Minor bug fixes
;	BJL Oct 3, 2002: Replace SPAWN with FILE_DELETE
;  SVN version 51
;	KDL 14 June 2006: added ability to do destretch w/o deblur,
;	set the no_limit keyword for ivm_xyoff2.  Mods mostly to make
;	things work better w/ Na data.  See updated batch_rivm.pro as well.
;	KDL 17 July 2006: decouple parasitic light correction from 
;	spurious modulation/extinction correction.  latter good for Na, 
;	former is not.   Also modified ivm_parasitic to set fraction 
;	to 0.0 if oblk wavelength is not 630.25, for calls within parfl2.
;	NOTE:  parasitic correction is always done for Iron AND Na
;	(set in batch_rivm and batch_rivm_na).
;       TRM 2006 Sep 12.  Added no_lim keyword.
;	KDL: added GEOM_REFLAT as part of the flat-fielding.
;  V5.0
;	ELW [20091112 forward to 2011]:
;	-Completely reorganized, in SVN repository.
;	-Scrapped quicklook*, and everything related to inversion code.
;	-Procedurized ivm_warp.
;	-Pointing in oblk fixed if needed.
;	-Mask used throughout, and warped by ivm_warp.
;	-clean_i used in place of prefilter/telluric code.
;	-Array I/O throughout.  As needed, convert between 3D(nx,ny,nimages) &
;		4D(nx,ny,nmod, nwavel*nrepeat).
;	-RAM usage:
;	   Unless 'fileoutput' is set, all computation from ivm_flatdata to
;	   the end is done in memory.  Amount of RAM necessary will be
;	   on the order of four (for geom+data I/O arrays) times the size
;	   of the arrays = 4 * (nx,ny,nmod,nwavel*nrepeat) * 4bytes/(float).
;	   Typically this will be around 1GB (4*512*512*4*30*2*4); but figure
; 	   about 30% more than that, maximum, for system/temp variables, and
;	   more for sodium data as well.
;	-Sodium keyword
;	-ivm_parasitic now has geometry output.  - ELW 20100511 
;	-shifter_2d now uses shrefidx (not shrefimg) for ref index. ELW20110406
;	-Removed dl option; rivm doesn't do inversions anymore.  ELW20110618
;	ELW Feb 2012 - journaling repairs (+telluric logic in clean_i)
;	ELW Apr 2012 - Fits file output, containing mask,noise,oblk extensions
;	ELW Aug 2012 - Mask computed by ivm_flatdata, and ANDed with raw mask
;	ELW Sep 2012 - Clean_i: if failure, ivm_telluric called separately
;	ELW Oct 2013 - Avoid overzealous mask cropping, byte scale gif image
;		       consistently (-.5,.5) to allow for better comparisons
;-

function rivm, calfile, datfile, save_binarr=save_binarr, $
	d0fil=d0fil, f0fil=f0fil, npoints=npoints, nfus=nfus,$
	redcal=redcal, redflat=redflat, $
	shut=shut, old_data=old_data, no_reg=no_reg, jump=jump, scat=scat,$
	fourstep=fourstep, calsmooth=calsmooth, $
	keep_int=keep_int, FileIO=FileIO, unblur=unblur,unbias=unbias,$
	warp=warp, matched=matched, spurious=spurious, parasitic=parasitic, $
	regparm=regparm, detrend=detrend, sodium=sodium, $
	seeing=seeing, mask=mask, nderiv=nderiv, dayflat=dayflat, $
	destretchonly=destretchonly, keep512=keep512, binstokes=binstokes, $
	no_lim=no_lim, shrefidx=shrefidx, NO_GEOMREFLAT=NO_GEOMREFLAT, $
	noise_calc=noise_calc, verbose=verbose, $   ; ELW-Feb-2010 additions
	GIF=GIF

version='5.0'

;
;;;========================SETUP======================================
;
setenv,'SHELL=/bin/sh'

if(KEYWORD_SET(keep_int)) THEN begin	; KEEP_INT
    keepint=1 & fileoutput=1 & save_binarr=1
endif else begin
    keepint=0 & fileoutput=0
endelse

if (keyword_set(FileIO)) then begin	; FileIO
    fileinput=1 & fileoutput=1
endif else begin
    fileinput=0
endelse

if keyword_set(sodium) then begin	; SODIUM
    sodium=1 & notelluric=1 & unbias=0
    print, '  RIVM:  sodium keyword was set.  Will skip telluric + unbias'
endif else begin
    sodium=0 & notelluric=0
endelse

if (~keyword_set(binstokes)) then binstokes=0
if (~keyword_set(keep512))   then keep512=0

if (keyword_set(verbose)) then vb=verbose else vb=0

; Set up journal file
stime = strcompress( systime() )
stime = str_sep( stime, ' ' )
ttime = str_sep( stime[3], ':' )
stime = stime[4]+stime[1]+stime[2]+'_'+ttime[0]+ttime[1]+ttime[2]
if( !journal ne 0 ) then begin
    flush, !journal             ; Flush output stream, then
    journal                     ; close existing journal file
endif
if (sodium) then $		; and start a new one
   journal, 'nalog_' + datfile + '_' +stime[0] $  
else $
   journal, 'log_' + datfile + '_' +stime[0]

print,' '
print,'***** IVM DATA REDUCTION PACKAGE v'+version+' *****'
print,' '
print, 'RIVM start: ', SYSTIME()
print
print, 'File: ', datfile

; Handle older data formats
IF( KEYWORD_SET( old_data ) ) THEN BEGIN		; Old data format
    hd = '' & hf = '' & hc = '' &  hm = ''		; Get header sizes
    PRINT, 'Old Data keyword called:'
    PRINT, 'Please enter the header sizes for the following files:'
    PRINT, '[or <Return> for the indicated defaults]'
    PRINT, 'Dark file [256]'
    READ, hd
    IF( hd EQ '' ) THEN hd=256 & PRINT, hd
    PRINT, 'Flat file [256]'
    READ, hf
    IF( hf EQ '' ) THEN hf=256 & PRINT, hf
    PRINT, 'Cal file [2176]'
    READ, hc
    IF( hc EQ '' ) THEN hc=2176 & PRINT, hc
    PRINT, 'Data file [2176]'
    READ, hm
    IF( hm EQ '' ) THEN hm=2176 & PRINT, hm
ENDIF ELSE BEGIN 					; New data format
    hd = 256 & hf = 256
    hc = 2176 & hm = 2176
ENDELSE

; Get the data header
oblk = IVM_RD_HEAD( datfile, MAGIC=magic )

;
;;; ============================= Dark files ============================= 
;
IF( KEYWORD_SET(d0fil) ) THEN begin
    PRINT, 'Using specified dark ',d0fil
ENDIF ELSE BEGIN
    ; Get the dark filename from oblk; make sure it exists.
    d0fil = STRING(oblk.dark0file)
    PRINT, '    Data file shows dark is ', d0fil
    check=find_file(d0fil, count=nf)    	; Look for specified dark file
    if (nf lt 1) || (check[0] eq '') then begin
        dkfiles=find_file('IDK*', count=nf)    	; Look for any dark file
	if (nf lt 1) then begin
	    print, "WARNING! No dark file could be found for this date.  Returning.'
	    return, -1
	endif
	print, "WARNING! Dark file specified in oblk "+ d0fil+" not found.  Using first one: "
	d0fil=dkfiles[0]
	print, "              ", d0fil
    endif
ENDELSE
d1fil = d0fil            	; Dark, geometry camera
j = STRPOS(d0fil,'IDKR') 	; See if it's a raw sequence
IF (j LT 0) THEN BEGIN   	; Yes, shuffle names and flag split need
    STRPUT, d1fil, '1', STRPOS(d1fil,'0')
END

;
;;; ============================= Flat files ============================= 
;
if (vb ge 1) then PRINT, '    Data file shows flat is ', STRING(oblk.flat0file)

if( keyword_set(f0fil)) then begin
    print, '    Using specified flat ',f0fil
endif else begin
    f0fil = STRING(oblk.flat0file)
endelse

j = STRPOS(f0fil,'IFSR')        ; See if it's a raw sequence
IF (j GE 0) THEN BEGIN          ; Yes, shuffle names and flag split need
    frfil = f0fil               ; save raw file name
    fs0fil = frfil
    STRPUT,fs0fil,'IFS0',j
    fs1fil = frfil
    STRPUT,fs1fil,'IFS1',j
    STRPUT,f0fil,'IFA0',j
    rawflatseq = 1
ENDIF ELSE BEGIN
    rawflatseq = 0
ENDELSE

; Set filenames
f1fil = f0fil                     ; Flat, geometry camera
STRPUT, f1fil,'1', STRPOS(f1fil,'0')
flatfg = f1fil+'_inv'             ; Inverse flat, geometry camera
datf = datfile
j = strpos(datf,'IMGR') & if j eq 0 then strput,datf,'IMG0',j
geof = datfile+'_df'
STRPUT, geof, 'IMG1'              ; Flatted geometry
calm = calfile+'_instr'           ; Instrument matrix

; Oblk / Mask / Noise output filenames:
datstring = strmid(datfile,5,13)
oblkptfile='oblk_pointing_'+datstring	; File in which to store poing-corrected oblk at end
rawmaskfile='mask_raw_'+datstring 	; File in which to store raw mask

; Default output is (only) FITS, so set empty filenames for binary output
; files, unless save_binarr is set
if (keyword_set(save_binarr)) then begin
    oblkfile='oblk_'+datstring		; File in which to store oblk at end
    oblkbinfile='oblk_bn_'+datstring	; File in which to store binned oblk
    maskfile='mask_'+datstring 		; File in which to store final mask
    maskbinfile='mask_bn_'+datstring 	; File in which to store binned mask
    noisefile='noise_' + strmid(datfile, strpos(datfile,'_')+1, 13)
    noisebinfile='noise_bn_'+strmid(datfile, strpos(datfile,'_')+1, 13)
endif else begin
    save_binarr = 0
    oblkfile=''
    oblkbinfile=''
    maskfile=''
    maskbinfile=''
    noisefile=''
    noisebinfile=''
endelse

;
;;;=============== Extract/Print Header Values ========================
;

; Get sizes
nn = long(oblk.spix)			; Image size, pixels, ASSUMED SQUARE
nx = nn
ny = long(oblk.ppix)
nxny = nx*ny
if (nn ne 512) then begin		; ELW20150213
   print, 'WARNING:  non-512 data.  nn=',strtrim(nn,2),'.  Continue at your own risk with .c'
   stop
endif
target = oblk.noaa_no
IF (target LT 0) THEN target = oblk.target          ;string, in 1999+ oblk

ntimes = oblk.ntimes		; Repeats of spectrum
nrepeat = ntimes
IF( ntimes EQ 0 ) THEN ntimes = 1

; Handle old data with bad header information
IF( N_ELEMENTS(npoints) GT 0 ) THEN $
	PRINT, 'Overriding NPTS from data file (',oblk.n_image_sets,')' $
ELSE npts = oblk.n_image_sets
nwavel = npts

fourf = 0
IF( N_ELEMENTS(fourstep) NE 0 ) THEN BEGIN
    PRINT, 'Overriding NFRAMES from data file (',oblk.nframes,')'
    fourf = fourstep
    nframes = 4
ENDIF ELSE BEGIN
    IF( oblk.nframes EQ 4 ) THEN fourf = 1
    nframes = oblk.nframes
ENDELSE
IF( fourf NE 0 ) THEN PRINT, '    Using 4-frame demodulation'

nmod = nframes

; Total number of images in each camera
nimages = ntimes * npts * nframes

; Spectral line  --------------------------------------------
waveln = oblk.f_wavelength
if( sodium && ((waveln LT 500) || (waveln GT 600))) then begin
   print, 'rivm: /sodium speified, but not a sodium wavelength. Exiting.'
   return, -1
endif
if( ~sodium && ((waveln LT 600) || (waveln GT 660))) then begin
   print, 'rivm: Unrecognized wavelength.  Exiting.'
   return, -1
endif


; IF( (FIX(waveln) LT 500) || (FIX(waveln) GT 700)) THEN BEGIN
;	PRINT, 'OBLK wavelength ',waveln,' is ridiculous. Assuming 6302.5A'
;	waveln = 630.25
; ENDIF
; waveln = waveln*10.		;Convert to Angstroms
;
; ; Get Lande g value
; glande = LANDE_G( waveln )
; The number of line wavelengths, NOT including the continuum wavelength
;_ nwave = npts - 1		; Wavelengths in line profile,

print,'    Wavelength is', oblk.f_wavelength, 'nm [in raw header].'
print,'    Image size is ', nn
print,'    Data are for NOAA',target,'   ',string(oblk.nss),'   ',string(oblk.ews)
print,'    Spectrum is repeated ', strtrim(ntimes,2),' times.'
print,'    Number of wavelength points:', npts
print,'    Total number of images in each camera:',nimages
print,'    Verbose level:  ', strtrim(vb,2), vb le 0 ? ' [NONE]' : $
			vb eq 1 ? ' [LOW]' : vb eq 2 ? ' [MEDIUM]' : $
			vb eq 3 ? ' [HIGH]' : ' [EXTREME]'

;
;;;=================== Processing options ===============================
;
; Spatial smoothing of calibration images
IF KEYWORD_SET(calsmooth) THEN csmth = calsmooth ELSE csmth = nn/16

;
;;;=================DARK, FLAT, CALIB=================================
;

;
; First, deal with flatting the flat, i.e. calculate the flat-field
; corrections for each wavelength.
;
IF(~KEYWORD_SET(redflat)) THEN BEGIN
    PRINT, 'RIVM V'+version+': USE RIVMCAL TO PROCESS FLATS.'
    IF( rawflatseq EQ 1) THEN BEGIN
	SPLIT_FIL,frfil, fs0fil, fs1fil
	IVM_MD_FLAT,fs0fil,f0fil
	IVM_MD_FLAT,fs1fil,f1fil
    ENDIF
    PARFL2, f0fil, d0fil, oblk, SCAT=scat, PARASITIC=parasitic
    IVM_DAY_FLAT, f0fil, flatfile
ENDIF ELSE BEGIN
    ; Get the daily flat filename
    IF( KEYWORD_SET(dayflat) ) THEN BEGIN
	flatfile = dayflat
    ENDIF ELSE BEGIN
	; Guess that the file or link to it is in the current directory
	flatfile = FINDFILE( 'IFI*' )
    ENDELSE
ENDELSE

;
; Second, the Calibrations - Reduce calibration sequence data
; Get the data dark for calibration
;
IF (~KEYWORD_SET(redcal)) THEN BEGIN
    IF( STRPOS(calfile,'ICAR') EQ 0 ) THEN $
	SPLIT_FIL, calfile, cf $	;separate out data images.
    ELSE cf = calfile
    OPENR, lun1, d0fil, /get_lun
    d0 = intarr(nn,nn)			; raw dark file
    readu, lun1, d0 			; get dark, ch 0
    CLOSE, lun1 & FREE_LUN, lun1
    d0 = REBIN( REBIN(d0,1,1), nn, nn) 	; Use mean

    calf=calfile+'_df'
    cals=calfile+'_df_st'
    ; 'calm' filename defined above for input -> bigmat (crosstalk)

    ; ELW parasitic will always be set
    IF KEYWORD_SET(parasitic) THEN BEGIN
        flatfh = f0fil+'_pl_inv'    ; Inverse flat, data camera
    ENDIF ELSE BEGIN
        flatfh = f0fil+'_inv'
    ENDELSE

    IVM_CALCAL, cf, d0, flatfh, bigmat, nn, npts, calm,$
	        SCAT=scat, KEEP=keep_int, FOURF=fourf, CSMTH=csmth
ENDIF
print, '    Calfile is: ', calfile

;
;;;=======================SPLIT DATA FILE=============================
;
gflag=1
;
;Look at the data file name. If it's IMGR..., call split_fil to
;	separate data and geometry frames. If IMG0..., assume already
;	split. (could use magic header, maybe should.)
;

j = STRPOS(datfile,'IMGR')
; 12/97 KDL - rearranged this so if you're jumping, don't have to split files
IF( j GE 0 ) THEN BEGIN
    IF (~keyword_set(jump)) THEN BEGIN
	if (vb ge 1) then print,'Splitting Data/Geometry camera images'
	SPLIT_FIL, datfile, daff, geff
    ENDIF ELSE BEGIN
	Print,'Jumping to ',jump,' - Not Splitting Files'
	daff=datfile
	geff=datfile
	strput,daff,'0',strpos(datfile,'R')
	strput,geff,'1',strpos(datfile,'R')
    ENDELSE   
ENDIF ELSE BEGIN
    j = STRPOS(datfile,'IMG0')
    IF( j GE 0 ) THEN BEGIN
	daff = datfile
	geff = daff
	STRPUT, geff, 'IMG1', j
    ENDIF ELSE BEGIN
	PRINT, 'Strange filename: ', datfile, '.'
	PRINT, 'Continuing, but no geometry file.'
	gflag = 0
	daff = datfile
    ENDELSE
ENDELSE

; Look at keyword 'jump' to skip into this program a ways
IF( KEYWORD_SET(jump) ) THEN CASE jump OF
    'shift':	GOTO, DOSHIFT
    'warp':     GOTO, DOWARP
    'para':     GOTO, DOPARA
    'spur':     GOTO, DOPARA
    'scat':	GOTO, DOSCAT
    'blur':	GOTO, DOBLUR
    'seeing':   GOTO, DOSEE
    'avg':	GOTO, DOAVG
    'demod':	GOTO, DODEMOD
    'calib':	GOTO, DOCAL
    'unbias':	GOTO, DOUNBIAS
;   'shutsub':  GOTO, DOSHUTSUB
ELSE:
ENDCASE

;
;;;===========================CORRECT POINTING========================
;
; Moved pointing here - ELW 20100621

; Create limb/puka mask from raw data [ALWAYS swap endian] and shifts array
if (vb ge 1) then print, 'ivm_raw_mask:  ', systime()
; ELW-20100323 - don't bother writing mask now... will write after warp at the end
raw_mask = ivm_raw_mask(datfile, raw_mask, swap_endian=1)
; Ignore repeat, if any, for now - ELW20100621
DEMOD, nn,nn,nwavel, 1,0,temp, rawdemod, rawfile=daff, fourf=fourf, /no_flt, /head, verbose=vb
temp=0	; Ignore input array for now - might be nice to pass to ivm_flatdata later, though - ELW20100621

if (vb ge 1) then print, 'Check/correct pointing: ', SYSTIME()

; Calculate quick & dirty Blong (and check line center), and get Icont

; Estimate Blong_qnd from Stokes demod data, and check line center
; for extreme cases.		ELW20110831
; Bl_QnD = total(rawdemod[*,*,3,0:12], 4)	; Old way
Bl_QnD = ivm_blong_qnd(rawdemod, nn, nframes, nwavel, lc=line_center, err=err)

if (err eq 0) then begin
    Icont  = rawdemod[*,*, 0, nwavel-1]
    rawdemod=0
    ; tv,bytscl(total(idats(*,*,3,0:12)),-.5,.5)
    ; tv,bytscl(total(idats(*,*,3,0:12)),-.05,.05)

    m=where(raw_mask eq 0, n_maskzeros) & m=0
    masked_percent = round(100.*n_maskzeros/(float(nx)*ny))
    if (vb ge 1) then begin
        print, 'Percentage of image masked (sky+pukas+edge):', masked_percent
    endif
    if (masked_percent ge 45) then print, 'WARNING!  Extreme limb may cause unreliable pointing.'

    ; Use Bl_QnD to attempt coalignment with MDI magnetogram
    ;    This will use the full image multiplied by mask (/full).
    ivm_fix_pointing, Bl_QnD, oblk, newoblk, raw_mask, verbose=vb, sign=sign, $
	icont=Icont, suffix='_rivm_raw', /full

    if (sign ne 0) then begin
        oblk = newoblk 			; Save updated oblk
        print,'    Pointing changed.'
        if (sign eq -1) then print, '    Original data has sign flip.'
        ; h = ivm_oblk2raw(oblk, out=oblkptfile)
    endif else begin
        print, '    Pointing was not changed.'
    endelse
endif else begin
    sign = 0
    print,'    Unable to compute Bl_Qnd.  Pointing was not changed.'
endelse
if (vb ge 1) then print, '    Computed line center:  ', strtrim(line_center,2)

;
;;;=======================FLAT-FIELD CORRECTION=======================
;

datf = datf+'_df'                 ; Flatted data

;; DATA camera   	[output is datd = fltarr(nx,ny,nimage)]
IVM_FLATDATA, daff, d0fil, flatfile[0], datf, datd, oblk, fileoutput, flatmask
print, '    Flatfile [data] is: ', flatfile[0]
IF( ~keepint && (daff NE datfile) ) THEN BEGIN
    PRINT, 'Removing ',daff ; , ' + ', flatfile[0]
    FILE_DELETE, daff	; , flatfile[0]
ENDIF

;; GEOMETRY camera   	[output is geod = fltarr(nx,ny,nimage)]
IF( (gflag NE 0) && (~KEYWORD_SET(no_reg)) ) THEN BEGIN
    geod=fltarr(nn,nn, nimages)
    print, '    Flatfile [geom] is: ', flatfile[1]
    IVM_FLATDATA, geff, d1fil, flatfile[1], geof, geod, oblk, fileoutput, flatmask2
    IF( ~keepint ) THEN BEGIN
	PRINT, 'Removing ', geff ;, ' + ', flatfile[1]
	FILE_DELETE, geff ;, flatfile[1]
    ENDIF
    IF (~KEYWORD_SET(no_geomreflat)) THEN BEGIN
	geomin=temporary(geod)
	geofin = geof+'_fldayonly'
	if (fileoutput) then FILE_MOVE,geof,geofin,/overwrite
	IVM_GEOMREFLAT,geomin,geod, oblk, fileinput, fileoutput, $
		geofin=geofin, geofout=geof, verbose=vb
    ENDIF
    flatmask = flatmask AND flatmask2
ENDIF

;
;;;==================TEMPORAL SEQUENCE REGISTRATION===================
;

DOSHIFT:
IF keyword_set(jump) THEN IF (jump EQ 'shift') THEN BEGIN
    PRINT, 'Need to set geof = name of flattded geom file'
    PRINT, 'and datf = name of flattded data file.'
    STOP, 'Stopping. Define geof, datf, then .continue'
ENDIF
IF( (gflag NE 0) && ( ~KEYWORD_SET(no_reg) ) ) THEN BEGIN
    IF( N_ELEMENTS(regparm) EQ 5 ) THEN BEGIN
	bx=regparm[0]
	by=regparm[1]
	hx=regparm[2]
	hy=regparm[3]
	no_lim=regparm[4]
    ENDIF

    ; Measure the spatial shifts
    SHIFTER_2D, nn, nimages, shifts, vb, fileoutput, geod=geod, geof=geof,$
	BX=bx,BY=by, HX=hx,HY=hy,DETREND=detrend,no_lim=no_lim,shrefidx=shrefidx

    ; Create limb/puka mask from raw data [ALWAYS swap endian] and shifts array
    if (vb ge 1) then print, 'ivm_raw_mask:  ', systime()
    ; ELW-20100323 - don't bother writing mask now... will write after warp at the end
    if (keepint) then $
	raw_mask = ivm_raw_mask(datfile, raw_mask, shifts_array=shifts, swap_endian=1, out=rawmaskfile) $
    else $
	raw_mask = ivm_raw_mask(datfile, raw_mask, shifts_array=shifts, swap_endian=1)

    ; Fracshift implemented as an alternative to shift_rigid. 	ELW-20100301
    use_fracshift=1
    if (use_fracshift) then begin
        datr = datf+'_rf'    		; Registered data using fracshift
        geor = geof+'_rf'		; Registered geo using fracshift

        if (vb ge 1) then print, 'ivm_fracshift_data [dat]: ', systime()
	ivm_fracshift_data, oblk, datd, datreg, shifts, fileinput,fileoutput,$
		fdatin=datf, fregout=datr

        if (vb ge 1) then print, 'ivm_fracshift_data [geo]: ', systime()
	ivm_fracshift_data, oblk, geod, datgeor, shifts, fileinput,fileoutput,$
		fdatin=geof, fregout=geor

    endif else begin
        datr = datf+'_re'                 ; Registered data using shift_rigid
        geor = geof+'_re'

	; Remove the integral pixel shifts to get close for deblur software    
	SHIFT_RIGID, datd, datreg, nn, nimages, shifts, fileoutput, regf=datr
	SHIFT_RIGID, geod, datgeor, nn, nimages, shifts, fileoutput, regf=geor
    endelse
    datd=0 & geod=0
    if (~keepint && fileoutput) then file_delete, datf, geof
    IF( keepint ) THEN BEGIN
	SAVE, shifts, FILE='shifts_' + datfile + '.sav'
    ENDIF

ENDIF ELSE BEGIN
    datr=datf & geor=geof 				; File
    datreg=temporary(datd) & datgeor=temporary(geod)	; Array
ENDELSE

; Make sure mask has been accounted for	prior to warp - ELW20120809
rawmask_set = (n_elements(raw_mask) eq nxny)
flatmask_set = (n_elements(flatmask) eq nxny)
if (flatmask_set) then begin
    raw_mask =  (rawmask_set) ? raw_mask AND flatmask : flatmask
endif else if (~rawmask_set) then $
    raw_mask = 0

;
;;;==================REFORM DATA FROM 3D to 4D========================
;
; Used 3D (nx,ny,nimages) up to this point; change dimensions to 4D, same data.
; After this point, any code that uses 3D for computation will accept 4D input,
; and output 4D (nx,ny,nmod,nwavel*nrepeat).
;
r1=data_3Dto4D(datreg,  nx, ny, nmod, nwavel*nrepeat)
r2=data_3Dto4D(datgeor, nx, ny, nmod, nwavel*nrepeat)
if (r1 ne 1) || (r2 ne 1) then begin
    print, 'rivm: Error converting data to 4D with dimensions: ',$
    		nx, ny, nmod, nwavel*nrepeat
    STOP, 'Stopping.  Fix datgeor/datreg arrays, or use file inputs for ivm_warp.'
endif

;
;;;==================DATA-GEOMETRY REGISTRATION=======================
;
DOWARP:
;
;; The data camera defines the spatial coordinates.
;; One transformation matrix is computed, then each of the geometry images
;; is transformed (along with the mask) using this matrix.
; define datr, geor
IF keyword_set(jump) && (jump EQ 'warp') THEN BEGIN
    PRINT, 'Need to set datgeor = registered (Shifter_2d) geom data array,'
    PRINT, 'and datreg = registered data array.  These must be 4D, of the'
    PRINT, 'form fltarr(nx, ny, nmod, nw*nrep).  Also need raw_mask,'
    PRINT, 'which can be read in directly from the mask file, to an array'
    PRINT, 'of the form:  bytarr(nx, ny) .  "igeomr" is output var.'
    STOP, 'Stopping. Define raw_mask, datgeor, datreg, then .continue'
ENDIF
IF (KEYWORD_SET(warp) || KEYWORD_SET(unblur)) THEN BEGIN
    georb = geor+'_wa'

    IVM_WARP, igeomr, oblk, fileinput, fileoutput, vb, $
    		geord=datgeor,datrd=datreg,georf=geor,datrf=datr,outfile=georb,$
		matched=matched, mask=raw_mask

    datgeor=0	; Free up input geometry array
    if (~keepint && fileoutput) then file_delete, geor

ENDIF ELSE BEGIN
    igeomr=temporary(datgeor)
    georb = geor
ENDELSE

;
;;;======= EXTINCTION, PARASITIC LIGHT & SPURIOUS MODULATION  CORRECTIONS ====
;
DOPARA:
IF keyword_set(jump) THEN IF (jump EQ 'para' || jump eq 'spur') THEN BEGIN
    PRINT, 'Need to set georb = name of registered, warped geom file'
    PRINT, 'and datr = name of registered data file.'
    STOP, 'Stopping. Define georb, datr, then .continue'
ENDIF
IF KEYWORD_SET(spurious) THEN BEGIN
    geomp = georb+'_ic'
    datp = datr+'_ic'

    IVM_SPUR_MOD, oblk, fileinput, fileoutput, vb, $
    			geomr=igeomr, datr=datreg, $	; Input arrays
			geomc=igeomp, datc=idatp,  $	; Output arrays
			fgeomr=georb, fdatr=datr,  $	; Input files
			fgeomc=geomp, fdatc=datp  	; Output files

    igeomr=0 & datreg=0		; Free up input arrays
    if (~keepint && fileoutput) then file_delete, georb, datr
ENDIF ELSE BEGIN
    igeomp=temporary(igeomr) & idatp=temporary(datreg) 	; Arrays
    geomp = georb & datp = datr				; Files
ENDELSE

IF KEYWORD_SET(parasitic) THEN BEGIN
    datrp = datp+'_pl' 		; 'datrp' file written only if 'fileoutput' set
    georp = geomp+'_pl' 	; 'georp' file written only if 'fileoutput' set

    ivm_parasitic, oblk, fileinput, fileoutput, vb, mask=raw_mask, $
	       datain=idatp, dataout=idatrp, dataf=datp, datoutf=datrp, $ ; DATA
	       geomin=igeomp, geomout=igeorp, geomf=geomp, geooutf=georp  ; GEOM

    idatp=0 & igeomp=0		; Free up input arrays
    if (~keepint && fileoutput) then file_delete, geomp, datp
ENDIF ELSE BEGIN
    idatrp=temporary(idatp) & igeorp=temporary(igeomp)	; Arrays
    datrp = datp & georp = geomp			; Files
ENDELSE

;
;;;==================SCATTERED LIGHT CORRECTION=======================
;
DOSCAT:
IF KEYWORD_SET(jump) THEN IF (jump EQ 'scat') THEN BEGIN
	PRINT, 'Need to set geomp = name of registered, warped, para geom file,'
	PRINT, ' and datrp = name of registered, para-corrected data file.'
	STOP, 'Stopping. Define georp and datrp, then .continue'
ENDIF
IF KEYWORD_SET(scat) THEN BEGIN
    datps = datrp+'_sc'
    geops = georp+'_sc'

    if (vb ge 2) then m0=memory(/current)

    descatter, oblk, fileinput, fileoutput, vb, nn,npts,ntimes, $
		scat=scat, fourstep=fourstep, 	$
		gin=igeorp, din=idatrp,		$	; Input arrays
		gout=igeops, dout=idatps, 	$  	; Output arrays
		fgin=georp, fdin=datrp,   	$	; Input files
		fgout=geops,fdout=datps   		; Output files

    if (vb ge 2) then begin
	mh=memory(/high) & m1=memory(/current) 
	print, 'DESCATTER Peak memory:',mh,' Current:', m1, ' Difference: ', m1-m0
    endif
    idatrp=0 & datrp=0
    if (~keepint && fileoutput) then file_delete, georp, datrp
ENDIF ELSE BEGIN
    datps = datrp & geops = georp	; Files
    idatps = temporary(idatrp)		; Data array
    igeops = temporary(igeorp)		; Data array
ENDELSE

;
;;;==================DEBLUR AND DESTRETCH=============================
;
DOBLUR:
IF( KEYWORD_SET(jump) ) THEN IF( jump EQ 'blur' ) THEN BEGIN
    PRINT, 'Need to set geops = name of scatter-corrected geom file,'
    PRINT, ' and datps = name of scatter-corrected data file.'
    STOP, 'Stopping. Define geops and datps, then .continue'
ENDIF

seef=''
if (save_binarr || keep_int) then $
    seef='seeing_'+STRMID(geof,0,STRPOS(geof,'_',5))
IF( KEYWORD_SET(unblur) ) THEN BEGIN
    datw=datps+'_ds'
    geow=geops+'_ds'

    ivm_desee, oblk, fileinput, fileoutput,vb, 		  $
	nodeblur=destretchonly, rmssee=rmssee, seef=seef, $
	gin=igeops, din=idatps, gout=igeow, dout=idatw,   $	; ARRAY
	fgin=geops, fdin=datps, fgout=geow, fdout=datw		; FILE

    idatps=0 & igeops=0
    if (~keepint && fileoutput) then file_delete, geops, datps, geow
ENDIF ELSE BEGIN
    datw=datps & geow=geops

    ivm_desee, oblk, fileinput, fileoutput,vb, gin=igeops, fgin=geops, $
	rmssee=rmssee, seef=seef, /seeingout  ; only compute/write seeing file

    idatw=temporary(idatps) & igeow=temporary(igeops)
    if (~keepint && fileoutput) then file_delete, geops, datps
ENDELSE

igeow=0 	; Finished with geometry, so free up the memory  ELW-20091222

;
;;;=======================IMAGE SUBSTITUTION====================
;
; Image-substitution, if seeing is variable, substitute
; bad-images for good; only works on multiple-spectral repeasts.  
; added 1/98 KDL 
DOSEE:
IF( KEYWORD_SET( jump ) ) THEN  IF( jump EQ 'seeing' ) THEN BEGIN
	PRINT, 'Need to set datw = name of deblurred data file.'
	PRINT, 'Make sure seeing_'+datfile,'exists!'
	STOP, 'Stopping. Define datw, then .continue'
ENDIF
IF(KEYWORD_SET(seeing)) THEN BEGIN
	gg = WHERE(rmssee LT (MEDIAN(rmssee) - STDEV(rmssee)))
	IF ((N_ELEMENTS(gg) GE (nimages/ntimes)) || (ntimes eq 1)) THEN BEGIN 
	    PRINT,'WARNING: Too many bad images or not enough repeats - skipping IVM_SUB'
	    idatws=temporary(idatw)
	    datws=datw
	ENDIF ELSE BEGIN 
	    datws=datw+'_sub'

	    if (vb ge 2) then m0=memory(/current)

	    ivm_sub, oblk, rmssee, fileinput, fileoutput, verbose=vb,  	$
	    		indat=idatw, outdat=idatws, outfile=datws, infile=datw
	    if (vb ge 2) then begin
	        mh=memory(/high) & m1=memory(/current)
	        print, 'ivm_sub Peak memory:',mh,' Current:', m1, ' Difference: ', m1-m0
	    endif
	    idatw=0
	    if (~keepint && fileoutput) then file_delete, datw
	ENDELSE
ENDIF ELSE BEGIN
	idatws=temporary(idatw)
	datws=datw
ENDELSE

;
;;;=======================AVERAGE SPECTRUM REPEATS====================
;
DOAVG:
IF( KEYWORD_SET( jump ) ) THEN  IF( jump EQ 'avg' ) THEN BEGIN
	PRINT, 'Need to set datws = name of deblurred data file.'
	STOP, 'Stopping. Define datws, then .continue'
	RESTORE, 'intensities_'+datfile
ENDIF

;********** need to get solar acceleration, jerk images before averaging
;
; Average spectrum-repeats if any; data will become (nx,ny,nmod,nwavel)
;
IF( ntimes GT 1 ) THEN BEGIN
    datra=datws+'_av'
    if (vb ge 2) then m0=memory(/current)
    if (vb ge 1) then print,'Average spectra: ',systime()
    SPEC_AVE, nn, nmod, nwavel, nrepeat, fileinput, fileoutput, /sum, $
   		datin=idatws, datout=idatra, $		; Array I/O
   		fdatin=datws, fdatout=datra		; File I/O
    if (vb ge 2) then begin
        mh=memory(/high) & m1=memory(/current)
        print,'SPEC_AVE memory, Start:',m0,'Peak:',mh,' Curr:',m1,' Diff:',m1-m0
    endif

    ; Commented out the following line, per discussion with KD:
    ;   Since all rivm output (and inversion) will only have one repeat,
    ;   keeping the original value in the oblk is useful going forward, so
    ;   we know the gram's history.		ELW-20110721
    ; oblk.ntimes=1L		; Now only one wavelength scan (aka repeat)

    idatws=0
    if (~keepint && fileoutput) then file_delete, datws
ENDIF ELSE BEGIN
    idatra=temporary(idatws)
    datra=datws
ENDELSE

;
;;;=======================STOKES DEMODULATION=========================
;
DODEMOD:
IF KEYWORD_SET(jump) THEN IF jump eq 'demod' THEN BEGIN
    PRINT,'Need to set datra = name of spectrum-averaged data file.'
    STOP,'Stopping. Define datra, then .continue'
ENDIF

dats=datra+'_st' 	; demodulated data file

; Depending on 'fourf', takes 4 or 6 frames input, outputs 4 frame uncal Stokes
; If data are sodium, skip the normalization for QUV.
demod, nn,nn,nwavel, fileinput, fileoutput, idatra, idats, $
			rawfile=datra, stkfile=dats, fourf=fourf, nonorm=sodium

oblk.nframes=4L		; Output definitively has 4 polarizations
idatra=0		; Free up input array memory
if (~keepint && fileoutput) then file_delete, datra

;
;;;=======================APPLY CALIBRATION MATRIX====================
;
docal:
if keyword_set(jump) then if jump eq 'calib' then begin
    print,'Need to set dats = name of demodulated data file.'
    stop,'Stopping. Define dats, then .continue'
endif
;	apply calibration matrix to data
datc=dats+'_ca'

; If calmat already exists, read it in
if (keyword_set(redcal)) then begin
    ; ELW 20091229 - Moved here from  'Second, the Calibrations'
    bigmat = FLTARR(nn,nn,4,4)	; for calib matrix
    openr, lun1, calm, /get_lun
    readu, lun1, bigmat
    close, lun1 & FREE_LUN, lun1
    if (vb ge 1) then begin
        print, '  Bigmat was read from: ' + calm
    endif
endif

;==========================================================================
; Determine whether to flip the sign in ivm_calstk, either to correct a
; v-sign error in the data, or to correct an error in the cal matrix.
;
; If the MDI pointing correction indicates that the v-sign is correct,
; and the [3,3] element of the calmatsum(bigmat,nn,nn) is negative, then
; there is a problem with the calibration matrix, and it needs to be
; post-multiplied by 'fixmat' in ivm_calstk.  In order to do this, the
; 'signflip' keyword is set to one, and a note is made.  ELW-20110616
;==========================================================================
bigmatsum = calmatsum(bigmat,nn,nn)
if (vb ge 1) then begin
    print, 'calmatsum(bigmat,nn,nn):'
    print, bigmatsum
endif

if (sign eq -1) then begin
    signflip = 1
endif else if (sign eq 0) then begin
    signflip = 0
endif else if (bigmatsum[3,3] < 0) then begin
    signflip = 1
    print, 'WARNING: bigmatsum[3,3] is ',strtrim(bigmatsum[3,3],2),' (< 0) but v-sign is fine.'

    ; Log this date to a file
    logfile = '/data/SOLDAT/IVM/NA-SURVEY/RightSignWrongCal.txt'
    ivm_tfile_log, logfile, datstring
    print, '  Added ', datstring, ' to list of dates where this is the case, in file:'
    print, '   ',logfile
    print, '  Post-multiplying calibration matrix to prevent extra v-sign flip.'
endif else begin
    signflip = 0
endelse

;
; Code to get wavelength dependence of calibration variation coefficients from
; a file was here.  Removed because only one such entry in a cal file exists,
; in ivm_calvar_db, which is for 20000929. 
; The contents are contained in this array:       [ELW-20100831]
;
calvar_20000929 = $
	[[[-0.163853,0.765526], [-0.338141,20.7501], [-0.135592,-23.0043]], $
	 [[-0.354664,-8.58510], [0.863920,-14.3640], [4.33266,-22.3554]], $
	 [[0.781244,-0.199167], [-3.17875,49.2536],  [-0.425414,-11.3019]]]

ivm_calstk, oblk, fileinput, fileoutput,vb,bigmat, idatc, iss=idats, $
	outf=datc, sfile=dats, caldbarr=calvar_20000929, signflip=signflip
			; Use cfile=calm to get bigmat from file
			; Use caldbfile=dbfile to read calvar from file

idats=0		; Free up memory held by input array
if (~keepint && fileoutput) then file_delete, dats

;
;;;=======================SHUTTER LEAKAGE SUBTRACTION====================
;
; Shuttersub code rediscovered, and added to repository 20100301.
; Untested.  Either use file I/O or modify it to take arrays.
; Only needed for OLD (before mid-90's w/LCD shutters) data.  ELW-20091229
;
;DOSHUTSUB:
;IF keyword_set(jump) then if jump eq 'shutsub' THEN BEGIN 
;    print,'Need to set datc = name of calibrated Stokes data file.'
;    STOP,'Stopping. Define datc, then .continue'
;ENDIF
;IF (keyword_set(shut)) THEN BEGIN
;    datsh=datc+'_sh'
;    shuttersub,datc,shut,datsh
;    IF ~keepint THEN BEGIN
;        print,'Removing ',datc
;        FILE_DELETE, datc
;    ENDIF
;ENDIF ELSE BEGIN
;    datsh = datc
;ENDELSE

;
;;;=======================CHECK AND ERODE MASK======================
;

; Create mask with ivm_sky_mask, if raw_mask was not set by ivm_raw_mask - ELW20100301
if (size(raw_mask, /n_dim) eq 0) then begin
    print, 'Obtaining mask from ivm_sky_mask: ', SYSTIME()
    mask=ivm_sky_mask(idatc[*,*,0,nwavel-1],oblk,limbdark,out=maskfile)
endif else begin
    ss=intarr(3,3)+1
    mask=erode(byte(raw_mask),ss)
    if (save_binarr) then begin
	if (vb ge 1) then print, 'Writing final mask to file:  ', maskfile
	openw, lun1, maskfile, /get_lun
	writeu,lun1, mask & close, lun1 & free_lun,lun1
    endif
endelse

;
;;;=======================PRE-CLEAN STOKES SPECTRA====================
;

; Revision 33 onward of Eric branch of rivm: removed prefilter/telluric code
; and committing to clean_i. 

; For sodium, "notelluric" should be set to 1.

tim = reform(idatc[*,*,0,*])
tim = clean_i(tim, oblk, verbose=verbose, mask=mask, notelluric=notelluric, err=err)

if (err ne 0) then begin
    print, 'clean_i failed.'
    datc = datc + '_t'
    if (err ne 2 && ~notelluric) then begin
        print, 'Telluric line removal was skipped due to clean_i failure.  Calling ivm_telluric separately.'
        tim = IVM_TELLURIC(tim, oblk); NOTE, scale=1.25 may be good here, don't want to hard-wire it, though
    endif
endif else begin
    if (sodium) then datc = datc+'_c' else datc = datc+'_ct'
endelse
idatc[*,*,0,*]=tim

; If unbias keyword not set, this will be the final output
if (save_binarr) && (~keyword_set(unbias) || keepint) then begin
    openw, lun, datc, /get_lun
    writeu,lun, idatc
    close, lun & free_lun, lun
endif

;
;;;=======================UNBIAS LINE POLARIZATION====================
;
DOUNBIAS:
IF keyword_set(jump) then if jump eq 'unbias' THEN BEGIN 
    print,'Need to set datc = name of calibrated Stokes data file to unbias'
    STOP,'Stopping. Define datc, then .continue'
ENDIF
IF (keyword_set(unbias)) THEN BEGIN
    datu=datc+'_unb'
    fout = (fileoutput || save_binarr)

    ; Call ivm_unbias, force file output if this is the last step
    ivm_unbias, nn,npts, fileinput, fout, idatc, idatu, verbose=vb, $
			infil=datc, outfil=datu

    if (vb ge 1) then print,'ivm_unbias wrote file: ', datu

    idatc=0	; Free up memory held by input array
    if (~keepint && fileoutput) then file_delete, datc
ENDIF ELSE BEGIN
    datu = datc
    idatu = temporary(idatc)
ENDELSE

;
;;;=======================SAVE OBLK HEADER=========================
;
if (binstokes eq 0 || keepint) then $
    rawoblk = ivm_oblk2raw(oblk, outfile=oblkfile)	; Write to 'oblkfile'

;
;;;===================SAVE REPRESENTATIVE GIFS=====================
;
if (keyword_set(gif)) then begin		; ELW-20110330

    ; Get a cropped mask, and only use it if it is not overzealous
    maskcr = ivm_mask_crop(mask, xr, yr)	; Compute cropped mask
    maskcrpercent = 100.*(xr[1]-xr[0]+1)*(yr[1]-yr[0]+1)/double(nn*nn)
    if (maskcrpercent le 50) then maskcr=mask

    ; Calculate Blong_qnd and save gif images
    Bl_QnD = ivm_blong_qnd(idatu, nn, nframes, nwavel, lc=line_center, err=err)
    img = bytscl(Bl_QnD * maskcr, -0.5,0.5)	; Mask, then byte-scale image
    nn_thum = 128				; Dimension of thumbnail image
    img_t = rebin(img, nn_thum, nn_thum)	; Create small thumbnail image
    gif_file = 'IMGS_'+strmid(datu,5,13) + '.gif'	; Gif filename
    gif_thum = 'IMGS_'+strmid(datu,5,13) + '_t.gif'	; Thumbnail gif filename
    write_gif, gif_file, img			; Write full-size gif
    write_gif, gif_thum, img_t			; Write thumbnail
endif

;
;;;===================BIN IMAGES, MODIFY OBLK====================
;
; If binstokes eq 1, bin down images by a factor of two by default
; File output is forced; filename will have _bn suffix.
; Mask and oblk are similarly binned and written out.
;

IF (binstokes ne 0) THEN BEGIN
    datb = datu + '_bn'
    IF (binstokes gt 1) THEN dbin=binstokes ELSE dbin = 2
    if (vb ge 1) then print, 'IVM_BINDOWN start:  ',systime()

    ivm_bindown, idatu,idatb,oblk, nn,nwavel,dbin,fileinput,1, $
				in=datu,out=datb,noblk=oblkbin
    if (vb ge 1) then print,'Wrote binned data to file: ', datb

    ; Bin the mask, too
    ivm_bindown_mask, mask, binmask, nn, dbin, 1, outfile=maskbinfile

    ; Write the new oblk corresponding to binned data
    rawoblk = ivm_oblk2raw(oblkbin, out=oblkbinfile)  ; Write to 'oblkbinfile'

    nn = nn/dbin
    nx = nx/dbin
    ny = ny/dbin

ENDIF ELSE BEGIN
    datb = datu
    idatb = temporary(idatu)
ENDELSE

;
;;;=======================NOISE CALCULATION======================
;
; Calculate noise level, write to file (noise_YYYYMMDD.HHMM) if
; save_binarr is set; otherwise, save array for FITS output.

; Noise calculation, to be written to noise file
if (vb ge 1) then print, 'IVM_NOISE_CALC: ', SYSTIME()

if (binstokes eq 0) then begin
    noise = ivm_noise_calc(idatb, oblk, noise, outfile=noisefile)
endif else begin
    if (keepint) then $			; Compute noise for unbinned
        noise = ivm_noise_calc(idatu, oblk, noise, outfile=noisefile)
    noise = ivm_noise_calc(idatb, oblkbin, noise, outfile=noisebinfile)
endelse

; If binning was done, reset mask and oblk to binned versions
if (binstokes ge 1) then begin
    mask = binmask
    oblk = oblkbin
endif

;
;;;====================LINE CENTER CALCULATION=====================
;
if (vb ge 1) then print, 'Calculating line center at each pixel: ', SYSTIME()

; Calculate line center at each unmasked pixel   - ELW20120820
lcarr=fltarr(nx,ny)
for j=0,ny-1 do begin
    for i=0,nx-1 do begin
	if (mask[i,j] ne 0) then begin
            line = reform(idatb[i, j, 0, *])
            lcarr[i,j] = (lc_find(line, 5, nwavel - 5 - 1, 9))[0]
	endif
    endfor
endfor

; Fix pathological cases		[per KDL-20131008]
oklc=where((mask eq 1) AND (lcarr gt 0) AND (lcarr lt nwavel-1), count)
if (count gt 0) then begin
   lcmed=median(lcarr(oklc))
   lcstdev=stdev(lcarr(where(mask eq 1)))
   badlc=where(lcarr lt 0 OR $
	       lcarr gt nwavel - 1 OR $
	       lcarr lt (lcmed - 10*lcstdev) OR $
	       lcarr gt (lcmed + 10*lcstdev), count)
   if (count gt 0) then begin
      if (verbose gt 1) then $
         print,'Inserting median LC value (', $
      		  strtrim(lcmed,2),') for ', strtrim(count,2), ' pixels.'
      lcarr(badlc) = lcmed
   endif
endif else begin
   print, 'WARNING!  Major LC calculation problem at all pixels.'
endelse

;
;;;=======================FITS FILE OUTPUT=======================
;

rivm_wrfits, oblk, idatb, datb, mask, noise, lcarr, vb

idatu=0 & idatb=0	; Done with the data arrays

;
;;;====================CLOSE JOURNAL & RETURN====================
;
print
PRINT, 'RIVM end: ', SYSTIME()
print
IF (!JOURNAL NE 0) THEN BEGIN
    FLUSH,!JOURNAL              ;Flush output stream, then
    JOURNAL                     ;close existing journal file
;   stime = STRCOMPRESS( SYSTIME() )	; ELW2012.02.14 commented out
;   stime = STR_SEP( stime, ' ' )
;   ttime = STR_SEP( stime[3], ':' )
;   stime = stime[4]+stime[1]+stime[2]+'_'+ttime[0]+ttime[1]+ttime[2]
;   JOURNAL,'log_'+stime[0]     ;and start a new one.
ENDIF

RETURN,0
END