function mk_bknpower,x,norm,gamma1,gamma2,pbreak,pder,Xo=Xo,$nobreak=nobreak,verbose=verbose, _extra=e ;+ ;function mk_bknpower ; returns a broken Power-law function, ; p(X) = NORM*(X/Xo)^gamma1 {X.le.PBREAK} ; = NORM*(PBREAK/Xo)^(gamma1-gamma2)*(X/Xo)^gamma2 {X.ge.PBREAK} ; ;syntax ; p=mk_bknpower(x,norm,gamma1,gamma2,pbreak,pder,Xo=Xo,/nobreak,verbose=verbose) ; ;parameters ; X [INPUT; required] where p(X) must be computed ; norm [INPUT; required] normalization ; gamma1 [INPUT; required] power-law index 1 (valid for 0PBREAK) ; pbreak [INPUT; necessary] break point ; * GAMMA2 and PBREAK are ignored if NOBREAK is set ; pder [OUTPUT; optional] partial derivatives of model wrt parameters ; at each X; calculated only if 6 parameters are supplied in call. ; * array of size [N(X),4], with columns containing the partial ; derivatives wrt NORM, GAMMA1, GAMMA2, and PBREAK respectively ; * to have it return PDER as [N(X),2] for a non-broken power-law ; distribution, set NOBREAK (must have GAMMA2 and PBREAK in the ; call, but they will be ignored) ; ;keywords ; Xo [INPUT] X-value at which which normalization is defined ; * default is 1.0 ; nobreak [INPUT] if set, ignores GAMMA2 and PBREAK, essentially assuming ; that PBREAK=$\infty$; verbose [INPUT] controls chatter ; _extra [JUNK] here only to prevent crashing the program ; ;example ; x=(findgen(100)+1)/10. & p=mk_bknpower(x,1.,-1,-3,1.5) & plot,x,p,/xl,/yl ; p2=mk_bknpower(x,1.,-1,-3,1.5,/nobreak) & oplot,x,p2,line=2,thick=2 ; ;history ; Vinay Kashyap (Apr2008) ;- ok='ok' & np=n_params() & nx=n_elements(x) nn=n_elements(norm) & ng1=n_elements(gamma1) yesb=1 & if keyword_set(nobreak) then yesb=0 ng2=n_elements(gamma2) & nb=n_elements(pbreak) if np eq 0 then ok='Insufficient parameters' else$ if np lt 3 then ok='Insufficient parameters' else $if keyword_set(yesb) and np lt 5 then ok='Insufficient parameters' else$ if nx eq 0 then ok='X is undefined' else $if nn eq 0 then ok='NORM is undefined' else$ if ng1 eq 0 then ok='GAMMA1 is undefined' else $if keyword_set(yesb) and ng2 eq 0 then ok='GAMMA2 is undefined' else$ if keyword_set(yesb) and nb eq 0 then ok='PBREAK is undefined' else $if nn gt 1 then ok='NORM must be scalar' else$ if ng1 gt 1 then ok='GAMMA1 must be scalar' else $if keyword_set(yesb) and ng2 gt 1 then ok='GAMMA2 must be scalar' else$ if keyword_set(yesb) and nb gt 1 then ok='PBREAK must be scalar' if ok ne 'ok' then begin print, 'Usage: p=mk_bknpower(x,norm,gamma1,gamma2,pbreak,pder,Xo=Xo) print, ' generates a broken power-law p(X)' if np ne 0 then message,ok,/informational return,-1L endif ; check inputs xx=x y0=norm[0] & g1=gamma1[0] g2=1.0 & xbrk=max([xx])+1. if keyword_set(yesb) then begin g2=gamma2[0] & xbrk=pbreak[0] endif ; keywords vv=0L & if keyword_set(verbose) then vv=long(verbose[0])>1 x0=1.0 & if keyword_set(Xo) then x0=Xo[0]*1.0 pp=fltarr(nx) o1=where([xx] le xbrk,mo1,complement=o2,ncomplement=mo2) ; do not multiply by NORM yet, in case it happens to be 0 if mo1 gt 0 then pp[o1]=(xx[o1]/x0)^(g1) if mo2 gt 0 then pp[o2]=((xbrk/x0)^(g1-g2))*((xx[o2]/x0)^g2) ; catch any NaNs and infinities and set them to 0 oy=where(finite(pp) eq 0,moy) if moy gt 0 then pp[oy]=0. ;compute partial dervatives ; PDER Order: wrt norm, gamma1, gamma2, pbreak if np ge 6 then begin if keyword_set(yesb) then pder=fltarr(nx,4) else pder=fltarr(nx,2) ;NORM pder[*,0]=pp ;GAMMA1 if mo1 gt 0 then pder[o1,1]=y0*pp[o1]*alog(xx[o1]/x0) if mo2 gt 0 then pder[o2,1]=y0*pp[o2]*alog(xbrk/x0) if keyword_set(yesb) then begin ;GAMMA2 if mo1 gt 0 then pder[o1,2]=0. if mo2 gt 0 then pder[o2,2]=-y0*pp[o2]*alog(xbrk/x0)+y0*pp[o2]*alog(xx[o2]/x0) ;PBREAK if mo1 gt 0 then begin dbrk=0.01*xbrk & if dbrk eq 0 then dbrk=0.001 pp1=mk_bknpower(xx[o1],y0,g1,g2,xbrk+dbrk) pp2=mk_bknpower(xx[o1],y0,g1,g2,xbrk-dbrk) pder[o1,3]=(pp1-pp2)/(dbrk*2.) endif if mo2 gt 0 then pder[o2,3]=y0*pp[o2]*(g1-g2)*(x0/xbrk) endif endif return,y0*pp ;multiply by NORM and exit end