pro sp_prep, $
sp4dfiles, $ ; array of file pathnames for one operation
l1index, l1data, $ ; return index(header) and data(*,*,4), for last file processed

pause=pause, $ ; pause after drift determination
display=display, _extra=_extra, version=version, $
verbose=verbose, $ ; if set, show filenames as processed
noinit=noinit, $ ; if set, skip pass-1 and restore thermdrift file instead
outdir=outdir, $ ; if set, output fits files are saved in this directory
; if not set, result fits are saved in level-1 tree
rapidtherm=rapidtherm,	$ ;  if set, shortens the smoothing width for
tryharder=tryharder	;  if set, does more thorough treatment to extract
;  data in the presence of packet loss

;  the thermal variations along the slit to accommodate the rapid
;  thermal variations when present


common sp_prep_com, $
scoef,degn,sltpos,nftot, vdeg,vigcoef, calpos,skx, $ 
;  fast, standard gain tables for default ROI
f_mflat, f_shftsl,f_specrfit,f_specshft,f_svrr, $; gainfile
s_mflat,s_shftsl,s_specrfit,s_specshft,s_svrr, $
;  darks for early mission
darkavg_std,darkavg_fast, $ 
;  eclipse data fast, standard gain tables for full ROI
sr_mflat,sr_shftsl,sr_specrfit,sr_specshft,sr_svrr, $
fr_mflat,fr_shftsl,fr_specrfit,fr_specshft,fr_svrr, $
;  eclipse darks, from full ROI
full_darkavg_std,full_darkavg_fast


