OV-GSO-DC / dustem-wrapper_idl

Commit 2b6fd350afab8f4331a05a28bdf76dba251093d1

Authored by Jean-Philippe Bernard
1 parent f8c3f95f
Exists in master

added pre-computing of sed total value

Showing 2 changed files with 110 additions and 12 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_define_grid.pro
  Diff comments   View file @2b6fd35
1   -PRO dustem_define_grid,table_name,filters
  1 +PRO dustem_define_grid,table_name,filters,help=help
  2 +
  3 +;+
  4 +; NAME:
  5 +; dustem_define_grid
  6 +; PURPOSE:
  7 +; loads an SED Grid table into !dustem_grid
  8 +; CATEGORY:
  9 +; Dustem
  10 +; CALLING SEQUENCE:
  11 +; dustem_define_grid,table_name,filters
  12 +; INPUTS:
  13 +; table_name : fits file name containing the grid
  14 +; filters : name of dustemwrap filters to be included in the GRID
  15 +; OPTIONAL INPUT PARAMETERS:
  16 +; NONE
  17 +; OUTPUTS:
  18 +; None
  19 +; OPTIONAL OUTPUT PARAMETERS:
  20 +; None
  21 +; ACCEPTED KEY-WORDS:
  22 +; help = If set, print this help
  23 +; COMMON BLOCKS:
  24 +; None
  25 +; SIDE EFFECTS:
  26 +; !dustem_grid is updated
  27 +; RESTRICTIONS:
  28 +; The DustEM fortran code must be installed
  29 +; The DustEMWrap IDL code must be installed
  30 +; PROCEDURE:
  31 +; None
  32 +; EXAMPLES
  33 +; dustem_make_sed_table,'DBP90',['(*!dustem_params).G0','(*!dustem_params).grains(0).mdust_o_mh'], $
  34 +; [0.1,1.e-3],[100.,1.e-1],[10,3],filename='/tmp/IRAS_GO.fits',log=[1,0],/show_seds, $
  35 +; filters=[dustem_instru2filters('IRAC'),dustem_instru2filters('SPIRE')]
  36 +; dustem_define_grid,/tmp/IRAS_GO.fits',['IRAC3','SPIRE2']
  37 +; help,!dustem_grid
  38 +; MODIFICATION HISTORY:
  39 +; Written by J.-Ph. Bernard (2023)
  40 +; Evolution details on the DustEMWrap gitlab.
  41 +; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
  42 +;-
  43 +
  44 +IF keyword_set(help) THEN BEGIN
  45 + doc_library,'dustem_define_grid'
  46 + goto,the_end
  47 +ENDIF
2 48  
3 49 ;=== read the grid table
4 50 message,'reading grid data from '+table_name,/continue
... ... @@ -25,14 +71,32 @@ table_filters=table_filters[1:*]
25 71 Ntfilters=n_elements(table_filters)
26 72  
27 73 table_params=['']
  74 +table_pmins=[0.d0]
  75 +table_pmaxs=[0.d0]
  76 +table_pnvals=[0L]
  77 +table_plogs=[0]
28 78 count=1
29 79 i=0L
30 80 WHILE count NE 0 DO BEGIN
31 81 pname=sxpar(h,'PNAME'+strtrim(i+1,2),count=count)
32   - IF count NE 0 THEN table_params=[table_params,strtrim(pname,2)]
  82 + IF count NE 0 THEN BEGIN
  83 + table_params=[table_params,strtrim(pname,2)]
  84 + pmin=sxpar(h,'PMIN'+strtrim(i+1,2),count=count)
  85 + pmax=sxpar(h,'PMAX'+strtrim(i+1,2),count=count)
  86 + pnval=sxpar(h,'PNVAL'+strtrim(i+1,2),count=count)
  87 + plog=sxpar(h,'PLOG'+strtrim(i+1,2),count=count)
  88 + table_pmins=[table_pmins,pmin]
  89 + table_pmaxs=[table_pmaxs,pmax]
  90 + table_pnvals=[table_pnvals,pnval]
  91 + table_plogs=[table_plogs,plog]
  92 + ENDIF
