FUNCTION eis_shift_spec, windata, offset, linear=linear, $ quiet=quiet ;+ ; NAME ; ; EIS_SHIFT_SPEC() ; ; PROJECT ; ; Hinode/EIS ; ; EXPLANATION ; ; When comparing two spectra within the WINDATA structure, there ; will be a wavelength offset between them that arises from the ; slit tilt and orbit variation. The wavelength offsets are stored ; in the WINDATA.WAVE_CORR tag (and also can be manually input ; through the OFFSET keyword). EIS_SHIFT_SPEC takes ; a WINDATA structure and the wavelength offset array and shifts each ; spectrum within WINDATA to lie on a reference wavelength ; scale. This requires data interpolation for both the intensity ; and error arrays. ; ; The user is *strongly* recommended to use the /refill option for ; eis_getwindata before using eis_shift_spec. ; ; Please check the programming notes section below for information ; about options for interpolation. Generally cubic interpolation ; (the default) is recommended but if there are a lot of missing ; pixels it may be safer to use linear interpolation. ; ; Note that, after running eis_shift_spec, the output WINDATA will ; have WINDATA.WAVE_CORR_SET=0 (normally this will be 1). This ; allows the user to check if a WINDATA structure has already been ; processed with EIS_SHIFT_SPEC. ; ; ; INPUTS ; ; WD A windata structure produced by EIS_GETWINDATA. It is ; recommended that the /REFILL keyword be set in the call ; to EIS_GETWINDATA. ; ; OPTIONAL INPUTS ; ; OFFSET An array of wavelength offsets of same size as the spatial ; dimensions of WD. The routine will take this array from ; the WD.WAVE_CORR by default, and so OFFSET can be used to ; over-ride the default value. ; ; KEYWORDS ; ; LINEAR If set, then linear interpolation is used (cubic ; interpolation is the default). ; ; QUIET If set, then no information messages will be printed to ; the screen. ; ; OUTPUTS ; ; A windata structure is returned that is identical to the input ; WD, except that the following tags are changed: ; .wvl Contains the wavelength array of the original windata. ; .int Contains the new intensity array, aligned to the ; wavelength scale of the reference pixel. ; .err Contains the new intensity array, aligned to the ; wavelength scale of the reference pixel. ; .wave_corr_set is set to 0. ; ; Note that the spectra are interpolated onto the windata.wvl ; wavelength array. ; ; EXAMPLES ; ; One example for using this routine is in conjunction with ; EIS_MASK_SPECTRUM which averages spectra over a ; user-specified pixel mask. If the spatial extent of the mask is ; large then wavelength offsets will tend to smear the ; spectra, making the lines appear broader. By first using ; EIS_SHIFT_SPEC to correct for the offsets, then the smearing can ; be prevented. ; ; PROGRAMMING NOTES ; ; Cubic interpolation is implemented through the interpolate.pro ; built-in IDL routine. The option cubic=-0.5 is used, as ; recommended in the IDL user guide. ; ; Note that cubic interpolation interpolates over 4 pixels to ; derive the new value for a particular pixel. If one pixel (of the ; 4) has a ; missing value, then this missing value is treated as real data ; and goes into the interpolation. If the /refill keyword has been ; used, then the number of missing pixels is very small and ; so this isn't a big problem. For this reason I print out ; the fraction of missing pixels in the data-set. ; ; For linear interpolation, only two pixels are used to derive the ; interpolated value, and I correctly keep track of missing pixels, ; so if the interpolated intensity has been derived partly from a ; missing pixel, then I just set the new intensity to missing (see ; the outflag array in the code). ; ; HISTORY ; ; Ver.1, 14-Dec-2009, Peter Young ; Ver.2, 17-Dec-2009, Peter Young ; default is now to use spline interpolation; linear ; interpolation is selected with the /linear keyword ; Ver.3, 19-Mar-2010, Peter Young ; updated header ; Ver.4, 17-May-2011, Peter Young ; routine now checks if the WD.WAVE_CORR tag exists. If yes, then ; the offset array is taken from this unless it has been ; given by the user. ; Ver.5, 20-May-2011, Peter Young ; fixed bug introduced in Ver. 4 ; Ver.6, 28-Jun-2011, Peter Young & Kuen Ko ; now deals correctly with data-sets that have more than one exposure ; per raster position (nexp_prp>1); optional input REFPIX has ; been made obsolete. ; Ver.7, 20-Aug-2012, Peter Young ; now sets the tag 'wave_corr_set' to 0 (should wave_corr also be ; set to 0's?). ; Ver.8, 17-Oct-2013, Peter Young ; the routine has been completely re-written in order to use the ; IDL interpolate routine which makes the code much faster. ;- IF n_params() LT 1 THEN BEGIN print,'Use: IDL> wdout = eis_shift_spec( wd [, offset, /linear, /quiet ])' print,'' print,' This routine uses cubic interpolation by default. If the windata structure has a lot of missing data, ' print,' then it may be better to use linear interpolation (/linear). Note that linear interpolation tends to' print,' lead to a slight broadening of the emission lines.' return,-1 ENDIF IF n_elements(offset) EQ 0 THEN BEGIN IF tag_exist(windata,'wave_corr') EQ 1 THEN BEGIN offset=windata.wave_corr ENDIF ELSE BEGIN print,'%EIS_SHIFT_SPEC: the windata structure does not contain the WAVE_CORR tag. Please either re-generate the ' print,'windata structure with eis_getwindata, or directly specify a wavelength offset array through the OFFSET input.' print,'Returning...' return,windata ENDELSE ENDIF n=n_elements(windata.int) k=where(windata.err EQ windata.missing,nk) ; frac_miss=double(nk)/double(n) perc_str=string(format='(f8.4,"%")',frac_miss*100.) IF NOT keyword_set(quiet) THEN print,'% EIS_SHIFT_SPEC: fraction of missing pixels is '+trim(perc_str) wd2=windata nl=windata.nl nx=windata.nx ny=windata.ny np=windata.hdr.nexp_prp ; ; work out average size of wavelength bin ; dlambda=mean(windata.wvl[1:*]-windata.wvl[0:nl-2]) ; ; ix and iy are 2D index arrays giving the pixel indices for the ; interpolated points. The spectra are assumed to be distributed on a ; wavelength scale of fixed spacing (dlambda), thus the wavelength ; pixels can be assigned integers 0, 1, 2, etc. For a given spectrum, ; if the wave_corr value is 0.05 A, for example, then this is ; converted to a pixel shift (dp) by dividing by dlambda. Thus the new ; spectrum is desired for pixel values of 0+dp, 1+dp, 2+dp, etc. ; ; The structures below are simply a device for efficiently creating ; the 2D ix and iy arrays. ; s={ ix: findgen(nl) } s2=replicate(s,ny) ix0=s2.ix s={ iy: findgen(ny) } s2=replicate(s,nl) iy=transpose(s2.iy) FOR j=0,np-1 DO BEGIN FOR i=0,nx-1 DO BEGIN img=reform(windata.int[*,i,*,j]) err=reform(windata.err[*,i,*,j]) flag=fltarr(nl,ny) k=where(err EQ windata.missing,nk) IF nk NE 0 THEN flag[k]=1.0 ; s={dp: offset[i,*]/dlambda} s2=replicate(s,nl) dp=transpose(s2.dp) ix=ix0+dp ; ; Note how the outflag array is used to keep track of missing data ; with the linear option. For the cubic option it's not ; possible to do this so any missing data will be included in the ; interpolation. ; IF keyword_set(linear) THEN BEGIN outimg=interpolate(img,ix,iy,missing=windata.missing) outflag=interpolate(flag,ix,iy) outerr=interpolate(err,ix,iy,missing=windata.missing) k=where(outflag NE 0.) IF nk NE 0 THEN BEGIN outimg[k]=windata.missing outerr[k]=windata.missing ENDIF ENDIF ELSE BEGIN outimg=interpolate(img,ix,iy,missing=windata.missing,cubic=-0.5) outerr=interpolate(err,ix,iy,missing=windata.missing,cubic=-0.5) ENDELSE ; wd2.int[*,i,*,j]=outimg wd2.err[*,i,*,j]=outerr ; ENDFOR ENDFOR ; ; This "un-sets" the wave_corr flag so that other routines can see ; that the wave_corr tag is no longer needed. ; wd2.wave_corr_set=0b return,wd2 END