PRO eis_auto_fit_old, windata, fitdata, refwvl=refwvl, wvlpix=wvlpix, $ phot_limit=phot_limit, xrange=xrange, yrange=yrange, $ uniform_backg=uniform_backg, pix_all=pix_all, $ interp_miss=interp_miss, miss_map=miss_map, $ pry_interp=pry_interp, keep=keep, frac=frac ;+ ; NAME ; ; EIS_AUTO_FIT_OLD ; ; PROJECT ; ; HINODE/EIS ; ; EXPLANATION ; ; **This is the old version of eis_auto_fit. It was renamed ; **eis_auto_fit_old.pro in Jan 2010. ; ; This routine takes the output from EIS_GETWINDATA and performs a fit ; over the intensity array, returning values for the intensity, line width ; and centroid. ; ; By default the routine ignores all missing pixels. To use the ; interpolated intensities for the missing pixels in the fitting, ; see the keywords /INTERP_MISS and /PRY_INTERP. Generally ; /PRY_INTERP gives better results. ; ; The routine goes through the image array, pixel-by-pixel, fitting ; each profile. This is a slow process in IDL, and it becomes even ; slower if the data are poor quality or the initial parameters are ; not good. To help the routine run quickly through the fits, the ; following features are implemented: ; - For each spatial pixel, intensity must be GE 1.5*error for at ; least one wavelength pixel. ; - For each spatial pixel, at least half of the pixels being fit ; must not be missing pixels. ; - The initial (Gaussian) width is set to 0.03 angstroms, ; corresponding to FWHM=0.07 angstroms. This is a typical width ; for EIS lines. ; - For each spatial pixel, the initial line peak is set to the ; difference between the minimum and maximum intensities in the ; spectrum (ignoring missing data). ; - For each spatial pixel, the initial line centroid is set to the ; wavelength pixel where the initial line peak is found (see ; above). ; - For each spatial pixel, the initial guess to the background is ; set to the minimum value in the spectrum (ignoring missing ; data). ; ; After fitting has been performed, the parameters can be rejected ; if any one of the following conditions holds: ; - The derived line centroid lies outside of the wavelength range ; for the fit. ; - The derived Gaussian width is > 2.0 angstroms (corresponding to ; FWHM > 4.7 angstroms). ; - The derived Gaussian peak is > factor 3 larger than the initial ; guess for the Gaussian peak. ; ; If a spatial pixel is ignored because of any one of the conditions ; described above, then the fit parameters are set to the value of ; FITDATA.MISSING. ; ; INPUTS ; ; WINDATA Structure output by EIS_GETWINDATA. ; ; OUTPUTS ; ; FITDATA A structure containing the fit details. The tags are: ; .int Intensity array. Units: erg/cm^2/s/sr. ; .interr 1-sigma error on intensity. Units: erg/cm^2/s/sr. ; .wid Gaussian width. Units: angstroms. Note: multiply by ; 2.*sqrt(2.*alog(2.)) to get FWHM. ; .widerr 1-sigma error on Gaussian width. Units: ; angstroms. ; .cen Gaussian centroid. Units: angstroms. ; .cenerr 1-sigma error on centroid. Units: angstroms. ; .vel Line velocity relative to REFWVL. Units: km/s. ; .velerr 1-sigma error on velocity. Units: km/s. ; .chi2 Chi^2 value for fit. ; .exp_start_times The start time of each exposure relative ; to the first exposure. Units: ; seconds. This is the same as the .time ; tag of WINDATA. ; .missing The 'missing' value for the data. Copied from ; WINDATA. ; .refwvl Contains the input REFWVL. Units: angstroms. ; .pix Two element array containing the pixel indices of ; the minimum and maximum wavelengths included in ; the fit. ; .pix_all Contains the input PIX_ALL, if specified. If ; not, then it is an array of same size as ; WINDATA.WVL, with elements set to -1. ; .bg1 The fitted background level at wavelength pixel ; FITDATA.PIX[0]. If /uniform_backg has been set, ; then BG1 contains the uniform background ; value. Units: erg/cm2/s/sr/angstrom. ; .bg2 The fitted background level at wavelength pixel ; FITDATA.PIX[1]. If /uniform_backg has been set, ; then BG2 is set to MISSING. Units: ; erg/cm2/s/sr/angstrom. ; .xrange Contains the contents of the input XRANGE. ; .yrange Contains the contents of the input YRANGE. ; .offset This is an empty array that will be filled by ; the orbit and slit tilt correction routines. ; .uniform_backg Set to 1 or 0 depending if /uniform_backg ; has been set or not. ; .interp_miss Set to 1 or 0 to indicate if /interp_miss ; has been set. ; ; ; OPTIONAL INPUTS ; ; WVLPIX By default eis_auto_fit assumes the entire wavelength range in ; the data array will be used for fitting. The range can be reduced ; by specifying the 2 element array WVLPIX. E.g., setting ; WVLPIX=[194.9,195.4] will use only the wavlength range 194.9 ; to 195.4 angstroms (if the Fe XII 195 line is being fit). See ; also the PIX_ALL optional input. ; ; PIX_ALL This allows a more general method of reducing the number of ; wavelength pixels involved in the fit compared to the WVLPIX ; input described above. PIX_ALL must be a vector of same size ; as the number of wavelengths in the input array. The values ; of the vector elements should be 0 or 1. 1 indicates the pixel ; will be included in the fit; 0 that the pixel will not be ; included. PIX_ALL can be useful if suitable areas of background ; are not directly adjacent to the emission line. ; ; REFWVL The reference wavelength of the line being fitted. This allows ; the velocity shift to be computed. ; ; XRANGE Allows a sub-range of the X-values to be fitted. E.g., setting ; equal to [10,30] means that only X-values from 10 to 30 ; (inclusive) are fit. ; ; YRANGE Allows a sub-range of the Y-values to be fitted. E.g., setting ; equal to [10,30] means that only Y-values from 10 to 30 ; (inclusive) are fit. ; ; OPTIONAL OUTPUTS ; ; MISS_MAP A 2D array of same size as FITDATA.INT that contains ; the number of missing pixels in the spectrum used for ; that fit at each spatial pixel. If PIX_ALL or WVLPIX have ; been set, then the number of missing pixels is those ; within the reduced spectrum. ; ; KEYWORDS ; ; UNIFORM_BACKG By default, the routine performs a linear fit to the ; spectrum background. By setting this keyword the background ; is assumed flat. ; ; INTERP_MISS By default, eis_auto_fit ignores missing data when ; performing fits. By setting /INTERP_MISS the ; interpolated intensity values inserted by eis_prep ; at the missing pixels will be used in the fit. The ; error bars on these interpolated intensity values ; are obtained by performing a linear fit to the ; intensity vs. error^2 plot for good pixels, and then ; using the interpolated intensity value to give the ; error estimate. ; ; PRY_INTERP Use the routine EIS_REPLACE_MISSING to interpolate ; missing data points. ; ; KEEP By default, the routine modifies the INT and ERR ; tags of WINDATA, replacing missing data with ; interpolated intensities and derived errors. If you ; prefer not to modify WINDATA, then use /KEEP. ; ; RESTRICTIONS ; ; Only performs a single Gaussian fit. ; ; CALLS ; ; GAUSS_SG, LINE_SG, MPFITEXPR, LAMB2V ; ; HISTORY ; ; 0.1, 16-Jan-2007, Peter Young ; 0.2, 18-Feb-2007, Peter Young ; corrected input to windata.fit tag ; 0.3, 22-Feb-2007, Peter Young ; I've changed the arrays to floats to save on memory ; 0.4, 10-Mar-2007, Peter Young ; Added XRANGE= and YRANGE= keywords. ; 0.5, 4-Apr-2007, Peter Young ; Changed check on photon counts to check on error (thus making photon ; array redundant). ; 0.6, 13-Jun-2007, Peter Young & Durgesh Tripathi ; Make compatible with output from eis_getwindata; added bg1 and bg2 ; output tags. ; 0.7, 19-Jun-2007, Peter Young ; Corrected minor bug; routine now assumes intensity in per angstrom ; units. ; 0.8, 10-Jul-2007, Peter Young ; Corrected bug in format statement. ; 0.9, 18-Oct-2007, Peter Young ; Following suggestion by Tom Van Doorsselaere, have given the output ; structure a name (allows structures to be concatenated easily). ; 1.0, 27-Nov-2007, Peter Young ; Added /uniform_backg keyword. ; 1.1, 16-Jan-2008, Peter Young ; Added PIX_ALL optional input ; 1.2, 30-Jan-2008, Peter Young ; Corrected bug: the intensity wasn't being calculated correctly. ; Instead of multiplying by sqrt(2*!pi), the factor used had been ; 1.06. ; 1.3, 7-Jan-2008, Peter Young ; Removed the structure name (v.0.9) as it was causing problems; ; also forced refwvl to be float when creating structure. ; 1.4, 13-Feb-2008, Peter Young ; Corrected calculation of pix when pix_all has been set. ; 1.5, 6-Nov-2008, Peter Young ; Added OFFSET tag to eis_auto_fit (not actually used by ; eis_auto_fit). ; 1.6, 20-Nov-2008, Peter Young ; OFFSET was being set to missing by default, I've now ; changed this to zero. ; 1.7, 11-Feb-2009, Peter Young ; Expanded header. Routine unchanged. ; 1.8, 25-Feb-2009, Peter Young ; Added /INTERP_MISS and MISS_MAP= keywords. The method for ; deriving the error bars for interpolated data is from ; H.P. Warren. Added interp_miss tag to FITDATA. ; 1.9, 10-Mar-2009, Peter Young ; Added new interpolation method (call to eis_replace_missing) ; and added /KEEP keyword. ; 2.0, 2-Feb-2010, Peter Young ; Renamed eis_auto_fit_old ;- IF n_params() LT 2 THEN BEGIN print,'Use: IDL> eis_auto_fit_old, windata, fitdata, pix_all=, refwvl=, wvlpix=' print,' xrange=, yrange=, /uniform_backg, /interp_miss' print,' miss_map=, /keep' return ENDIF print,'' print,' EIS_AUTO_FIT_OLD was written by Peter Young (NRL/GMU).' print,' Please report any errors to pyoung@ssd5.nrl.navy.mil' print,'' nx=windata.nx ny=windata.ny wvl=windata.wvl int=windata.int err=windata.err miss=windata.missing ; ; If /INTERP_MISS is set, then the error values for missing data ; (which are set by eis_prep to have the value WD.MISSING) must be set ; to a realistic value. ; ; The realistic value is obtained by plotting err^2 vs. int for the ; good data, and fitting a straight line. Using this line, one can ; take the interpolated intensity values for the missing pixels, and ; convert them to an error estimate. ; ; If /INTERP_MISS is not set, then the interpolated intensity values ; need to be replaced with missing values. ; CASE 1 OF keyword_set(pry_interp): BEGIN IF keyword_set(keep) THEN replace=0 ELSE replace=1 eis_replace_missing,windata,replace=replace,wdout=int,errout=err,frac=frac,/quiet END keyword_set(interp_miss): BEGIN i=where(err NE miss) x=int[i] y=err[i]^2 f=poly_fit(x,y,1) ; ; replace missing data with error estimate ; i=where(err EQ miss) err[i]=sqrt(int[i]*f[1]+f[0]) END ELSE: BEGIN i=where(err EQ miss,ni) IF ni NE 0 THEN int[i]=miss interp_miss=0 END ENDCASE nwvl=n_elements(wvl) ; ; DEFINE PIX ; ---------- IF n_elements(wvlpix) NE 0 THEN BEGIN getmin=min(abs(wvl-wvlpix[0]),i0) getmin=min(abs(wvl-wvlpix[1]),i1) pix=[i0,i1] ENDIF ; IF n_elements(pix) EQ 0 THEN pix=[0,nwvl-1] exp_start_times=windata.time IF n_elements(xrange) NE 0 THEN BEGIN nx=xrange[1]-xrange[0]+1 exp_start_times=exp_start_times[xrange[0]:xrange[1]] ENDIF ELSE BEGIN xrange=[0,nx-1] ENDELSE IF n_elements(yrange) NE 0 THEN ny=yrange[1]-yrange[0]+1 ELSE yrange=[0,ny-1] IF n_elements(refwvl) EQ 0 THEN BEGIN refwvl=average(wvl[pix[0]:pix[1]]) print,'%EIS_AUTO_FIT: reference wavelength not specified.' print,' Velocities derived using refwvl='+ $ trim(string(refwvl,format='(f10.3)')) ENDIF IF NOT keyword_set(uniform_backg) THEN uniform_backg=0 ; ; DEFINE OUTPUT STRUCTURE, FITDATA ; -------------------------------- darr=fltarr(nx,ny)+miss fitdata={ int: darr, $ interr: darr, $ cen: darr, $ cenerr: darr, $ wid: darr, $ widerr: darr, $ vel: darr, $ velerr: darr, $ chi2: darr, $ exp_start_times: exp_start_times, $ missing: miss, $ refwvl: float(refwvl), $ pix: pix, $ pix_all: intarr(nwvl), $ bg1: darr, $ bg2: darr, $ xrange: float(xrange), $ yrange: float(yrange), $ offset: fltarr(nx,ny), $ uniform_backg: uniform_backg} miss_map=intarr(nx,ny) ; ; CHECK PIX_ALL INPUT ; pix_all is assumed to be a vector of same size as number of wavelengths ; in WINDATA. The values of each element are either 0 or 1. 0 indicates ; the pixel will not be included in fit. 1 indicates the pixel will be ; included in fit. ; ; If pix_all was not specified as input, then each element is set to -1. ; IF n_elements(pix_all) EQ 0 THEN BEGIN fitdata.pix_all=-1 pix_all=fitdata.pix_all ENDIF ELSE BEGIN IF n_elements(pix_all) NE nwvl THEN BEGIN print,'%EIS_AUTO_FIT: number of elements of PIX_ALL input does not match the wavelength' print,' dimension of data.' print,' Check the routine header for more details.' print,' No. elements of pix_all: '+trim(n_elements(pix_all)) print,' Wavelength dimension of windata: '+trim(nwvl) return ENDIF ELSE BEGIN fitdata.pix_all=pix_all ipix=where(pix_all EQ 1) fitdata.pix=[min(ipix),max(ipix)] ENDELSE ENDELSE FOR i=yrange[0],yrange[1] DO BEGIN FOR k=xrange[0],xrange[1] DO BEGIN ; ; Load up XX, YY and EE arrays that will be used for fitting. ; There are two cases depending on whether PIX_ALL has been set. ; IF pix_all[0] EQ -1 THEN BEGIN xx=wvl[pix[0]:pix[1]] yy=int[pix[0]:pix[1], k, i] ee=err[pix[0]:pix[1], k, i] ENDIF ELSE BEGIN ipix=where(pix_all EQ 1) xx=wvl[ipix] yy=int[ipix,k,i] ee=err[ipix,k,i] ENDELSE IF max(yy) GT 1.5*max(ee) THEN BEGIN ; j=where(yy NE miss,nj) ; ; If more than half the data-points are missing, then don't try a fit. ; IF nj/float(n_elements(yy)) GT 0.5 THEN BEGIN ; ; Calculate initial fit parameters ; -------------------------------- getmax=max(yy[j],imax) init=[max(yy[j])-min(yy[j]),xx[j[imax]],0.03,min(yy[j]),min(yy[j])] IF keyword_set(uniform_backg) THEN init=init[0:3] ; ; Define fit function ; ------------------- IF keyword_set(uniform_backg) THEN BEGIN expr='gauss_sg(x,p[0:2]) + p[3]' ENDIF ELSE BEGIN expr='gauss_sg(x,p[0:2]) +line_sg(x,p[3:4])' ENDELSE ; ; Perform fit ; ----------- aa = MPFITEXPR(expr, xx[j], yy[j], ee[j], init, $ perr=sigmaa, /quiet, bestnorm=bestnorm,yfit=yfit) ; ; ; The following restrictions aim to prevent bad fits going into the ; output structure ; IF aa[1] LE max(xx[j]) AND aa[1] GE min(xx[j]) AND aa[2] LE 2.0 AND aa[0] LE 3.*(max(yy[j])-min(yy[j])) THEN BEGIN ; ; Calculate intensity ; ------------------- ; k0=k-xrange[0] i0=i-yrange[0] fitdata.int[k0,i0]=aa[0]*aa[2]* sqrt(2.*!pi) fitdata.interr[k0,i0]=max([aa[0]*sigmaa[2],aa[2]*sigmaa[0]])*sqrt(2.*!pi) ; fitdata.cen[k0,i0]=aa[1] fitdata.cenerr[k0,i0]=sigmaa[1] fitdata.wid[k0,i0]=abs(aa[2]) fitdata.widerr[k0,i0]=sigmaa[2] fitdata.chi2[k0,i0]=bestnorm/(n_elements(j)-n_elements(aa)) fitdata.bg1[k0,i0]=aa[3] IF NOT keyword_set(uniform_backg) THEN fitdata.bg2[k0,i0]=aa[4] ; IF n_elements(refwvl) NE 0 THEN BEGIN fitdata.vel[k0,i0]=lamb2v(aa[1]-refwvl,refwvl) fitdata.velerr[k0,i0]=lamb2v(sigmaa[1],refwvl) ENDIF ; j=where(yy EQ miss,n_miss) miss_map[k0,i0]=n_miss ENDIF ENDIF ENDIF ENDFOR ENDFOR END