dustem_create_rfield.pro
5.4 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
FUNCTION DUSTEM_CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=chi, fname=fname, help=help
if N_PARAMS() LT 1 or keyword_set(help) then begin
print,'------------------------------------------------------------------------------------------------------'
print,'FUNCTION DUSTEM_CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=chi, fname=fname'
print,'------------------------------------------------------------------------------------------------------'
print,''
print,' generates the radiation field for DUSTEM (ISRF.DAT)'
print,' in erg/cm2/s/Hz (4*!pi*I_nu)'
print,' RFIELD = ISRF + WDIL*PI*BB(TEMP) '
print,' ISRF from Mathis et al. 1983, A&A 128, 212'
print,' NB: if blackbody is isotropic set WDIL to 4 (to get 4*!pi factor)'
print,''
print,' TEMP (I): blackbody temperature (can be an array). If 0 only ISRF.'
print,' X (I): X-grid for the ISRF in microns. Default is 200 pts over 0.01-10^5 microns.'
print,' (including WCUT point). You want to include WCUT for accurate edges.'
print,' ISRF (I): if set to 0 no ISRF is added, 1 is Mathis (default), 2 is Habing'
print,' WDIL (I): blackbody dilution factor (can be an array)'
print,' FNAME(I): filename for ISRF.DAT'
print,' WCUT (I): for wave < wcut radiation field is 0 '
print,' G0 (O): factor flux(6-13.6eV) wrt. Mathis field '
print,' CHI (O): scaling factor at 100nm wrt. Mathis field'
print,''
print,'Example: tt = dustem_create_rfield([2d4,5d4],wdil=[1.d-14,1.d-16],x=x,fname=''ISRF.DAT'')'
print,''
print,' Created Aug. 2009, L. Verstraete, IAS'
print,' Force the WCUT point, May 2011, LV'
print,''
print,'------------------------------------------------------------------------------------------------------'
RETURN,0
endif
; inits
unit='W'
ly_limit = 9.11267101235d-2
xi = [ ly_limit, 2.07d-1 ] ; for G0 6-13.6 eV in microns
x1 = 1.d-1 ; for CHI
if n_elements(WCUT) EQ 0 then wcut=ly_limit
if n_elements(x) EQ 0 then begin
nx = 199
xb = [ 1.d-2, 1.d5 ] ; wave boundaries in microns
dx = (alog(xb(1))-alog(xb(0))) / (nx-1)
x = alog(xb(0)) + dindgen(nx)*dx
x = exp(x)
x = [x, wcut] ; add wcut to avoid coarse edge
x = x(UNIQ(x,SORT(x)))
nx = n_elements(x)
if x(nx-1) NE xb(1) then x(nx-1)=xb(1) ; check rounding
endif else begin
print,'(W) DUSTEM_CREATE-RFIELD : using input wave x'
x = double(x)
nx = n_elements(x)
ix = WHERE( ABS(x-wcut)/x LE 0.01, cx )
if cx EQ 0 then begin
print,'(W) DUSTEM_CREATE_RFIELD: your x-grid does not contain wcut.'
print,' Should be included if radiation is 0 below wcut.'
endif
endelse
rfield = dblarr(nx)
; get Habing for normalization
rhabing = DUSTEM_HABING_FIELD(x,unit=unit)
rhabing = 1.d3*x*rhabing/3.d14 ; erg/cm2/s/Hz
; get isrf
rmathis = DUSTEM_MATHIS_FIELD(x,unit=unit)
rmathis = 1.d3*x*rmathis/3.d14 ; erg/cm2/s/Hz
if n_elements(ISRF) EQ 0 then begin
print,'(W) DUSTEM_CREATE_RFIELD : ISRF keyword not set, using default Mathis'
ISRF=1
endif
if ISRF EQ 0 then begin
print,'(W) DUSTEM_CREATE_RFIELD : not adding diffuse ISM ISRF component'
endif else if ISRF EQ 1 then begin
print,'(W) DUSTEM_CREATE_RFIELD : adding Mathis ISRF'
rfield = rmathis
endif else if ISRF EQ 2 then begin
print,'(W) DUSTEM_CREATE_RFIELD : adding Habing ISRF'
rfield = rhabing
endif
; get blackbody
ntemp = n_elements(TEMP)
if ntemp GT 0 then begin
bb = DUSTEM_BLACKBODY( x, temp, unit=unit, wdil=wdil )
bb = bb * 1.d3 * !pi * x^2 / 3.d14 ; erg/cm2/s/Hz
print,'(W) DUSTEM_CREATE_RFIELD : adding BB with T= ',temp, format='(A38,10(1E10.4,1x))'
print,' dilution factor wdil= ',wdil, format='(A38,10(1E10.4,1x))'
endif else bb=0.D0
rfield = rfield + bb
; apply cut
ix = WHERE( x LT WCUT, cx )
if cx GT 0 then rfield(ix)=0.d0
; get G0
ig = where( x GE xi(0) AND x LE xi(1), cg)
if cg GT 0 then begin
xx=x(ig)
yy=rfield(ig)
rr=rmathis(ig)
s1 = TOTAL( (yy(1:cg-1)+yy(0:cg-2))*(xx(1:cg-1)-xx(0:cg-2))/2.D0 )
s2 = TOTAL( (rr(1:cg-1)+rr(0:cg-2))*(xx(1:cg-1)-xx(0:cg-2))/2.D0 )
g0 = s1/s2
endif else g0 = 0.d0
print,'(W) DUSTEM_CREATE_RFIELD : G0 =',g0
; 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 = rfield(ig) / rmathis(ig)
endif else chi = 0
print,'(W) DUSTEM_CREATE_RFIELD : chi =',chi
; write file
if n_elements(FNAME) NE 0 then begin
OPENW, iu, fname, /get_lun
printf, iu, '# DUSTEM: exciting radiation field featuring '
if ISRF EQ 1 then printf, iu, '# Mathis ISRF'
if NTEMP GT 0 then begin
a = '# Blackbody with T='
a1 = '# dilution factor wdil='
for i = 0, ntemp - 1 do begin
a = a + string( format='(2x,1E11.4)', temp(i) )
a1 = a1 + string( format='(2x,1E11.4)', wdil(i) )
endfor
printf, iu, a
printf, iu, a1
endif
printf, iu, '# Nbr of points'
printf, iu, '# wave (microns), 4*pi*Inu (erg/cm2/s/Hz)'
printf, iu, nx, format='(i4)'
for i=0,nx-1 do begin
printf, iu, x(i), rfield(i), format='(2(1E13.6,2x))'
endfor
FREE_LUN, iu
print,'(W) DUSTEM_CREATE_RFIELD: radiation field written in ', strtrim(fname,2)
ENDIF
RETURN, RFIELD
END