; Name: sp_prep
;
; Purpose: Calibrate SOT/SP Level0 -> Level1
; This routine makes two passes through the data. The first pass
; determines the thermal shifts in the spectral and slit dimensions,
; and the second pass calibrates, merges the two CCDSIDEs, and
; applies the drift corrections.
;
; Input Parmeters:
; sp4dfiles - list of one or more SOT SP IQUV 4D files
;
; Output Paramters:
; lindex, ldata - Level1 Calibrated 'index,data' (see Restrictions)
;
; Keyword Parameters:
; _extra - keyword inherit to calib_sbsp.pro
; show_params (switch) - if set, show params->tty (verbose!)
; param_history (switch) - if set, show all params -> .HISTORY
; display=display (switch) - if set, do display
;
; Calling Sequence:
; IDL> sp_prep,<SP4D QUIV>, l1index, l1data [,/display] [,/show_params]
;
;
; Side Effects:
; If number of input files gt 1 or no output parameters
; specified, then output file generation is assumed
;
; Restrictions:
;
; History:
; 17-Oct-2006 - S.L.Freeland - Wrapper for SSW-version of
; Bruce Lites SP calibration suite
; Lots of .HISTORY info added
; 23-Jan-2007 K.Ichimoto
; handle sp4dfiles(*) as a set, write fits in outdir with binary extension
; add 'noinit' keyword, and common block to store calib data
; 25-jun-2007 - S.L.Freeland - define OUTDIR if not present, concat_dir
version=0.11 ; 2007.1.24 BL/KI
version=0.12 ; 2007.1.29 KI, add SPWLSFT0 and SPSLSFT0 keywords
; 24-oct-2007 - S.L.Freeland, B.W. Lites - new version in preparation
;		for bulk processing.  Includes:
;		-- fix for allowing processing of short (1 or greater)
;		   sequences of files
;		-- set smoothing width for thermal drifts as a fixed 
;		   time interval, not a fixed number of files
;		-- update gapsmth_sbsp.pro for the above smoothing 
;		-- make plotting of thermal drift only when /display set
;		-- update getdark_sbsp.pro to account for long-term drifts
;		   of dark levels, adjust the offsets if CAMDACA, CAMDACB
;		   are 5 instead of the usual 7 used so far
;		-- includes new method of correcting for flat field
;		   with 2-dimensional variation of the vignetting in the
;		   slit scan image
version=1.00 ; 2007.10.24  (above changes)
;  16-nov-2007 - B. Lites
;		-- refined gapsmth_sbsp.pro to better treat ends of 
;		   segments near gaps: it smoothly transitions from the 
;		   smooth function with /edge_truncate to a linear function
;		   fitted to the smoothing width at the end of each segment
;		-- in thermd_sbsp.pro, insert tests for bad data: stops the 
;		   processing when encountered (data must be fixed and sp_prep 
;		   run manually)
;		-- in thermd_sbsp.pro, and calib_sbsp.pro, insert checks to be 
;		   sure polynomial fits have adequate number of samples, 
;		   otherwise no polynomial fits are attempted
;		-- correct an error that previously halved the smoothing width:
;		   removed definition of smoothing width in gapsmth_sbsp.pro
;		-- in sp_prep.pro, use gapsmth_sbsp for smoothing segements of
;		   spectral drift (previously not done)
version=1.01 ; 2007.11.16  (above changes)
;  03-ded-2007 - B. Lites
;		-- Using values computed by R. Centeno Elliott from a comparison
;		   with G-band, update the header values for xscale,yscale,
;		   xcen,ycen, and add the keyword FGSPROT: rotation angle 
;		   between G-band and SP slit
;		-- update stksimages_sbsp.pro to calculate images of pointing
;		   information based on the header pointing values, and
;		   to add scattered light profile to output.  Add or modify
;		   routines to do this: julian1200_sbsp.pro, solareph_sbsp.pro,
;		   solcor_sbsp.pro, sphtri0_sbsp.pro, sqrtfn_sbsp.pro, 
;		   bapp_sbsp.pro, sunframe_sbsp.pro
version=1.02 ; 2007.12.03  (above changes)
;		-- Correct updates of xscale,yscale keywords to reflect the
;		   presence of binning (i.e. fast maps)
;		-- update the pvpvcomp_sbsp.pro, stksimages_sbsp.pro to better
;		   treat the behavior of pv near the noise threshold (according
;		   to suggestion by Stenflo)
;		-- remove from stksimages IDL-save file output some unneeded 
;		   parameters from fitting the thermal drifts.  These parameters
;		   remain in the output file thermdrift.save
;		-- add verbose keyword to suppress printing while running
;		   (added to thermd_sbsp.pro and calib_sbsp.pro)
;		-- change the output FITS file names to SP3D, in 
;		   stksimages_sbsp.pro find FITS files containing only
;		   SP3D for constructing output images
;		-- alter bapp_sbsp.pro to have a softer median filter, and one
;		   that operates along the slit instead of in the scan direction.
;		   This allows stksimages_sbsp.pro to work with single FITS files
version=1.03 ; 2008.01.30  (above changes)
;		-- Change calib_sbsp.pro and thermd_sbsp.pro to allow for only
;		   one CCDSIDE digitized
;		-- Correct thermd_sbsp.pro to exit if it encounters files with
;		   all zero values, and print out the offending file names.
;		-- Correct gapsmth_sbsp.pro to allow for segments <3 steps long
;		-- Correct sp_prep.pro for correct action when keyword noinit is set
version=1.04 ; 2008.02.22  (above changes)
;		-- add mods to correct the sign of the dispersion (CDELT1)
;		   and to adjust CRVAL1 to its actual value based on the
;		   spectral ROI
;		-- introduces a new keyword rapidtherm that, when set, does a
;		   much narrower boxcar smoothing on the thermal drift along the
;		   slit, but still does the normal smoothing of the wavelength
;		   drifts
version=1.05 ; 2008.04.18  (above changes)
;		-- handle errors in thermd gracefully (AdW)
version=1.06 ; 2008.05.09  (above changes)
;		-- Due to X-band failure, frequent packet loss in data in 2008.
;                  Adapted to handle frequent packet loss.  Standard option:
;                  rejects frames with any packet loss. Those files are considered 
;		   irrecovably lost, and are not retained in the level1 data set.
;                  /tryharder keyword: when set, attempts to make
;                  use of available good images in frame with some packet loss
version = 1.04 ; 4-Sep-2008 - B. Lites (above changes)
;		-- New version of gapsmth_sbsp.pro that treats ends of segments
;		   between gaps better: ends are fit with linear trend to avoid
;		   end effects of smooth.pro.  Better estimates of time steps in 
;		   gapsmth_sbsp.pro.  Writes more variables to thermdrift.save
;		-- Modifications to correct current versions to restart after pass1
;		-- Modification to treat gaps near ends of map
;	        -- Modified gain/dark to treat different ROI, new gain/dark data
version = 1.07 ; 27-Jan-2010 - B. Lites (above changes)
;		-- Simplified the calculation of gains/darks and put all of the
;		   selection of such into recovergndk_sbsp.pro.  Required changing
;		   the common area sp_prep_com (much simpler now)
;		-- Updated the residual I->QUV crosstalk computations with data from
;		   several years, now selects closest data for this correction.
;		   The relevant data files are placed in the SSW SP CALIBRATION data
;		   base.  NOTE: the sp_prep.pro code must be updated with 
;		   hardwired file names and years as newer crosstalk data becomes
;		   available, usually with eclipse season dark/flat measurements
;		-- Sort the list of input file names to be in chronological order
version = 1.08 ; 27-Jul-2011 - B. Lites (above changes)
;		-- update call to calib_sbsp to pass the year of observation so that
;		   routine may use the new residcross2_sbsp.pro that works differently
;		   for 2013 onward
version = 1.09 ; 12-Nov-2013 - B. Lites (above changes)
;		-- update yresid and xtalkfiles for 2014 season flat/xtalk data
version = 1.10 ; 6-Apr-2015 - B. Lites (above changes)
;
;      6-Dec-2019 - R. Centeno: Scan FITS files for corrupt headers
;                   that break the Lev1 processing. All changes are in
;                   sections that start with --- By RCE and end with
;                   --- END By RCE. These changes propagate to
;                   thermd_sbsp.pro, getdark_sbsp.pro and the
;                   automated telemetry dropout detector.
;
;-
;


