pro planckcal, B, T, spectrum, interval=interval, system=system, terse=terse
;+
; planck, B, T, spectrum, interval=interval, system=system, terse=terse
;     
;    calculates planck function in standard
;    physical units.
; 
; B = planck function (output float array)
; T = temperature (float)
; spectrum = wavelength/frequency/wavenumber values (float array)
; interval = spectrum units (either 'wavelength', 'frequency', 'wavenumber')
; system   = units (either 'cgs' or 'mks')
; terse  0 = announce units of B (default), 1 = do not announce units of B
;
;-

lowerlimit = 1e-30
n = n_elements(spectrum)
if not keyword_set(system)   then system   = 'cgs'
if not keyword_set(interval) then interval = 'wavelength'
if not keyword_set(terse)    then terse	   = 0


case system of
'mks':	    begin
	    h  = 6.626e-34        ; J s
	    c  = 2.998e+08        ; m s-1
	    kb = 1.381e-23        ; J K-1
	    bunits =               'J s-1 m-2 sr-1 '
	    end
'cgs':	    begin
	    h  = 6.626e-27        ; erg s
	    c  = 2.998e+10        ; cm s-1
	    kb = 1.381e-16        ; erg K-1
	    bunits =		   'erg s-1 cm-2 sr-1 '
	    end
endcase


case interval of
'wavelength': begin
	      B = double((2* h* c^2 / spectrum^5) / (exp(h* c/ (kb* T* spectrum)) - 1))
	      if system eq 'mks' then bunits = bunits + 'm-1' $
	      else bunits = bunits + 'cm-1'
	      end
'frequency':  begin
	      B = double((2* h* spectrum^3 / c^2) / (exp(h* spectrum / (kb* T)) - 1))
	      bunits = bunits + 'Hz-1'
	      end
'wavenumber': begin
	      B = double((2* h* c^2 * spectrum^3) / (exp(h* c* spectrum / (kb* T)) - 1))
	      if system eq 'mks' then bunits = bunits + 'm' $
	      else bunits = bunits + 'cm'
	      end
endcase

if not terse then print, 'Planck function returned in units of '+bunits+'. '

return
end