function mk_lorentz,x,mean,core,peak,pder,betap=betap,normflx=normflx,$missing=missing,fwhm=fwhm, _extra=e ;+ ;function mk_lorentz ; returns a modified Lorentzian, Lb(X) ; ;syntax ; l=mk_lorentz(x,mean,core,peak,pder,/betap,/normflx,missing=m,/fwhm) ; ;parameters ; x [INPUT array; required] where G(X) must be computed ; mean [INPUT; default: mid_point(X)] position of peak ; core [INPUT; default: 0.1*range(X)] core width ; peak [INPUT; default: 1] value of Lb(X=MEAN) ; pder [OUTPUT; optional] partial derivatives of model wrt ; parameters at each X; calculated only if 5 parameters ; are supplied in call. ; * array of size [N(X),4], with columns containing the partial ; derivatives wrt MEAN, CORE, PEAK, and BETAP respectively ; ;keywords ; betap [INPUT; default=1] the index of the beta-profile. default ; is regular Lorentzian. ; * if NORMFLX is set and BETAP.le.0.5, the default is used ; normflx [INPUT] if set, {\int_{-\infty}^{\infty} dX Lb(X) = PEAK} ; missing [INPUT] 3 element array to populate missing values of ; MEAN, CORE, and PEAK ; fwhm [INPUT] if set, assumes that the input CORE is actually ; given as full-width at half-max, and converts to the true ; core width prior to calculation, and then converts back to ; full-width at half-max post-calc. ; _extra [JUNK] here only to prevent crashing ; ;description ; The Lorentzian is ; L(X) = PEAK / ( 1 + ((X-MEAN)/CORE)^2 ) ; The modified Lorentzian is ; Lb(X) = PEAK / ( 1 + ((X-MEAN)/CORE)^2 ) ^ BETAP ; ; When integrated over the real axis (Gradshteyn & Ryzhik, 3.251,2), ; \int dX Lb = PEAK*CORE*B(1/2, BETAP-1/2), ; where B(x,y) is the Beta-function, with BETAP > 1/2 ; ; hence, if NORMFLX is set, ; Lb(X) = (PEAK/CORE/B(1/2,beta-1/2))/(1+((X-MEAN)/CORE)^2)^BETAP ; ; The hwhm is defined as that value of X=W such that ; Lb(W)=0.5*Lb(X=MEAN) ; i.e., Lb(W)=(stuff)/(1+((W-MEAN)/CORE)^2)^(BETAP)=0.5*(stuff) ; hence, hwhm = W-MEAN = CORE*sqrt(2^(1/BETAP)-1) ; or, fwhm = CORE*2*sqrt(2^(1/BETAP)-1) ; ;usage summary ; * call as a function ; * generates modified Lorentzian model only at specified points X ; * needs MEAN, CORE, PEAK for complete specification ; ;examples ; x=findgen(1000)*0.01 & peasecolr & b=2.5 ; L1=mk_lorentz(x,5,1,1) & L2=mk_lorentz(x,5,1,1,/fwhm) ; Lb1=mk_lorentz(x,5,1,1,betap=b) & Lb2=mk_lorentz(x,5,1,1,betap=b,/fwhm) ; Lf1=mk_lorentz(x,5,1,1,/normflx) & Lf2=mk_lorentz(x,5,1,1,betap=b,/norm) ; plot,x,L1 & oplot,x,L2,col=2 ; oplot,x,Lb1,thick=2 & oplot,x,Lb2,thick=2,col=2 ; oplot,x,Lf1,thick=2,line=1 & oplot,x,Lf2,thick=2,col=3,line=1 ; ;subroutines ; NONE ; ;history ; vinay kashyap (Oct98) ; changed keyword BETA to BETAP; corrected normalization (VK; Dec98) ; recomputed partial derivatives (VK; MarMM) ; now works even if X are integers (VK; Jul01) ; what if PEAK=0? (VK; Apr02) ; converted array notation to IDL 5 (VK; Apr02) ; added partial derivative wrt BETAP to PDER output array (VK; Jun02) ; changed keyword NORM to NORMFLX (VK; Oct02) ; added keyword FWHM (VK; Apr03) ;- np=n_params() if np lt 1 then begin print, 'Usage: Lb=mk_lorentz(x,mean,core,peak,pder,betap=b,missing=m,/normflx,/fwhm)' print, ' generates a modified Lorentzian Lb(x;beta)' return,[-1L] endif ;initialize nx=n_elements(x) & x0=x[nx/2] & mxx=max(x,min=mnx) ; figure out the defaults if not keyword_set(betap) then b=1. else b=betap[0] if keyword_set(normflx) and b lt 0.5 then begin message,'Normalization becomes infinite! Resetting BETAP to 0.5',/info b=0.5 endif if not keyword_set(missing) then missing=[x0,0.1*(mxx-mnx),1.] if np lt 4 then p=missing[2] else p=peak[0] if np lt 3 then c=missing[1] else c=core[0] if c lt 0 then c=abs(c) & if c eq 0 then c=missing[1] if np lt 2 then m=missing[0] else m=mean[0] ; fwhm? if keyword_set(fwhm) then c2f=2.*sqrt(2.^(1./b)-1.) else c2f=1. c=c/c2f ;if input is FWHM, convert to CORE ; renorm if keyword_set(normflx) then begin if b gt 0.5 then bnorm=1./beta(0.5,b-0.5)/c else bnorm=0. endif else bnorm=1. p=p*bnorm ; compute function z=(x-m+0.0)/c & z=alog10(1.+z^2) & z=-b*z Lb=make_array(size=size(0.*x)) oz=where(z gt -30,moz) & if moz gt 0 then Lb[oz]=10.^(z[oz]) ; compute partial derivatives if np ge 5 then begin pder = fltarr(nx,4) z=(x-m+0.0)/c & zz=(1.+z^2) ; partial wrt MEAN Lbm=p*(zz^(-b-1))*(2.*b/c^2)*(x-m+0.0) pder[*,0] = Lbm[*] ; partial wrt CORE Lbc=p*(zz^(-b-1))*2.*b*(x-m+0.0)^2/c^3 if keyword_set(normflx) then Lbc=Lbc-p*Lb/c pder[*,1] = Lbc[*] ; partial wrt PEAK Lbp=Lb*bnorm pder[*,2]=Lbp[*] ; partial wrt BETAP if not keyword_set(normflx) then begin Lbb=Lb*p & oo=where(Lbb gt 0,moo) if moo gt 0 then pder[oo,3]=-Lbb[oo]*alog(Lbb[oo]) endif else begin ; because this is easier than differentiating the beta function delb=0.01 & b2=b+delb & b1=(b-delb)>0.5 tmp1=mk_lorentz(x,m,c,p,betap=b1,/normflx) tmp2=mk_lorentz(x,m,c,p,betap=b2,/normflx) pder[*,3]=(tmp2-tmp1)/delb/2. endelse endif Lb=Lb*p ; fwhm? c=c*c2f ;if input was FWHM, convert back from CORE ;if np ge 5 then begin ; pder = fltarr(nx,3) ; z=(x-m)/c ; ;z=(1.+z^2) & z=-p*(b+1)/z^(b+1) & oz=where(abs(z) gt 1e-10,moz) ; z2=(1.+z^2) & z3=p*(-b)/z2^(b+1) & oz=where(abs(z) gt 1e-10,moz) ; ; partial wrt MEAN ; Lbm=0*Lb & if moz gt 0 then Lbm(oz)=z3(oz)*(1./c^2)*(2.*(x-m)*(-1)) ; pder(*,0) = Lbm(*) ; ; partial wrt CORE ; Lbc=0*Lb & if moz gt 0 then Lbc(oz)=z(oz)*(x-m)^2*(-2./c^3) ; if keyword_set(normflx) then$ ; if moz gt 0 then Lbc(oz)=Lbc(oz)-Lb(oz)/c ; pder(*,1) = Lbc(*) ; ; partial wrt PEAK ; Lbp=Lb*bnorm/p ; pder(*,2) = Lbp(*) ;endif return,Lb end