;+
; NAME:
;  planck_nu
; PURPOSE: (one line)
;  Return Planck function (brightness of a blackbody) in erg/cm^2/s/Hz
; DESCRIPTION:
;  Return Planck function (brightness of a blackbody) in erg/cm^2/s/Hz
; CATEGORY:
;  Astronomy
; CALLING SEQUENCE:
;  b = planck_nu(t, nu)
; INPUTS:
;  t - temperature in K 
;  nu - frequency (Hz)
; OPTIONAL INPUT PARAMETERS:
;  none
; KEYWORD INPUT PARAMETERS:
;  none
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;  Planck function
; COMMON BLOCKS:
;  uses share_physconst
; SIDE EFFECTS:
; RESTRICTIONS:
;  Must run phys_const first
;  For now, T and nu are scalars
; PROCEDURE:
; MODIFICATION HISTORY:
;  Written 2000 April, by Leslie Young, SwRI
;  Modified Jan 2001 LAY SwRI
;-

function planck_nu, t, nu

;common share_physconst
; c = speed_of_light & h = planck_constant & k = boltzman_constant
 physconstants
 c = !phys.c & h = !phys.h & k = !phys.k
 e = h*nu ; energy of a photon
 
 b = 2. * (e) * (nu/c)^2. / expm1( e / (k*t) )
 
 return, b
end

pro planck_nu_test
;common share_physconst
;c = speed_of_light
c = !phys.c
rs = 695265.e5 ; radius of sun in cm
au = 1.496e13 ; AU in cm
  t = 5770.
  lam_mu = findgen(1000)*0.01 + 0.1
  invcm = 1.e4/lam_mu
  nu = c * invcm
  b = fltarr(1000)
  for i=0,999 do b(i) = planck_nu(t, nu(i))
  plot, lam_mu, b * (nu/lam_mu) * !pi * (rs/au)^2. * 1.e-3, xr=[0,4]
end