OV-GSO-DC / dustem-wrapper_idl

Commit e2ce1dee049fa0065eec7deb9eba846b1a750811

Authored by Jean-Philippe Bernard
1 parent 54f8eed5
Exists in master

defined grid structure

Showing 3 changed files with 74 additions and 98 deletions   Show diff stats
Inline Side-by-side
LabTools/IRAP/JPB/dustem_get_all_filter_names.pro
  Diff comments   View file @e2ce1de
1   -FUNCTION dustem_get_all_filter_names
2   -;filters=dustem_get_all_filter_names()
  1 +FUNCTION dustem_get_all_filter_names,all_instrus=all_instrus
3 2  
  3 +;all_filters=dustem_get_all_filter_names(all_instrus=all_instrus)
  4 +
  5 +;=== This was the structure to get filters out of *!dustem_filters.
  6 +;=== However some of them are not in filtes.xcat, which later makes dustem_filter2wav to fail.
  7 +;=== so for now, I do everything from instrument_description.xcat only
  8 +goto,from_instrument_description
4 9 instrus=tag_names(*!dustem_filters)
5 10 Ninstrus=n_elements(instrus)
6 11  
7 12 all_filters=['NONE']
  13 +all_instrus=['NONE']
8 14 FOR i=0L,Ninstrus-1 DO BEGIN
9   - all_filters=[all_filters,(*!dustem_filters).(i).filter_names]
  15 + these_filters=(*!dustem_filters).(i).filter_names
  16 + Nfilt=n_elements(these_filters)
  17 + all_filters=[all_filters,these_filters]
  18 + all_instrus=[all_instrus,replicate(instrus[i],Nfilt)]
10 19 ENDFOR
11 20 all_filters=all_filters[1:*]
  21 +all_instrus=all_instrus[1:*]
  22 +
  23 +wavs=dustem_filter2wav(all_filters)
  24 +
  25 +ind=where(wavs NE la_undef(),count)
  26 +all_filters=all_filters[ind]
  27 +all_instrus=all_instrus[ind]
  28 +
  29 +goto,the_end
  30 +;stop
  31 +
  32 +from_instrument_description:
  33 +file=!dustem_wrap_soft_dir+'instrument_description.xcat'
  34 +st=read_xcat(file,/silent)
  35 +all_filters=st.filter
  36 +all_instrus=st.instru
  37 +
  38 +the_end:
12 39  
13 40 RETURN,all_filters
14 41  
... ...
LabTools/IRAP/JPB/make_sed_phangs_tables.pro
  Diff comments   View file @e2ce1de
1 1 PRO make_sed_phangs_tables,test=test,show_seds=show_seds,grid_type=grid_type,isrf_class=isrf_class
2 2  
3   -;make_sed_phangs_tables ;That's a test
  3 +;make_sed_phangs_tables,/show ;That's a test
4 4 ;make_sed_phangs_tables,/test,grid_type=2,/show_seds,isrf_class=15
5 5 ;make_sed_phangs_tables,/test,grid_type=3,/show_seds
6 6  
... ... @@ -29,15 +29,15 @@ CASE use_grid_type OF
29 29 iv_min = [0.1]
30 30 iv_max = [10]
31 31 iv_Nvalues=[2]
32   - log=[1]
  32 + plog=[1]
33 33 fpd=['(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
34 34 '(*!dustem_params).grains(1).mdust_o_mh']
35 35 fiv=[1.e-3,1.e-2]
36 36 table_name=dir+'TEST_'+model+'_G0.fits'
37 37 dustem_init,model=model
38   - ;=== select filters to be used for the grid
39   - filters=[(*!dustem_filters).nircam.filter_names,(*!dustem_filters).miri.filter_names]
40   - dustem_make_sed_table,model,pd,iv_min,iv_max,iv_Nvalues,fpd=fpd,fiv=fiv,filename=table_name,log=log,show_seds=show_seds,/print_params,filters=filters
  38 + ;=== select filters to be used for the grid (default=ALL)
  39 + ;filters=[(*!dustem_filters).nircam.filter_names,(*!dustem_filters).miri.filter_names]
  40 + dustem_make_sed_table,model,pd,iv_min,iv_max,iv_Nvalues,fpd=fpd,fiv=fiv,filename=table_name,plog=plog,show_seds=show_seds,/print_params,filters=filters