; set path to calibration files
scal=get_logenv('SOT_SP_CALIBRATION')
if scal eq '' then set_logenv,'SOT_SP_CALIBRATION',$
   concat_dir('$SSW_SOT','calibration')

;----------- include all the seldom-changed fixed parameters ------------------------
params_file=file_search('$SOT_SP_CALIBRATION','params_sbsp*.com')
calparams_file=file_search('$SOT_SP_CALIBRATION','calparams_sbsp*.com')
sotcal=[last_nelem(params_file),last_nelem(calparams_file)] ; << may want close
sotcal=str_replace(sotcal,'//','/')
npf=n_elements(sotcal)
case 1 of
   total(file_exist(param_files)) eq npf: pdata=rd_tfiles(param_files)
   total(file_exist(sotcal)) eq npf: pdata=rd_tfiles(sotcal)
   else: begin
      box_message,['Cannot find parameter files under $SOT_SP_CALIBRATION',$
                       'returning...']
      return
   endcase
endcase
pedata=strnocomment(strarrcompress(pdata),comment=';')
if keyword_set(show_params) then box_message,pedata
pstat=1
for i=0,n_elements(pedata)-1 do begin
   estat=execute(pedata(i))
   pstat=pstat and estat
endfor
if not pstat then box_message,'At least one error setting parameters??'
  
;  keyword tryharder activation
tryhard = 0
if keyword_set(tryharder) then tryhard = 1

;**** pass-1 ***************************************************************
;------- restore calibration data ----------------------------------

; read in residual crosstalk vs slit scan coefficients
;  define years for which residual crosstalk files are present
yresid = [2006.,2007.,2009.,2010.,2011.,2012.,2013.,2014.,2016.,2017., $
	2018.,2019.]
