FUNCTION eis_update_fitdata, fitdata, line=line, yrange=yrange, offset=offset, $ fit_error=fit_error, refpix=refpix, fitdiff=fitdiff, $ quiet=quiet, median=median ;+ ; NAME ; ; EIS_UPDATE_FITDATA() ; ; PROJECT ; ; Hinode/EIS ; ; EXPLANATION ; ; This routine updates the orbit correction stored in the OFFSET ; tag of FITDATA by using the line fits within FITDATA. ; ; A typical example of the use of EIS_UPDATE_FITDATA is if a first ; approximation to the orbit variation has been derived with ; EIS_WAVE_CORR in order to perform the Gaussian fit. After the ; fit has been performed the line fits themselves are used to ; derive the orbit variation of the line centroid. By specifying ; YRANGE a sub-set of the Y-pixels can be used for the new orbit ; variation. ; ; INPUTS ; ; FITDATA An EIS fit structure produced by EIS_AUTO_FIT. ; ; OPTIONAL INPUTS ; ; YRANGE Specifies the Y-pixels to be used to determine the ; orbit variation. Can be a two element array, indicating ; a range, or a single element, indicating a single ; pixel. E.g., YRANGE=[0,19] implies Y-pixels 0 to 19 ; (inclusive) will be used; YRANGE=19 implies only ; Y-pixel 19 will be used. If not specified, then the ; entire Y-range is used. ; ; LINE If a multi-Gaussian fit has been performed then LINE= ; is used to specify which of the lines within FITDATA is ; to be used to update the orbit correction. If not ; specified then LINE=0. ; ; OPTIONAL OUTPUTS ; ; OFFSET The offset tag within the FITDATA structure (which is ; the one input to EIS_AUTO_FIT) is updated by ; EIS_UPDATE_FITDATA. Now the FITDATA array sizes can be ; different to those in the original WINDATA structure if ; EIS_BIN_WINDATA, or the XRANGE and YRANGE keywords have ; been used in the call to EIS_AUTO_FIT. By specifying ; the optional output OFFSET, an OFFSET array of the same ; size as the original WINDATA arrays is returned. If you ; have used XRANGE in the call to EIS_AUTO_FIT then it is ; not possible to obtain the orbit correction for the ; entire original X-range of WINDATA, so the OFFSET ; values for X-values outside this range are set to ; FITDATA.MISSING. It is recommended that you do not use ; this optional output if you have used XRANGE or XBIN. ; ; FIT_ERROR This is the error in the fit to the centroid ; variation in angstroms. It is computed as the ; standard deviation of (y-yfit) where y is the ; variation of the centroid in time, and yfit is the ; fit to the variation. ; ; REFPIX Contains the reference pixel. This is needed if the ; error due to the uncertainty in the EIS slit tilt ; needs to be computed. ; ; KEYWORDS ; ; QUIET If set, then a plot is not made and no messages are ; sent to the terminal. ; ; MEDIAN By default the routine averages the centroids along ; the specified Y-range. Setting /MEDIAN makes the ; routine take the median instead. This can be useful ; if strong Doppler shifts in small region are ; affecting the centroid variation curve. ; ; OUTPUTS ; ; A new FITDATA structure is returned, with the following tags ; modified: ; .refwvl This is updated with the average value of the line ; centroid over YRANGE. For multi-Gaussian fits the ; difference between the new refwvl value and the old ; one for the line specified by LINE= is used to update ; the refwvl values for the remaining lines. ; .offset A new offset array is computed by using the emission ; line fits to obtain a new orbit correction. The slit ; tilt will be exactly the same as the original array ; (assuming you have used eis_wave_corr). ; .aa The fit parameters for the line centroids are ; corrected for the new offset array. ; .x_bg1 The wavelength value that defines the lower edge of the ; wavelength window is adjusted for the new offset. ; .x_bg2 The wavelength value that defines the upper edge of the ; wavelength window is adjusted for the new offset. ; All other tags remain the same as the original FITDATA structure. ; ; EXAMPLE ; ; The following first performs an auto-fit using an approximate ; offset (derived using eis_wave_corr, for example), then uses the ; fitdata structure to derive a new offset which is then used to ; update the centroid information stored in FITDATA. The keyword ; YRANGE= is used to specify the Y-pixels to be used to determine ; the orbit variation. In this case pixels 0 to 19 are used. ; ; IDL> eis_wave_corr, windata, offset ; IDL> eis_auto_fit, windata, fitdata, offset=offset ; IDL> newfitdata=eis_update_fitdata(fitdata,yrange=[0,19]) ; ; CALLS ; ; EIS_GET_FITDATA, EIS_CUBIC_SPLINE, EIS_SLIT_TILT_ARRAY ; ; HISTORY ; ; Version 1, 2-Feb-2010, Peter Young ; Version 2, 4-Feb-2010, Peter Young ; updated header ; Version 3, 13-Apr-2010, Peter Young ; added fit_error= and refpix= optional output ; Version 4, 15-Apr-2010, Peter Young ; missing data in OFFSET are now set correctly ; Version 5, 23-Apr-2010, Peter Young ; if FITDATA contains multiple lines and the user has not ; specified LINE=, then the routine now exits and prints an ; information message; now prints out wavelength of line used on ; plot; exposures with zero exposure time are set to be missing ; in OFFSET ; Version 6, 8-Jan-2013, Peter Young ; added /QUIET and /MEDIAN keywords. ; Version 7, 11-Feb-2020, Peter Young ; added information message if raster duration is too short for ; the orbit correction. ;- ;------- ; NOTE: FITDATA remains unchanged throughout the code. OUTDATA is a ; copy of FITDATA that will be modified and returned to the user as ; the output. ; outdata=fitdata ; ; The input LINE determines which of the emission lines (in the case ; of a multi-Gaussian fit) will be used to determine the orbit ; variation. The default is the first line. ; ng=fitdata.ngauss IF n_elements(line) EQ 0 AND ng GT 1 THEN BEGIN print,'' print,'%EIS_UPDATE_FITDATA: The fit structure contains '+trim(ng)+' emission lines. Please select' print,' which line to use by specifying the keyword LINE=.' print,' (LINE takes values 0, 1, 2, ....)' return,-1 ENDIF ; ; For the single Gaussian case set line=0 ; IF n_elements(line) EQ 0 THEN line=0 ; ; We want to use the centroid data in FITDATA to derive a new orbit ; correction. However the centroids have already been corrected for the orbit ; correction contained in FITDATA.OFFSET, so we need to 'uncorrect' ; this first. ; int=eis_get_fitdata(fitdata,line=line) cen=eis_get_fitdata(fitdata,/cen,line=line) imiss=where(cen EQ fitdata.missing,nmiss) cen=cen+fitdata.offset IF nmiss GT 0 THEN cen[imiss]=fitdata.missing tt=fitdata.exp_start_times ; ; The slit tilt depends on the wavelength channel, so determine which ; channel data belongs to ; IF mean(fitdata.refwvl) GE 230. THEN long=1 ELSE long=0 missing=fitdata.missing nx=fitdata.nx ny=fitdata.ny slit_ind=fitdata.slit_ind ybin=fitdata.ybin IF n_elements(yrange) EQ 0 THEN yrange=[0,ny-1] ; ; 'refpix' is the Y-pixel where the tilt correction is set to be ; zero. It is derived from YRANGE. Note that this 'refpix' is ; different to the one used later for the full size offset array. ; IF n_elements(yrange) EQ 1 THEN refpix=yrange ELSE refpix=round(mean(yrange)) ; ; Apply the tilt correction to the CEN array. Giving REFPIX as an ; input means that the tilt correction at REFPIX will be zero. ; tilt_corr=eis_slit_tilt_array(nx,ny,ybin,slit_ind,fitdata.yip,fitdata.date_obs, $ refpix=refpix, $ ystart=fitdata.start_pix[1], long=long) cen=cen-tilt_corr IF nmiss GT 0 THEN cen[imiss]=missing IF nx EQ 1 THEN BEGIN print,'%EIS_UPDATE_FITDATA: there is only one exposure in this data-set, so the orbit' print,' correction can not be performed. Returning...' return,-1 ENDIF t=tt/60. ; convert times to minutes error=0 ; ; Do not proceed if the raster duration is less than 10 minutes ; IF max(t)-min(t) LT 10 THEN BEGIN error=1 dt=trim(string(format='(f6.1)',max(t)-min(t))) print,'%EIS_UPDATE_FITDATA: the raster duration is '+dt+' mins, which is too short to perform the ' print,' orbit correction. Returning...' return,-1 ENDIF ; ; Collapse the 3D array along the slit (2nd dimension) ; IF n_elements(yrange) EQ 1 THEN BEGIN cenx=reform(cen[*,yrange]) ENDIF ELSE BEGIN nyx=yrange[1]-yrange[0]+1 IF keyword_set(median) THEN BEGIN cenx=fmedian(cen[*,yrange[0]:yrange[1]],1,nyx,missing=missing) ENDIF ELSE BEGIN cenx=average(cen[*,yrange[0]:yrange[1]],2,missing=missing) ENDELSE ENDELSE ; ; Check for T<1440 to filter out points with anomalously high time ; values. ; i_notmiss=where(cenx NE missing AND t LE 1440.,n_notmiss) IF n_notmiss NE 0 THEN BEGIN scenx=cenx[i_notmiss] tx=t[i_notmiss] ENDIF ELSE BEGIN tx=t ENDELSE ; ; Plot the non-missing centroid values (CENX) as a function of time ; IF NOT keyword_set(quiet) THEN BEGIN plot,tx,scenx,/xsty,xra=[min(tx),max(tx)], /ynozero, $ xtit='Time / minutes', $ ytit='Centroid variation / '+string(197b), psym=psym, symsiz=1.5, $ charsiz=1.4 ENDIF ; ; Get spline fit to line centroid and overplot the spline fit; first ; need to make sure that TX is in ascending time order (this is a ; requirement of eis_cubic_spline). ; ; Note: this is the only place in the code where I have to worry about ; the order of time. The rest of the code is independent of the time ; ordering. ; IF fitdata.raster EQ 1 THEN BEGIN tx=reverse(tx) scenx=reverse(scenx) ENDIF y_spline=eis_cubic_spline(tx,scenx,nodes,yfit=yfit) ; fit_error=stdev(yfit-scenx) IF NOT keyword_set(quiet) THEN BEGIN oplot,tx,yfit,th=2 xr=!x.crange yr=!y.crange errstr=trim(string(format='(f10.4)',fit_error))+' '+string(197b) wvlstr=trim(string(format='(f10.3)',fitdata.refwvl[line]))+' '+string(197b) xyouts,0.95*xr[0]+0.05*xr[1],0.07*yr[0]+0.93*yr[1],'Line = '+wvlstr,charsiz=1.4 xyouts,0.95*xr[0]+0.05*xr[1],0.10*yr[0]+0.90*yr[1],'Error = '+errstr,charsiz=1.4 ENDIF ; ; For raster data need to reverse YFIT, since 1st element corresponds ; to time=0 ; IF fitdata.raster EQ 1 THEN yfit=reverse(yfit) ; ; The new 'refwvl' value in OUTDATA is just the average of the spline ; fits YFIT. ; ; Note that, for multi-Gauss fits, each of the 'refwvl' values is ; corrected by the correction found for 'line' ; refwvl=mean(yfit) refwvl_corr=refwvl-fitdata.refwvl[line] outdata.refwvl=fitdata.refwvl+refwvl_corr IF NOT keyword_set(quiet) THEN print,'%EIS_UPDATE_FITDATA: refwvl['+trim(line)+'] updated with value '+string(format='(f7.3)',refwvl) ; ; Now derive the orbit variation for all times (including ; missing). Note that for exposures with bad time values yfit_all will ; probably take a crazy value. ; y2=spl_init(nodes,y_spline) yfit_all=spl_interp(nodes,y_spline,y2,t) ; ; Create the new offset array that goes into FITDATA. Note that it is ; defined such that the offset should that centroids for the row of ; pixels identified by REFPIX should be close to REFWVL. ; offset=fltarr(nx,ny) FOR i=0,ny-1 DO offset[*,i]=tilt_corr[*,i]+yfit_all-refwvl ; ; Replace OFFSET in OUTDATA and update the centroid information stored ; in FITDATA.AA. The wavelengths for the background definition also ; need to be modified. ; outdata.offset=offset ngauss=fitdata.ngauss ; this is number of lines in fit FOR i=0,ngauss-1 DO outdata.aa[i*3+1,*,*]=outdata.aa[i*3+1,*,*]+fitdata.offset-offset outdata.x_bg1=fitdata.x_bg1+fitdata.offset-offset outdata.x_bg2=fitdata.x_bg2+fitdata.offset-offset ; ; The above OFFSET applies to FITDATA. If the user has used XBIN, ; YBIN, XRANGE or YRANGE to modify the original size of the raster, ; then OFFSET will be a different size to the data in the original ; WINDATA structure. The code below uses the newly-derived orbit ; correction to create a new offset array that matches the size of the ; original raster (this is referred to as FULL_OFFSET). Since the ; Y-variation of offset is determined ; entirely by the slit tilt then extending to the full length of the ; slit is trivial. For the X-direction, however, if XRANGE has been ; set then it is not possible to recover the orbit variation outside ; of XRANGE. The offsets in this region are simply set to MISSING. ; IF fitdata.nx NE fitdata.wd_nx OR fitdata.ny NE fitdata.wd_ny THEN BEGIN nx=fitdata.wd_nx ny=fitdata.wd_ny ; ; Initially set the full offset array to missing. Some of these ; missing values will remain if a sub-range in X was specified by the ; user in the call to eis_auto_fit ; full_offset=fltarr(nx,ny)+fitdata.missing ; ; ; Specifying refpix here is crucial for ensuring that the wavelength ; offsets in the full array will yield zero velocity in the sub-array ; specified by FITDATA. Note that I'm specifying ; fitdata.start_pix[1] - this is the YSTART input to eis_bin_windata, ; which is zero if eis_bin_windata was note used. ; ; In the code below, be careful with YRANGE (optional input to this ; routine) and FITDATA.YRANGE which was an optional input to ; eis_auto_fit and is different. ; IF n_elements(yrange) EQ 1 THEN BEGIN refpix=(mean([yrange,yrange+1]) + fitdata.yrange[0])*ybin+fitdata.start_pix[1] ENDIF ELSE BEGIN refpix=(mean([yrange[0],yrange[1]+1]) + fitdata.yrange[0])*ybin+fitdata.start_pix[1] ENDELSE refpix=round(refpix) print,'**REFPIX for full_offset: ',refpix ; tilt_corr=eis_slit_tilt_array(nx,ny,1,slit_ind,fitdata.yip,fitdata.date_obs, $ refpix=refpix,long=long) ; xr=fitdata.xrange FOR i=0,ny-1 DO full_offset[xr[0]:xr[1],i]=tilt_corr[xr[0]:xr[1],i]+yfit_all-refwvl ; ; The following checks for exposures with missing exposure ; times. These are flagged by anomalously high exp_start_times values ; k=where(fitdata.exp_start_times GT 43200.,nk) IF nk NE 0 THEN BEGIN FOR ik=0,nk-1 DO BEGIN j=k[ik] full_offset[j,*]=fitdata.missing ENDFOR ENDIF ; ; ; The full size offset is returned as the OFFSET optional output ; offset=full_offset ENDIF return,outdata END