dustem_plugin_dl07_isrf_model.pro
3.81 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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
FUNCTION dustem_plugin_dl07_isrf_model, key=key, val=val, scope=scope, paramtag=paramtag, help=help, Umean=Umean,st=st
;+
; NAME:
; dustem_plugin_dl07_isrf_model
;
; PURPOSE:
; Returns the dust spectrum after heating by a Draine-like ISRF
; with Umin, Umax, gamma, alpha parameters described by the val keyword
;
; CATEGORY:
; DUSTEM Wrapper, Plugin
;
; CALLING SEQUENCE:
; sed=dustem_plugin_dl07_isrf_model(key=,val=,Umean=,/paramtag][,/scope][,/help])
;
; INPUTS:
; None
;
; OPTIONAL INPUT PARAMETERS:
; key = input parameter number.
; 1: gamma
; 2: alpha
; 3: Umin
; 4: Umax
; val = input parameter values for the above.
;
; OUTPUTS:
; None
;
; OPTIONAL OUTPUT PARAMETERS:
; Umean = the average radiation field
;
; ACCEPTED KEY-WORDS:
; help = if set, print this help
;
; COMMON BLOCKS:
; None
;
; SIDE EFFECTS:
; None
;
; RESTRICTIONS:
; The DustEM fortran code must be installed
; The DustEMWrap IDL code must be installed
; The path of the ISRF must be assigned to one of the DustEMWrap free pointers.
; The user-defined ISRF needs to be on the same dustem ISRF grid.
;
; PROCEDURE:
;
; EXAMPLES
; res=dustem_dl07_isrf_model()
;
; MODIFICATION HISTORY:
; Evolution details on the DustEMWrap gitlab.
; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
;-
IF keyword_set(help) THEN BEGIN
doc_library,'dustem_dl07_isrf_model'
goto,the_end
ENDIF
IF keyword_set(scope) THEN BEGIN
out=0
goto, the_end
ENDIF
IF keyword_set(paramtag) THEN BEGIN
out=0
goto, the_end
ENDIF
;dustem_print_params
;stop
spec=0.
;==== default parameter values
Umin=1.
Umax=1.e6
alpha=2.
gamma=0. ;default is no contribution from U<>1.
IF keyword_set(key) THEN BEGIN
ind=where(key EQ 1,count) ;gamma parameter as in DL07
IF count NE 0 then begin
gamma=val[ind[0]]
ENDIF
ind=where(key EQ 2,count) ;gamma parameter as in DL07
IF count NE 0 then begin
alpha=val[ind[0]]
ENDIF
ind=where(key EQ 3,count) ;gamma parameter as in DL07
IF count NE 0 then begin
Umin=val[ind[0]]
ENDIF
ind=where(key EQ 4,count) ;gamma parameter as in DL07
IF count NE 0 then begin
Umax=val[ind[0]]
ENDIF
ENDIF
;This is the average radiation field in DL07
Umean=(1.-gamma)*Umin+gamma*Umin*alog(Umax/Umin)/(1.-Umin/Umax)
lambir=dustem_get_wavelengths()
Nwavs=n_elements(lambir)
spec=fltarr(Nwavs,3)
spec2=fltarr(Nwavs,3)
;==== Run dustem with G0(a.k.a. U)=1.
param_desc=['(*!dustem_params).G0']
param_values=[Umin]
dustem_init_params,!dustem_model,param_desc,param_values
;stop
st=dustem_run(param_values)
spec1=st.sed.em_tot
spec=(1.-gamma)*spec1
;==== add contribution of dust seeing a power law of U
spec2[*]=0.
;NUs=3L
NUs=10L ;This seems sufficient
;NUs=50L
;U values over which integration will be performed
Us=range_gen(NUs,[Umin,Umax],/log)
specs=fltarr(NUs,Nwavs)
;stop
ffacts=fltarr(NUs)
IF gamma NE 0. THEN BEGIN
fact=gamma*(alpha-1.)/(Umin^(1.-alpha)-Umax^(1.-alpha))
FOR i=0L,NUs-1 DO BEGIN
param_values[0]=Us[i];^(-1.*alpha)
;stop
st=dustem_run(param_values)
ffacts[i]=Us[i]^(-1.*alpha)*fact
specs[i,*]=ffacts[i]*st.sed.em_tot
ENDFOR
;integrate over U
;stop
FOR j=0L,Nwavs-1 DO BEGIN
integ=integral(Us,reform(specs[*,j]),Umin,Umax)
spec[j]=spec[j]+integ[0]
ENDFOR
ENDIF
;stop
;This is to update keyword st with final spectrum. Useless, as there is no st keyword !
;Not necessary here, but to be used in plugins with scope 'REPLACE_FORTRAN' which actually do not call the Fortran at all.
;Ngrains=3
;st=dustem_make_empty_fortran_st(Ngrains)
st.sed.em_tot=spec
;dustem_print_params
;stop
the_end:
scope='REPLACE_FORTRAN'
;parameter tags.
paramtag=['gamma','alpha','Umin','Umax']
RETURN,spec
END