dustem_planck_function.pro
2.01 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
function dustem_planck_function,temp,wave,help=help
;+
; 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.
; ACCEPTED KEYWORDS
; help = If set, print this help
;
; 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 keyword_set(help) THEN BEGIN
doc_library,'dustem_planck_function'
bbflux=-1.
goto,the_end
ENDIF
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
the_end:
; return, bbflux*1.E-8 ; Convert to ergs/cm2/s/A
RETURN,bbflux
end