dustem_init.pro
12.6 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
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
PRO dustem_init,dir=dir,wraptest=wraptest,plot_it=plot_it,model=model,help=help,noerase=noerase,st_model=st_model,pol=pol,verbose=verbose
;+
; NAME:
; dustem_init
; CALLING SEQUENCE:
; dustem_init
; PURPOSE:
; Defines system variables to run the DUSTEM code
; INPUTS:
; model = Selects one of the dust mixture used by dustem
; 'COMPIEGNE_ETAL2010' from Compiegne et al 2010 (default)
; 'DBP90' from Desert et al 1990
; 'DL01' from Draine & Li 2001
; 'DL07' from Draine & Li 2007
; dir = overrides default directory where to read data
; wraptest = test running dustem f90 through the wrapper
; plot_it = plots run result (only if wraptest)
; help = If set print this help
; OPTIONAL INPUT:
; None
; OUTPUTS:
; None
; PROCEDURE AND SUBROUTINE USED
; 1/ The folowing system variables are initialized to null values
; or empty pointers:
; !dustem_fit
; !dustem_data
; !dustem_inputs
; !dustem_params
; !run_ionfrac
; !dustem_idl_continuum
; !dustem_filters
; !dustem_verbose
; !dustem_show_plot
; 2/ optionally runs the fortran code (if wraptest)
; This step calls the Fortran code to run a standard model
; 3/ !dustem_filters is initialized calling dustem_filter_init.pro.
; RESTRICTIONS:
; The dustem fortran code must be installed
; The dustem idl wrapper must be installed
; EXAMPLES
; dustem_init,/wrap,/plot_it
; MODIFICATION HISTORY:
; See cvs logs
; Written by JPB,NF,DP Jan-2007
; see evolution details on the dustem cvs maintained at CESR
;-
IF keyword_set(help) THEN BEGIN
doc_library,'dustem_init'
goto,the_end
ENDIF
defsysv,'!dustem_verbose',0 ;1=verbose
IF keyword_set(verbose) THEN !dustem_verbose=1
defsysv,'!run_pol',0 ;0=no polar
IF keyword_set(pol) THEN !run_pol=1
;=== color correction system variables
defsysv,'!dustem_previous_cc',ptr_new() ;store previous color correction values
defsysv,'!dustem_do_cc',1 ;indicates to dustem if color corrections are needed. Will be changed when fitting
defsysv,'!dustem_never_do_cc',0 ;set to 1 to never do color corrections (not recommanded, but faster).
;This is for polarization
defsysv, '!run_circ', 0.
defsysv, '!run_anis', 0.
defsysv, '!run_rrf', 0.
defsysv, '!run_univ', 0.
defsysv, '!run_lin', 0.
defsysv, '!dustem_redshift', 0. ;in case the SED needs to be shifted for a non zero redshift
;==== define a stellar modyfier to the ISRF
toto=dustem_define_isrf_stellar_modifyer_variable()
file=!dustem_wrap_soft_dir+'instrument_description.xcat'
st=read_xcat(file,/silent)
defsysv,'!dustem_instrument_description',st
defsysv,'!indef', exists=exists_indef
IF NOT exists_indef THEN defsysv,'!indef',-32768.
;define !dustem_fit
dustem_fit_st={data:ptr_new(),wavelength:ptr_new(), $
param_descs:ptr_new(),param_init_values:ptr_new(),param_func:ptr_new(), $
fixed_param_descs:ptr_new(),fixed_param_init_values:ptr_new(), $
chi2:0.D0,rchi2:0.D0, $
current_param_values:ptr_new(), current_param_errors:ptr_new() $
}
defsysv, '!dustem_fit', ptr_new(dustem_fit_st) ;Data to fit
if keyword_set(pol) then begin
defsysv, '!dustem_data', { $ ;Data to fit
sed: ptr_new(), $
ext: ptr_new(), $
polext: ptr_new(), $
polsed: ptr_new(), $
polfrac: ptr_new(), $ ; after this line you can put the new variables to fit
qsed: ptr_new(), $
used: ptr_new(), $
qext: ptr_new(), $
uext: ptr_new() $
}
; rchi2_weight = {sed: 0.,ext: 0.,polext: 0.,polsed: 0.,polfrac:0.} ;old line
rchi2_weight = {sed: 0.,ext: 0.,polext: 0.,polsed: 0.,polfrac:0.,qsed:0.,used:0.,qext:0.,uext:0.}
endif else begin
defsysv, '!dustem_data', { $ ;Data to fit
sed: ptr_new(), $
ext: ptr_new() $
}
rchi2_weight = {sed: 0.,ext: 0.}
endelse
tagnames = tag_names(!dustem_data)
instr="defsysv, '!fit_rchi2_weight', {"
FOR i = 0, n_elements(tagnames)-1 DO BEGIN
instr+=tagnames(i)+': rchi2_weight.'+tagnames(i)
IF i NE n_elements(tagnames)-1 THEN instr+=','
ENDFOR
instr+='}'
toto=execute(instr)
;IDL> help,*!dustem_data,/str
;** Structure <2173b20>, 4 tags, length=448, data length=448, refs=1:
; INSTRU_NAMES STRING Array[14]
; FILT_NAMES STRING Array[14]
; VALUES FLOAT Array[14]
; SIGMA FLOAT Array[14]
defsysv,'!dustem_inputs',ptr_new() ;dustem input files
;=== This is not a system variable anymore
;defsysv, '!dustem_ext', ptr_new() ;Contains the values of extinction
;=== This is not a system variable anymore
;defsysv, '!dustem_specem', ptr_new() ;Contains the values of emission
defsysv, '!dustem_params', ptr_new() ;Contains the values of all Desert Model parameters (except extinction !)
;IDL> help,*!dustem_params,/str
;** Structure <2926e10>, 7 tags, length=13312, data length=13282, refs=1:
; GRAINS STRUCT -> <Anonymous> Array[3]
; GEMISSIV STRUCT -> <Anonymous> Array[300]
; ISRF STRUCT -> <Anonymous> Array[64]
; QABS STRUCT -> <Anonymous> Array[64]
; CALOR STRUCT -> <Anonymous> Array[5]
; SPECEM STRUCT -> <Anonymous> Array[1]
; IONFRAC STRUCT -> <Anonymous> Array[10]
;defsysv, '!dustem_fit_sed',ptr_new() ;Best fit SED
;The following is not used in the WEB3p8 version of the code
;Only in the VERSTRAETE version of the code.
;IF getenv('DUSTEM_WHICH') EQ 'VERSTRAETE' THEN BEGIN
defsysv, '!run_ionfrac', 0. ; 0:no call to dustem_create_ionfrac UNLESS it is called as a param
; 1:call to dustem_create_ionfrac with the MIX_xxx.DAT files
; 2:call to dustem_create_ionfrac with the SIZE_xxx.DAT files
;ENDIF
;WE DO NOT NEED THESE VARIABLES ANYMORE BUT HAVEN'T TESTED YET
defsysv, '!dustem_idl_continuum', 0.
defsysv, '!dustem_idl_freefree', 0.
defsysv, '!dustem_idl_synchrotron', 0.
;
;
; ;a new plugin will have to be added to this list
; stt={synchrotron:ptr_new(),free_free:ptr_new(), nir_continuum:ptr_new(), grey_body:ptr_new()}
; defsysv, '!dustem_plugins', ptr_new(stt)
defsysv,'!dustem_filters',ptr_new() ;filter information (filled by dustem_filter_init.pro)
;; IDL> help,*!dustem_filters,/str
;; ** Structure <1b95608>, 17 tags, length=8232, data length=8118, refs=1:
;; IRAC STRUCT -> <Anonymous> Array[1]
;; MIPS STRUCT -> <Anonymous> Array[1]
;; MSX STRUCT -> <Anonymous> Array[1]
;; IRAS STRUCT -> <Anonymous> Array[1]
;; DIRBE STRUCT -> <Anonymous> Array[1]
;; SPM STRUCT -> <Anonymous> Array[1]
;; ISOCAM STRUCT -> <Anonymous> Array[1]
;; ISOPHOTP STRUCT -> <Anonymous> Array[1]
;; ISOPHOTC STRUCT -> <Anonymous> Array[1]
;; EFIRAS STRUCT -> <Anonymous> Array[1]
;; ARCHEOPS STRUCT -> <Anonymous> Array[1]
;; HFI STRUCT -> <Anonymous> Array[1]
;; LFI STRUCT -> <Anonymous> Array[1]
;; WMAP STRUCT -> <Anonymous> Array[1]
;; SPIRE STRUCT -> <Anonymous> Array[1]
;; PACS STRUCT -> <Anonymous> Array[1]
;; AKARI STRUCT -> <Anonymous> Array[1]
defsysv,'!dustem_verbose',0 ;1=verbose
defsysv,'!dustem_show_plot',1 ;0=show no plot
;apparently not needed anymore. Was removed.
;defsysv,'!ki',0. ;apparently needed by mpfit, but what is it ?
;== This is for the default mode
dustem_inputs={GRAIN:'GRAIN.DAT',ISRF:'ISRF.DAT',ALIGN:'ALIGN.DAT'}
!dustem_inputs=ptr_new(dustem_inputs)
;===This is the Desert option
IF keyword_set(model) THEN BEGIN
CASE model OF
'COMPIEGNE_ETAL2010':BEGIN
(*!dustem_inputs).grain='GRAIN_MC10.DAT'
END
'DBP90':BEGIN
(*!dustem_inputs).grain='GRAIN_DBP90.DAT'
END
'DL01':BEGIN
(*!dustem_inputs).grain='GRAIN_DL01.DAT'
END
'DL07':BEGIN
(*!dustem_inputs).grain='GRAIN_DL07.DAT'
END
'AJ13':BEGIN
(*!dustem_inputs).grain='GRAIN_J13.DAT'
END
'G17_MODELA':BEGIN
(*!dustem_inputs).grain='GRAIN_G17_MODELA.DAT'
(*!dustem_inputs).align='ALIGN_G17_MODELA.DAT'
END
'G17_MODELB':BEGIN
(*!dustem_inputs).grain='GRAIN_G17_MODELB.DAT'
(*!dustem_inputs).align='ALIGN_G17_MODELB.DAT'
END
'G17_MODELC':BEGIN
(*!dustem_inputs).grain='GRAIN_G17_MODELC.DAT'
(*!dustem_inputs).align='ALIGN_G17_MODELC.DAT'
END
'G17_MODELD':BEGIN
(*!dustem_inputs).grain='GRAIN_G17_MODELD.DAT'
(*!dustem_inputs).align='ALIGN_G17_MODELD.DAT'
END
ELSE:BEGIN
message,'model '+model+' unknown',/continue
message,'Known models are COMPIEGNE_ETAL2010,DBP90,DL01,DL07,AJ13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue
stop
END
ENDCASE
ENDIF
;Initialize !dustem_params structure
IF not keyword_set(dir) THEN BEGIN
; CASE getenv('DUSTEM_WHICH') OF
CASE !dustem_which OF
'DESERT':dir_in=!dustem_wrap_soft_dir+'/Data/DESERT_POST_ISO_MORE/'
'COMPIEGNE': dir_in=!dustem_wrap_soft_dir+'/Data/MC_DAT/'
; 'VERSTRAETE': dir_in=!dustem_wrap_soft_dir+'/Data/d_3.5/'
'VERSTRAETE': dir_in=!dustem_wrap_soft_dir+'/Data/LV_DAT/'
'WEB3p8': dir_in=!dustem_soft_dir
;'WEB4p3': dir_in=!dustem_soft_dir ; I ADDED THIS
;'WEB4p3': dir_in='/Users/ilyeschoubani/Soft_Librairies/DUSTEM/Data/les_DAT/' ;dir_in=!dustem_soft_dir+ '/Data/' ; I added this
ELSE: dir_in=!dustem_wrap_soft_dir+'/Data/les_DAT/';!dustem_wrap_soft_dir+'/Data/les_DAT/';I CHANGED THIS!!!!
ENDCASE
ENDIF ELSE BEGIN
dir_in=dir
ENDELSE
;stop
;dir_in=getenv('DUSTEM_SOFT_DIR')+'/src/'+dir
;dir_in=!dustem_soft_dir+'/src/'+dir
st_model=dustem_read_all(dir_in,/silent)
!dustem_params=ptr_new(st_model)
;=== create dynamical storage if needed
IF not file_test(!dustem_dat,/dir) THEN BEGIN
spawn,'mkdir '+!dustem_dat
ENDIF
IF not file_test(!dustem_res,/dir) THEN BEGIN
spawn,'mkdir '+!dustem_dat
ENDIF
;=== remove dynamical storage directories
;=== This is needed for iterative use of DUSTEM from IDL
if not keyword_set(noerase) then begin
CASE !dustem_which OF
'VERSTRAETE':BEGIN
IF file_test(!dustem_dat,/dir) THEN BEGIN
spawn,'rm -rf '+!dustem_dat+'/les_DAT'
spawn,'rm -rf '+!dustem_dat+'/les_QABS'
spawn,'rm -rf '+!dustem_dat+'/les_CAPA'
ENDIF
IF file_test(!dustem_res,/dir) THEN BEGIN
spawn,'rm -rf '+!dustem_res+'/les_RES'
ENDIF
END
'WEB3p8':BEGIN ; I CHANGED THIS
IF file_test(!dustem_dat,/dir) THEN BEGIN
spawn,'rm -rf '+!dustem_dat+'/data'
spawn,'rm -rf '+!dustem_dat+'/hcap'
spawn,'rm -rf '+!dustem_dat+'/oprop'
spawn,'rm -rf '+!dustem_dat+'/out'
ENDIF
END
ELSE: BEGIN
END
ENDCASE
;=== create storage directories
;=== This is needed for iterative use of DUSTEM from IDL
CASE !dustem_which OF
'VERSTRAETE':BEGIN
IF file_test(!dustem_dat,/dir) THEN BEGIN
spawn,'mkdir '+!dustem_dat+'/les_DAT'
spawn,'mkdir '+!dustem_dat+'/les_QABS'
spawn,'mkdir '+!dustem_dat+'/les_CAPA'
ENDIF
IF file_test(!dustem_res,/dir) THEN BEGIN
spawn,'mkdir '+!dustem_res+'/les_RES'
ENDIF
END
'WEB3p8':BEGIN ; IN CHANGED THIS AS WELL
IF file_test(!dustem_dat,/dir) THEN BEGIN
spawn,'mkdir '+!dustem_dat+'/data'
spawn,'mkdir '+!dustem_dat+'/hcap'
spawn,'mkdir '+!dustem_dat+'/oprop'
spawn,'mkdir '+!dustem_dat+'/out'
ENDIF
END
ELSE: BEGIN
END
ENDCASE
endif ; /noerase
;=== If needed, test runing the model through the wrapper
IF keyword_set(wraptest) THEN BEGIN
dir_out=!dustem_dat
dustem_write_all,st_model,dir_out
; stop
st=dustem_run()
IF keyword_set(plot_it) THEN BEGIN
loadct,13
dustem_plot_nuinu_em,st,yr=[1e-14,1.e-6],/ysty,xr=[1,5e3],/xsty
ENDIF
ENDIF
;stop
;Initialize filters
dustem_filter_init,dir=dir
;=== test plugins
plugins_desc_file=getenv('DUSTEM_WRAP_SOFT_DIR')+'/plugins_description.xcat' ;added a slash and changed DUSTEM_SOFT_DIR to DUSTEM_WRAP_SOFT_DIR
plugins_st=read_xcat(plugins_desc_file,/silent)
Nplugins=n_elements(plugins_st)
FOR i=0L,Nplugins-1 DO BEGIN
str='res='+'dustem_'+plugins_st[i].PLUGIN_NAME+'()'
IF !dustem_verbose NE 0 THEN BEGIN
message,'Executing:'+str,/continue
ENDIF
toto=execute(str)
ENDFOR
;stop
the_end:
END