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,, 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