dustem_read_all.pro
9.28 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
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
FUNCTION dustem_read_all,dir_in,silent=silent,help=help
;+
; NAME:
; dustem_read_all
; PURPOSE:
; Manages reading DustEM input data about properties of grain mixture
; CATEGORY:
; DustEM, Distributed, High-Level, Initialization
; CALLING SEQUENCE:
; st=dustem_read_all(dir_in,silent=silent,help=help)
; INPUTS:
; dir_in: input directory
; OPTIONAL INPUT PARAMETERS:
; None
; OUTPUTS:
; Structure containing information about properties of grain mixture
; OPTIONAL OUTPUT PARAMETERS:
; None
; ACCEPTED KEY-WORDS:
; silent = If set, function is silent
; 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
; PROCEDURE:
; None
; EXAMPLES
; st=dustem_read_all(!dustem_wrap_soft_dir)
; MODIFICATION HISTORY:
; Written by JP Bernard 2007
; 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_read_all'
st=0.
goto,the_end
ENDIF
;== PUBLIC RELEASE VERSION
IF !dustem_which EQ 'WEB3p8' or !dustem_which EQ 'RELEASE' THEN BEGIN
st=dustem_read_all_release(dir_in,silent=silent)
;goto,ionfrac
goto, the_end
ENDIF
;== IAS VERSION -- TO BE DEPRECATED IN FUTURE DISTRIB
IF !dustem_which EQ 'VERSTRAETE' THEN BEGIN
st=dustem_read_all_lv(dir_in,silent=silent)
goto,ionfrac
ENDIF
;== IRAP/IPAG VERSIONS -- TO BE DEPRECATED IN FUTURE DISTRIB
;== note: the output keywords are set only for LV version
file=dir_in+'GRAIN.DAT'
st_grains=dustem_read_grain(file,silent=silent)
file=dir_in+'TAILLE.DAT' ;I put this back to its original form
st_taille=dustem_read_taille(file,silent=silent)
file=dir_in+'PROPMASS.DAT'
st_propmass=dustem_read_propmass(file,silent=silent)
file=dir_in+'SPECEM.DAT'
st_specem=dustem_read_specem(file,silent=silent)
file=dir_in+'CALOR.DAT'
st_calor=dustem_read_calor(file,silent=silent)
file=dir_in+'ISRF.DAT'
st_isrf=dustem_read_isrf(file,silent=silent,Nisrf=Nisrf)
IF !dustem_which EQ 'DESERT' THEN BEGIN
file=dir_in+'GEMISSIV.DAT'
st_gemissiv=dustem_read_gemissiv(file,silent=silent)
file=dir_in+'QABS.DAT'
st_qabs=dustem_read_qabs_desert(file,silent=silent)
file=dir_in+'IONFRAC.DAT'
st_ionfrac=dustem_read_ionfrac(file,silent=silent)
ENDIF
IF !dustem_which EQ 'COMPIEGNE' THEN BEGIN
Qabs_files=findfile(dir_in+'QABS_*.DAT',count=Nfiles)
st_qabs=ptrarr(Nfiles)
FOR i=0,Nfiles-1 DO BEGIN
st=dustem_read_qabs(Qabs_files(i),silent=silent)
st_qabs(i)=ptr_new(st)
ENDFOR
file=dir_in+'LAMBDA.DAT'
st_lambda=dustem_read_lambda(file,silent=silent)
ENDIF
Ngrains=n_elements(st_grains)
one_stt_grains={ $
type:st_grains(0).type, $
albmax:st_grains(0).albmax, $
densgr:st_grains(0).densgr, $
ngem:st_grains(0).ngem, $
hydro:st_grains(0).hydro, $
ioni:st_grains(0).ioni, $
alpha:st_taille(0).alpha, $
tmin:st_taille(0).tmin, $
tmax:st_taille(0).tmax, $
ndiscr:st_taille(0).ndiscr, $
propmass:st_propmass(0).propmass}
sst_grains=replicate(one_stt_grains,Ngrains)
sst_grains.type=st_grains.type
sst_grains.albmax=st_grains.albmax
sst_grains.densgr=st_grains.densgr
sst_grains.ngem=st_grains.ngem
sst_grains.hydro=st_grains.hydro
sst_grains.ioni=st_grains.ioni
sst_grains.alpha=st_taille.alpha
sst_grains.tmin=st_taille.tmin
sst_grains.tmax=st_taille.tmax
sst_grains.ndiscr=st_taille.ndiscr
sst_grains.propmass=st_propmass.propmass
IF !dustem_which EQ 'DESERT' THEN BEGIN
st={Ngrains:Ngrains,grains:sst_grains, $
gemissiv:st_gemissiv,isrf:st_isrf,qabs:st_qabs,calor:st_calor,specem:st_specem,ionfrac:st_ionfrac}
ENDIF
IF !dustem_which EQ 'COMPIEGNE' THEN BEGIN
st={Ngrains:Ngrains,grains:sst_grains, $
isrf:st_isrf,qabs:st_qabs,calor:st_calor,specem:st_specem,lambda:st_lambda}
ENDIF
; according to logic, should will never reach these lines
;IF !dustem_which EQ 'VERSTRAETE' THEN BEGIN
;
;st={Ngrains:Ngrains,grains:sst_grains, $
; isrf:st_isrf,qabs:st_qabs,calor:st_calor,specem:st_specem,lambda:st_lambda}
;ENDIF
ionfrac:
IF !dustem_which EQ 'VERSTRAETE' THEN BEGIN
IF !dustem_verbose NE 0 THEN message, "============= IONization FRACtion of the PAHs =============", /info
;Get master-keywords
master_keywords = strsplit(st.keywords,' ',/regex,/extr)
;check presence of IONFRACPAH in master-keywords
tmp = where(master_keywords eq 'IONFRACPAH',test_ionfracpah)
;check presence of dustem_create_ionfrac in the *!param_desc
;JPB: added test on validity of param_desc
IF ptr_valid((*!dustem_fit).param_descs) THEN BEGIN
tmp2 = total(stregex(*(*!dustem_fit).param_descs,'ionfrac',/bool,/fold_case))
ENDIF ELSE BEGIN
tmp2=0
ENDELSE
test_dustemcreateionfrac = (tmp2 GT 0) ? 1 : 0
;if IONFRACPAH or "dustem_create_ionfrac" is present ...
IF test_ionfracpah or test_dustemcreateionfrac THEN BEGIN
IF !dustem_verbose NE 0 THEN message, "You want the IONization FRACtion of the PAHs to be computed by the model", /info
; Check if there are both PAH1 and PAH0 provided
pah0 = (where(strmid(st.grains.type,0,4) eq 'PAH0', count0))(0)
pah1 = (where(strmid(st.grains.type,0,4) eq 'PAH1', count1))(0)
IF count0 eq count1 and count0 ne 0 THEN BEGIN
IF !dustem_verbose NE 0 THEN BEGIN
message, "and you have provided me with both PAH0 and PAH1 components.", /info
message, "I will calculate the fractions.", /info
message, "I will assume you have forced them to keep the same abundance (initial values + tied).", /info
ENDIF
ENDIF ELSE BEGIN
IF !dustem_verbose NE 0 THEN BEGIN
message, "but you have NOT provided me with both PAH0 and PAH1 components.", /info
message, "Please correct and rerun the code."
ENDIF
ENDELSE
; Check presence of SIZE in master-keywords
tmp = where(master_keywords eq 'SIZE',test_size)
; Check if MIX is a keyword provided in GRAIN.DAT for PAH0 and PAH1
pah0_keyword = st.grains(pah0).flag
pah1_keyword = st.grains(pah1).flag
tmp = where(pah0_keyword eq 'MIX',test_mix_pah0)
tmp = where(pah1_keyword eq 'MIX',test_mix_pah1)
; if SIZE is not present ...
IF not test_size THEN BEGIN
; if MIX is not present for both PAH components ...
IF not test_mix_pah0 or not test_mix_pah1 THEN BEGIN
IF !dustem_verbose NE 0 THEN BEGIN
message, "You did not provide me with either SIZE master-keyword nor MIX keyword for both PAHs components.", /info
message, "I will add the MIX keyword for both PAHs components and generate 50%-50% MIX_xxx.DAT files.", /info
message, "I will assume you have forced them to have the same size distributions.", /info
ENDIF
; Add MIX keyword to both PAH components (had to concatenate them
; to order the F90 code to read them as 1 column)
if not test_mix_pah0 then $
st.grains(pah0).flag = st.grains(pah0).flag+'MIX'
if not test_mix_pah1 then $
st.grains(pah1).flag = st.grains(pah1).flag+'MIX'
; Update structure and generate MIX_xxx.DAT files
; full_st0 = ptr_new({file:getenv('DUSTEM_DAT')+'/MIX_'+st.grains(pah0).type+'.DAT',fmix:fltarr(st.grains(pah0).ndiscr)+0.5})
; full_st1 = ptr_new({file:getenv('DUSTEM_DAT')+'/MIX_'+st.grains(pah1).type+'.DAT',fmix:fltarr(st.grains(pah1).ndiscr)+0.5})
full_st0 = ptr_new({file:!dustem_dat+'/MIX_'+st.grains(pah0).type+'.DAT',fmix:fltarr(st.grains(pah0).ndiscr)+0.5})
full_st1 = ptr_new({file:!dustem_dat+'/MIX_'+st.grains(pah1).type+'.DAT',fmix:fltarr(st.grains(pah1).ndiscr)+0.5})
full_st = ptrarr(2) & full_st(0)=full_st0 & full_st(1)=full_st1
; dustem_write_mix, getenv('DUSTEM_DAT'), full_st
dustem_write_mix,!dustem_dat, full_st
st.mix(pah0) = ptr_new(dustem_read_mix((*full_st0).file))
st.mix(pah1) = ptr_new(dustem_read_mix((*full_st1).file))
!run_ionfrac = 1.
ENDIF ELSE BEGIN
IF !dustem_verbose NE 0 THEN BEGIN
message, "You provided me with the MIX keyword for both PAHs components.", /info
message, "I will use your MIX_xxx.DAT files to update the fractions.", /info
ENDIF
!run_ionfrac = 1.
ENDELSE
; if SIZE is present
ENDIF ELSE BEGIN
IF !dustem_verbose NE 0 THEN BEGIN
message, "You provided me with the SIZE master-keyword.", /info
message, "I will use your SIZE_xxx.DAT files to update the fractions.", /info
message, "I will assume you have forced them to have the same size distributions.", /info
ENDIF
!run_ionfrac = 2.
ENDELSE
;if IONFRACPAH is NOT present ...
ENDIF ELSE BEGIN
IF !dustem_verbose NE 0 THEN BEGIN
message, "You do not want the IONization FRACtion of the PAHs to be computed by the model", /info
ENDIF
!run_ionfrac = 0.
ENDELSE
IF test_dustemcreateionfrac THEN BEGIN
!run_ionfrac = -1. * !run_ionfrac
ENDIF
IF !dustem_verbose NE 0 THEN BEGIN
message, "============= IONization FRACtion of the PAHs =============", /info
wait, 5 ;might be a bit long, but I want to be sure the user sees it
ENDIF
;JPB: The above is stupid, especially before I added the !dustem_verbose test
ENDIF
the_end:
RETURN,st
END