dustem_write_isrf_lv.pro
6.57 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
PRO dustem_write_isrf_lv,file,st,help=help
;+
; NAME:
; dustem_write_isrf_lv
; PURPOSE:
; Writes file ISRF.DAT used by DustemWrapper
; CATEGORY:
; Dustem
; CALLING SEQUENCE:
; dustem_write_isrf_lv,file,st[,/help]
; INPUTS:
; file : file name
; st : input ISRF structure
; OPTIONAL INPUT PARAMETERS:
; None
; 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 dustem idl wrapper must be installed
; PROCEDURE:
; Write ISRF for the DustEmWrap
;-
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).spec) 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 G0_mathis = (*(*!dustem_fit).param_init_values)[ind_G0]
;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 G0_mathis = (*(*!dustem_fit).param_init_values)[ind_GASG0]
;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=4 ;default number of comments in the .DAT file
c=strarr(Ncomments)
c(0)='# DUSTEM: exciting radiation field featuring'
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 NE 0 THEN BEGIN ;storing the mathis isrf in variable mathis_isrf if the user wants to use it
;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)
stop
final_isrf=final_isrf/G0_mathis[0]+ma_isrf.isrf
Ncomments+=1
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+=1
c=strarr(Ncomments)
c(0)='# DUSTEM: exciting radiation field featuring'
c(1) = '#STELLAR POPULATION ISRF'
c(2) = '#MATHIS ISRF'
c(3) ='#USER-DEFINED ISRF'
ENDIF
c(Ncomments-3) = '#MATHIS ISRF'
ENDIF
;Will this enable us to have the current 'updated' ISRF through the !dustem_params system variable?
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