dustem_plugin_dl07_isrf_model.pro
3.39 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
FUNCTION dustem_plugin_dl07_isrf_model, key=key, val=val, scope=scope, paramtag=paramtag,help=help,spec=spec,Umean=Umean
;+
; NAME:
; dustem_plugin_dl07_isrf_model
;
; PURPOSE:
;
;
; CATEGORY:
; DUSTEM Wrapper
;
; CALLING SEQUENCE:
; sed=dustem_plugin_dl07_isrf_model(key=,val=,/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:
; None
;
; 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
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
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)
st=dustem_run(param_values)
ffacts[i]=Us[i]^(-1.*alpha)*fact
specs[i,*]=ffacts[i]*st.sed.em_tot
ENDFOR
;integrate over U
FOR j=0L,Nwavs-1 DO BEGIN
integ=integral(Us,reform(specs[*,j]),Umin,Umax)
spec[j]=spec[j]+integ[0]
ENDFOR
;FOR i=0L,Nus-1 DO BEGIN
; print,Us[i],ffacts[i]
;ENDFOR
ENDIF
;This is to update keyword st with final spectrum
st.sed.em_tot=spec
the_end:
scope='REPLACE_FORTRAN'
;parameter tags.
paramtag=['gamma','alpha','Umin','Umax']
RETURN,spec
END