; ; voi and voi2 are auxiliary functions used by getew_xy.pro ; pro voi,x,b,f common cte,center,height,continuum f=b(0)*voigt(b(1),(x-center)*center^2/b(2)) + continuum end ; pro voi2,x,b,f common cte,center,height,continuum f=b(0)*voigt(b(1),(x-center)*center^2/b(2)) + b(3) end ; pro getew_xy,x,y,ew,a=a,b=b,noplot=noplot,cont=cont,restrict=restrict,\$ yerror=yerror,_extra=e ;+ ; ; Automatic determination of the EW of a spectral line by fitting ; Gaussian and Voigt profiles. The fits are forced to match the continuum ; at 1.0 and, in the case of a Voigt profile, the center of the line at ; is fixed to the position determined by the Gaussian fit. ; ; IN: x - flt/dblarr Vector with wavelengths ; y - flt/dblarr Vector with normalized flux ; ; OUT: ew - dblarr(2) Equivalent width determined from the Gaussian ; and the Voigt fits (angstroms). A zero value is returned ; when the measurement could not be performed. ; ; KEYWORDS: a float Gaussian parameters ; f=a(0)*exp(-((x-a(1))/a(2))^2./2) + a(3) ; b flt array(2) Voigt profile parameters ; f=b(0)*voigt(b(2),(x-b(1))*b(1)^2/b(3))+b(4) ; ; cont when switched on, the continuum location ; becomes a fitting parameter (otherwise ; fixed to 1.0) ; ; restrict when on, the initial search for the ; line center (mideshift_xy) is restricted ; to 'restrict' from the mean of x ; (same units as in x) ; ; yerror std. deviation between data and model in the area ; used for least-squares fitting. As ew, this is a 2-element ; array, with the first element for the Gaussian fit and the ; second for the Voigt fit ; ; Extra keywords are passed along to plot (e.g. yrange, charsize, etc.) ; ; ; Note: the numerical integration of the fits goes on only out to the ; limits of the input wavelength window. The only points taken ; into account in the fit are those from the line center up to ; the first place where the slope changes sign (to avoid blends). ; ; The routine curvefit.pro, which is used by this procedure, ; will sometimes crash when there are not enough points to fit. ; A modified curvefit.pro (see callende library) that will exit ; more graciously is recommended for batch applications. ; ; Carlos Allende Prieto, UT, Oct 2001 ; , UT, April 2006, Voigt fitting enabled, 'height' ; and 'continuum' are included as fitting ; parameters, and so is 'center' for the ; Gaussian case. Keywords 'cont' and ; 'restrict' added ; , UT, February 2009, yerror keyword added ; , IAC, July 2011, accepted extra keywords for 'plot' ; ;- itmax=100 common cte,center,height,continuum npar = n_params() if (npar eq 0) then begin print,'use -- getew_xy,x,y,ew[,a=a,b=b,cont=cont,restrict=restrict,yerror=yerror,_extra=e]' return endif ew=dblarr(2) if keyword_set(restrict) then w=where(abs(x-mean(x)) le restrict) else \$ w=indgen(n_elements(x)) mideshift_xy,x[w],y[w],xmin,ymin,xsig,/noplot,/int,/nois if (xmin eq -1000) then begin print,'% GETEW_XY: no minimum could be measured by mideshift' return endif height=-(1.-ymin) & center=xmin & continuum=1.00 if not keyword_set(noplot) then begin plotsym,0,0.6,/fill plot,x,y,xr=[min(x),max(x)],yr=[min(y)*.9,1.03],ystyl=1,xstyle=1,psy=8,\$ xtitle='!3Wavelength ('+string("305B)+')',ytitle='Normalized flux',charsiz=1.4,_extra=e xyouts,0.7,0.03,systime(0),/norm,charsize=1. endif ; cleaning - from the line center we go up only until the slopes change ; sign (leaving a margin from the center) margin=4e-6 blue_limit=max(where(x lt xmin*(1.-margin) and deriv(x,y) ge 0.0)) red_limit=min(where(x gt xmin*(1.+margin) and deriv(x,y) le 0.0)) if (blue_limit(0) eq -1) then blue_limit=0 if (red_limit(0) eq -1) then red_limit=n_elements(x)-1 sy=y(blue_limit:red_limit) sx=x(blue_limit:red_limit) ;oplot,replicate(xmin*(1.-margin),2),[0,2] ;oplot,replicate(xmin*(1.+margin),2),[0,2] if keyword_set(cont) then begin f=gaussfit(sx,sy,a,nter=4,yerror=yerror1) gaussianmia,x,a,f ew(0)=int_tabulated(x,a[3]-f) endif else begin f=gaussfit(sx,1.-sy,a,nter=3,yerror=yerror1) gaussianmia,x,[a,0.0],f ew(0)=int_tabulated(x,f) f=1.-f a=[-a(0),a(1:2),continuum] endelse ;update value of the center, which will be FIXED ;for the Voigt fitting center=a(1) if not keyword_set(noplot) then begin oplot,sx,sy,thick=7,col=30 oplot,x,f,col=80,thick=3 endif weights=replicate(1.0d0,n_elements(sx)) if keyword_set(cont) then begin b=dblarr(4) b(0)=height & b(1)=0.034237377 & b(2)=2311136.7 & b(3)=1.00 r=curvefit(sx,sy,weights,b,sigma,function_name='voi2',/noder,itmax=itmax,iter=it,yerror=yerror2) voi2,x,b,f ew(1)=int_tabulated(x,b(3)-f) b=[b(0),center,b(1:3)] endif else begin b=dblarr(3) b(0)=height & b(1)=0.034237377 & b(2)=2311136.7 & continuum=1.00 r=curvefit(sx,sy,weights,b,sigma,function_name='voi',/noder,itmax=itmax,iter=it,yerror=yerror2) voi,x,b,f ew(1)=int_tabulated(x,1.-f) b=[b(0),center,b(1:2),continuum] endelse if (n_elements(it) lt 1) then begin ew(1)=0.0d0 print,'% GETEW_XY: Voigt fit did not succeed' return endif if (it gt itmax) then begin ew(1)=0.0d0 print,'% GETEW_XY: Voigt fit did not converge after ',itmax,' iterations' return endif if not keyword_set(noplot) then begin oplot,x,f,col=140,thick=3 endif yerror=[yerror1,yerror2] end