33 93 i=i+1
34 94 ENDWHILE
35 95 table_params=table_params[1:*]
  96 +table_pmins=table_pmins[1:*]
  97 +table_pmaxs=table_pmaxs[1:*]
  98 +table_pnvals=table_pnvals[1:*]
  99 +table_plogs=table_plogs[1:*]
36 100 Ntparams=n_elements(table_params)
37 101  
38 102 ;=== select only requested filters
... ... @@ -47,21 +111,32 @@ ENDFOR
47 111  
48 112 ;stop
49 113  
50   -;sst=st
  114 +;totals=st.total
51 115 FOR i=0L,Ntfilters-1 DO BEGIN
52 116 IF keep_tfilters[i] EQ 0 THEN BEGIN
53 117 col_name=tags[i+Ntparams]
54 118 st=structure_remove_column(st,col_name)
55 119 ENDIF
56 120 ENDFOR
57   -
  121 +;==== recompute totals
  122 +Nlines=n_elements(st)
58 123 Nsed=n_elements(st)
  124 +;==== recompute totals (needed if some columns have been removed)
  125 +FOR i=0L,Nsed-1 DO BEGIN
  126 + this_total=0.d0
  127 + FOR j=0L,Nfilters-1 DO BEGIN
  128 + this_total=this_total+(st[i]).(j+Ntparams)
  129 + ENDFOR
  130 + st[i].total=this_total
  131 +ENDFOR
59 132  
60 133 ;stop
61 134  
62   -dustem_grid={Nsed:Nsed,Nfilters:Nfilters,Nparams:Ntparams,filters:filters,params:table_params,seds:st}
  135 +dustem_grid={Nsed:Nsed,Nfilters:Nfilters,Nparams:Ntparams,filters:filters,params:table_params,seds:st,pmin_values:table_pmins,pmax_values:table_pmaxs}
63 136  
64 137 defsysv,'!dustem_grid',dustem_grid
65 138 message,'Defined !dustem_grid',/continue
66 139  
  140 +the_end:
  141 +
67 142 END
68 143 \ No newline at end of file
... ...
src/idl/dustem_make_sed_table.pro
  Diff comments   View file @2b6fd35
1   -PRO dustem_make_sed_table,model,parameters_description,iv_min,iv_max,iv_Nvalues, $
2   - fpd=fpd,fiv=fiv,filename=filename,filters=filters,log=log,show_seds=show_seds,help=help
  1 +PRO dustem_make_sed_table,model, $
  2 + parameters_description, $
  3 + iv_min, $
  4 + iv_max, $
  5 + iv_Nvalues, $
  6 + fpd=fpd, $
  7 + fiv=fiv, $
  8 + filename=filename, $
  9 + filters=filters, $
  10 + log=log, $
  11 + show_seds=show_seds, $
  12 + help=help
3 13  
4 14 ;+
5 15 ; NAME:
... ... @@ -16,7 +26,7 @@ PRO dustem_make_sed_table,model,parameters_description,iv_min,iv_max,iv_Nvalues,
16 26 ; parameters_description : dustemwrap parameter description array
17 27 ; parvalues_min : minimum values for each parameter
18 28 ; parvalues_max : maximum values for each parameter
19   -; parvalues_Nvalues : number of parameter values
  29 +; parvalues_Nvalues : number of parameter values for each parameter