;  define the corresponding yearly residual calibration crosstalk files
;  that are resident in the Hinode SOT/SP SSW calibration data base
;xtalkfiles = ['residxtalk_20061026.geny', 'residxtalk_20071015.geny', $
;	'residxtalk_20090731.geny','residxtalk_20100519.geny', $
;	'residxtalk_20110723.geny' $
;	]
;  new processing, variation along slit only
;xtalkfiles = ['residxtalkb_20061026.geny', 'residxtalkb_20071015.geny', $
	;'residxtalkb_20090731.geny','residxtalkb_20100519.geny', $
	;'residxtalkb_20110723.geny','residxtalk_20120702_114140.geny', $
;  progressively replace residual crosstalk data with newer analysis
;  These newer files are placed in the calibration data volume on SSW
xtalkfiles = ['residxtalk_20061026_230257.geny', 'residxtalkb_20071015.geny', $
	'residxtalk_20090731_140100.geny','residxtalk_20110723_125235.geny', $
	'residxtalk_20110723_125235.geny','residxtalk_20120702_114140b.geny', $
	'residxtalk_20130713_053711.geny','residxtalk_20140628_044121.geny', $
	'residxtalk_20160524_014233.geny','residxtalk_20170716_094445.geny', $
	'residxtalk_20180729_112219.geny','residxtalk_20190616_024408.geny' $
	]
;  residual crosstalk will be read in below once year of data being
;  processed is established

; read in slit scan vignetting function: cubic fit coefficients
vignet_func_files=file_search('$SOT_SP_CALIBRATION','slitvig*')
vignetfunc=last_nelem(vignet_func_files) 
restgenx,file=vignetfunc,slitvig


; read polarization coeffs
polar_coeff_files=file_search('$SOT_SP_CALIBRATION','cal_fit_coef_sb*')
pcoef=last_nelem(polar_coeff_files) ; may want different logic in future?
; ; ie, closest or most recent?
restgenx,file=pcoef,calpos,skx

;;;;;;;;;;;;;;GAIN DARK TABLES;;;;;;;;;;;;;;;;;;;;;;

;  determine the spectral ROI for the present data set
;  First read in the data and the FITS header
        hdr=headfits(sp4dfiles(0))
;  Extract some header parameters relating to CCD region of interest
;  spectral pixel range
        spc0 = sxpar(hdr,'SPCCDIY0')
        spc1 = sxpar(hdr,'SPCCDIY1')
;  Extract the date of the observations
	dateobs = ' '
	dateobs = sxpar(hdr,'DATE_OBS')
	year = strmid(dateobs,0,4)
	fxyear = fix(year)
	month = strmid(dateobs,5,2)

;  recover the gain, dark data based on year of observation
recovergndk_sbsp,fxyear


;  for non-standard ROI, use the eclipse season full ROI darks, flats
if spc0 ne 56 or spc1 ne 167 then begin
	print,'using SP non-standard ROI flatfield data'
	s_mflat = sr_mflat(spc0:spc1,*,*)
	s_shftsl = sr_shftsl(spc0:spc1,*)
	s_spcrfit = sr_specrfit
	s_specshft = sr_specshft
	s_svrr = sr_svrr
	f_mflat = fr_mflat(spc0:spc1,*,*)
	f_shftsl = fr_shftsl(spc0:spc1,*)
	f_spcrfit = fr_specrfit
	f_specshft = fr_specshft
	f_svrr = fr_svrr
	darkavg_std = full_darkavg_std(spc0:spc1,*,*)
	darkavg_fast = full_darkavg_fast(spc0:spc1,*,*)
endif


;------ evaluate drift -------------------------------------------------------
if not keyword_set(noinit) then begin	; skip the first pass for noinit
; first pass through data to get the estimates of drift along the slit
; get list of path names for this directory
nfile = n_elements(sp4dfiles)
print,'number of files to process = ',nfile
;  resort the files in chronological order for the rare instances where
;  the file list sp4dfiles is not.  Use header times to sort
;  File clock time, take some liberty in terms of fractional years
;  for the purose of sorting: assume months all have 31 days
dbltime = dblarr(nfile)

;--- By RCE: define ssm here, instead of below. Fill with values from
;    fits headers in the following loop (rather than inside
;    thermd_sbsp.pro). Define a corrupt_header flag
corrupt_hdr_flag = boolarr(nfile)
ssm = intarr(nfile) ; CAMSSUM value, summing along the slit (1 or 2)
;--- END By RCE


for kk = 0,nfile-1 do begin
;  read the FITS header
	hdr=headfits(sp4dfiles(kk))
	dateobs = sxpar(hdr,'DATE_OBS')
        ayr = double(strmid(dateobs,0,4))

       ; --- By RCE:  Get camera parameters to scan for corrupt header values
        camdaca = SXPAR(hdr, 'CAMDACA') ; typical parameters with corrupt values
        camdacb = SXPAR(hdr, 'CAMDACB') 
        camgain = SXPAR(hdr, 'CAMGAIN')
        camamp = SXPAR(hdr, 'CAMAMP')
                                ; CAMSSUM often time takes value 1
                                ; when it should take value 2. This
                                ; breaks in thermd_sbsp.pro
        ssm[kk] = SXPAR(hdr, 'CAMSSUM') ; typical parameters with corrupt values
        IF NOT(camdaca EQ 5 or camdaca EQ 7) THEN corrupt_hdr_flag[kk] = 1
        IF NOT(camdacb EQ 5 or camdacb EQ 7) THEN corrupt_hdr_flag[kk] = 1
        IF (camamp NE 1) THEN corrupt_hdr_flag[kk] = 1
        IF (camgain NE 2) THEN corrupt_hdr_flag[kk] = 1
        ;--- END By RCE
        
;  save year of observation for first file
	if kk eq 0 then yearobs = ayr
	amnth = double(strmid(dateobs,5,2))
	aday = double(strmid(dateobs,5,2))
	strtime = sxpar(hdr,'TIME-OBS')
	ahr = double(strmid(strtime,0,2))
	amn = double(strmid(strtime,3,2))
	asec = double(strmid(strtime,6,6))
;  now compute fractional years (FOR SORTING ONLY, NOT EXACT!)
;  This fudge presumes each month has 31 days
	dbltime(kk) = ayr + double(amnth/12.) + double(aday/(31.*12.))  $
		+ double(ahr + amn/60. + asec/3600.)/double(24.*31.*12.)
endfor
;  sort the files in terms of increasing time
isrt = sort(dbltime)
sp4dfiles = sp4dfiles(isrt)

;--- By RCE: sort the new variables created to search for corrupt headers
corrupt_hdr_flag = corrupt_hdr_flag(isrt) ; By RCE: sort corrupt flag
ssm = ssm(isrt)
;--- END By RCE

;  with year of data now established, select and read in the residual 
;  crosstalk data based on the year of this observation
nyresid = n_elements(yresid)
for kyr = 0,nyresid-1 do begin
	deltyear = abs(yearobs - yresid)
	mmnn = min(deltyear,iyr)
endfor
xtalk_file=file_search('$SOT_SP_CALIBRATION',xtalkfiles(iyr))
;  do this for 2007, 2009, 2010, 2011, 2012 as well
;if yearobs lt 2007. then restgenx,file=xtalk_file(0),scoef,degn,sltpos,nftot
;if yearobs ge 2007. then restgenx,file=xtalk_file(0),xtalkimg,sltpos,nftot
;  can now use new processing for all years since calibration files have
;  been updated
restgenx,file=xtalk_file(0),xtalkimg,sltpos,nftot


; set some arrays for saving the thermal drift in the slit direction
wdelw = fltarr(nfile) ; result of cross-correlation
ftime = fltarr(nfile) ; time in fractional hours from first file in directory

;--- By RCE: define ssm above, prior to reading headers, and fill in
;    the values
;ssm = intarr(nfile) ; CAMSSUM value, summing along the slit (1 or 2)
; --- END By RCE
fitww = fltarr(nfile) ; Smoothed thermal drift in the slit direction
fitsp = fltarr(nfile) ; Smoothed thermal drift in the spectral direction
sltdr = fltarr(nfile) ; Ichimoto prediction of thermal drift along slit
spcdr = fltarr(nfile) ; Ichimoto/Kubo prediction of thermal spectral drift
dthrm = fltarr(nfile) ; total of orbital and thermal spectral drift, pixels
dopcv = fltarr(nfile) ; orbital Doppler shift DOP_CVR, pixels
avctr = fltarr(nfile) ; slit-average min. intensity position of 6301.5
fitav = fltarr(nfile) ; smoothed slit-average min. intensity position of 6301.5
; ADJUSTED SO AS TO PLACE THE AVERAGE LINE CENTER AT FINAL
; SPECTRAL PIXEL 29

kntr = 0


; FIRST DERIVE THE THERMAL DRIFT ALONG THE SLIT USING CROSS-CORRELATION
; first pass through data: derive estimates of thermal drift along slit
; loop over files in the directory
erfiles = [-1]
stksi_ok = bytarr(nfile)
; file index of 1st good frame: determined in thermd_sbsp
	firstgood = -1
;  initialize data for thermal drifts
	avgctr=0.
;  begin loop over files in map
for ks = 0,nfile-1 do begin
	erflag = 0
        if keyword_set(verbose) then print,' starting file ',ks
        fits_input = sp4dfiles(ks)
        ;  set previous values
        if kntr ge 1 then begin
           wdelw_prev = wdelw(kntr-1)
           avgctr_prev = avgctr
        endif
        
        ;--- By RCE: if camera header values are corrupt, set error flag
        ; and check that CAMSSUM has correct value
        IF corrupt_hdr_flag[ks] THEN BEGIN
           erflag = 1
                                ; Temporarily correct the values of
                                ; ssm for the files with corrupt headers
           IF ks GT 0 THEN BEGIN ; if the corrupt file is not the first one in the map
              if ssm[ks] NE ssm[ks-1] THEN ssm[ks]=ssm[ks-1]
           ENDIF ELSE BEGIN ; if the corrupt file is the first one in the map
              if ssm[ks] NE ssm[ks+1] THEN ssm[ks]=ssm[ks+1]
           ENDELSE
        ENDIF
                                ; bad files are
                                ; blacklisted automatically by
                                ; checkdropout_sbsp.pro
        
        ;--- END By RCE

	thermd_sbsp, fits_input, $
	wdelw,kntr,ssm,slitdrft,specdrft,doprcv,ftime, $
	avgctr,erflag, $
    	darkavg_std,darkavg_fast, $
    	s_mflat,s_shftsl,s_specrfit,s_specshft,s_svrr, $
    	f_mflat,f_shftsl,f_specrfit,f_specshft,f_svrr, $
   	wdelw_prev,avgctr_prev,tryhard,firstgood,sidebad,stksiok, $
	_extra=_extra, verbose=verbose

;  check for single-sided data, adjust the size of side_bad accordingly
;  set packet loss save arrays
	if kntr eq 0 then begin
		nsides = sizeof_sbsp(sidebad,1)
		side_bad = fltarr(nsides,4,nfile)
		stksi_ok = bytarr(nfile)
	endif
;  store data on bad packets
	stksi_ok(ks) = stksiok
	side_bad(*,*,ks) = sidebad

; load the projected thermal drifts into output arrays
	sltdr(ks) = slitdrft
	spcdr(ks) = specdrft
	dopcv(ks) = doprcv
	dthrm(ks) = doprcv+specdrft
	avctr(ks) = avgctr
; increment the file counter
	kntr = kntr+1
	if erflag ne 0 then erfiles = [erfiles,ks]
endfor

;  reset the array erfiles to only those bad files, then remove them
;  from the processing
nerr = n_elements(erfiles)
if nerr gt 1 then begin
	erfiles = erfiles(1:*)
	nerr = nerr-1
endif else begin
	nerr = 0
endelse
;  set array of bad files for pass2 processing
badindx = intarr(nfile) & badindx(*) = 0 
if nerr gt 0 then begin
	badindx([erfiles]) = 1
	goodfiles = indgen(nfile)
	goodfiles([erfiles]) = -10
	whr = where(goodfiles ge 0)
	goodfiles = goodfiles(whr)
;  reset data from thermd to exclude error files
	wdelw = wdelw(goodfiles)
	ftime = ftime(goodfiles)
	sltdr = sltdr(goodfiles)
	spcdr = spcdr(goodfiles)
	dopcv = dopcv(goodfiles)
	dthrm = dthrm(goodfiles)
	avctr = avctr(goodfiles)
	fitww = fitww(goodfiles)
	fitsp = fitsp(goodfiles)
	fitav = fitav(goodfiles)
endif


;  if error flag is set, make note of bad files skipped over
if nerr gt 0 then begin
	for ks=0,nerr-1 do print, 'error encountered in file ', $
		sp4dfiles[erfiles[ks]]
endif

;BWL- updated drift conditioning along slit
; CONDITION THE DERIVED DRIFT ALONG THE SLIT
; set the time from start of sequence
ftime = ftime - ftime(0)
; smooth the output drifts with a median filter to remove some outlying
; data values. Ignore portions of map where observations are entirely off
; the limb, as those points are set to -1000 in thermd_sbsp
;whrdsk = where(wdelw gt -100.,complement=whrout,whcount)
;  get points off disk in both line center, intensity variables)
whrdsk = where(wdelw gt -100. and avctr gt 0.,complement=whrout,whcount)
if whcount gt 0 then begin
        wdelon = wdelw(whrdsk)