41 41 END
42 42 1: BEGIN ;======= This is to produce a grid with the DBP90 model + Mathis field for G0 and PAH abundance
43 43 !quiet=1
... ... @@ -49,7 +49,7 @@ CASE use_grid_type OF
49 49 '(*!dustem_params).grains(1).mdust_o_mh'] ;VSG mass fraction
50 50 iv_min = [0.1,1.e-4,1.e-4]
51 51 iv_max = [100,1.e-1,1.e-1]
52   - log=[1,1,1]
  52 + plog=[1,1,1]
53 53 fpd=[]
54 54 fiv=[]
55 55 IF keyword_set(test) THEN BEGIN
... ... @@ -63,7 +63,7 @@ CASE use_grid_type OF
63 63 dustem_init,model=model
64 64 ;=== select filters to be used for the grid
65 65 filters=[(*!dustem_filters).nircam.filter_names,(*!dustem_filters).miri.filter_names]
66   - dustem_make_sed_table,model,pd,iv_min,iv_max,iv_Nvalues,fpd=fpd,fiv=fiv,filename=table_name,filters=filters,log=log,show_seds=show_seds
  66 + dustem_make_sed_table,model,pd,iv_min,iv_max,iv_Nvalues,fpd=fpd,fiv=fiv,filename=table_name,filters=filters,plog=plog,show_seds=show_seds
67 67 END
68 68 2: BEGIN ;This is DPB90 model with PHANGS ISRF classes
69 69 !quiet=1
... ... @@ -76,7 +76,7 @@ CASE use_grid_type OF
76 76 '(*!dustem_params).grains(1).mdust_o_mh'] ;VSG mass fraction
77 77 iv_min = [0.1,1.e-4,1.e-4]
78 78 iv_max = [100,1.e-1,1.e-1]
79   - log=[1,1,1]
  79 + plog=[1,1,1]
80 80 Nclass=31L
81 81 fpd=['(*!dustem_params).gas.g0','(*!dustem_params).g0','dustem_plugin_phangs_class_isrf_1'] ;ISRF class to be used
82 82 dustem_init,model=model
... ... @@ -101,7 +101,7 @@ CASE use_grid_type OF
101 101 iv_Nvalues=[50,20,20]
102 102 table_name=dir+model+'_MuseISRF_JWST_G0_YPAH_YVSG_4Phangs'+isrf_class_str+'.fits'
103 103 ENDELSE
104   - dustem_make_sed_table,model,pd,iv_min,iv_max,iv_Nvalues,fpd=fpd,fiv=fiv,filename=table_name,filters=filters,log=log,show_seds=show_seds
  104 + dustem_make_sed_table,model,pd,iv_min,iv_max,iv_Nvalues,fpd=fpd,fiv=fiv,filename=table_name,filters=filters,plog=plog,show_seds=show_seds
105 105 END
106 106 3: BEGIN ;This is DPB90 model with DL07 ISRF prescription (a la Chastanet)
107 107 !quiet=1
... ... @@ -117,7 +117,7 @@ CASE use_grid_type OF
117 117 ;iv_min = [0.1,0. ,1.e-4,1.e-4] ;does not run for values =0., because of use of pmin (division by initial value)
118 118 iv_min = [0.1,1.e-10 ,1.e-4,1.e-4]
119 119 iv_max = [100,0.02,1.e-1,1.e-1]
120   - log=[1,0,1,1]
  120 + plog=[1,0,1,1]
