pro flatgen3_sbsp,input_dir,xctr,yctr,darkpath,goodfiles,roistart,skewname ; ;+ ; Name: flatgen3_sbsp ; ; ; Purpose: Generate flat-field data including flat images, variation of the ; intensity along the slit, spectral curvature, skew, spectral offsets, etc ; from the average spectra of a sequence of quiet Sun maps taken during ; eclipse season with attendant darks during S/C night ; ; The selection of valid flat field files is done external to this routine ; ; This routine does four processes in sequence to arrive at a HINODE SP flat ; correction. It was modified from flatgen2_sbsp.pro. ; ; The prior version of flatgen_sbsp.pro skipped PROCESS1 and PROCESS2 below ; since the data had already been calibrated and thermal drift files ; (thermdrift.save) were available. For non-standard spectral ROI positions, ; standard calibration is not yet available, so the full procedure must be ; carried out. This procedure should also work for standard ROI data acquired ; from here on out. ; ; PROCESS 1: for each data directory, dark corrects the data, extracts ; the variation of the continuum intensity along the slit for each file, ; uses pre-flight slit variation to determine the thermal drift along the ; slit for each file, smooths the slit variation in time, shifts and ; averages the continuum intensity variations along the slit for all ; files in the directory, saves each the continuum intensity variations ; along the slit for each CCDSIDE and file in an IDL-save file ; This procedure is primarily used to extract thermal drifts. Older ; versions of the slit variation are still used in calibration. ; ; PROCESS 2: average and normalize the continuum intensity variations ; over all directories ; ; PROCESS 3: pass through all the data again to remove the slit intensity ; variations from each frame, then average the frames to get an average ; flat image corrected for the slit intensity variations. This requires ; shifting the mean slit intensity pattern along the slit according to the ; time-smoothed thermal drift, and correcting for spectral skew of the ; slit intensity pattern ; ; PROCESS 4: call routine genflt3_sbsp.pro (modified from genflt_sbsp.pro ; to input a specified name for output save file) to generate the gain correction ; data. This passes the slit intensity variations through to be written ; out as ancellary data accompanying the gain correction images. ; This last process, genflt3_sbsp.pro, does all the output ; ; Inputs: ; input_dir = array of paths to directories containing the disk center ; observations to be averaged and processed. Note that ; it is best to have directories separated for each map ; so that there are no large time jumps in the data sets. ; This is because the thermal drifts will be smoothed. ; xctr,yctr = [x,y] centers of map. This is used in calculation of the ; limb darkening function for each map. These are ; input in the input_dir ascii file and passed to this ; routine via the run macro ; darkpath = path to idl save file containing averaged darks taken at ; a time close to that of the flat observations ; goodfiles = name of IDL save file containing the list of good file names ; roistart = starting pixel of data ; ; Output: (See genflt3_sbsp.pro) ; skewname = pathname to most recent output of skew parameters for ROI56 ; this will resemble shftavg_20100516_1721_roi56.save ; ; History: ; 6-Jan-2007 - Bruce Lites, HAO. ; 25-Jan-2008 - Bruce Lites, HAO: use concat_dir to construct restored ; path names ; 3-Dec-2009 - Bruce Lites, HAO, adapted for individual flat/dark observations ; during eclipse season ;- ; version=1.0 ; first specify the base for output file names ; get the string representation of the start pixel of the ROI rstrt = strcompress(string(roistart),/remove_all) ; get the time ID for the filename ppos = strpos(input_dir,'/20') plgth = strlen(input_dir)-ppos - 2 timid = strmid(input_dir,ppos+1,plgth) ; generate the output filename base for the selected files for averaging outfilbase = '_' + timid + '_roi' + rstrt display=keyword_set(display) or n_params() eq 1 ; include all the seldom-changed fixed parameters scal=get_logenv('SOT_SP_CALIBRATION') if scal eq '' then set_logenv,'SOT_SP_CALIBRATION',$ concat_dir('$SSW_SOT','calibration') params_file=file_search(get_logenv('$SOT_SP_CALIBRATION'),'params_sbsp*.com') calparams_file=file_search(get_logenv('$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??' ; Now get updated one for fast/normal map, comment this out ;; CALIBRATION DATA THAT IS INDEPENDENT OF SLIT SCAN POSITION ;; read standard slit intensity variation ; open slit variation files slit_inten_var_files=file_search(get_logenv('$SOT_SP_CALIBRATION'),'sltvar*') slitinten=last_nelem(slit_inten_var_files) restgenx,file=slitinten,sltvar ; get the stored dark images ; Below is for the special processing for selected spectral ROI ; stored as a standard IDL-save file. print,' Using special spectral ROI dark data: ',darkpath restore,darkpath ; APPLY THE SLIT SCAN VIGNETTING FUNCTION TO ; THE FLATS FIRST TO TAKE CARE OF RESIDUAL SLOW VARIATIONS ALONG THE ; SLIT. WILL ALSO TAKE OUT LIMB-DARKENING ; new 2-D slit scan vignetting function vignet_func_files=file_search(get_logenv('$SOT_SP_CALIBRATION'),'slitvig*') vignetfunc=last_nelem(vignet_func_files) restgenx,file=vignetfunc,slitvig ; PROCESS 1: GET THE AVERAGE VARIATION OF INTENSITY ALONG THE SLIT FOR EACH ; CCDSIDE AND WRITE OUT FOR EACH DIRECTORY SPECIFIED ABOVE ndir = 1 for idir = 0,ndir-1 do begin ; find the file names in the current directory ; previously determined good files restore,goodfiles nfile = n_elements(derf) ftime = fltarr(nfile) wdelw = fltarr(nfile,2) fitww = fltarr(nfile,2) ; for fast map must bypass the slitshftnp procssing datt = readfits(derf(0),hdr) ssum = sxpar(hdr,'CAMSSUM') if ssum eq 2 then goto,jump443 ; For standard ROI runs (ROI + 56) determine the spectral skew ; for each of the raw images, average the results, then save them ; for output for the flat field data ; This ROI must be done first, and results saved if rstrt eq '56' then begin ; begin loop over files in this directory for kntr = 0,nfile-1 do begin sp4fits = derf(kntr) slitshftnp_sbsp,sp4fits,cc1,cc2,sp1,sp2,sl1,sl2,ssum, $ ssl1,ssl2,nsearch,shftsl,delsh ; set some arrays if kntr eq 0 then begin shftavg = shftsl shftavg(*,*) = 0. delshavg = 0. charr = ' ' nxx = sizeof_sbsp(shftsl,1) endif ; average the slit variations shftavg = shftavg + shftsl delshavg = delshavg + delsh endfor ; renormalize shftavg = shftavg/float(nfile) delshavg = delshavg/float(nfile) ; save temporary files with the good shifts skewname = 'shftavg'+outfilbase+'.save' save,filename=skewname,shftavg,delshavg,ssl1,ssl2,sl1,sl2 ; end procesing of skew results endif jump443: ; begin loop over files in this directory for kntr = 0,nfile-1 do begin sp4fits = derf(kntr) ; First read in the data and the FITS header print,' Begin processing file: ',sp4fits hdr=headfits(sp4fits) read_sot,sp4fits,index,dat dat = float(dat) ; define arrays etc if kntr eq 0 then begin ; IMAGE DIMENSIONS ; x is slit dimension, y is spectral dimension ; Fortuitously, the X-direction in the naming convention for SP header ; keywords is in the dimension along the slit, and the Y-dimension is ; in the spectral direction. Instead of retaining that convention, which as ; it happens was the naming convention for DLSP reduction software ; also, to avoid confusion we use nw as the number of pixels in the wavelength (Y) ; dimension, and nsl as the number of pixels along the slit. nw = sizeof_sbsp(dat,1) & nsl = sizeof_sbsp(dat,2) nsl1 = nsl-1 & nwm1 = nw-1 nsl2 = 2*nsl ; twice the slit length for apodizing ; continuum intensity vs slit length, each side avgcont = fltarr(nsl,2,nfile) ; arrays independent of directory if idir eq 0 then begin ; define intermediate arrays (spectral,spatial,stokes) sttmp = fltarr(nw,nsl,4) stt = sttmp stt_zero=stt ; use inside loop instead of stt(*,*,*)=0 contin = fltarr(nsl) ; polynomial degree, temporary array for fitting residual shifts xx = findgen(nsl) ; continuum intensity vs slit length, each side contm = fltarr(nsl,2) ; array for averaged observed slit variation for each directory svr = fltarr(nsl,2) svv = fltarr(nsl) ; gain correction image for slit width variations sltgn = fltarr(nw,nsl) endif endif ; EXTRACT PARAMETERS FROM HEADER ; Extract the number of pixels binned along the slit (along serial direction) ssum = sxpar(hdr,'CAMSSUM') ; Extract some header parameters relating to CCD region of interest ; spectral pixel range spccdiy0 = sxpar(hdr,'SPCCDIY0') spccdiy1 = sxpar(hdr,'SPCCDIY1') ; TEMP KLUDGE TO CORRECT THE OCCASIONAL ERRONEOUS VALUE NEAR HOUR CHANGE if spccdiy1 ne (spccdiy0 + 112 -1) then spccdiy1 = spccdiy0+112-1 nnw = spccdiy1-spccdiy0+1 ; check to see that image size does not change during this map if nnw ne nw then begin print,'spectral array size (',nnw,') not consistent with first file' stop endif ; Range of pixels along the slit ; Note that the pixel indices along the slit direction are reversed from ; those during pre-launch testing. The images are now reversed in ; the vertical direction during reformatting. Also, pixel indices in ; the header like SPCCDIX0, SPCCDIX1 refer to the pre-launch CCD coordinates, ; not the coordinates of LEVEL0 data. spccdix0 = sxpar(hdr,'SPCCDIX0') spccdix1 = sxpar(hdr,'SPCCDIX1') ; define these pixel positions in terms of the LEVEL0 image coordinates ; reverse order so that low limit of active range is in spccd0 spccd0 = 1023-spccdix1 spccd1 = 1023-spccdix0 ; bit shifting for SP data bitshft = sxpar(hdr,'SPBSHFT') ; extract the index of the slit scan position slpos = sxpar(hdr,'SLITPOS') ; File clock time strtime = sxpar(hdr,'TIME-OBS') ahr = float(strmid(strtime,0,2)) amn = float(strmid(strtime,3,2)) asec = float(strmid(strtime,6,6)) ftime(kntr) = ahr + amn/60. + asec/3600. ; temp print times print,kntr,'hour,min,sec',ahr,amn,asec ; account for possibility of change across the day boundary if kntr gt 0 then begin if ahr lt fix(ftime(kntr-1)) then $ ftime(kntr) = ftime(kntr) + 24. endif ; ADJUST SOME LIMITS IF SUMMING, OR CUTOUT ALONG SLIT, OR BOTH ; if summing along slit, adjust appropriately sltvr = sltvar ; !!!TEMP for now,just rebin the slit variation if summing, later will ; get slit variation from fast map flats if ssum eq 2 then sltvr = rebin(sltvar,512) ; make a gain correction image for slit width variations ; also adjust for cutouts along slit xcl = spccd0/ssum & xcu = spccd1/ssum sltvrc = sltvr(xcl:xcu) for ii = 0,nwm1 do sltgn(ii,*) = mean_sbsp(sltvrc)/sltvrc ; Derive active range along slit for shifting in recorded pixels, not CCD pixels ssl1 = max([sl1(1),spccd0]) - spccd0 ; always truncate the limit in the case of summing along slit ssl2 = (min([sl2(1),spccd1]) - spccd0)/ssum if ssum eq 2 then begin if (ssl1 mod 2) eq 1 then begin ssl1 = (ssl1+1)/ssum endif else begin ssl1 = ssl1/ssum endelse endif xxr = xx(ssl1:ssl2) ; pixel range for resid. spectral shift corr. ; (see Operation 7 below) ; check to see if data dimensions are those expected from the header if nw ne (spccdiy1-spccdiy0+1) then begin print,'Spectral data dimension (',nw,') not equal to ', $ 'SPCCDIY1-SPCCDIY0+1: (',spccdiy1-spccdiy0+1,')' stop endif if nsl ne (spccdix1-spccdix0+1)/ssum then begin print,'Dimension along slit (',nsl,') not equal to ', $ '(SPCCDIX1-SPCCDIX0+1)/CAMSSUM: (',(spccdix1-spccdix0+1)/ssum,')' stop endif ; EXTRACT DARK FIELDS APPROPRIATE FOR THIS FILE ; since flats should always be taken in a standard mode (normal, fast, ; dynamics), and since we read in a specialized dark for these flat data, ; just use the dark average data as-is ; dark = darkavg_std ; for specific dark/flat eclipse season data, just use the standard darks dark = darkavg ; if ssum eq 2 then dark = darkavg_fast ; DIGITAL WRAP, BIT SHIFTING CORRECTIONS ; correct for digital wrap-around of Stokes I, if wrap-around is present temp = dat(*,*,*,0) whr = where(temp lt 0., countwrap) if countwrap ne 0 then temp(whr) = temp(whr) + 65536. dat(0,0,0,0) = temp ; account for bit shifting of the data by multiplying the appropriate ; Stokes images by 2 switch bitshft of 3: dat(*,*,*,1:2) = 2.*dat(*,*,*,1:2) 2: dat(*,*,*,3) = 2.*dat(*,*,*,3) 1: dat(*,*,*,0) = 2.*dat(*,*,*,0) endswitch ;BWL- !!!THIS IS A TEMPORARY FIX FOR ERRORS IN THE 1ST FILE HEADER ;BWL- !!!THIS SHOULD BE REMOVED WHEN REFORMATTER FIXES THIS PROBLEM if spccdiy1 ne (spccdiy0 + 112 -1) then spccdiy1 = spccdiy0+112-1 ; FIND APPROXIMATE THERMAL FLEXURE SHIFT FROM VIGNETTING ; AT TOP OF IMAGE OF CCDSIDE1 IF THE DATA SAMPLE THAT REGION vigsh = -2000 if spccd1 eq 1023 then vigshift_sbsp,dat(*,*,1,0),ssum,vigsh if vigsh gt -2000 and ssum eq 2 then vigsh = vigsh - 1 ; !!!FOR ROI0, FIRST TWO COLUMNS BAD, REPLACE WITH THIRD COL ; !!!FOR ROI0, FIRST TWO COLUMNS BAD, REPLACE WITH EXTRAPOLATION if roistart eq 0 then begin ; !!!FOR ROI0, problem with first pixels of row 400 for iside = 0,1 do for istks = 0,3 do begin dat400 = mean(dat(6:10,400,iside,istks)) dat(0:5,400,iside,istks) = dat400 endfor ; for Stokes I, do quadratic fit for iside = 0,1 do begin for iyy = 0,nsl-1 do begin xxx = findgen(10) + 2. yyy = dat(3:12,iyy,iside,0) coeff = poly_fit(xxx,yyy,2,yfit=yfitt) for ii = 0,2 do begin fitt = 0. for kk = 0,2 do begin fitt = fitt + coeff(kk)*float(ii)^kk dat(ii,iyy,iside,0) = fitt endfor endfor ;dat(0:2,iyy,iside,0) = yfitt endfor endfor ; for QUV, do simple average ;for istks=0,3 do begin for istks=1,3 do begin for iside = 0,1 do begin for iyy = 0,nsl1 do begin for ixx = 0,2 do begin dat(ixx,iyy,iside,istks) = mean(dat(3:10,iyy,iside,istks)) endfor endfor endfor endfor endif ; BEGIN LOOP OVER CCDSIDEs for iside = 0,1 do begin stt=stt_zero ; fill the working array for this CCDSIDE, Stokes I only sttmp(0,0,0) = dat(*,*,iside,0) ; DARK/FLAT CORRECTION stt(0,0,0) = (sttmp(*,*,0) - dark(*,*,iside)) ; replicate the unused spectral columns to avoid discontinuities if msh1 gt 0 then begin repvec = stt(msh1,*,0) for ii = 0,msh1-1 do stt(ii,*,0) = repvec endif if msh2 lt nwm1 then begin repvec = stt(msh2,*,0) for ii = msh2+1,nwm1 do stt(ii,*,0) = repvec endif ; now, replicate the unused spatial rows to avoid discontinuities if ssl1 gt 0 then begin repvec = stt(*,ssl1,0) for jj = 0,ssl1-1 do stt(*,jj,0) = repvec endif if ssl2 lt nsl1 then begin repvec = stt(*,ssl2,0) for jj = ssl2+1,nsl1 do stt(*,jj,0) = repvec endif ; determine the continuum intensity as the average of spectral pixels for jj = 0,nsl1 do begin contm(jj,iside) = total(stt(1:11,jj,0))/11. avgcont(jj,iside,kntr) = contm(jj,iside) endfor ; for flat correction, no flat field as yet. Use just dark corrected ; data sout = stt ; SENSING THERMAL FLEXURE OF STRUCTURE, COMPENSATION ; cross-correlate the slit intensity pattern to find the thermal ; flexure shift along the slit. The intensity pattern is determined ; from the flat-field data and passed through to this routine ; select limits along slit to allow room for apodization. Allow at ; least 20 points total liml = max([ssl1,10]) limu = min([ssl2,nsl-11]) ; !!!this could be done once outside this routine, array is constant ; Ichimoto-san's 2nd derivative for reference variation, exclude ends tvvt = extend_sbsp(sltvrc,liml,limu) d2sl = (2.*tvvt-shift(tvvt,1)-shift(tvvt,-1)) ; get the measured continuum intensity for jj = 0,nsl1 do contin(jj) = mean_sbsp(sout(cc1:cc2,jj,0)) ; check for off-limb measurements. If mean continuum along the less than 1000 DN ; then set the shift to -1000 if mean_sbsp(contin) lt 1000. then begin delw = -1000. goto,jump5 endif ; Ichimoto-san's 2nd derivative for continuum, exclude ends tvvt = extend_sbsp(contin,liml,limu) ; normalize 2nd derivative by the local continuum but don't over-enhance ; the darkest points, hence don't normailze by less than 5000 DN tvvnorm = tvvt>5000. d2cc = (2.*tvvt-shift(tvvt,1)-shift(tvvt,-1))/tvvnorm ; cross-correlate, searching +/- 20 pixels. This should be ample even ; for data early in the mission. Otherwise, we can restrict this lag = findgen(41)-20 cccr = c_correlate(d2cc,d2sl,lag) ; find maximum, which is shift of d2sl wrt d2cc ccmax = max(cccr,iccmx) ; fit a parabola around maximum parabofit,xx(iccmx-1:iccmx+1),cccr(iccmx-1:iccmx+1),delw delw = delw-20 jump5: wdelw(kntr,iside) = delw idelw = round(delw) ; end loop over CCDSIDEs endfor ; end loop over files in directory endfor ; get time from start of sequence ftime = ftime - ftime(0) ; smooth the thermal drift ; smooth the output drifts for iside = 0,1 do begin ; median filter the curve to get rid of spurious points temp = wdelw(*,iside) medfilt_sbsp,temp,7,2,tempf ; go through again looking for discontinuities for kntr = 1,nfile-1 do begin if (abs(tempf(kntr-1) - tempf(kntr)) gt 3. and $ (ftime(kntr) - ftime(kntr-1)) lt 0.0333) then begin tempf(kntr) = tempf(kntr-1) endif endfor nsmthw = min([50,nfile/2]) ; use smoothing instead of polynomial fit fitws = smooth(tempf,nsmthw,/edge_truncate) fitww(*,iside) = fitws wdelw(*,iside) = tempf ; detect the transition to vignetting at top of image, and set the continuum ; intensity to its median value along the slit above that transition ; get array of top of image ; hardwire nvigtop to 36 to search for vignetting at top nvigtop = 36 sobtop = fltarr(nfile,nvigtop) for ii = 0,nfile-1 do sobtop(ii,*) = avgcont(nsl-nvigtop:nsl1,iside,ii) sobtop = sobel(sobtop) ; sobel filter the result medsob = median(sobtop) ; find points where sobel filtered image is greater than 5 times its median whr = where(sobtop gt 4.0*medsob,kcsob) sobtop(*,*) = 0. if kcsob gt 0 then sobtop(whr) = 1. ; fill in points in the top of the image with the last bright point below ;rawvig(*,*,iside) = 1. for ii = 0,nfile-1 do begin medslit = median(avgcont(*,iside,ii)) for jj = 2,nvigtop-1 do begin nslst = nsl-nvigtop+jj if sobtop(ii,jj) gt 0. then begin avgcont(nslst-1:*,iside,ii) = medslit ; rawvig(ii,jj:*,iside) = 0 goto,kilroy endif endfor kilroy: endfor ; end loop over sides endfor ; temporary to output avgcont display dispob = fltarr(nfile,nvigtop*2+2) dispob(*,*) = 0. for jj = 0,nvigtop-1 do begin for ii = 0,nfile-1 do begin dispob(ii,jj) = avgcont(nsl-nvigtop+jj,0,ii) dispob(ii,jj+nvigtop+1) = avgcont(nsl-nvigtop+jj,1,ii) endfor endfor ; now shift the slit intensity variation and average them svr(*,*) = 0. kaddv = fltarr(nsl) ; initialize the array for indicating adding these data into average kaddv(*) = 1. svv(*) = 0. for iside = 0,1 do begin for ifil = 0,nfile-1 do begin temp = avgcont(*,iside,ifil) svv = fshft_sbsp(temp, $ fitww(ifil,iside),ssl1,ssl2,nsl2) idelw = round(fitww(ifil,iside)) ; repeat wrapped rows with nearest actual variation. No repeat ; if shift rounds to zero if idelw lt 0 then begin for jj = nwm1-idelw,nwm1 do svv(jj) = $ svv(nwm1-idelw-1) endif if idelw gt 0 then begin for jj = 0,idelw-1 do svv(jj) = $ svv(idelw) endif svr(*,iside) = svr(*,iside) + svv endfor ; end loop over sides endfor svr = svr/float(nfile) pmax = max(wdelw) & pmin = min(wdelw) save,filename=input_dir(idir) + 'thermdrift.save',wdelw,ftime, $ fitww,avgcont,svr,nfile,nw,nsl,dispob ; end grand loop over directories endfor jump1: ; PROCESS 2: AVERAGE AND NORMALIZE THE SLIT VARIATIONS OVER ALL THE ; DATA ; recover the slit variations and average them for all observations ; to avoid large numbers, divide by 1000. here and near end of this ; routine where totf is also used for average flat images totf = 0. for idir = 0,ndir-1 do begin restore,concat_dir(input_dir(idir),'thermdrift.save') fnf = float(nfile)/1000. if idir eq 0 then begin nsl = sizeof_sbsp(svr,1) svrr = fltarr(nsl,2) svrr(*,*) = 0. endif ; average the files according to the number of observations in the file. ; scale to avoid very large numbers svrr = svrr + svr*fnf totf = totf +fnf endfor svrr = svrr/totf ; normalize the observed slit variations for kk = 0,1 do svrr(*,kk) = svrr(*,kk)/mean_sbsp(svrr(*,kk)) ; THIS ADDITION TO RECOVER THE GOOD SLIT VARIATION PREVIOUSLY DERIVED goto,jumpmm jump2: ; get the individual CCDSIDE slit intensity variations as used in the ; standard calibration ; use the SSW restore sswsvrr=file_search(get_logenv('$SOT_SP_CALIBRATION'),'ssw_svrr*geny') sswsvrr=last_nelem(sswsvrr) if not file_exist(sswsvrr) then begin box_message,'Cannot find svrr files under $SOT_SP_CALIBRATION' stop endif restgenx,file=sswsvrr(0),s_svrr,f_svrr jumpmm: ; PROCESS 3: GET THE AVERAGE FLAT IMAGES HAVING COMPENSATED FOR SLIT ; INTENSITY VARIATIONS ; construct an unshifted slit variation image for each CCDSIDE ; go back through all the data again, normalizing by the observed, shifted ; slit variations ; begin grand loop over maps of QS. Each of these needs to be processed totfile = 0. for idir = 0,ndir-1 do begin ; get the smoothed drift curve for this directory restore,concat_dir(input_dir(idir),'thermdrift.save') ; begin loop over files in this directory for kntr = 0,nfile-1 do begin sp4fits = derf(kntr) ; First read in the data and the FITS header print,' Begin processing file: ',kntr,sp4fits hdr=headfits(sp4fits) read_sot,sp4fits,index,dat dat = float(dat) ; DIGITAL WRAP, BIT SHIFTING CORRECTIONS ; correct for digital wrap-around of Stokes I, if wrap-around is present temp = dat(*,*,*,0) whr = where(temp lt 0., countwrap) if countwrap ne 0 then temp(whr) = temp(whr) + 65536. dat(0,0,0,0) = temp ; account for bit shifting of the data by multiplying the appropriate ; Stokes images by 2 switch bitshft of 3: dat(*,*,*,1:2) = 2.*dat(*,*,*,1:2) 2: dat(*,*,*,3) = 2.*dat(*,*,*,3) 1: dat(*,*,*,0) = 2.*dat(*,*,*,0) endswitch ;BWL- !!!THIS IS A TEMPORARY FIX FOR ERRORS IN THE 1ST FILE HEADER ;BWL- !!!THIS SHOULD BE REMOVED WHEN REFORMATTER FIXES THIS PROBLEM if spccdiy1 ne (spccdiy0 + 112 -1) then spccdiy1 = spccdiy0+112-1 ; !!!FOR ROI0, FIRST TWO COLUMNS BAD, REPLACE WITH EXTRAPOLATION if roistart eq 0 then begin ; !!!FOR ROI0, problem with first pixels of row 400 for iside = 0,1 do for istks = 0,3 do begin dat400 = mean(dat(6:10,400,iside,istks)) dat(0:5,400,iside,istks) = dat400 endfor ; for Stokes I, do quadratic fit for iside = 0,1 do begin for iyy = 0,nsl-1 do begin xxx = findgen(10) + 2. yyy = dat(3:12,iyy,iside,0) coeff = poly_fit(xxx,yyy,2,yfit=yfitt) for ii = 0,2 do begin fitt = 0. for kk = 0,2 do begin fitt = fitt + coeff(kk)*float(ii)^kk dat(ii,iyy,iside,0) = fitt endfor endfor ;dat(0:2,iyy,iside,0) = yfitt endfor endfor ; for QUV, do simple average ;for istks=0,3 do begin for istks=1,3 do begin for iside = 0,1 do begin for iyy = 0,nsl1 do begin for ixx = 0,2 do begin dat(ixx,iyy,iside,istks) = mean(dat(3:10,iyy,iside,istks)) endfor endfor endfor endfor endif ; IMAGE DIMENSIONS ; x is slit dimension, y is spectral dimension ; Fortuitously, the X-direction in the naming convention for SP header ; keywords is in the dimension along the slit, and the Y-dimension is ; in the spectral direction. Instead of retaining that convention, which as ; it happens was the naming convention for DLSP reduction software ; also, to avoid confusion we use nw as the number of pixels in the wavelength (Y) ; dimension, and nsl as the number of pixels along the slit. nw = sizeof_sbsp(dat,1) & nsl = sizeof_sbsp(dat,2) nsl1 = nsl-1 & nwm1 = nw-1 nsl2 = 2*nsl ; twice the slit length for apodizing if kntr eq 0 then begin ; Extract the number of pixels binned along the slit (along serial direction) ssum = sxpar(hdr,'CAMSSUM') ; initiate the summing array for multiplicative gain image if idir eq 0 then fltav = dblarr(nw,nsl,2) ; save the pathname for the singleframe image that will be used in the ; routine genflt3_sbsp.pro to get the relative shift between CCDSIDEs and ; the skew for each CCDSIDE if idir eq 0 then singleframe = derf(0) ; more elaborate determination of skew parameters for ROI56 restore,skewname shftsl = shftavg delsh = delshavg ; for fast map, adjust the previously determined values from ; normal map processing if ssum eq 2 then begin shftsl = shftsl/2. delsh = delsh/2. endif ; define intermediate arrays (spectral,spatial,stokes) sttmp = fltarr(nw,nsl,4) stt = sttmp stt_zero=stt ; use inside loop instead of stt(*,*,*)=0 contin = fltarr(nsl) vigg = fltarr(nsl) lmbdkslit = fltarr(nsl) ; polynomial degree, temporary array for fitting residual shifts xx = findgen(nsl) ; continuum intensity vs slit length, each side contm = fltarr(nsl,2) avgcont = fltarr(nsl,2,nfile) ; array for averaged observed slit variation for each directory svr = fltarr(nsl,2) svv = fltarr(nsl) ; gain correction image for slit width variations sltgn = fltarr(nw,nsl) ; get the solar radius by calling solar_eph.pro ; File clock time strtdate = ' ' strtdate = sxpar(hdr,'DATE_OBS') iyr = fix(strmid(strtdate,0,4)) imnth = fix(strmid(strtdate,5,2)) iday = fix(strmid(strtdate,8,2)) ahr = float(strmid(strtdate,11,2)) amn = float(strmid(strtdate,14,2)) asec = float(strmid(strtdate,17,6)) utimed = ahr + amn/60. + asec/3600. solareph_sbsp,iyr,imnth,iday,utimed,ras,dcs,gstd,b0ang, $ peeang,ctrlong,srad ; now get the theoretical limb darkening over the full range possible for a map ; The zero point of the arrays above is index 1049. This array will be zero at ; the same point. sltpos = findgen(2100) - 1049. ; calculate the limb darkening vs slit position for this map lmbdrk_sbsp,sltpos,srad,xctr(idir),yctr(idir),lmdrk ; resize slitvig, lmdrk if fast map if ssum eq 2 then lmdrk = rebin(lmdrk,2100,512) if idir eq 0 and ssum eq 2 then begin slitvig = rebin(slitvig,2100,512) endif endif ; EXTRACT PARAMETERS FROM HEADER ; Extract the number of pixels binned along the slit (along serial direction) ssum = sxpar(hdr,'CAMSSUM') ; Extract some header parameters relating to CCD region of interest ; spectral pixel range spccdiy0 = sxpar(hdr,'SPCCDIY0') spccdiy1 = sxpar(hdr,'SPCCDIY1') ; TEMP KLUDGE TO CORRECT THE OCCASIONAL ERRONEOUS VALUE NEAR HOUR CHANGE if spccdiy1 ne (spccdiy0 + 112 -1) then spccdiy1 = spccdiy0+112-1 nnw = spccdiy1-spccdiy0+1 ; check to see that image size does not change during this map if nnw ne nw then begin print,'spectral array size (',nnw,') not consistent with first file' stop endif ; Range of pixels along the slit ; Note that the pixel indices along the slit direction are reversed from ; those during pre-launch testing. The images are now reversed in ; the vertical direction during reformatting. Also, pixel indices in ; the header like SPCCDIX0, SPCCDIX1 refer to the pre-launch CCD coordinates, ; not the coordinates of LEVEL0 data. spccdix0 = sxpar(hdr,'SPCCDIX0') spccdix1 = sxpar(hdr,'SPCCDIX1') ; define these pixel positions in terms of the LEVEL0 image coordinates ; reverse order so that low limit of active range is in spccd0 spccd0 = 1023-spccdix1 spccd1 = 1023-spccdix0 ; bit shifting for SP data bitshft = sxpar(hdr,'SPBSHFT') ; extract the index of the slit scan position slpos = sxpar(hdr,'SLITPOS') ; obtain the slit scan vignetting function for this slit position sltindx = fix(slpos + 1049) ; check to see if beyond range, then set to limit if sltindx gt 2099 then begin print,' slit index for file ',kntr,' exceeds limit of 2099' print,' slit scan vignetting set to limit for this position' sltindx = 2099 endif ; for specific dark/flat eclipse season data, just use the standard darks dark = darkavg ; if ssum eq 2 then dark = darkavg_fast ; variation of vignetting function along slit for this scan position vigg(*) = slitvig(sltindx,*) ; variation of limb darkening function along slit for this scan pos. lmbdkslit(*) = lmdrk(sltindx,*) ; extract the Stokes I Q-table parameter qtabl1 = sxpar(hdr,'QTABLE1') ; THESE FLATS FOR QTABLE1 = 2 ONLY, STOP IF ENCOUNTER OTHER TABLE if qtabl1 ne 2 then begin print,'QTABLE1=',qtabl1,' for file ',sp4fits print,'THIS VALUE NOT ALLOWED, MUST HAVE ALL FLATS WITH ' print,'JPEG COMPRESSION VALUE 75' stop endif ; ADJUST SOME LIMITS IF SUMMING, OR CUTOUT ALONG SLIT, OR BOTH ; Derive active range along slit for shifting in recorded pixels, not CCD pixels ssl1 = max([sl1(1),spccd0]) - spccd0 ; always truncate the limit in the case of summing along slit ssl2 = (min([sl2(1),spccd1]) - spccd0)/ssum if ssum eq 2 then begin if (ssl1 mod 2) eq 1 then begin ssl1 = (ssl1+1)/ssum endif else begin ssl1 = ssl1/ssum endelse endif ; EXTRACT DARK FIELDS APPROPRIATE FOR THIS FILE ; for specific dark/flat eclipse season data, just use the standard darks dark = darkavg ; if ssum eq 2 then dark = darkavg_fast ; generate the slit scan vignetting spectral image for this scan position ; also include limb darkening sltvig = fltarr(nw,nsl) sltvig(*,*) = replicate(1., nw) # (vigg*lmbdkslit) ; BEGIN LOOP OVER CCDSIDEs for iside = 0,1 do begin stt=stt_zero ; fill the working array for this CCDSIDE, Stokes I only sttmp(0,0,0) = dat(*,*,iside,0) ; DARK/FLAT CORRECTION stt(0,0,0) = (sttmp(*,*,0) - dark(*,*,iside)) ; replicate the unused spectral columns to avoid discontinuities if msh1 gt 0 then begin repvec = stt(msh1,*,0) for ii = 0,msh1-1 do stt(ii,*,0) = repvec endif if msh2 lt nwm1 then begin repvec = stt(msh2,*,0) for ii = msh2+1,nwm1 do stt(ii,*,0) = repvec endif ; now, replicate the unused spatial rows to avoid discontinuities if ssl1 gt 0 then begin repvec = stt(*,ssl1,0) for jj = 0,ssl1-1 do stt(*,jj,0) = repvec endif if ssl2 lt nsl1 then begin repvec = stt(*,ssl2,0) for jj = ssl2+1,nsl1 do stt(*,jj,0) = repvec endif ; appropriately reverse shift the gain correction image to remove ; the slit pattern, taking into account the skew temp = svrr(*,iside) ; reference is average skew over same wavelength range as done above for ; deriving the slit variation of the continuum (pixels cc1:cc2 on each CCDSIDE) skew0 = mean_sbsp(shftsl(cc1:cc2,0)) for ii = 0,nwm1 do begin fskew = shftsl(ii,iside)-skew0-fitww(kntr,0) ;fskew = shftsl(ii,iside)-skew0-fitww(kntr,iside) sltgn(ii,*) = fshft_sbsp(temp,fskew,ssl1,ssl2,nsl2) endfor ; now also perform the slit scan vignetting correction ;sltgn = 1./sltgn sltgn = 1./(sltgn*sltvig) ; check to see that the intensity corresponds to quiet Sun. Reject ; low intensities if (mean(stt(*,*,0),/double) lt 1000.) then begin ; printf,luout,' bad frame at directory, file: ',idir,kntr goto,jumpss endif ; gain correct the image and sum, renormalizing to make reasonable numbers ; for grand sum fltav(*,*,iside) = fltav(*,*,iside) + sltgn*stt(*,*,0)/10000. ; incriment the counter totfile = totfile+1. jumpss: ; end loop over sides endfor ; end loop over files endfor ; end loop over directories endfor ;close,luout ; renormalize the average flat image, taking into account the factor ; of 1000 smaller for totf, and the factor of 10000 in the image summing ; new average fltav = 10000.*fltav/totfile save,filename='avgflatimg' + outfilbase + '.save',outfilbase,fltav,singleframe, $ shftsl,totfile,cc1,cc2,sp1,sp2,sl1,sl2,wc1,wc2,sr1,sr2,nxtnd,nsearch ; PROCESS 4: process the average flat image to get true gain table, and ; output all the results ;; use date of last data file in flat sequence to identify these flats ; modified routine for different spectral regions, processing genflt3_sbsp, outfilbase, fltav,singleframe,shftavg,delshavg,svrr, $ cc1,cc2,sp1,sp2,sl1,sl2,wc1,wc2,sr1,sr2,nxtnd,nsearch,/noplot return end