endif else begin
        print, ' error, found no on-disk data in this map'
        stop
endelse

;  perform median filtering for requisite number of points
tempf = wdelon
npser = n_elements(wdelon)
nptop = mfltrang
if npser lt mfltrang then nptop = npser
if (npser ge 2) then medfilt_sbsp,wdelon,nptop,1,tempf
; smooth the drifts along the slit
fitww(*) = 0.
xxx = ftime(whrdsk)
yyy = tempf  ; modified 26-jun-2007 per B.L. 

; get mean time in hours between the observations
nwdisk = n_elements(whrdsk)
dtt = (xxx(nwdisk-1)-xxx(0))/float(nwdisk)

;  smoothing should be done on a fixed time interval (4 min or less),
;  not a fixed number of points.  The given value is typical of a 
;  normal map
nsmoot = fix(nsmoothr*4.8/(3600.*dtt))

;  temporary fix for rapid changes: alters smoothing for thermal
;  drifts along slit only
nsmoot5 = nsmoot
rpdthrm = keyword_set(rapidtherm)

if rpdthrm then nsmoot5 = nsmoot/5

; new smoothing routine for drifts along the slit taking into
; account gaps in the data
gapsmth_sbsp,nsmoot5,xxx,yyy,fitwh
fitww(whrdsk) = fitwh
ftime(whrdsk) = xxx

