dustem_planck_function.pro 1.82 KB
function dustem_planck_function,temp,wave

;+
; NAME:
;	dustem_planck_function
;
; PURPOSE: 
;	To calculate the Planck function in units of MJy/sr 
;
; CATEGORY:
;    DustEMWrap, Distributed, Low-Level, User convenience, Analysis
;
; 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:
;	BBFLUX - Scalar or vector giving the Planck function at the specified
;		wavelength points.
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;    None
;
; EXAMPLES:
;
; RESTRICTIONS:
;	Values less than approximately 1E-24 are truncated to 0.
;;
; MODIFICATION HISTORY:
;    Written by JPB Apr-2011
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
;-

 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