function voigtfit,x,y,a,verbose=verbose,sigma=sigma,weight=w,plot=plot, $ derivative=derivative,emission=emission,guess=guess, $ vpguess=vpguess,doppguess=doppguess,quiet=quiet,chi2=chi2, $ goodfit=goodfit,nofitvp=nofitvp,nofitdopp=nofitdopp common voigt_fit_private,npoints,nterms,npder,one_over_sqrt_pi,nofita,voigtp,nofitd,doppl on_error,2 if keyword_set(nofitvp) then nofita=1 else nofita=0 if keyword_set(nofitdopp) then nofitd=1 else nofitd=0 if n_elements(vpguess) NE 1 then vpguess=0.1 one_over_sqrt_pi = 1.0d0/sqrt(!pi) nterms = 5L npoints = n_elements(x) npder = nterms*npoints if n_elements(y) NE npoints then $ message,'x and y must have the same dimension' if n_elements(guess) NE nterms then begin if keyword_set(verbose) and NOT keyword_set(quiet) then begin if keyword_set(emission) then $ message,/info,'Deriving initial guess ... assuming emission line' $ else $ message,/info,'Deriving initial guess ... assuming absorption line' endif miny = min(y,minplace) maxy = max(y,maxplace) guess = fltarr(nterms) if keyword_set(emission) then guess(0) = miny $ ; continuum else guess(0) = maxy if n_elements(doppguess) NE 1 then begin ; Doppler width is half the half width (assumes a full spectral line) if keyword_set(emission) then $ halfmax = where(y GT (miny+maxy)/2.0,nhalfmax) $ else $ halfmax = where(y LT (miny+maxy)/2.0,nhalfmax) if nhalfmax GE 1 then $ guess(3) = 0.5*abs(x(halfmax(1))-x(halfmax(nhalfmax-2))) $ else message,'Cant guess doppler width' endif else guess(3)=doppguess guess(2) = vpguess ; damping parameter ... first guess hardwired (yuk) voigt,guess(2),0.,h0,f0 if keyword_set(emission) then begin guess(1) = (guess(0)-maxy)/h0 ; strength of line guess(4) = x(maxplace) ; line center endif else begin guess(1) = (guess(0)-miny)/h0 ; strength of line guess(4) = x(minplace) ; line center endelse endif voigtp = guess(2) doppl = guess(3) if n_elements(w) NE npoints then w = fltarr(npoints)+1.0 a = guess q=!quiet ; Don't print the convergence failure messages if NOT keyword_set(verbose) then !quiet=1 fit = curvefit(x,y,w,a,sigma,function_name='voigt_funct') !quiet=q ; Check that we got a good fit ; Guess that errors are sqrt(y) to compute chi squared chi2 = total(w*(y-fit)^2/abs(y)) chi_limit = chi_sqr(0.001,npoints-nterms) goodfit = chi2 LE chi_limit if (NOT goodfit) AND (NOT keyword_set(quiet)) then $ message,/info,strcompress('Questionable fit '+string(chi2)+' > '+ $ string(chi_limit)) if keyword_set(derivative) then begin v = (x-a(4))/a(3) voigt,a(2),v,h,f dIdx = -2.0*a(1)*(a(2)*f-v*h)/a(3) endif if keyword_set(plot) then begin plot,x,y oplot,x,fit,thick=2 endif if keyword_set(verbose) and NOT keyword_set(quiet) then begin print,chi2,chi_limit,format="(' Chi**2:',f9.2,' Limit:',f8.2)" print,a(0),guess(0),format="('Continuum:',f9.2,' from:',f9.2)" print,a(1),guess(1),format="(' Strength:',f9.2,' from:',f9.2)" print,a(2),guess(2),format="('Damping P:',f9.2,' from:',f9.2)" print,a(3),guess(3),format="('Doppler W:',f9.2,' from:',f9.2)" print,a(4),guess(4),format="('L. Center:',f9.2,' from:',f9.2)" endif if keyword_set(derivative) then return,dIdx else return,fit end