;+
; NAME:
;	Black_Body
; PURPOSE:
;	Compute black body emission spectrum (Planck function)
;	in units of Janskys/arcsec^2 at either an array of temperatures
;	or array of wavelengths.
; CALLING:
;	BB_spectrum = Black_Body( Temperature, dBdT, WAVELENGTH= )
; INPUTS:
;	Temperature = scalar or array, degrees Kelvin.
; KEYWORDS:
;	WAVELENGTH = wavelength in microns, array or scalar,
;		(stored into common block for reusage).
; OUTPUTS:
;	dBdT = optional output, derivative of Planck function respect to Temp.
;
;	Function returns array of Planck spectrum in units of Janskys/arcsec^2.
;	(Jansky = 1e-26 Watts/m^2/Hz)
; COMMON BLOCKS:
;	common Black_Body, wavel	(keeps wavelength array for reuse).
; EXTERNAL CALLS:
;	Assumes system variable !CV is defined
;	(cgs constants defined in Varosi's idl_startup.pro).
; HISTORY:
;	Written: Frank Varosi, NASA/GSFC, 1996.
;-

function Black_Body, Temperature, dBdT, WAVELENGTH=wavelength, GET_DERIV=gderiv

   common Black_Body, wavel

	if N_elements( wavelength ) GT 0 then wavel = wavelength
	Nwv = N_elements( wavel )
	Ntemp = N_elements( Temperature )

	if (Nwv LE 0) OR (Ntemp LE 0) then begin
		print,"syntax:  BB_spec = Black_Body( T, WAVEL=wavelens )"
		return,0
	   endif

	if (Nwv gt 1) and (Ntemp gt 1) and (Nwv NE Ntemp) then begin
		message,"warning:",/INFO
		print," Temperature & wavelength arrays are not same size"
	   endif

	arcsec2 = ( !DTOR/60/60 )^2
	hc = 2 * !cv.h * (1e-7 * 1e26) * !cv.c * 1e-2 * arcsec2
	hcw = hc * (1.e6/wavel)^3

	xT = ( !cv.h * !cv.c * 1e4 / !cv.k ) / ( Temperature * wavel )
	w = where( xT LT 88.7, nok )
	Nout = N_elements( xT )
	BT = fltarr( Nout )

	if nok gt 0 then begin
		eT = EXP( xT(w) )
		if (Nwv gt 1) then BT(w)=hcw(w)/(eT-1) else BT(w)=hcw/(eT-1)
	   endif

	if (N_params() GE 2) or keyword_set( gderiv ) then begin
		exprat = replicate( 1.0, Nout )
		if nok gt 0 then exprat(w) = eT/(eT-1)
		dBdT = BT * (xT/Temperature) * exprat
	   endif

return, BT
end