function mk_bbang,x,norm,temp,pder,verbose=verbose, _extra=e ;+ ;function mk_bbang ; returns the Planck function B(lam) [ph/s/cm^2/Ang] as a ; function of lam [Ang] 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 lam=c/nu, ; B(lam)d(lam) = (R/d)^2 * (2*pi*c) * d(lam) * (1/lam^4)/(exp(hc/lam/k/T)-1) ; == A * (lam^-4/(exp(hc/lam/k/T)-1)) * d(lam) ; ;syntax ; b=mk_bbang(x,norm,temperature,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_bbang(x,norm,temperature,pder,verbose=verbose)' print,' compute and return a Planck spectrum B(lambda) [ph/s/cm^2/Ang]' 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(Ang;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] lamerg=hh*cc/(x*1d-8) ;[erg] xx=(lamerg/kB/tt) < 69. tmp=(exp(xx)-1.) > 1d-30 ;B(l)dl = (R/d)^2*(2*pi*c) * (1/l^4)/(exp(hc/lkT)-1) * dl Bx = exp( alog(aa)-4.*alog(x)-alog(tmp) ) if vv gt 50 then plot,x,Bx,xtitle='lambda [Ang]',ytitle='B(lambda)' ; compute partial derivatives if np gt 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