pro thermd_sbsp, sp4fits,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, $ show_params=show_params, display=display,verbose=verbose, $ incalpos=incalpos, param_history=param_history ; ;+ ; Name: thermd_sbsp ; ; Project: Hinode/SOT ; ; Purpose: Compute the thermal drift of image along the slit from Stokes I ; only. This will be smoothed over the duration of the measurements ; before executing the final full calibration pass. This routine does ; many of the functions of calib_sbsp.pro, but for Stokes I only. ; ; Inputs: ; sp4fits = SOT Level0 SP IQUV 4D FITS file ; darkavg_std, darkavg_fast = dark data for respective 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 ; kntr = file counter in series to enble access to thermal drift history ; l1index,l1data = fits-header structure and data of L1 data ; wdelw_prev = previous frame wdelw value ; avgctr_prev = previous avg. line center value ; tryhard - if set, then tries to find shifts from good portions of data ; firstgood = frame number of first good frame without packet loss ; ; ; Output: ; wdelw = raw thermal drift along slit (from cross-correlation) ; in this program ; ssm = summing along slit (1 or 2) ; slitdrft = Ichimoto prediction of thermal drift along slit ; specdrft = Ichimoto/Kubo prediction of thermal spectral drif ; doprcv = orbital Doppler shift DOP_CVR, pixels ; ftime = time in fractional hours from first file in directory ; sidebad = float index (#CCDSIDEs,4), when ne 0., indicates packet loss ; TEST CENTER FINDING ; avgctr = slit-average min. intensity position of 6301.5 ; ; 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> thermd_sbsp,,l1index,l1data, ; scoef,degn,sltpos,nftot,sltvar,vdeg,vigcoef,calpos,skx, $ ; darkavg_std,darkavg_fast,rmflat,rshftsl,rspecrfit,rspecshft, $ ; [,/show_params][,/display] ; (Generally called via 'sp_prep' wrapper ) ; ; History: ; 23-Dec-2006 - Bruce Lites - Adaptation from calib_sbsp.pro for calculating ; thermal drift only from Stokes I image. ; 14-Jan-2007 - Bruce Lites - extensive modifications ;version=0.11 ; 27-Jan-2007 - K.Ichimoto, center of gravity method for determining the drift along slit ; 24-May-2007 - S.L.Freeland - SOT_CALIBRATION -> SOT_SP_CALIBRATION ;version=1.0 ; 19-Feb-2008 - B. Lites adjust for possibility of NUM_SIDE=1 ; (one CCDSIDE) ; Pass error flag back to sp_prep if a file is encountered with ; all zeros, & print out bad file names ; ; ; ; -- remove depreciated calculation of thermal drift from vignetting on CCDSIDE1 ; -- just ignore the vignetting shift calculation. version = 1.01 ; 23-Apr-2008 - B. Lites (above change) ; -- processing set to skip over a file if at specified fradtion ; of the Stokes I image of either CCDSIDE = 0. ; This is the characteristic of telemetry dropouts after the ; X-band failure. Those files are considered irrecovably lost, ; and are not retained in the level1 data set. ; -- Processing continues to try to get shifts etc if tryhard ; is set, uses available good portions of frame, or substitutes ; prior frame values if information is not present version = 1.02 ; 2008.09.03 (above change) ; -- Adjust the algorithm for finding the drift of the spectral image ; along the slit to avoid regions off the limb in the cross-correlation ; with the standard slit intensity variation ; -- Change the algorithm for finding packet loss ; -- Changes in algorithm for finding dropouts (de Wijn) version = 1.03 ; 2010.01.28 (above changes) ; 6-Dec-2019 - R. Centeno - Changes to detect corrupt header ; values, particularly in camera parameters. Changes ; propagate from sp_prep.pro, and start with --- By ; RCE and end with --- END By RCE ; Files with bad headers are automatically ; blacklisted, but they still break the Level 1 ; processing both in Pass 1 and Pass 2, so they need ; to be handled separately in sp_prep, thermd_sbsp ; and getdark_sbsp. ; ;- ; display=keyword_set(display) or n_params() eq 1 ; include all the seldom-changed fixed parameters from the params_sbsp* file 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??' ; 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 dat = float(dat) l1index=index l1index.obs_type='Level 1: ' + index.obs_type update_history,l1index,/caller,version=version; report VERSION->.HISTORY 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) ww = findgen(nw) ; side-summed continuum intensity contm = fltarr(nsl) ; gain correction image for slit width variations sltgn = fltarr(nw,nsl) ; EXTRACT PARAMETERS FROM HEADER ; Extract the number of pixels binned along the slit (along serial direction) ; --- By RCE: CAMSSUM is read in ; sp_prep.pro and passed through the ; header to thermd_sbsp in ssm(Nfiles) ;ssum = sxpar(hdr,'CAMSSUM') ssum = ssm(kntr) ;--- END By RCE ; Extract exposure time, scan summing index exptime = sxpar(hdr,'EXPTIME') scnsum = sxpar(hdr,'SCN_SUM') ; 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') ; 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 ; 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. ; 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 ; determine if data are 1 or 2 CCDSIDEs numside = sxpar(hdr,'NUM_SIDE') ; get S/C pointings to test for N or S limbs ;XSUN = sxpar(hdr,'SC_ATTX') YSUN = sxpar(hdr,'SC_ATTY') ; ADJUST SOME LIMITS IF SUMMING, OR CUTOUT ALONG SLIT, OR BOTH ; use the slit variation derived from the detailed flat field ; for either fast map of full map case ssum of 1: sltvr = s_svrr(*,0) 2: sltvr = f_svrr(*,0) endcase ; adjust for cutouts along slit xcl = spccd0/ssum & xcu = spccd1/ssum sltvrc = sltvr(xcl:xcu) ; 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 GAIN AND DARK FIELDS APPROPRIATE FOR THIS FILE ; Call routine to find the appropriate dark images based upon the header information getdark_sbsp,hdr,darkavg_std,darkavg_fast,dark ; 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 endif else begin print,' invalid ssum = ',ssum,' in thermd_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) ; 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 badpacket = 0 ; flag set to 1 if this obs. has bad packet stksiok = 0 ; flag set to 1 if both CCDSIDE Stokes I is good ; Read dropout blacklist (falls back on heuristics) sidebad = checkblacklist_sbsp(index,dat) ; set flag badpacket if packet loss happens in any of the Stokes images if mean_sbsp(sidebad) ne 0. then badpacket = 1 ; if tryhard set, then examine some tests if tryhard ne 0 then begin erflag = 0 ; the next case will continue processing in this routine ; **determine if packet loss is not in either CCDSIDE of Stokes I if mean_sbsp(sidebad(*,0)) eq 0. then stksiok = 1 ; the next two cases will exit the routine and reject this frame ; **processing starts with first frame without any packet loss ; **also, the first good frame without packet loss is set if firstgood lt 0 and badpacket eq 0 $ then firstgood = kntr if firstgood lt 0 and badpacket ne 0 then begin erflag = 1 return endif ; **if both CCDSIDEs of any Stokes parameter IQUV are corrupted, then the ; **frame must be rejected for istks = 0,3 do begin if mean_sbsp(sidebad(*,istks)) eq 1. then begin erflag = 1 return endif endfor endif ; for standard processing, reject frame if any bad packets if badpacket ne 0 and tryhard eq 0 then begin erflag = 1 return endif ; BEGIN LOOP OVER CCDSIDEs for iside = 0,numside-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 ; subtract darks and correct for flat field ; This flat field intentionally does not correct for variations of ; intensity along the slit, since that is what we are seeking here stt(0,0,0) = (sttmp(*,*,0) $ - dark(*,*,iside))*mflat(*,*,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 ; 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, Stokes I only slshi_sbsp,stt,shftz,ssl1,ssl2 ; load results back into dat array dat(0,0,iside,0) = stt(*,*,0) ; end loop over CCDSIDEs endfor ; bypass merging if one CCDSIDE only, also if packet loss in ; CCDSIDE1 only if numside lt 2 then goto,oneside ; continue this processing if both CCDSIDEs of Stokes I are good. ; otherwise, just use thermd shift from previous frame. If this is ; the first frame of the series, we have abandoned it. if stksiok eq 0 and tryhard ne 0 then goto,useprevious ; MERGING CCDSIDES ; First find the precise spectral shifts ; 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 ; find umbrae (continuum intensity lt avoidumb* median along slit) contm = contm(ssl1:ssl2) ; isolate range along slit whrumb = where(contm gt avoidumb*median(contm),kountu ) ; 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),kountu ) endif if kountu le 0 then begin print,'error in thermd, probably all zero data at file ',kntr print,'file: ',sp4fits print,'Examine it. This file should probably be deleted.' erflag = 1 return endif ; smooth the residual shifts using 3rd-degree polynomial ; exclude very dark areas (umbrae) from fit 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 ; 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. ; work with Stokes I only temp0(0,0) = dat(*,*,0,0) for ii = 0,nsl1 do begin temp2(*) = dat(*,ii,1,0) ; 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 avg0 = mean_sbsp(temp0(msh1+3:msh2-3,ssl1:ssl2)) avg1 = mean_sbsp(temp1(msh1+3:msh2-3,ssl1:ssl2)) anorm = avg0/avg1 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,0) = 0.5*(temp0 + temp1) ; save the difference of the Stokes I images and display medint = median(temp1) diff = temp0 - temp1 oneside: if numside lt 2 then for kk = 0,3 do sout(*,*,kk) = dat(*,*,0,kk) ; CORRECTION FOR SPECTRAL CURVATURE, ORBITAL DOPPLER SHIFT ; correct for orbital Doppler shift, but not for smoothed drift ; from temperature variations because we don't yet have the history of ; that parameter to smooth ; Shift Stokes spectra to correct for spectral curvature for ii = 0,nsl1 do begin temp0(0,ii) = fshft_sbsp(sout(*,ii,0), $ ; at this point, 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. -spcrfit(ii)+doprcv,msh1, msh2,nxtnd) endfor sout(*,*,0) = temp0 ; endfor ; load output array, reversing the spectral direction to have increasing ; wavelength corresponging to increasing spectral pixel sout = reverse(sout,1) ; TEST CENTER FINDING ; LOCATE MINIMUM POSITION OF 6301.5 STOKES I, AVERAGE ALONG SLIT EXCLUDING ; SUNSPOTS, OFF-LIMB POSITIONS WITH LOW INTENSITY ; set some parameters del = round(500./disper) ; search limit around line,pixels lctr = fix(stdlnctr) ; standard line center position iw1 = lctr-del & iw2 = lctr+del ; pixel range for search ; Identify off-limb measurements as having continuum less than ; parameter ondisk, appropriately scaled for exposure time etc ; (ondisk usually 1000) ondsk = ondisk*exptime/4.8 if ssum eq 2 then ondsk = ondisk*scnsum*ssum*exptime/4.8 ; get the measured continuum intensity, transform it for jj = 0,nsl1 do contin(jj) = mean_sbsp(sout(cc1:cc2,jj,0)) ctrline = fltarr(nsl) ; line centers vs. slit length ; loop over length of slit to find line centers at each point for jj = 0,nsl1 do begin ; smooth the Stokes I profile smthi = smooth(sout(*,jj,0),5,/edge_truncate) ; if there are off-limb points identified by small intensities, don't ; try to fit. ; check for off-limb measurements. If mean continuum along the less than ; 5*ondisk DN (usually = 5000), then set the shift to -1000 to flag it ; for avoidance. This value is high enough that at the limb the ; 6301.5 line (but not 6302.5!!) still has an absorption core below ; the continuum if contin(jj) lt 5.*ondsk then begin ctrln = -1000. goto,jump3 endif prfmin = min(smthi(iw1:iw2),imdx) imdx = imdx + iw1 ; fitting parabola of three points is equivalent to inverting a 3x3 matrix ; to get the 2nd order polynomial coefficients yy = a + b x + c x^2 if (smthi(imdx-1) eq smthi(imdx)) and (smthi(imdx) eq smthi(imdx+1)) $ then begin ctrln = float(imdx) goto, jump3 endif parabofit,ww(imdx-1:imdx+1),smthi(imdx-1:imdx+1),ctrln jump3: ctrline(jj) = ctrln endfor ; get the average of all points along the slit with continuum intensity gt ; avoidpen*median of continuum intensity .and. ondsk. This will avoid sunspots. whr = where(ctrline gt 0.,kount) if kount gt 0 then begin contg = contin(whr) ctll = ctrline(whr) mdni = median(contg) whrc = where(contg gt ondsk and contg gt avoidpen*mdni) avgctr = mean_sbsp(ctll(whrc)) endif else begin avgctr = -1000. endelse ; 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,napod2dr]) limu = min([ssl2,nsl-napod2dr-1]) ; !!!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)) / tvvt ; /tvvt by KI 2007.1.29 ; 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 slit ; is less than ondisk DN (usually 1000), then set the shift to -1000 ; as a flag to indicate off-limb measurements if mean_sbsp(contin) lt ondsk then begin delw = -1000. goto,jump5 endif ; Ichimoto-san's 2nd derivative for continuum, exclude ends tvvt = extend_sbsp(contin,liml,limu) ; modification to avoid N,S Limbs in cross correlation, test on the ; S/C y-position if abs(ysun) gt 850. then begin maxcont = max(contin) if ysun lt 850. then begin ; south limb llimu = limu srtarr = sort(contin) srtcont = contin(srtarr) ; cont in ascending order whrsrt = where(srtcont gt 0.56*maxcont) lliml = whrsrt(0) ; pick first occurrence > 0.5*Imax endif else begin ; north limb lliml = liml srtarr = sort(contin) srtcont = contin(reverse(srtarr)) ; cont in descending order whrsrt = where(srtcont gt 0.56*maxcont,count) llimu = whrsrt(count-1) ; pick last occurrence > 0.5*Imax endelse tvvt = extend_sbsp(contin,lliml,llimu) endif else begin tvvt = extend_sbsp(contin,liml,limu) endelse ; 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>dsknorm 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 lagctr = nlag/2 lag = findgen(nlag)-lagctr cccr = c_correlate(d2cc,d2sl,lag) ; find maximum, which is shift of d2sl wrt d2cc ccmax = max(cccr,iccmx) ; fit a parabola around maximum ; iccmx = min([nlag-2,iccmx]) ; iccmx = max([1,iccmx]) ; parabofit,xx(iccmx-1:iccmx+1),cccr(iccmx-1:iccmx+1),delw ;-- center of gravity method to find a correlation peak, KI 2007.1.27 -- xa=findgen(nlag)-lagctr i1=max([0,iccmx-1]) & i2=min([iccmx+1,nlag-1]) xa1=xa(i1:i2) & cccr1=cccr(i1:i2) cccr1=(cccr1+0.35)>0 ; this formula is obtained after some trial and error to minimize ; the pixel boundary ripple in drift curve, 2007.1.29, KI delw=total(xa1*cccr1)/total(cccr1) ;;plot,cccr,psym=4,yrange=[-0.6,0.6] & oplot,delw*[1,1],[-1,1],line=1 jump5: ; load the final output numbers wdelw(kntr) = delw ;--- By RCE: ssm is obtained in ; sp_prep and passed through the ; header to thermd_sbsp so that I ; can correct corrupt values when ; header parameters are bad. ; ssm(kntr) = ssum ;--- END By RCE useprevious: if stksiok eq 0 and tryhard ne 0 then begin wdelw(kntr) = wdelw_prev avgctr = avgctr_prev endif ; if reached this point, then set error flag to 0 erflag = 0 return end