;+ ; NAME ; ; EIS_AUTO_FIT_GEN ; ; PROJECT ; ; HINODE/EIS ; ; EXPLANATION ; ; A routine for automatically fitting multiple Gaussians to the ; window structures produced by eis_getwindata. ; ; INPUTS ; ; WINDATA A structure output by EIS_GETWINDATA() ; ; INITDATA A fit template structure ; ; OPTIONAL INPUTS ; ; WVLPIX By default eis_auto_fit_gen 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. ; ; 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. ; ; 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. ; ; KEYWORDS ; ; QUIET If set, then no information is printed to the screen. ; ; OUTPUTS ; ; FITDATA A structure with the following tags: ; .aa The final fit parameters for each spatial pixel. ; .sigmaa Errors on the fit parameters. ; .int Intensity array for each Gaussian. ; .interr Intensity error array for each Gaussian. ; .chi2 Chi^2 value for each spatial pixel. ; .offset Wavelength offset value for each spatial ; pixel. (Set to zero by eis_auto_fit_gen.) ; ; CALLS ; ; EIS_READ_FIT_TEMPLATE ; ; HISTORY ; ; Beta 1, 6-Nov-2008, Peter Young ; Beta 2, 16-Dec-2008, Peter Young ; modified INIT for individual pixels ; Beta 3, 17-Apr-2009, Mike Marsh & Peter Young ; corrected bug when computing the intensity errors ; Beta 4, 26-Oct-2009, Peter Young ; added pix_all keyword ; Beta 5, 5-Nov-2009, Paul Bryans & Peter Young ; changed pix_all to be an intarr (rather than bytarr) ; Beta 6, 20-Nov-2009, Paul Bryans ; passed parinfo information to mpfit ; Beta 7, 13-Dec-2009, Peter Young ; now checks that number of data points for a given fit exceeds ; the number of free parameters. ;- ;------------------- PRO eis_auto_fit_gen, windata, initdata, fitdata, wvlpix=wvlpix, $ xrange=xrange, yrange=yrange, $ pix_all=pix_all IF n_params() LT 3 THEN BEGIN print,'Use: IDL> eis_auto_fit_gen, windata, initdata, fitdata [, wvlpix=, ' print,' xrange=, yrange= , /pix_all]' return ENDIF IF NOT keyword_set(quiet) THEN BEGIN print,'' print,' EIS_AUTO_FIT_GEN was written by Peter Young (NRL/GMU).' print,' Please report any errors to pyoung@ssd5.nrl.navy.mil' print,'' ENDIF ; ; The input INITDATA can either be the structure returned by ; eis_read_fit_template, or it can be the name of a fit template file ; (which will be read by eis_read_fit_template). ; IF n_tags(initdata) EQ 0 THEN BEGIN chck=file_exist(initdata) IF chck EQ 1 THEN BEGIN initdata=eis_read_fit_template(initdata) ENDIF ENDIF nx=windata.nx ny=windata.ny wvl=windata.wvl int=windata.int err=windata.err miss=windata.missing exp_start_times=windata.time ; ; Set any missing pixels (flagged in the error array) to be missing in ; the intensity array. ; i=where(err EQ miss,ni) IF ni NE 0 THEN int[i]=miss nwvl=n_elements(wvl) IF n_elements(wvlpix) EQ 0 THEN wvlpix=initdata.wrange getmin=min(abs(wvl-wvlpix[0]),i0) getmin=min(abs(wvl-wvlpix[1]),i1) pix=[i0,i1] np=i1-i0+1 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] darr=fltarr(nx,ny)+miss ngauss=n_elements(initdata.lines) nback=initdata.background.nback nparams=n_elements(initdata.init) refwvl=fltarr(ngauss) FOR i=0,ngauss-1 DO refwvl[i]=initdata.lines[i].centroid fitdata={exp_start_times: exp_start_times, $ missing: miss, $ pix: pix, $ xrange: xrange, $ yrange: yrange, $ ngauss: ngauss, $ nback: nback, $ refwvl: refwvl, $ nx: nx, $ ny: ny, $ offset: fltarr(nx,ny), $ aa: fltarr(nparams,nx,ny), $ sigmaa: fltarr(nparams,nx,ny), $ int: fltarr(ngauss,nx,ny), $ interr: fltarr(ngauss,nx,ny), $ chi2: fltarr(nx,ny), $ pix_all: intarr(nwvl)} ; ; 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_GEN: 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 expr='' FOR i=0,ngauss-1 DO BEGIN expr=expr+'gauss_sg(x,p['+trim(i*3)+':'+trim(i*3+2)+']) + ' ENDFOR IF nback EQ 2 THEN BEGIN expr=expr+'line_sg(x,p['+trim(ngauss*3)+':'+trim(ngauss*3+1)+'])' ENDIF ELSE BEGIN expr=expr+'p['+trim(ngauss*3)+']' ENDELSE FOR i=xrange[0],xrange[1] DO BEGIN FOR k=yrange[0],yrange[1] DO BEGIN IF pix_all[0] EQ -1 THEN BEGIN xx=wvl[pix[0]:pix[1]] yy=int[pix[0]:pix[1],i,k] ee=err[pix[0]:pix[1],i,k] ENDIF ELSE BEGIN ipix=where(pix_all EQ 1) xx=wvl[ipix] yy=int[ipix,i,k] ee=err[ipix,i,k] ENDELSE ;; IF n_elements(pix) NE 0 THEN BEGIN ;; ENDIF ELSE BEGIN ;; xx=wvl ;; yy=reform(int[*,i,k]) ;; ee=reform(err[*,i,k]) ;; ENDELSE ; IF max(yy) GT 3.*max(ee) THEN BEGIN j=where(yy NE miss,nj) IF nj GT nparams THEN BEGIN ; init=initdata.init ; ; Scale the peak values stored in INIT ; using the maximum of the current spectrum. ; i_init=indgen(ngauss)*3 peaks=init[i_init] init[i_init]=peaks/max(peaks)*max(yy) ; aa = MPFITEXPR(expr, xx[j], yy[j], ee[j], init, $ perr=sigmaa, /quiet, bestnorm=bestnorm,yfit=yfit, $ parinfo=initdata.parinfo) i0=i-xrange[0] k0=k-yrange[0] fitdata.aa[*,i0,k0]=aa fitdata.sigmaa[*,i0,k0]=sigmaa fitdata.chi2[i0,k0]=bestnorm/(n_elements(j) - n_elements(aa)) ENDIF ENDIF ENDFOR IF NOT keyword_set(quiet) THEN BEGIN ncol=xrange[1]-xrange[0]+1 print,'Completed column '+trim(i+1-xrange[0])+' of '+trim(ncol) ENDIF ENDFOR ; ; Calculate the intensity and intensity error using the peak and width ; fit parameters. ; FOR i=0,ngauss-1 DO BEGIN i1=i*3 ; index of line peak in aa i2=i*3+2 ; index of line width in aa ; fitdata.int[i,*,*]=fitdata.aa[i1,*,*]*fitdata.aa[i2,*,*]*sqrt(2.*!pi) fitdata.interr[i,*,*]=fitdata.int[i,*,*]*sqrt(0.5*((fitdata.sigmaa[i1,*,*]/fitdata.aa[i1,*,*])^2. + (fitdata.sigmaa[i2,*,*]/fitdata.aa[i2,*,*])^2. )) ;; fitdata.interr[i,*,*]=max( [fitdata.sigmaa[i1,*,*]*fitdata.aa[i2,*,*], $ ;; fitdata.aa[i1,*,*]*fitdata.sigmaa[i2,*,*] ] ) $ ;; * sqrt(2.*!pi) ENDFOR END