nsmthw = min([nsmoot5,whcount/2])

; set the off-limb values to the mean of the on-limb values
if whrout(0) ge 0 then begin
        mnwhr = mean_sbsp(fitww(whrdsk))
        wdelw(whrout) = mnwhr
        fitww(whrout) = mnwhr
endif
dispflag = keyword_set(display)
if dispflag then begin
	plot,ftime,wdelw,psym=3,title= $
	'history of thermal drift of image along slit, smoothed fit'
	oplot,ftime,fitww
endif




fitsp(*) = 0.
yyy = spcdr(whrdsk)
;  use gapsmth function for data with gaps
gapsmth_sbsp,nsmoot,xxx,yyy,fitaaa
fitsp(whrdsk) = fitaaa
; smooth the empirical spectral variation too
; line center wavelength pixel for points off the limb set to standard
; value of pixel 29 (params_sbsp value stdlnctr, but fitav is value
; of shift relative to this standard pixel center)
fitav(*) = 0.
; median filter time series to purge spurious values. These occur when
; there are e.g. dropout frames.  
;  Use disk determination from wdelw above, so comment the next line
;whrdsk = where(avctr gt 0.,complement=whrout,whcount)
yyy = avctr(whrdsk)

;  perform median filtering for requisite number of points
tempav = yyy
npdsk = n_elements(yyy)
nptop = mfltrang
if npdsk lt mfltrang then nptop = npdsk
if (npdsk ge 2) then medfilt_sbsp,yyy,nptop,1,tempav
;  use gapsmth function for data with gaps
gapsmth_sbsp,nsmoot,xxx,tempav,fitaaa
tempav = fitaaa
; fitav = the pixels to shift the data in order for average pixel of
; 6301.5 to be at pixel 29.0 in final data.
; Apply no shifts for points substantially outside
; the limb. There might be some artificial discontinuities at the very
; limb as a result. Oh well....
if whcount gt 0 then begin
        fitav(whrdsk) = stdlnctr-tempav
	if whrout(0) ge 0 then begin
		fitav(whrout) = 0.
		avctr(whrout) = stdlnctr
	endif
