function dustem_planck_function,temp,wave

;+
; NAME:
;	dustem_planck_function
; PURPOSE: 
;	To calculate the Planck function in units of MJy/sr 
;
; CALLING SEQUENCE: 
;	bbflux = dustem_planck_function(temp,wave) 
;
; INPUT PARAMETERS: 
;	WAVE   in microns
;	TEMP   Scalar giving the temperature of the planck function in degree K
;
; OUTPUT PARAMETERS:
;	BBFLUX - Scalar or vector giving the planck function at the specified
;		wavelength points.
;
; EXAMPLES:
;
; RESTRICTIONS:
;	Values less than approximately 1E-24 are truncated to 0.
;
; PROCEDURE:
;
; MODIFICATION HISTORY:
;-

 On_error,2

 if ( N_elements(wave) LT 1 ) then begin
     print,'Syntax - bbflux = dustem_planck_function(temp,wave)'
     print,'wave in microns, Temp in K, bbflux in MJy/sr'
     return,0
  endif    

;  if ( N_elements( temp ) NE 1 ) then $
;      read,'Enter a blackbody temperture',temp

  Nwav=n_elements(wave)
  Ntemp=n_elements(temp)
  Nel=Nwav*Ntemp
  bbflux = fltarr(Nel)
;  bbflux = wave*0.

; Gives the blackbody flux (i.e. PI*Intensity) ergs/cm2/s/a
  w = wave / 1.e4   ;microns to cm
;  w = wave / 1.E8                              ; Angstroms to cm    
  c1 =  3.74185E-5             ;constants appropriate to cgs units.
  C2 =  1.43883
  val =  c2/w/temp               
  good = where( val LT 88, Ngood )    ;Avoid floating underflow
  if ( Ngood GT 0 ) then  $
      bbflux[ good ] =  C1 / ( w[good]^5 * ( exp( val[good])-1. ) )

; fact is 1.e4/1.e7/2.99792e10*1e20/!pi
  fact=1.06177e+06
  bbflux=bbflux*w*w*fact ;MJy/sr

;  return, bbflux*1.E-8              ; Convert to ergs/cm2/s/A
  RETURN,bbflux

  end