pro calib_sbsp, sp4fits, 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,sidebad,stksiok, $ xtalkimg, yearobs, $ show_params=show_params, display=display, verbose=verbose, $ inskx=inskx, incalpos=incalpos, param_history=param_history ; ;+ ; Name: calib_sbsp ; ; Project: Hinode/SOT ; ; Purpose: generate calibrated, merged Stokes spectra from 4D SP IQUV ; ; Inputs: ; sp4fits = SOT Level0 SP IQUV 4D FITS file ; fitww = structure thermal drift along slit, smoothed, pixels ; fitav = structure thermal drift, spectral direction, smoothed, pixels ; kntr = file counter in series to enble access to thermal drift history ; in this program ; scoef,degn,sltpos,nftot = residual crosstalk coefficients, degree of ; polynomial fit, slit scan positions where taken, ; and number of samples ; slitvig = 2-D slit scan vignetting image, centered on x pixel 1049 ; (full resolution image, overscanned, [2100,1024] ) ; calpos,skx = polarization calibration data ; darkavg_std,darkavg_fast = dark data for standard and fast maps ; (s,f)_mflat = flat field multiplicative images (standard,fast maps) ; (s,f)_shftsl = spectral skew fits ; (s,f)_specrfit = fit to spectral curvature ; (s,f)_specshft = spectral shift between the two CCDSIDEs ; (s,f)_svrr = intensity pattern along slit due to slit width variations ; tryhard - if set, then tries to find shifts from good portions of data ; rsdsh_prev = residual shifts from previous frame ; anorm_prev = intensity CCDSIDE normalization from previous frame ; sidebad = float index (#CCDSIDEs,4), when ne 0., indicates packet loss ; stksiok = byte index = 1 if both Stokes I CCDSIDEs are good, 0 otherwise ; ; ; Output: ; l1index,l1data - fits-header structure and data of L1 data ; wlshft = pixel wavelength shift applied, for L1 KEYWORD SPWLSHFT ; slshft = pixel shift applied, for L1 KEYWORD SPSLSHFT ; ; Keyword Inputs: ; inskx - optional polarization calibration coeff [4,4,4,4,2,8] ; incalpos - optional slit scan positions where cal. data ; (defaults for skx&calpos now via $SOT_CALIABRATION files) ; show_params (switch) - if set, echoback parameter definitions ; verbose - if set, show filenames as processed ; display (switch) - if set, display to ; param_history - if set, add full parameter settings to ; .HISTORY ; default is paramter file names ; ; Calling Sequence: ; IDL> calib_sbsp,,l1index,l1data, ; scoef,degn,sltpos,nftot,sltvar,vdeg,vigcoef,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, $ ; [,/show_params][,/display] ; (Generally called via 'sp_prep' wrapper ) ; ; History: ; 6-Oct-2006 - Bruce Lites, HAO. ; 13-Oct-2006 - S.L.Freeland - ssw hooks, environmentals etc. ; parameter files via $SSW_CALIBRATION ; skx & calpos derived from cal files, not input by this routine ; 16-Oct-2006 - pick up.. add /DISPLAY (default is no display ; rename sizeof -> sizeof_sbsp (ssw-gen conflict) ; All restores -> restgenx; add date to param/geny fnames ; 17-Oct-2006 - S.L.Freeland - add some HISTORY; end-to-end testing ; 23-Nov-2006 - Bruce Lites - Major structural changes, reorganization, ; revisions for better calibration ; 23-Jan-2007 - Bruce Lites - many revisions throughout the code. ; ;- ; ;version=1.0 ;version = 1.01 ; correct sign of coordinate rotation angle in OP. 9 ; ;version = 1.02 ; 19-Feb-2008 - B. Lites - adapt to be able to ; handle NUM_SIDE=1 data (one CCDSIDE) ;version = 1.03 ; 23-Apr-2008 - B. Lites - remove depreciated calculation ; of thermal drift from vignetting on CCDSIDE1 ; -- just ignore the vignetting shift calculation. version = 1.04 ; 4-Sep-2008 - B. Lites ; -- adapted to handle frequent packet loss. Standard option ; rejects frames with any packet loss. tryhard option (set ; with /tryharder keyword in sp_prep.pro, attempts to make ; use of available good images in frame with some packet loss version = 1.05 ; 20-Jan-2009 - B. Lites ; -- slit scan range can exceed available data for slitvig ; function. Use limits of slitvig beyond the available range version = 1.06 ; 28-Jul-2011 - B. Lites ; -- add the residual crosstalk correction icrossv_sbsp.pro for ; Stokes I->V row-by-row correction. Does this for both one- ; and two-sided data, but only if the integrated Stokes V/I over ; the whole wavelength range is less than 0.02. Change the ; residual crosstalk correction for variations along slit, ; dont have variations along spectrum. See new residcross_sbsp.pro version = 1.07 ; 12-Nov-2013 - B. Lites ; -- update residual crosstalk correction procedure, add input ; variables xtalkimg, yearobs display=keyword_set(display) or n_params() eq 1 ; include all the seldom-changed fixed parameters ; 24-may-2007 - S.L.Freeland - SOT_CALIBRATION -> SOT_SP_CALIBRATION scal=get_logenv('SOT_SP_CALIBRATION') if scal eq '' then set_logenv,'SOT_SP_CALIBRATION',$ concat_dir('$SSW_SOT','calibration') 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??' ;BWL-READING OF POLARIZATION CALIBRATION COEFFICIENTS DONE ONCE IN THE ;BWL-SP_PREP ROUTINE. ARE USER OVERRIDES TO BE DONE HERE: if keyword_set(incalpos) then calpos=incalpos ; user override if keyword_set(inskx) then skx=inskx ; user override ; First read in the data and the FITS header hdr=headfits(sp4fits(0)) if keyword_set(verbose) then begin print,' Begin processing file: ',sp4fits(0) read_sot,sp4fits(0),index,dat endif else begin read_sot,sp4fits(0),index,dat,/silent endelse ;BWL-INPUT DATA IS INTEGER,MUST BE CONVERTED TO FLOATS FOR PROCESSING!! -BWL ;BWL-!!!NEXT LINE CAN BE OMITTED IF READ_SOT CONVERTS SP DATA TO FLOATS dat = float(dat) l1index=index l1index.obs_type='Level 1: ' + index.obs_type ; update the history for calib_sbsp.pro update_history,l1index,/caller,version=version; report VERSION->.HISTORY ; update history to include names of the current calibration files used update_history,l1index,str_replace(sotcal,scal,'$SOT_SP_CALIBRATION') if keyword_set(param_history) then $ update_history,l1index,str_replace(pdata,'; ',''),/noroutine update_history,l1index,str_replace(pcoef,scal,'$SOT_SP_CALIBRATION') ; 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 ; 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) ; prepare output merged Stokes spectra images temp0 = fltarr(nw,nsl) temp1 = fltarr(nw,nsl) temp2 = fltarr(nw) sout = fltarr(nw,nsl,4) ; polynomial degree, temporary array for fitting residual shifts ndegree = 3 rsdsh = fltarr(nsl) xx = findgen(nsl) ; side-summed continuum intensity contm = fltarr(nsl) ; gain correction image for slit width variations sltgn = fltarr(nw,nsl) tmpsv = fltarr(nw,nsl,2) ; 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') ;BWL- !!!THIS IS A TEMPORARY FIX FOR ERRORS IN THE 1ST FILE HEADER OR NEAR ; HOUR CHANGE ;BWL- !!!THIS SHOULD BE REMOVED WHEN REFORMATTER FIXES THIS PROBLEM 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) ; avoid exceeding slit scan range of available slitvig data if slpos lt -1049 then sltindx = 0 if slpos gt 1050 then sltindx = 2099 ; variation of vignetting function along slit for this scan position vigg = slitvig(sltindx,*) ;vigg = 0. ;for kv = 0,vdeg do vigg = vigg + vigcoef(kv)*slpos^kv ; spacecraft-sun relative velocity (m/sec, positive = redshift), ; converted to SP pixels. HARDWIRED for SP. doprcv = sxpar(hdr,'DOP_RCV')*wavl*1000./(cccc*disper) ; get the CCD temperatures for estimating wavelength, spatial drifts tfgccd = sxpar(hdr,'T_FGCCD') tspccd = sxpar(hdr,'T_SPCCD') ; use Ichimoto-san's linear fits to get the thermal offsets in both ; directions, in pixels. Note that this formula for specdrft has ; redshift positive! specdrft = 0.57342*tfgccd - 0.03594*tspccd + 21.86 slitdrft = -0.00359*tfgccd + 0.44022*tspccd + 18.41 ; rotation angle between solar north and the SP CCD crota2 = sxpar(hdr,'CROTA2') rotn = crota2*!dtor ; determine if data are 1 or 2 CCDSIDEs numside = sxpar(hdr,'NUM_SIDE') ; 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 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 ; OPERATION 1: EXTRACT GAIN AND DARK FIELDS APPROPRIATE FOR THIS FILE ; Call routine to find the appropriate dark images based upon the header information ; the dark routine already selects the proper region along the slit and ; adjusts for summing along the slit, exposure times, etc. getdark_sbsp,hdr,darkavg_std,darkavg_fast,dark,l1index=l1index ; select gain table for either fast map or standard map if ssum eq 1 then begin mflat = s_mflat shftsl = s_shftsl spcrfit = s_specrfit spcshft = s_specshft svrr = s_svrr endif else if ssum eq 2 then begin mflat = f_mflat shftsl = f_shftsl spcrfit = f_specrfit spcshft = f_specshft svrr = f_svrr ; for slit scan vignetting, rebin the function to 512 if fast map vigg = rebin(vigg,512) endif else begin print,' invalid ssum = ',ssum,' in calib_sbsp' stop endelse ; select these quantities for the portion of the slit observed (cutout) ssr0 = spccd0/ssum & ssr1 = spccd1/ssum mflat = mflat(*,ssr0:ssr1,*) spcrfit = spcrfit(ssr0:ssr1) spcshft = spcshft(ssr0:ssr1) ; for slit scan vignetting function vigg = vigg(ssr0:ssr1) ; OPERATION 2: 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 ;; OPERATION 2.1: FIND APPROXIMATE THERMAL FLEXURE SHIFT FROM VIGNETTING ;; AT TOP OF IMAGE OF CCDSIDE1 IF THE DATA SAMPLE THAT REGION ; vigsh = -2000 ; if (numside eq 2) and (spccd1 eq 1023) then $ ; vigshift_sbsp,dat(*,*,1,0),ssum,vigsh ; OPERATION 3: DERIVE INTERPOLATED INVERSE X-MATRIX ; selected for this slit position and spectral ROI terp_xmat_sbsp,slpos,calpos,skx,spccdiy0,spccdiy1,xinv ; reverse the inverse matrix elements along the slit because the flight ; data reformatter now flips the images top-to-bottom relative to the GSE ; reformatter that the cal data were based upon xinv = reverse(xinv,2) ; for matrix inverse, isolate the region of interest along the slit xinv = xinv(*,spccd0:spccd1,*,*,*) ; rebin matrix inverse in the case of summing along the slit if ssum eq 2 then xinv = rebin(xinv,nw,nsl,4,4,2) ; BEGIN LOOP OVER CCDSIDEs for iside = 0,numside-1 do begin stt=stt_zero ; fill the working array for this CCDSIDE for kk = 0,3 do sttmp(0,0,kk) = dat(*,*,iside,kk) ; OPERATION 4: DARK/FLAT CORRECTION ; subtract darks and correct for flat field sttmp(0,0,0) = (sttmp(*,*,0) $ - dark(*,*,iside))*mflat(*,*,iside) for kk = 1,3 do sttmp(0,0,kk) = sttmp(*,*,kk) $ *mflat(*,*,iside) ; now correct for shifted slit intensity variation using pre-determined ; empirical thermal drifts from first pass through data (thermd_sbsp.pro) ; to shift the pattern, and build a multiplicative gain correction image ; for the slit intensity variation that encorporates the skew svrg = svrr(ssr0:ssr1,iside) ; reference is average skew over same wavelength range as done above for ; deriving the slit variation of the continuum (pixels 1:11 on each CCDSIDE) ; here since we are using the shift determined from CCDSIDE0 for both ; CCDSIDEs, slit variation gain image is calculated using the full ; skew determination skew0 = mean_sbsp(shftsl(0:10,0)) for ii = 0,nwm1 do begin fskew = shftsl(ii,iside)-skew0-fitww(kntr) sltgn(ii,*) = fshft_sbsp(svrg,fskew,ssl1,ssl2,nsl2) endfor tmpsv(*,*,iside) = sltgn sltgn = 1./sltgn ; gain correct the Stokes images for this slit intensity pattern for kk = 0,3 do sttmp(0,0,kk) = sltgn*sttmp(*,*,kk) ; OPERATION 5: POLARIZATION CALIBRATION WITH SUN TEST X-MATRIX for jj = 0,3 do begin for ii = 0,3 do stt(0,0,jj) = stt(*,*,jj) + $ xinv(*,*,jj,ii,iside)*sttmp(*,*,ii) endfor ; replicate the unused spectral columns to avoid discontinuities if msh1 gt 0 then begin for kk = 0,3 do begin repvec = stt(msh1,*,kk) for ii = 0,msh1-1 do stt(ii,*,kk) = repvec endfor endif if msh2 lt nwm1 then begin for kk = 0,3 do begin repvec = stt(msh2,*,kk) for ii = msh2+1,nwm1 do stt(ii,*,kk) = repvec endfor endif ; now, replicate the unused spatial rows to avoid discontinuities if ssl1 gt 0 then begin for kk = 0,3 do begin repvec = stt(*,ssl1,kk) for jj = 0,ssl1-1 do stt(*,jj,kk) = repvec endfor endif if ssl2 lt nsl1 then begin for kk = 0,3 do begin repvec = stt(*,ssl2,kk) for jj = ssl2+1,nsl1 do stt(*,jj,kk) = repvec endfor endif ; OPERATION 6: SKEW CORRECTION, COMPENSATION FOR SHIFTS ALONG SLIT BETWEEN SIDES ; remove skew of each CCDSIDE, which also rectifies the realtive shift along ; the slit between the two CCDSIDEs shftz = shftsl(*,iside) ; perform shifting slsh_sbsp,stt,shftz,ssl1,ssl2 ; load results back into dat array for kk = 0,3 do dat(0,0,iside,kk) = stt(*,*,kk) ; end loop over CCDSIDEs endfor ; bypass merging if one CCDSIDE only if numside lt 2 then goto,oneside ; OPERATION 7: MERGING CCDSIDES ; First find the precise spectral shifts ; This can only be done if both CCDSIDE Stokes I images are good. ; If one Stokes I CCDSIDE image is corrupted, use previous frame shifts if tryhard eq 0 or stksiok eq 1 then begin ; Do preliminary spectral shift of Stokes I between the two CCDSIDES ; using the spectral shifts computed from the flat fields. ; Shift CCDSIDE1 to match CCDSIDE0. temp0(0,0) = dat(*,*,0,0) for ii = 0,nsl1 do begin temp2(*) = dat(*,ii,1,0) temp2 = fshft_sbsp(temp2,spcshft(ii),ssh1,ssh2,nxtnd) temp1(*,ii) = temp2 endfor ; now find the residual spectral shift of these data along the slit ; find only very local maximum, searching 5 pixels. This will ; be added to the preliminary shift to get the total shift for jj = ssl1, ssl2 do begin tempr = temp0(wc1:wc2,jj) temps = temp1(wc1:wc2,jj) rsdsh(jj) = corshft_sbsp(temps,tempr,nsrch=5,/noslope) ; get the spectral median continuum intensity contm(jj) = median(tempr+temps) endfor ; avoid umbrae (continuum intensity lt avoidumb of median along slit) contm = contm(ssl1:ssl2) ; isolate range along slit whrumb = where(contm ge avoidumb*median(contm),countumb ) ; FOR OFF-LIMB, CHOOSE MEAN OF CONTM INSTEAD OF MEDIAN if median(contm) lt 0.05*max(contm) then begin whrumb = where(contm gt avoidumb*mean_sbsp(contm),countumb ) endif if countumb le 0 then begin print,' warning in calib_sbsp.pro: no non-umbral points found' endif ; smooth the residual shifts using 3rd-degree polynomial ; exclude very dark areas (umbrae) from fit, perform smoothing if meaningful if countumb gt (ndegree+1) then begin rrr = rsdsh(ssl1:ssl2) coef = poly_fit(xxr(whrumb),rrr(whrumb),ndegree) rsdsh(*) = coef(0) for idg = 1,ndegree do rsdsh = rsdsh + coef(idg)*xx^idg endif rsdsh_prev = rsdsh endif ; if using tryharder option and one of Stokes I images has packet loss, then ; use spectral shifts and intensity normalizations from previous frame if tryhard eq 1 and stksiok eq 0 then begin rsdsh = rsdsh_prev endif ; now apply spectral shifts and sum the CCDSIDEs ; loop over Stokes parameters, scaling and shifting them spectrally, ; then adding the two sides to compute merged images. Shift CCDSIDE1 to ; match CCDSIDE0. ; begin loop over Stokes parameters for kk = 0,3 do begin temp0(0,0) = dat(*,*,0,kk) for ii = 0,nsl1 do begin temp2(*) = dat(*,ii,1,kk) ; use the computed value for shift from gain, plus the residual correction temp2 = fshft_sbsp(temp2,spcshft(ii)+rsdsh(ii), $ ssh1,ssh2,nxtnd) temp1(*,ii) = temp2 endfor ; Rectify intensity differences between the CCDSIDEs ; there can be a slight residual intensity offset between the two sides ; no apparent difference in slope of continuum however, so just adjust the ; overall difference in intensity of CCDSIDE1 to that of CCDSIDE1 ; find mean of Stokes I images of each CCDSIDE excluding first and last ; three columns if kk eq 0 then begin avg0 = mean_sbsp(temp0(msh1+3:msh2-3,ssl1:ssl2)) avg1 = mean_sbsp(temp1(msh1+3:msh2-3,ssl1:ssl2)) anorm = avg0/avg1 if tryhard eq 1 then begin if stksiok eq 1 then begin anorm_prev = anorm endif else begin anorm = anorm_prev endelse endif endif temp1 = temp1*anorm ; After polarization calibration, Stokes parameters are sum of two CCDSIDEs ; Use the mean of the sides to avoid digital overflow on conversion to integers ;sout(0,0,kk) = 0.5*(temp0 + temp1) sidegood = 1.-sidebad sout(0,0,kk) = (sidegood(0,kk)*temp0 + $ sidegood(1,kk)*temp1)/total(sidegood(*,kk)) ; save the difference of the Stokes I images and display if kk eq 0 and display then begin medint = median(temp1) diff = temp0 - temp1 if display then begin delvarx,xx wdef,xx,image=diff,/ur tvscl,(diff/medint < 0.03)>(-0.03) endif endif ; end loop over Stokes parameters endfor oneside: if numside lt 2 then begin for kk = 0,3 do sout(*,*,kk) = dat(*,*,0,kk) if display then begin delvarx,xx wdef,xx,image=sout(*,*,3),/ur endif endif ; OPERATION 8: CORRECTION FOR SPECTRAL CURVATURE, ORBITAL DOPPLER SHIFT ; Shift Stokes spectra to correct for spectral curvature for kk = 0,3 do begin ; use empirical correction to put average pixel of 6301.5 at 29.0 thrmdr = doprcv-fitav(kntr) for ii = 0,nsl1 do begin ; also correct for spacecraft-Sun Dopplershift ; presuming the convention for the header information is POSITIVE = REDSHIFT ; note at this point in the analysis, redshift is toward lower spectral ; pixel value, so the applied shift to correct for a reshift is positive. ; The spectra are spectrally reversed farther down in this ; code in order to make redshift correspond to increasing pixel value. temp0(0,ii) = fshft_sbsp(sout(*,ii,kk), $ -spcrfit(ii)+thrmdr,msh1, msh2,nxtnd) endfor sout(*,*,kk) = temp0 endfor ; OPERATION 9: POLARIZATION REFERENCE FRAME CORRECTION, AND ; REVERSAL OF SPECTRAL DIMENSION ; convert the Stokes vectors from the reference frame of the FPP to that of the ; Sun. crota2 = sxpar(hdr,'CROTA2') rotn = crota2*!dtor ; CROTA2 is the angle, in degrees, of the slit orientation measured ; positive clockwise from solar north. ; The rotation angle in the coordinate rotation is positive in the ; CCL direction. Thus if data are taken at a positive CROTA2, the ; coordinate frame must be rotated by a positive amount (+CROTA2) to ; arrive in the solar reference frame with +Q along solar E-W. filtq = sout(*,*,1) filtu = sout(*,*,2) crot = cos(2.*rotn) & srot = sin(2.*rotn) sout(0,0,1) = crot*filtq + srot*filtu sout(0,0,2) = -srot*filtq + crot*filtu ; load output array, reversing the spectral direction to have increasing ; wavelength corresponging to increasing spectral pixel sout = reverse(sout,1) ; OPERATION 10: RESIDUAL I->QUV CROSSTALK COMPENSATION ; to get true residual, bypass all corrections ;goto,jump443 ; get pre-determined residual crosstalk correction arrays ; if yearobs lt 2013. then residcross_sbsp,slpos,scoef,degn,sltpos,nftot,nw,nsl,resid ; if yearobs ge 2013. then residcross2_sbsp,slpos,xtalkimg,sltpos,nftot,nw,nsl,resid ; can now use this new process for all years since calibration files are updated residcross2_sbsp,slpos,xtalkimg,sltpos,nftot,nw,nsl,resid ; now use the fitted crosstalk from the quiet Sun map. Will ; do linear interpolation in the slit scan range. This ; does not remove the row-by-row crosstalk, but it removes ; the overall offset and trends. stksi = sout(*,*,0) for istks = 1,3 do sout(0,0,istks) = sout(*,*,istks) - $ resid(*,*,istks-1)*stksi ; Since there is some residual feed-through of the granulation pattern into ; Stokes Q,U, but not V, correct those spectra on a spectrum-by-spectrum ; basis ;!!!TEMP BYPASS ICROSS ;goto,jump443 icrossqu_sbsp,sout,ssl1,ssl2 ; Generally for 2-sided data, no V crosstalk correction is needed. For ; 1-sided data, however, it is quite apparent, so apply a correction ; apply to all data ;if numside lt 2 then icrossv_sbsp,sout,ssl1,ssl2 icrossv_sbsp,sout,ssl1,ssl2 jump443: ; OPERATION 11: SHIFT SPECTRA UP OR DOWN TO CORRECT FOR THERMAL FLEXURE ; OF STRUCTURE, THEN CORRECT FOR SLIT SCAN VIGNETTING ; Sensing was done in first pass by thermd_sbsp.pro, now just do correction idelw = round(fitww(kntr)) ; generate the slit scan vignetting spectral image for this scan position sltvig = fltarr(nw,nsl) sltvig(*,*) = replicate(1., nw) # vigg ; now shift the data ; begin loop over Stokes parameters for istks = 0,3 do begin for ii = 0,nwm1 do begin sout(ii,*,istks) = fshft_sbsp(sout(ii,*,istks), $ fitww(kntr),ssl1,ssl2,nsl2) endfor ; 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 sout(*,jj,istks) = $ sout(*,nwm1-idelw-1,istks) endif if idelw gt 0 then begin for jj = 0,idelw-1 do sout(*,jj,istks) = $ sout(*,idelw,istks) endif ; now correct only for slit scan vignetting, width variations are now ; compensated in the gain correction for each CCDSIDE sout(0,0,istks) = sout(*,*,istks)/sltvig ; end loop over Stokes parameters endfor ; OPERATION 12: BACK-CONVERSION TO INTEGER, OUTPUT TO FITS ; convert to integer, rescale the data and bit-shift it to the ; original raw data. Ensure that there are no negative intensities ; that might occur in off-limb observations temp = sout(*,*,0) whr = where(temp lt 0.,count) if count gt 0 then temp(whr) = 0. sout(*,*,0) = temp ; bit shift the data back again to original scaling in order to avoid ; any potential overflows, especially in Stokes I. The final data will ; need to have bit shifting restored! switch bitshft of 3: sout(*,*,1:2) = 0.5*sout(*,*,1:2) 2: sout(*,*,3) = 0.5*sout(*,*,3) 1: sout(*,*,0) = 0.5*sout(*,*,0) endswitch ; final data output, alter header l1data = fix(round(sout)) l1index.naxis=3 l1index.naxis1=data_chk(sout,/nx) l1index.naxis2=data_chk(sout,/ny) l1index.naxis3=data_chk(sout,/nimage) ; four Stokes parameters ;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 slshft = fitww(kntr) wlshft = fitav(kntr) ; OPERATION 13: FINAL DISPLAY OF DATA (OPTIONAL) ; setup final result display if display then begin bstks = bytarr(nw*4,nsl) medi = 0.005*median(sout(*,*,0)) for kk = 0,3 do begin xstrt = kk*nw if kk eq 0 then begin bstks(xstrt:xstrt+nwm1,*) = $ bytscl(sout(*,*,kk)) endif else begin bstks(xstrt:xstrt+nwm1,*) = $ bytscl(sout(*,*,kk),max=medi,min=-medi) endelse endfor wdelete,xx delvarx,xx wdef,image=bstks tvscl,bstks endif return end