function mk_bbkev,x,norm,temp,pder,verbose=verbose, _extra=e ;+ ;function mk_bbkev ; returns the Planck function B(E) [ph/s/cm^2/keV] as a ; function of E [keV] at a given temperature T. ; ; In general, ; B(nu)d(nu) = (R/d)^2*(2*pi/c^2) * d(nu)*nu^2/(exp(h*nu/kT)-1) ; for E=h*nu, ; B(E)dE = (R/d)^2 * (2*pi/c^2/h^3) * dE*E^2/(exp(E/kT)-1) ; == A * (E^2/(exp(E/kT)-1)) * dE ; ;syntax ; b=mk_bbkev(x,norm,temp,pder,verbose=verbose) ; ;parameters ; x [INPUT array; required] where B(X) must be computed ; norm [INPUT; required] normalization for B(X) ; temp [INPUT; required] temperature in [K] ; pder [OUTPUT; optional] partial derivatives of model wrt parameters ; at each X; calculated only if 4 parameters are supplied in call. ; ;keywords ; verbose [INPUT] controls chatter ; _extra [JUNK] here only to prevent crashing ; ;usage summary ; * call as a function ; ;subroutines ; NONE ; ;history ; vinay kashyap (Aug08; based loosely on 1993 blackbody.pro) ;- ; usage ok='ok' & np=n_params() nx=n_elements(x) & nn=n_elements(norm) & nt=n_elements(temp) if np lt 3 then ok='Insufficient parameters' else \$ if nx eq 0 then ok='X is not defined' else \$ if nn eq 0 then ok='NORM is not defined' else \$ if nt eq 0 then ok='TEMPERATURE is not defined' else \$ if nn gt 1 then ok='NORM should not be an array' else \$ if nn ne nt then ok='TEMPERATURE incompatible with NORM' if ok ne 'ok' then begin print,'Usage: B=mk_bbkev(x,norm,temperature,pder,verbose=verbose)' print,' compute and return a Planck spectrum B(E) [ph/s/cm^2/keV]' if np ne 0 then message,ok,/info return,-1L endif ; inputs vv=0 & if keyword_set(verbose) then vv=long(verbose[0])>1 ; aa=norm[0] tt=temp[0] if vv gt 20 then begin ct=strtrim(tt,2) cc='B(keV;T='+ct+')' message,'Generating a Planck curve '+cc,/informational endif ; compute function hh=6.6261760d-27 ;[erg s] cc=2.9979246d+10 ;[cm/s] kB=1.3806620e-16 ;[erg/K] ee=1.6021892e-19 ;[C] onekeV=ee*1e7*1e3 ;[erg/keV] Eerg=x*onekeV ;[keV]*[erg/keV]=[erg] xx=(Eerg/kB/tt) < 69. tmp=(exp(xx)-1.) > 1d-30 ;B(E)dE = (R/d)^2*(2*pi/c^2/h^3) * E^2/(exp(E/kT)-1) * dE Bx = exp( alog(aa)+2.*alog(x)-alog(tmp) ) if vv gt 50 then plot,x,Bx,xtitle='E [keV]',ytitle='B(E)' ; compute partial derivatives if np ge 3 then begin pder=fltarr(nx,2) ; partial wrt norm pder[*,0]=Bx/aa ; partial wrt T tmp2=exp(alog(Bx)+alog(xx)-alog(tmp)+xx) pder[*,1]=tmp2 endif return,Bx end