dustem_blackbody.pro
4.79 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
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
FUNCTION DUSTEM_BLACKBODY, x, temp, unit=unit, wdil=wdil, check=check, g0=g0, chi=chi, $
bb_to_hab=bb_to_hab, NORM=norm
if n_params() EQ 0 then begin
print,'------------------------------------------------------------------------------------------'
print,'FUNCTION DUSTEM_BLACKBODY, x, temp, unit=unit, wdil=wdil, NORM=norm, check=check, g0=g0, chi=chi'
print,'------------------------------------------------------------------------------------------'
print,' computes a (or a sum of) blackbody brightness(es) B(x,T) in W/m2/(x unit)/sr '
print,' X (I): blackbody abscissa in cm-1, GHz, microns or eV'
print,' TEMP (I): blackbody temperature (can be an array)'
print,' WDIL (I): dilution factor for each TEMP. [Default=1].'
print,' UNIT (I): blackbody x unit ''K'' is B_nu (W/m2/cm-1/sr), x in cm-1 '
print,' ''F'' is B_nu (W/m2/Hz/sr), x in GHz [default]'
print,' ''W'' is B_lambda (W/m2/um/sr), x in microns.'
print,' ''E'' is B_E (W/m2/eV/sr), x in eV'
print,' NORM (I): overall normalizing factor, returns NORM*BB/MAX(BB)'
print,' CHECK (O): check energy conservation for the BB (power is sigma*T^4)'
print,' CHI (O): scaling factor @ 1000 A in ISRF unit (1.6e-3 erg/cm2/s)'
print,' G0 (O): scaling factor for integrated 5-13.6 eV flux in ISRF unit.'
print,''
print,' Created 2001, L Verstraete IAS'
print,' More units and checks, Oct 2010 LV'
print,'---------------------------------------------------------------------------------------'
return,0
endif
nx = n_elements( X )
nt = n_elements( TEMP )
nw = n_elements( WDIL )
if nw EQ 0 then begin
wdil = dblarr(nt) + 1.
endif else begin
if nw LT nt then begin
print,'DUSTEM_BLACKBODY (F): WDIL and TEMP have different sizes'
return, 0
endif
endelse
if n_elements( UNIT ) EQ 0 then unit='F' $
else begin
unit = strupcase(strcompress(unit,/rem))
if unit NE 'K' AND unit NE 'F' AND unit NE 'W' AND unit NE 'E' then begin
print,'DUSTEM_BLACKBODY (W): not a valid BB unit --> make a B_nu'
unit = 'F'
endif
endelse
if n_elements( CHECK ) GT 1 then check = 0
x = double(x)
temp = double(temp)
wdil = double(wdil)
xpi = 3.1415926535897932385d
clight = 2.9979246d08 ; m/s
kb = 1.3806503d-23
hp = 6.62606876d-34
hck = hp*clight / kb
ze = 1.60217653d-19 ; electron charge
hev = hp/ze ; h in eV unit
xly = 12.4
isrf0 = 1.6d-6 ; 4*pi*nu*I_nu (W/m2) in solar neighborhood
;
; loop on temperatures
;
bbt = dblarr(nx)
acheck = 0.
FOR j = 0, nt-1 do begin
xi = [6., 13.6]
CASE UNIT OF
'K': begin
be = exp( -(1d2*x) * hck / temp(j) )
bb = be * (1d2*x)^3 / (1.d0 - be)
bb = bb * 2d0 * hp * clight^2
bb = bb * 1d2 ; m-1 to cm-1
xi = xi / hev/clight/1d2
x1 = xly / hev/clight/1d2
end
'F': begin
be = exp( -(1d9*x) * hp/kb/ temp(j) )
bb = be * (1d9*x)^3 / (1.d0 - be)
bb = bb * 2d0 * hp / clight^2
xi = xi / hev / 1d9
x1 = xly / hev / 1d9
end
'W': begin
be = exp( -hp*clight/kb/ (1d-6*x) / temp(j) )
bb = be / (1d-6*x)^5 / (1.d0 - be)
bb = bb * 2d0 * hp * clight^2
bb = 1d-6 * bb ; m-1 to micron-1
xi = hev*clight*1d6 / xi
x1 = hev*clight*1d6 / xly
end
'E': begin
be = exp( -(ze*x)/kb / temp(j) )
bb = be * (ze*x)^3 / (1.d0 - be)
bb = bb * 2d0 / hp^3 / clight^2
bb = bb * ze ; J-1 to eV-1
xi = xi
x1 = xly
end
ENDCASE
ck = 0.
for jx = 1L,nx-1 do $
ck = ck + ( bb(jx)+bb(jx-1) )*( x(jx)-x(jx-1) ) / 2.d0
ck = xpi* ck / 5.6697d-8 / temp(j)^4
if n_elements( CHECK ) NE 0 then begin
if check EQ 1 then $
print,'DUSTEM_BLACKBODY (W): integrated BB('+strtrim(temp(j),2)+$
' K) is '+ strtrim(ck,2)+ ' times Sigma*T^4'
endif
acheck = [ acheck, ck ]
bbt = bbt + bb*wdil(j)
ENDFOR
check = acheck(1:*)
; get G0
ig = where( x GE MIN(xi) AND x LE MAX(xi), cg)
if cg GT 0 then begin
xx=x(ig)
yy=xpi*x(ig)*bbt(ig)
s1 = TOTAL( (yy(1:cg-1)+yy(0:cg-2))*(xx(1:cg-1)-xx(0:cg-2))/2.D0 )
g0 = s1/isrf0
endif else g0 = 0.d0
; get chi
ig = where( abs(2.*(x-x1)/(x+x1)) LE 1.d-1, cg)
if cg GT 0 then begin
tmp = min( abs(x(ig)-x1), i_min )
ig = ig(i_min)
ig = ig(0)
chi = xpi * x(ig) * bbt(ig) / isrf0
endif else chi = 0
if n_elements( NORM ) NE 0 then bbt = norm * bbt/max(bbt)
the_end:
RETURN, BBT
END ;FUNCTION BLACKBODY