dustem_write_isrf_release.pro
5.78 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
PRO dustem_write_isrf_release,file,st,help=help
;+
; NAME:
; dustem_write_isrf_release
;
; PURPOSE:
; Writes file ISRF.DAT used by DustemWrap
;
; CATEGORY:
; DustEMWrap, Distributed, Low-level, Fortran
;
; CALLING SEQUENCE:
; dustem_write_isrf_release,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_release'
goto,the_end
ENDIF
;LOCATING THE ISRF PLUGINS
;IF PLUGINS AREN'T SET. WRITE THE DEFAULT DUSTEM ISRF.DAT FILE
Nplugins=n_elements(*!dustem_plugin)
scopes = strarr(Nplugins)
IF (*!dustem_plugin)[0].name EQ 'NONE' THEN GOTO,write_isrf_file
;IF PLUGINS ARE SET, LOOP OVER THEIR SCOPES.
;Getting the actual scopes
FOR i=0L, Nplugins-1 DO BEGIN
;IF isa((*!dustem_plugin).(i).scope) THEN BEGIN
scopes[i] = (*!dustem_plugin)[i].scope
;ENDIF ELSE goto, the_after
ENDFOR
Ncomments=3 ;default number of comments in the .DAT file (minus the ISRF used)
c=strarr(Ncomments)
c[*]='#'
c[0]='# DUSTEM_WRAP: Interstellar radiation field '
;Testing for the two current ISRF plugin scopes
ind_replace_isrf = where(scopes EQ 'REPLACE_ISRF',count_replace_isrf)
ind_add_isrf = where(scopes EQ 'ADD_ISRF',count_add_isrf)
;stop
IF count_add_isrf NE 0 OR count_replace_isrf NE 0 THEN BEGIN
lambir=dustem_get_wavelengths()
;IF ctpop NE 0 or ctusrisrf NE 0 THEN BEGIN
;NB: Since we're using G0_mathis=1 by default now. We will test on the !dustem_params sysvar instead of the parameter description vectors
G0_mathis = 1. ;Default value
IF ((*!dustem_params).g0) NE 1. THEN BEGIN
G0_mathis = ((*!dustem_params).g0)
ENDIF
; Fortran now uses negative numbers to specify Gas.dat G0
; If positive G0 present in Gas.dat, it takes precedence over the G0 value in GRAIN.DAT
; IF ((*!dustem_params).gas.g0) NE 1. then G0_mathis = float((*!dustem_params).gas.g0)
IF ((*!dustem_params).gas.g0) gt 0. THEN BEGIN
G0_mathis = float((*!dustem_params).gas.g0)
ENDIF
;Getting the ISRF default wavelengths
;using this structure
;st=((*!dustem_params).isrf)
final_isrf = fltarr(n_elements(st)) ;Initializing the ISRF vector
;=== ADD to the ISRF
;IF ctpop NE 0. THEN final_isrf=final_isrf+(*(*!dustem_plugin)[tstpop].spec)
;JPB: This line below can have an undefined plugin ISRF if the REPLACE_ISRF plugin has not been called earlier except with /scope or paramtags keywords
IF count_replace_isrf NE 0. THEN BEGIN
IF ptr_valid((*!dustem_plugin)[ind_replace_isrf].spec) THEN BEGIN
c[1]='# DUSTEM_WRAP: replaced by radiation field from plugin '+(*!dustem_plugin)[ind_replace_isrf].name
;final_isrf=final_isrf+(*(*!dustem_plugin)[tstusrisrf].spec)
;final_isrf=interpol(*(*!dustem_plugin)[ind_replace_isrf].spec,lambir,st.lambisrf)
final_isrf=*(*!dustem_plugin)[ind_replace_isrf].spec
ENDIF ELSE BEGIN
message,'Replacement ISRF not found. Using Mathis',/continue
;stop
mathis_isrf_datfile=!dustem_soft_dir+'data/ISRF_MATHIS.DAT'
sst=dustem_read_isrf(mathis_isrf_datfile)
final_isrf=sst.isrf
;stop
;LEFT blank intensionally
ENDELSE
ENDIF
IF count_add_isrf NE 0 THEN BEGIN
IF ptr_valid((*!dustem_plugin)[ind_add_isrf].spec) THEN BEGIN
c[2]='# DUSTEM_WRAP: added radiation field from plugin '+(*!dustem_plugin)[ind_replace_isrf].name
final_isrf=final_isrf+(*(*!dustem_plugin)[ind_add_isrf].spec)
;final_isrf=final_isrf+interpol(*(*!dustem_plugin)[ind_add_isrf].spec,lambir,st.lambisrf)
ENDIF ELSE BEGIN
;stop
;LEFT blank intensionally
ENDELSE
ENDIF
;JPB: I don't think this loop is in fact needed, as the default value (without plugins) is the Mathis ISRF and ADD_ISRF will just add to this
;IF G0_mathis GT 1E-10 THEN BEGIN ;discussed with JP
; ;This is reading the mathis isrf from the fortran release
; mathis_isrf_datfile=!dustem_soft_dir+'data/ISRF_MATHIS.DAT'
; mathis_isrf=dustem_read_isrf(mathis_isrf_datfile)
; final_isrf=final_isrf/G0_mathis+mathis_isrf.isrf ;This is because Fortran uses ISRF*G0 which is now final_isrf+G0*mathis_isrf
; c=[c,'# Mathis ISRF']
; Ncomments=Ncomments+1
;ENDIF ELSE BEGIN
; ;final_isrf = final_isrf/G0_mathis ;This is because Fortran uses ISRF*G0
;ENDELSE
;=== set the ISRF in the st structure
st.isrf = final_isrf
ENDIF ELSE BEGIN ;This is Mathis field only
;==== default comment lines for ISRF.DAT
Ncomments=4
c=strarr(Ncomments)
c[1]='# Mathis ISRF'
ENDELSE
c[Ncomments-2]='# Nbr of points'
c[Ncomments-1]='# wave (microns), 4*pi*Inu (erg/cm2/s/Hz)'
;==== Write the ISRF file
write_isrf_file:
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
;print,(*!dustem_params).g0
;print,(*!dustem_params).gas.g0
;cgplot,st.lambisrf,st.isrf,/xlog,/ylog
;stop
the_end:
END