121 121 ;fpd=['(*!dustem_params).gas.g0', $ ;gas.G0=-1 (not used)
122 122 ; '(*!dustem_params).g0', $ ;g0=1 (not used)
123 123 ; 'dustem_plugin_dl07_isrf_model_2', $ ;alpha=1
... ... @@ -146,7 +146,7 @@ CASE use_grid_type OF
146 146 iv_Nvalues=[6,5,20,20]
147 147 table_name=dir+model+'_DL07ISRF_YPAH_YVSG_Umin_gamma_4Phangs'+'.fits'
148 148 ENDELSE
149   - dustem_make_sed_table,model,pd,iv_min,iv_max,iv_Nvalues,fpd=fpd,fiv=fiv,filename=table_name,filters=filters,log=log,show_seds=show_seds
  149 + dustem_make_sed_table,model,pd,iv_min,iv_max,iv_Nvalues,fpd=fpd,fiv=fiv,filename=table_name,filters=filters,plog=plog,show_seds=show_seds
150 150 END
151 151  
152 152 ENDCASE
... ...
src/idl/dustem_make_sed_table.pro
  Diff comments   View file @e2ce1de
... ... @@ -7,7 +7,7 @@ PRO dustem_make_sed_table,model, $
7 7 fiv=fiv, $
8 8 filename=filename, $
9 9 filters=filters, $
10   - log=log, $
  10 + plog=plog, $
11 11 show_seds=show_seds, $
12 12 print_params=print_params, $
13 13 help=help
... ... @@ -33,7 +33,7 @@ PRO dustem_make_sed_table,model, $
33 33 ; fpd : fixed parameter description
34 34 ; fiv : fixed parameter values
35 35 ; filename : output fits file name for the table (default ='/tmp/dustem_seds.fits')
36   -; log : array indicating if a given parameter is to be sampled in linear (0) or log scale (1)
  36 +; plog : array indicating if a given parameter is to be sampled in linear (0) or log scale (1)
37 37 ; show_seds : if set, plots SEDs as they are computed.
38 38 ; OUTPUTS:
39 39 ; None
... ... @@ -52,7 +52,7 @@ PRO dustem_make_sed_table,model, $
52 52 ; Filters in the table are ordered by increasing wavelength
53 53 ; EXAMPLES
54 54 ; dustem_make_sed_table,'DBP90',['(*!dustem_params).G0','(*!dustem_params).grains(0).mdust_o_mh'], $
55   -; [0.1,1.e-3],[100.,1.e-1],[10,3],filename='/tmp/IRAS_GO.fits',log=[1,0],/show_seds, $
  55 +; [0.1,1.e-3],[100.,1.e-1],[10,3],filename='/tmp/IRAS_GO.fits',plog=[1,0],/show_seds, $
56 56 ; filters=[dustem_instru2filters('IRAC'),dustem_instru2filters('SPIRE')]
57 57 ; MODIFICATION HISTORY:
58 58 ; Written by J.-Ph. Bernard (2023)
... ... @@ -76,15 +76,23 @@ dustem_init,model=model,pol=use_polarisation
76 76 !dustem_which='RELEASE'
77 77  
78 78 ;=== define default filters
79   -use_filters=dustem_get_all_filter_names()
  79 +use_filters=dustem_get_all_filter_names(all_instrus=all_instrus)
80 80 IF keyword_set(filters) THEN use_filters=filters
  81 +ind=where(use_filters NE 'SPECTRUM',count)
  82 +IF count NE 0 THEN BEGIN
  83 + use_filters=use_filters[ind]
  84 +ENDIF ELSE BEGIN
  85 + message,'No usable filter found',/info
  86 + stop
  87 +ENDELSE
  88 +
81 89 wavs=dustem_filter2wav(use_filters)
82 90 order=sort(wavs)
83 91 use_filters=use_filters[order]
84 92 Nfilters=n_elements(use_filters)
85 93  
86 94 ;==== compute all combinations of model parameter values
87   -parameter_values=dustem_param_range2param_values(iv_min,iv_max,iv_Nvalues,Nc=Nc,log=log)
  95 +parameter_values=dustem_param_range2param_values(iv_min,iv_max,iv_Nvalues,Nc=Nc,log=plog)
88 96  
89 97 ;=== initialize filters for the SED
90 98 Nfilt=n_elements(use_filters)
... ... @@ -144,13 +152,15 @@ pd=parameters_description
144 152 ;=== user to (i) change default values and (2) see the adopted values in the graphical output windows
145 153 ;fpd=parameters_description
146 154 seds=ptrarr(Nc)
  155 +em_spectra=ptrarr(Nc)
  156 +ext_spectra=ptrarr(Nc)
147 157 yrange=[1.e-5,1.e5]
148 158 colors=long(range_gen(Nc,[0,255]))
149 159 IF keyword_set(show_seds) THEN loadct,13
150 160 FOR i=0L,Nc-1 DO BEGIN
151 161 ;dustem_init,model=model,pol=use_polarisation
152 162 ;dustem_set_data,m_fit=sed,m_show=sed,x_fit=ext,x_show=ext
153   - dummy=0 ;otherwise model is not recomputed !!
  163 + st=0 ;otherwise model is not recomputed !!
154 164  
155 165 ;fpval = *parameter_values[i]
156 166 pval = *parameter_values[i]
... ... @@ -163,14 +173,15 @@ FOR i=0L,Nc-1 DO BEGIN
163 173 ;=== initialise all this information into dustemwrap's brain
164 174 dustem_init_params,model,pd,pval,fpd=fpd,fiv=fiv,pol=use_polarisation,/force_reset
165 175  
166   - IF keyword_set(print_params) THEN BEGIN
167   - dustem_print_params,0
168   - ;stop
169   - ENDIF
  176 + ;IF keyword_set(print_params) THEN BEGIN
  177 + ; dustem_print_params,0
  178 + ; ;stop
  179 + ; ENDIF
170 180  
171 181 ;=== compute the predicted SED and extinction curve for the dust-model
172   - ;=== and any plugins
173   - dustem_Ised=dustem_compute_sed(pval,st=dummy)
  182 + ;=== and any plugins
  183 + dustem_Ised=dustem_compute_sed(pval,st=st) ;is in MJy/sr/1e20 H/cm2
  184 +
174 185 dustem_Qsed=dustem_Ised*0.
175 186 dustem_Used=dustem_Ised*0.
176 187  
... ... @@ -190,93 +201,31 @@ FOR i=0L,Nc-1 DO BEGIN
190 201 ;==== fill in dependent columns of the SED.
191 202 sed=dustem_fill_sed_dependent_columns(sed)
192 203 seds[i]=ptr_new(sed)
  204 + em_spectra[i]=ptr_new(st.sed)
  205 + ext_spectra[i]=ptr_new(st.ext)
193 206 IF keyword_set(show_seds) THEN BEGIN
  207 + xrange=[0.1,1.e3]
194 208 IF i EQ 0 THEN BEGIN
195 209 tit='Grid SEDs model '+model
196 210 xtit='Wavelength [mic]'
197 211 ytit='SED [MJy/sr for NH=1.e20 H/cm^2]'
198   - cgplot,(*seds[i]).wave,(*seds[i]).stokesI,psym=-4,/ylog,yr=yrange,xtit=xtit,ytit=ytit,tit=tit,/xlog
  212 + wavs=(*seds[i]).wave
  213 + Is=(*seds[i]).stokesI
  214 + cgplot,(*seds[i]).wave,(*seds[i]).stokesI,psym=-4,/ylog,yr=yrange,xtit=xtit,ytit=ytit,tit=tit,/xlog,xrange=xrange,/xsty
199 215 ENDIF ELSE BEGIN
200 216 cgoplot,(*seds[i]).wave,(*seds[i]).stokesI,psym=-4,color=colors[i]
201 217 ENDELSE
  218 + ;stop
202 219 ENDIF
203 220 ;stop
204 221 ENDFOR
205 222 ;stop
206 223  
207   -str='one_sed={'
208   -FOR i=0L,Nparams-1 DO BEGIN
209   - str=str+'P'+strtrim(i+1,2)+':0.d0,'
210   -ENDFOR
211   -FOR i=0L,Nfilters-1 DO BEGIN
212   - ;istr=strtrim(i+1,2)
213   - istr=strtrim(use_filters[i],2)
214   - str=str+'I'+istr+':0.d0'
215   - IF use_polarization EQ 1 THEN str=str+',Q'+istr+':0.d0,U'+istr+':0.d0'
216   - ;IF i NE Nfilters-1 THEN str=str+',' ELSE str=str+'}'
217   - str=str+','
218   -ENDFOR
219   -str=str+'total:0.d0}'
220   -toto=execute(str)
221   -
222   -stop
223   -Noffset=1
224   -IF use_polarization EQ 1 THEN Noffset=3
225   -seds_array=replicate(one_sed,Nc)
226   -FOR i=0L,Nc-1 DO BEGIN
227   - FOR j=0L,Nparams-1 DO BEGIN
228   - seds_array[i].(j)=(*parameter_values[i])[j]
229   - ENDFOR
230   - sum=0.d0
231   - FOR j=0L,Nfilters-1 DO BEGIN
232   - pos=Nparams+j*Noffset
233   - seds_array[i].(pos)=((*seds[i]).stokesI)[j]
234   - sum=sum+((*seds[i]).stokesI)[j] ;======== compute the sum of SEDs
235   - ENDFOR
236   - seds_array[i].total=sum
237   -ENDFOR
238   -
239   -;======== save SED tables
240   -
241   -fits_file='/tmp/dustem_seds.fits'
242   -IF keyword_set(filename) THEN fits_file=filename
243   -
244   -mwrfits,seds_array,fits_file,header,/create
245   -
246   -st=mrdfits(fits_file,0,h)
247   -st1=mrdfits(fits_file,1,h1)
248   -
249   -sxaddpar,h1,'MODEL',model,'dustemwrap model used',after='TFORM'+strtrim(Nparams+Nfilters,2)
250   -sxaddpar,h1,'POLAR',use_polarization,'dustemwrap polarization ?',after='MODEL'
251   -FOR j=0L,Nparams-1 DO BEGIN
252   - value=sxpar(h1,'TTYPE'+strtrim(j+1,2))
253   - sxaddpar,h1,'TTYPE'+strtrim(j+1,2),value,'dustemwrap '+parameters_description[j]
254   - sxaddpar,h1,'PNAME'+strtrim(j+1,2),parameters_description[j],'dustemwrap param name'
255   - ;iv_min,iv_max,iv_Nvalues,Nc=Nc,log=log
256   - sxaddpar,h1,'PMIN'+strtrim(j+1,2),iv_min[j],'minimum param value'
257   - sxaddpar,h1,'PMAX'+strtrim(j+1,2),iv_max[j],'maximum param value'
258   - sxaddpar,h1,'PNVAL'+strtrim(j+1,2),iv_Nvalues[j],'Number of param values'
259   - sxaddpar,h1,'PLOG'+strtrim(j+1,2),log[j],'1 if log scale sampling'
260   - last_pname='PLOG'+strtrim(j+1,2)
261   - IF j EQ 0 THEN first_pname='PNAME'+strtrim(j+1,2)
262   -ENDFOR
263   -j0=j
264   -sxaddpar,h1,'COMMENT','**** Dustemwrap PARAMETERS ***',before=first_pname
265   -sxaddpar,h1,'COMMENT','**** Dustemwrap FILTERS ***',after=last_pname
266   -FOR j=0L,Nfilters-1 DO BEGIN
267   - jj=j+j0
268   - value=sxpar(h1,'TTYPE'+strtrim(jj+1,2))
269   - sxaddpar,h1,'TTYPE'+strtrim(jj+1,2),value,'Intensity in '+use_filters[j]
270   - sxaddpar,h1,'FNAME'+strtrim(j+1,2),use_filters[j],'dustemwrap filter name'
271   -ENDFOR
272   -
273   -;stop
274   -
275   -mwrfits,seds_array,fits_file,h1,/create
276   -message,'Wrote '+fits_file,/continue
  224 +dustem_write_sed_grid,model,parameters_description,iv_min,iv_max,iv_Nvalues,plog $
  225 + ,parameter_values,seds,em_spectra,ext_spectra,fpd=fpd,fiv=fiv,filename=filename,use_polarisation=use_polarisation
277 226  
278   -toto=mrdfits(fits_file,1,header)
279 227  
  228 +;toto=mrdfits(fits_file,1,header)
280 229 ;hprint,header
281 230  
282 231 ;stop
... ...