function mk_slant,x,mean,core,peak,pder,angle=angle,betap=betap,$normflx=normflx,missing=missing, _extra=e ;+ ;function mk_slant ; returns the product of a slanted line with a modified ; Lorentzian, Lb(X) ; ; the idea is to model the asymmetrical line profiles seen in ; CXO grating data. ; ;syntax ; Lbs=mk_slant(x,mean,core,peak,pder,angle=angle,betap=betap,$ ; /normflx,missing=missing) ; ;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),5], with columns containing the partial ; derivatives wrt MEAN, CORE, PEAK, BETAP, and ANGLE respectively ; ;keywords ; angle [INPUT; default=0] the atan(slope) of the straight line ; that makes the Lorentzian asymmetric ; 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 Lbs(X) = PEAK} ; missing [INPUT] 3 element array to populate missing values of ; MEAN, CORE, and PEAK ; _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 ; The asymmetric modified Lorentzian is ; Lbs(X) = (PEAK + PEAK*SLOPE*(X-MEAN)) / $; ( 1 + ((X-MEAN)/CORE)^2 ) ^ BETAP ; ; When integrated over the real axis (Gradshteyn & Ryzhik, 3.251,2), ; the second terms drops out and the integral is identical to that ; of the modified Lorentzian itself, i.e., ; \int dX Lbs = CORE*PEAK*B(1/2,BETAP-1/2) ; where B(x,y) is the Beta-function, with BETAP > 1/2 ; ; if NORMFLX is set, ; Lbs(X) = ((PEAK+PEAK*SLOPE*(X-MEAN))/CORE/B(1/2,beta-1/2)) /$ ; (1+((X-MEAN)/CORE)^2)^BETAP ; ;usage summary ; * call as a function ; * generates modified Lorentzian model only at specified points X ; * needs MEAN, CORE, PEAK for complete specification ; ;subroutines ; NONE ; ;history ; vinay kashyap (MarMM, based on MK_LORENTZ) ; now works even if X are integers (VK; Jul01) ; converted array notation to IDL 5 (VK; Apr02) ; accounted for case PEAK=0, added partial derivatives of BETAP and ; ANGLE to PDER (VK; Jun02) ; changed keyword NORM to NORMFLX (VK; Oct02) ;- np=n_params() if np lt 1 then begin print, 'Usage: Lbs=mk_slant(x,mean,core,peak,pder,angle=a,betap=b,missing=m,/normflx)' print, ' generates an asymmetric modified Lorentzian Lbs(x;beta,angle)' 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=float(betap[0]) if not keyword_set(angle) then ang=0. else ang=float(angle[0]) while ang gt 90 do ang=ang-90. & if ang eq 90. then ang=ang-1e-3 while ang lt -90 do ang=ang+90. & if ang eq -90. then ang=ang+1e-3 slope=tan(ang*!pi/180.) 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] ; 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 Lbs=make_array(size=size(0.0*x)) oz=where(z gt -30,moz) & if moz gt 0 then Lbs[oz]=10.^(z[oz]) Lbs=Lbs*(1.+slope*(x-m)) > 0. ; compute partial derivatives if np ge 5 then begin pder = fltarr(nx,5) z=(x-m+0.0)/c & zz=(1.+z^2) & pp=p*(1.+slope*(x-m+0.0)) ; partial wrt MEAN Lbsm=pp*(zz^(-b-1))*(2.*b/c^2)*(x-m+0.0)-(zz^(-b))*p*slope pder[*,0] = Lbsm[*] ; partial wrt CORE Lbsc=pp*(zz^(-b-1))*2.*b*(x-m+0.0)^2/c^3 if keyword_set(normflx) then Lbsc=Lbsc-Lbs/c pder[*,1] = Lbsc[*] ; partial wrt PEAK Lbsp=Lbs*bnorm pder[*,2] = Lbsp[*] ; partial wrt BETAP if not keyword_set(normflx) then begin Lbsb=Lbs*p & oo=where(Lbsb gt 0,moo) if moo gt 0 then pder[oo,3]=-Lbsb[oo]*alog(Lbsb[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_slant(x,m,c,p,angle=ang,betap=b1,/normflx) tmp2=mk_slant(x,m,c,p,angle=ang,betap=b2,/normflx) pder[*,3]=(tmp2-tmp1)/delb/2. endelse ; partial wrt ANGLE dela=0.01 & a1=ang-dela & a2=ang+dela tmp1=mk_slant(x,m,c,p,angle=a1,betap=b,normflx=normflx) tmp2=mk_slant(x,m,c,p,angle=a2,betap=b,normflx=normflx) pder[*,4]=(tmp2-tmp1)/dela/2. endif Lbs=Lbs*p return,Lbs end