endif else begin
        fitav = stdlnctr-tempav
endelse



; SMOOTH THE PREDICTED THERMAL DRIFT IN SPECTRAL, SPATIAL DIRECTIONS
; the thermal drift shows the noise in the temperature measurements, which
; cause variations in the line center position. Smooth these the same as
; the inferred thermal drift parameters
; output the thermal drift parameters:
; wdelw = result of cross-correlation
; ftime = time in fractional hours from first file in directory
; tempf = median-filtered, discontinuity-scrubbed cross-correlation result
; sltdr = Ichimoto prediction of thermal drift along slit
; spcdr = Ichimoto/Kubo prediction of thermal drift in spectral direction
; fitsp = smoothed predicted spectral drift
; dopcv = orbital Doppler shift prediction from DOP_CVR, in pixels
; dthrm = total of orbital and thermal spectral drift, pixels
; avctr = slit-average min. intensity position of 6301.5
; fitav = smoothed slit-average min. intensity position of 6301.5
; badindx = array of indices of bad files
; ADJUSTED SO AS TO PLACE THE AVERAGE LINE CENTER AT FINAL
; SPECTRAL PIXEL 29
;save these results to an IDL save file
save,filename=concat_dir(outdir,'thermdrift.save'),wdelw,ftime,tempf,fitww,sltdr, $
spcdr,dopcv,dthrm,fitsp,avctr,fitav,badindx,side_bad,stksi_ok
;print,'!!!EXAMINE THERMAL DRIFT HISTORY FOR ANOMALIES!!!'

; if noinit is set, then skipping pass1 to here
endif else restore,concat_dir(outdir,'thermdrift.save')

if keyword_set(pause) then begin ; 2007.1.27 k.i.
print,'drift fit complete...'
print,'type .c to continue'
stop
endif


;***** pass-2 *********************************************************
nsp=n_elements(sp4dfiles)

; this variable set to 2 identifies when F-P exception occurs
;!except=2

if nsp eq 0 then begin
   box_message,['Need SOT/SP4D IQUV file input...,returning']
   return
endif

fexist=total(file_exist(sp4dfiles))
if nsp ne fexist then begin
   box_message,['Not all Leve0 SP files are found, returning...']
   return
endif


;BWL- START LOOP OVER FILES TO BE PROCESSED HERE
kntr = 0
nfile=n_elements(sp4dfiles)
for ks = 0,nfile-1 do begin

;  bypass processing if file has error found in thermd_sbsp
	if badindx(ks) ne 0 then goto, skipass2

;fits_input = derf(ks)
fits_input = sp4dfiles(ks) ;; k.i. 2007.1.18

; Call Calibration Suite (w/ssw mods...)
calib_sbsp, fits_input,l1index,l1data, $
fitww,fitav,kntr, $
    scoef,degn,sltpos,nftot, $
    slitvig, $
    calpos,skx, $
    darkavg_std,darkavg_fast, $
    s_mflat,s_shftsl,s_specrfit,s_specshft,s_svrr, $
    f_mflat,f_shftsl,f_specrfit,f_specshft,f_svrr, $
    wlshft,slshft, $
    tryhard,rsdsh_prev,anorm_prev,side_bad(*,*,ks), $
    stksi_ok(ks), xtalkimg, yearobs, $
display=display,_extra=_extra, verbose=verbose
;_extra=_extra