20 30 ; OPTIONAL INPUT PARAMETERS:
21 31 ; filters : name of dustemwrap filters to be included in the SED calculations (default = IRAS filters)
22 32 ; fpd : fixed parameter description
... ... @@ -40,8 +50,9 @@ PRO dustem_make_sed_table,model,parameters_description,iv_min,iv_max,iv_Nvalues,
40 50 ; PROCEDURE:
41 51 ; None
42 52 ; EXAMPLES
43   -; dustem_make_sed_table,'DBP90',['(*!dustem_params).G0'],[0.1],[100.],10,filename='/tmp/IRAS_GO.fits',/log,/show_seds, $
44   -; filters=[dustem_instru2filters('SPIRE'),dustem_instru2filters('SPIRE')]
  53 +; dustem_make_sed_table,'DBP90',['(*!dustem_params).G0','(*!dustem_params).grains(0).mdust_o_mh'], $
  54 +; [0.1,1.e-3],[100.,1.e-1],[10,3],filename='/tmp/IRAS_GO.fits',log=[1,0],/show_seds, $
  55 +; filters=[dustem_instru2filters('IRAC'),dustem_instru2filters('SPIRE')]
45 56 ; MODIFICATION HISTORY:
46 57 ; Written by J.-Ph. Bernard (2023)
47 58 ; Evolution details on the DustEMWrap gitlab.
... ... @@ -208,8 +219,10 @@ FOR i=0L,Nfilters-1 DO BEGIN
208 219 istr=strtrim(use_filters[i],2)
209 220 str=str+'I'+istr+':0.d0'
210 221 IF use_polarization EQ 1 THEN str=str+',Q'+istr+':0.d0,U'+istr+':0.d0'
211   - IF i NE Nfilters-1 THEN str=str+',' ELSE str=str+'}'
  222 + ;IF i NE Nfilters-1 THEN str=str+',' ELSE str=str+'}'
  223 + str=str+','
212 224 ENDFOR
  225 +str=str+'total:0.d0}'
213 226 toto=execute(str)
214 227  
215 228 ;stop
... ... @@ -220,10 +233,13 @@ FOR i=0L,Nc-1 DO BEGIN
220 233 FOR j=0L,Nparams-1 DO BEGIN
221 234 seds_array[i].(j)=(*parameter_values[i])[j]
222 235 ENDFOR
  236 + sum=0.d0
223 237 FOR j=0L,Nfilters-1 DO BEGIN
224 238 pos=Nparams+j*Noffset
225 239 seds_array[i].(pos)=((*seds[i]).stokesI)[j]
  240 + sum=sum+((*seds[i]).stokesI)[j] ;======== compute the sum of SEDs
226 241 ENDFOR
  242 + seds_array[i].total=sum
227 243 ENDFOR
228 244  
229 245 ;======== save SED tables
... ... @@ -245,7 +261,12 @@ FOR j=0L,Nparams-1 DO BEGIN
245 261 value=sxpar(h1,'TTYPE'+strtrim(j+1,2))
246 262 sxaddpar,h1,'TTYPE'+strtrim(j+1,2),value,'dustemwrap '+parameters_description[j]
247 263 sxaddpar,h1,'PNAME'+strtrim(j+1,2),parameters_description[j],'dustemwrap param name'
248   - last_pname='PNAME'+strtrim(j+1,2)
  264 + ;iv_min,iv_max,iv_Nvalues,Nc=Nc,log=log
  265 + sxaddpar,h1,'PMIN'+strtrim(j+1,2),iv_min[j],'minimum param value'
  266 + sxaddpar,h1,'PMAX'+strtrim(j+1,2),iv_max[j],'maximum param value'
  267 + sxaddpar,h1,'PNVAL'+strtrim(j+1,2),iv_Nvalues[j],'Number of param values'
  268 + sxaddpar,h1,'PLOG'+strtrim(j+1,2),log[j],'1 if log scale sampling'
  269 + last_pname='PLOG'+strtrim(j+1,2)
249 270 IF j EQ 0 THEN first_pname='PNAME'+strtrim(j+1,2)
250 271 ENDFOR
251 272 j0=j
... ... @@ -258,6 +279,8 @@ FOR j=0L,Nfilters-1 DO BEGIN
258 279 sxaddpar,h1,'FNAME'+strtrim(j+1,2),use_filters[j],'dustemwrap filter name'
259 280 ENDFOR
260 281  
  282 +;stop
  283 +
261 284 mwrfits,seds_array,fits_file,h1,/create
262 285 message,'Wrote '+fits_file,/continue
263 286  
... ...