dustem_write_isrf_lv.pro
7.07 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
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
PRO dustem_write_isrf_lv,file,st,help=help
;+
; NAME:
; dustem_write_isrf_lv
;
; PURPOSE:
; Writes file ISRF.DAT used by DustemWrap
;
; CATEGORY:
; DustEMWrap, Distributed, Low-level, Fortran
;
; CALLING SEQUENCE:
; dustem_write_isrf_lv,file,st[,/help]
;
; INPUTS:
; file : file name
; st : input structure containing ISRF information (will be modified)
;
; OPTIONAL INPUT PARAMETERS:
; None
;
; OUTPUTS:
; None
;
; OPTIONAL OUTPUT PARAMETERS:
; None
;
; ACCEPTED KEY-WORDS:
; help = If set print this help
;
; COMMON BLOCKS:
; None
;
; SIDE EFFECTS:
; Writes a file
;
; RESTRICTIONS:
; The DustEM fortran code must be installed
; The DustEMWrap IDL code must be installed
;
; PROCEDURE AND SUBROUTINES
;
; EXAMPLES
;
; MODIFICATION HISTORY:
; Written by Ilyes Choubani 2022
; 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_write_isrf_lv'
goto,the_end
ENDIF
;Replacing the old block with an iteration over the plugins' scopes
scopes = strarr(n_tags(*!dustem_plugin))
IF (tag_names(*!dustem_plugin))[0] EQ 'NONE' then goto, the_after
;Getting the actual scopes
FOR i=0L, n_tags(*!dustem_plugin)-1 DO BEGIN
if isa((*!dustem_plugin).(i).scope) THEN BEGIN
scopes[i] = (*(*!dustem_plugin).(i).scope)
ENDIF ELSE goto, the_after
ENDFOR
tstpop = where(scopes EQ 'STELLAR_POPULATION',ctpop)
tstusrisrf = where(scopes EQ 'USER_ISRF',ctusrisrf)
IF ctpop NE 0 or ctusrisrf NE 0 THEN BEGIN
G0_mathis = 0. ;Just an initialization
;The initial test was going to be on ((*!dustem_params).g0) and ((*!dustem_params).gas.g0) but it isn't enough...
;this is because the case where mathis ISRF is used with G0=1 oughtn't be forgotten and should be distinguishable
;A solution is to instead test on the parameter vectors (fixed and thawed alike)
;Locating G0 in the free parameters vector
ind_G0 = where(strupcase(strmid((*(*!dustem_fit).param_descs),5,/reverse)) EQ 'MS).G0', countg0)
IF countg0 NE 0 THEN BEGIN
IF isa((*!dustem_fit).current_param_values) THEN G0_mathis = (*(*!dustem_fit).current_param_values)[ind_G0] $
ELSE G0_mathis = (*(*!dustem_fit).param_init_values)[ind_G0]
ENDIF
;Locating G0 in the fixed parameters vector
;under the assumption that the user will always have a free parameter vector set but not necessarily a fixed one
IF isa((*!dustem_fit).fixed_param_descs) THEN BEGIN
ind_G0 = where(strupcase(strmid((*(*!dustem_fit).fixed_param_descs),5,/reverse)) EQ 'MS).G0', countg0)
IF countg0 NE 0 THEN G0_mathis = (*(*!dustem_fit).fixed_param_init_values)[ind_G0]
ENDIF
;Locating gas.G0 in the free parameters vector
ind_GASG0 = where(strupcase(strmid((*(*!dustem_fit).param_descs),5,/reverse)) EQ 'GAS.G0', countgasg0)
IF countgasg0 THEN BEGIN
IF isa((*!dustem_fit).current_param_values) THEN G0_mathis = (*(*!dustem_fit).current_param_values)[ind_GASG0] $
ELSE G0_mathis = (*(*!dustem_fit).param_init_values)[ind_GASG0]
ENDIF
;Locating gas.G0 in the fixed parameters vector
;under the assumption that the user will always have a free parameter vector set but not necessarily a fixed one
IF isa((*!dustem_fit).fixed_param_descs) THEN BEGIN
ind_GASG0 = where(strupcase(strmid((*(*!dustem_fit).fixed_param_descs),5,/reverse)) EQ 'GAS.G0', countgasg0)
IF countgasg0 THEN G0_mathis = (*(*!dustem_fit).fixed_param_init_values)[ind_GASG0]
ENDIF
;Getting the ISRF default wavelengths
;using this structure
st=((*!dustem_params).isrf)
final_isrf = fltarr(n_elements(st)) ;Initializing the ISRF vector
Ncomments=3 ;default number of comments in the .DAT file (minus the ISRF used)
c=strarr(Ncomments)
c(0)='# DUSTEM: exciting radiation field featuring'
;stop
IF ctpop NE 0 THEN final_isrf=final_isrf+(*(*!dustem_plugin).(tstpop).spec)
IF ctusrisrf NE 0 THEN final_isrf=final_isrf+(*(*!dustem_plugin).(tstusrisrf).spec)
IF ctpop NE 0 and ctusrisrf EQ 0 THEN c(1) = '# STELLAR POPULATION ISRF'
IF ctpop EQ 0 and ctusrisrf NE 0 THEN c(1) = '# USER-DEFINED ISRF'
IF ctpop NE 0 and ctusrisrf NE 0 THEN BEGIN
Ncomments+=1
c=strarr(Ncomments)
c(0)='# DUSTEM: exciting radiation field featuring'
c(1) = '# STELLAR POPULATION ISRF'
c(2) = '# USER-DEFINED ISRF'
ENDIF
IF G0_mathis EQ 0 THEN G0_mathis=1.
;the below condition will always be valid
IF G0_mathis NE 0 THEN BEGIN
;This is reading the mathis isrf from the fortran release
;Let's hope they keep releasing it otherwise we'll have to include it ourselves
ma_isrf_dir=!dustem_soft_dir+'data/ISRF_MATHIS.DAT'
ma_isrf=dustem_read_isrf(ma_isrf_dir)
final_isrf=final_isrf/G0_mathis[0]+ma_isrf.isrf
Ncomments=4
c=strarr(Ncomments)
c(0)='# DUSTEM: exciting radiation field featuring'
IF ctpop NE 0 and ctusrisrf EQ 0 THEN c(1) = '# STELLAR POPULATION ISRF'
IF ctpop EQ 0 and ctusrisrf NE 0 THEN c(1) = '# USER-DEFINED ISRF'
IF ctpop NE 0 and ctusrisrf NE 0 THEN BEGIN
Ncomments+=2
c=strarr(Ncomments)
c(0)='# DUSTEM: exciting radiation field featuring'
c(1) = '# STELLAR POPULATION ISRF'
c(2) ='# USER-DEFINED ISRF'
c(3) = '# MATHIS ISRF'
ENDIF ELSE c(Ncomments-3) = '# MATHIS ISRF'
ENDIF
st.isrf = final_isrf
;Last .DAT header comments lines
c(Ncomments-2)='# Nbr of points'
c(Ncomments-1)='# wave (microns), 4*pi*Inu (erg/cm2/s/Hz)'
;internally displaying the 'G0' of the new ISRF (integral(ISRF)/integral(MATHIS_ISRF))
; un=uniq(ma_isrf.isrf)
;
; int_mathis=INT_TABULATED((ma_isrf.lambisrf)[un],(ma_isrf.isrf)[un])
;
;
; ;integrating the composite isrf
; int_isrf=INT_TABULATED((ma_isrf.lambisrf)[un],(st.isrf)[un])
; print, 'G0_stellar_population='
; print, int_isrf/(int_mathis)
;WRITING THE NEW ISRF
file=!dustem_dat+'data/ISRF.DAT'
openw,unit,file,/get_lun
FOR i=0,Ncomments-1 DO BEGIN
printf,unit,c(i)
ENDFOR
n_waves=n_elements(st)
printf,unit,n_waves
FOR i=0L,n_waves-1 DO BEGIN
printf,unit,st(i).lambisrf,st(i).isrf
ENDFOR
close,unit
free_lun,unit
goto, the_end
ENDIF
the_after:
;Defaut ISRF writing (MATHIS)
;This was changed because 6 comments aren't needed
;Ncomments=6
Ncomments=4
c=strarr(Ncomments)
;First lines of the ISRF.DAT file
c(0)='# DUSTEM: exciting radiation field featuring'
c(1)='# Mathis ISRF'
;c(2)='# Blackbody with T= 0.0000E+00'
;c(3)='# dilution factor wdil= 1.0000E+00'
c(2)='# Nbr of points'
c(3)='# wave (microns), 4*pi*Inu (erg/cm2/s/Hz)'
openw,unit,file,/get_lun
FOR i=0,Ncomments-1 DO BEGIN
printf,unit,c(i)
ENDFOR
n_waves=n_elements(st)
printf,unit,n_waves
FOR i=0L,n_waves-1 DO BEGIN
printf,unit,st(i).lambisrf,st(i).isrf
ENDFOR
close,unit
free_lun,unit
the_end:
END