;BWL- ADD NEW KEYWORDS FOR LEVEL1 DATA HERE:
;BWL- SPSLSHFT = slshft = PIXEL SHIFT APPLIED ALONG SLIT
;BWL- SPWLSHFT = wlshft = PIXEL SHIFT APPLIED IN SPECTRAL DIRECTION
;BWL- FGSPROT = fgsprot = rotation angle clockwise from from N of FG pixel column to SP slit
l1index=boost_tag(l1index,wlshft,'SPWLSHFT') ; wl shift after smoothing in pix, KI 2007.1.24
l1index=boost_tag(l1index,slshft,'SPSLSHFT') ; drift correction
l1index=boost_tag(l1index,dthrm(kntr),'SPWLSFT0') ; KI 2007.1.27
l1index=boost_tag(l1index,wdelw(kntr),'SPSLSFT0')

;BWL- 3-dec-2007 CHANGE AND ADD HEADER VALUES FOR POINTING - FROM COALIGNMENT WITH G-BAND
;BWL- SEE WRITEUP BY R. CENTENO ELLIOTT
;BWL- FGSPROT = fgsprot = rotation angle clockwise from from N of FG pixel column to SP slit
;BWL- fgsprot is hardwired to value determined by R. Centeno Elliott = 0.2636 degrees
fgsprot = 0.2636
l1index=boost_tag(l1index,fgsprot,'FGSPROT') ; rotation of SP slit relative to G-band
;BWL- update the header values
;SPscan = l1index.XSCALE
SPscan = 0.1476
; New values for XSCALE and YSCALE, scaled appropriately for summing along slit
; and in scan direction
l1index.XSCALE = 0.14857*l1index.scn_sum
l1index.YSCALE = 0.15999*l1index.camssum

; Read SLITPOS, and original positioning SC_ATTX and SC_ATTY
slpos = l1index.SLITPOS
XSUN =  l1index.SC_ATTX
YSUN =  l1index.SC_ATTY

; Compute new values for XCEN and YCEN according to the corrections we
; found
;---- XCEN:
a = -4.33d
b = 649.8d
c = 0.551d
d = 0.0066d
xoff = -34.57d
SPXCEN = XSUN + SPscan * (a * cos(2d0*!PI/b*slpos+c) $
        + (1d0+d) * slpos + xoff)
;---- YCEN 
a = -0.00234
yoff = 49.35
SPYCEN = YSUN + SPscan * (a * slpos + yoff)
;  correct ycen for the offset cutout along the slit, using the current
;  pixel size along the slit
SPYCEN = SPYCEN +  $
(l1index.naxis2/2+0.5-l1index.crpix2)*l1index.YSCALE
; change header values
l1index.XCEN = SPXCEN
l1index.YCEN = SPYCEN


;  update the wavelength dispersion to be positive, update CRVAL1
cdelt1 = abs(l1index.CDELT1)
l1index.CDELT1 = cdelt1
spccdiy0 = l1index.SPCCDIY0
spccdiy1 = l1index.SPCCDIY1
;  reference pixel for 6301
crval1 = 0.5*(float(spccdiy0)+float(spccdiy1))
;  for usual spectral ROI, the 6301.5091 line is forced to fall at 
;  pixel 29 in the calibrated (wavelength-reversed, so red is 
;  higher pixel values) data.  Spectral ROI for the raw data
;  is usually [56,167], thus 6301.5091 falls at (167-29) = 138.
crval1 = 6301.5091 + cdelt1*(138. - crval1)
l1index.CRVAL1 = crval1









update_history,l1index,/caller,version=version ; add sp_prep version

;BWL- UPDATE FITS HEADER HERE

; --- fits output, 2007.01.19 K.I.
;  Change SP4D to SP3D, insert C in file name
p1=strpos(fits_input,'/',/reverse_search)
fname=strmid(fits_input,p1+1,strlen(fits_input)-p1)
p2 = strpos(fname,'SP4D')
if p2 gt 0 then begin
	fnpre = strmid(fname,0,p2)
	fname = fnpre + 'SP3D' + strmid(fname,p2+4)
endif else begin
	fname = 'SP3D' + strmid(fname,p2+4)
endelse
fname=strmid(fname,0,strpos(fname,'.fits'))
outfile=concat_dir(outdir,fname+'C.fits')
l1hdr = struct2fitshead(l1index)
fits_write,outfile, l1data, l1hdr
dat=readfits(fits_input,exhdr,/silent,exten=1)
writefits,outfile,dat,exhdr,/append

; increment the file counter
kntr = kntr+1


;BWL- CLOSE LOOP OVER FILES TO BE PROCESSED
skipass2:
endfor

;  run the stksimages_sbsp program
stksimages_sbsp,outdir,outdir=outdir

end