OV-GSO-DC / dustem-wrapper_idl

Commit 7a1d28bfcef8f1d02584f754edc2ddd9c23a3138

Authored by Jean-Philippe Bernard
2 parents 7fd64468 7d2d66e7
Exists in master

Merge branch 'master' of https://gitlab.irap.omp.eu/OV-GSO-DC/dustem-wrapper_idl

Showing 9 changed files with 180 additions and 78 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_activate_plugins.pro
  Diff comments   View file @7a1d28b
... ... @@ -83,13 +83,18 @@ FOR i=0L,n_elements(p_min)-1 DO BEGIN
83 83 if strmid(ftn,0,13) eq 'dustem_plugin' then begin
84 84  
85 85  
86   - k=where(strmid(tag_names(*!dustem_scope),0,4) eq strmid(strupcase(strmid(ftn,14)),0,4),counte) ; Selecting a plugin through matching the string name of the plugin form the scope system variable with the one read from the parameter description vector
  86 + k=where(strmid(tag_names(*!dustem_scope),0,8) eq strmid(strupcase(strmid(ftn,14)),0,8),counte) ; Selecting a plugin through matching the string name of the plugin form the scope system variable with the one read from the parameter description vector
87 87  
88 88 ;Dry run of the plugins to obtain their scopes and run them accordingly (an advanced user might want to add their own lines here)
89 89 str='toto='+ftn+'(scope=scope)'
90 90 toto=execute(str)
  91 +
91 92 (*!dustem_scope).(k)=ptr_new(scope)
92 93  
  94 +
  95 + ;print,
  96 +
  97 +
93 98  
94 99 str='((*!dustem_plugin).('+strtrim(k,2)+'))=ptr_new('+ftn+'(key=index,val=value)'+')' & str=str(0)
95 100 toto=execute(str) & IF !dustem_verbose NE 0 THEN message,strupcase(strmid(ftn,7)),/info
... ...
src/idl/dustem_compute_polsed.pro
  Diff comments   View file @7a1d28b
... ... @@ -58,6 +58,16 @@ ENDFOR
58 58 ;-----------------------------------------
59 59  
60 60  
  61 +;MODIFYING PLUGIN----------------
  62 +FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
  63 + IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'REPLACE_POLSED') THEN BEGIN
  64 + P=(*(*!dustem_plugin).(i))[*,3]
  65 + P+=(*(*!dustem_plugin).(i))[*,3] ;3 is for the polarized portion of the plugin because this is polsed
  66 + ENDIF
  67 +ENDFOR
  68 +;-----------------------------------------
  69 +
  70 +
61 71 dustem_polsed = (*!dustem_data.sed).values * 0.
62 72  
63 73  
... ...
src/idl/dustem_compute_stokes.pro
  Diff comments   View file @7a1d28b
... ... @@ -60,7 +60,7 @@ frac=P/stI
60 60 tes=where(finite(frac) eq 0)
61 61 frac(tes)=0.
62 62  
63   -
  63 +;if Q and U have already been computed elsewhere
64 64 FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
65 65  
66 66 IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'MODIFY_DUST_POLAR_ANGLE') THEN BEGIN
... ... @@ -70,11 +70,27 @@ FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
70 70  
71 71 ENDIF
72 72 ENDFOR
  73 +
  74 +
  75 +FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
  76 +
  77 + IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'REPLACE_QSED') THEN BEGIN
  78 +
  79 + Q_spec=(*(*!dustem_plugin).(i))[*,1]
  80 + ENDIF
  81 +
  82 + IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'REPLACE_USED') THEN BEGIN
  83 + U_spec=(*(*!dustem_plugin).(i))[*,2]
  84 + ENDIF
  85 +
  86 +ENDFOR
  87 +
  88 +
73 89  
74 90 IF ~isa(Q_spec) && ~isa(U_spec) THEN polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(0.,Nwaves)
75 91  
76 92 FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
77   -
  93 +
78 94 IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_QSED') THEN BEGIN
79 95 Q_spec+=(*(*!dustem_plugin).(i))[*,1]
80 96 ENDIF
... ... @@ -135,14 +151,6 @@ IF count_spec NE 0 THEN BEGIN
135 151 ENDIF
136 152 out_st=stp
137 153  
138   -;out=fltarr(Nwaves)
139   -
140   -; IF keyword_set(dustem_qsed) THEN out=dustem_qsed
141   -; IF keyword_set(dustem_used) THEN out=dustem_used
142   -; IF keyword_set(Q_spec) THEN out=Q_spec
143   -; IF keyword_set(U_spec) THEN out=U_spec
144   -
145   -
146 154  
147 155 ;clean pointers
148 156 heap_gc
... ...
src/idl/dustem_fit_sed_polsed_readme.pro
  Diff comments   View file @7a1d28b
... ... @@ -93,14 +93,18 @@ st=dustem_set_data(sed=sed,rchi2_weight=rchi2_weight,f_HI=f_HI)
93 93 pd = [ $
94 94 '(*!dustem_params).G0', $ ;G0
95 95 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
96   - '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
97   - '(*!dustem_params).grains(2).mdust_o_mh',$ ;amCBEx
  96 + '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH0 mass fraction
  97 + '(*!dustem_params).grains(2).mdust_o_mh',$ ;PAH1 mass fraction
  98 + '(*!dustem_params).grains(3).mdust_o_mh',$ ;amCBEx
98 99 'dustem_plugin_synchrotron_2', $ ;Synchrotron amplitude
99 100 'dustem_plugin_synchrotron_4', $ ;Synchrotron polarization angle
100   - 'dustem_plugin_modify_dust_polar_angle_1' $ ;Polarization angle applied to the Fortran dust model via a core plugin
  101 + 'dustem_plugin_modify_dust_polar_angle_1', $ ;Polarization angle applied to the Fortran dust model via a core plugin
  102 + 'dustem_plugin_modify_spinning_sed_1', $
  103 + 'dustem_plugin_modify_spinning_sed_2' $
101 104 ]
102 105  
103   -p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04,5e-3,20.,28.]
  106 +p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04,1.6500E-04,5e-3,20.,68.,1e-3,20]
  107 +;p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04,1.6500E-04,5e-3,20.,0.2,20]
104 108  
105 109 ;=== set initial value of parameters you want to fit
106 110 ;iv=p_truth ;exact solution
... ... @@ -108,8 +112,8 @@ p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04,5e-3,20.,28.]
108 112  
109 113 ;UNCOMMENT THIS BLOCK AFTER THE TESTING
110 114  
111   - iv=p_truth+[0.,8.e-4,8.e-4,8.e-4,8e-3,-8.,-10.] ;shifted from solution
112   -
  115 +iv=p_truth+[0.,8.e-4,8.e-4,8.e-4,8.e-4,8e-3,-8.,-10.,1.2,10] ;shifted from solution
  116 +;iv=p_truth+[0.,8.e-4,8.e-4,8.e-4,8.e-4,8e-3,-8.,1.2,10] ;shifted from solution
113 117  
114 118 ;TESTING THE ISRF PLUGIN
115 119  
... ... @@ -120,6 +124,9 @@ llims=replicate(0,Npar)
120 124  
121 125 ;== SET THE FITTED PARAMETERS (now done in dustem_init_parinfo)
122 126 dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
  127 +(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg-spin'
  128 +(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg-spin'
  129 +
123 130  
124 131 dustem_init_plugins,pd
125 132  
... ...
src/idl/dustem_fit_sed_readme.pro
  Diff comments   View file @7a1d28b
... ... @@ -76,7 +76,7 @@ IF keyword_set(png) THEN BEGIN
76 76 force_mkdir,dir_png
77 77 ENDIF
78 78  
79   -polarization=0 ;default is no polarization in models
  79 +polarization=0. ;default is no polarization in models
80 80  
81 81 ;COMMENT: for some reson polarization=1 doesn't work for ALL of the models except for the two last A and B models
82 82  
... ... @@ -90,7 +90,7 @@ CASE use_model OF
90 90 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
91 91 '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
92 92 '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx
93   - 'dustem_create_continuum_2'] ;Intensity of NIR continuum
  93 + 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum
94 94 iv = [1.0, 4.3e-4, 4.7e-4,6.4e-3,0.001]
95 95 Npar=n_elements(pd)
96 96 ulimed=replicate(0,Npar)
... ... @@ -105,7 +105,7 @@ CASE use_model OF
105 105 '(*!dustem_params).grains(2).mdust_o_mh', $ ;Gra
106 106 '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra
107 107 '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil
108   - 'dustem_create_continuum_2'] ;Intensity of NIR continuum
  108 + 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum
109 109 iv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001]
110 110 Npar=n_elements(pd)
111 111 ulimed=replicate(0,Npar)
... ... @@ -119,9 +119,9 @@ CASE use_model OF
119 119 '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
120 120 '(*!dustem_params).grains(2).mdust_o_mh', $ ;Gra
121 121 '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra
122   - '(*!dustem_params).grains(4).mdust_o_mh' $ ;aSil
123   - ;'dustem_create_continuum_2'] ;Intensity of NIR continuum
124   - ]
  122 + '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil
  123 + 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum
  124 + ;]
125 125  
126 126 ; pd = [ $
127 127 ; '(*!dustem_params).gas.G0', $ ;G0
... ... @@ -137,7 +137,7 @@ CASE use_model OF
137 137 ; 'dustem_create_stellar_population_B45'$ ; number of stars of the B4 stellar population
138 138 ; ]
139 139  
140   - iv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3];,0.001];,10,1.,10.,1.]
  140 + iv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001];,10,1.,10.,1.]
141 141  
142 142  
143 143 Npar=n_elements(pd)
... ... @@ -155,8 +155,6 @@ CASE use_model OF
155 155 '(*!dustem_params).gas.G0', $ ;G0
156 156 'dustem_plugin_continuum_2', $ ;intensity of NIR continuum
157 157 'dustem_plugin_synchrotron_2',$
158   - ;'!dustem_isrf_star_add[0].distance', $ ;distance of first stellar contribution to ISRF
159   - ;'!dustem_isrf_star_add[0].amplitude', $ ;amplitude of first stellar contribution to ISRF
160 158 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
161 159 '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
162 160 '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx
... ... @@ -167,9 +165,7 @@ CASE use_model OF
167 165 iv = [ $
168 166 1.0, $
169 167 0.002, $ ;intensity of NIR continuum
170   - 0.01,$
171   - ;0.003, $ ;distance of first stellar contribution to ISRF
172   - ;1.e-3, $ ;amplitude of first stellar contribution to ISRF
  168 + 0.01,$
173 169 7.8e-4, $ ;mass fraction of PAH0
174 170 7.8e-4, $ ;mass fraction of PAH1
175 171 1.65e-4, $ ;mass fraction of amCBEx
... ... @@ -305,7 +301,7 @@ CASE use_model OF
305 301 '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx
306 302 '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx
307 303 '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil
308   - 'dustem_create_continuum_2'] ;Intensity of NIR continuum
  304 + 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum
309 305 iv = [1.0, 7.8e-4, 7.8e-4,1.65e-4,1.45e-3,7.8e-3,0.001]
310 306 Npar=n_elements(pd)
311 307 ulimed=replicate(0,Npar)
... ... @@ -322,7 +318,7 @@ CASE use_model OF
322 318 '(*!dustem_params).grains(2).mdust_o_mh', $ ;Gra
323 319 '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra
324 320 '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil
325   - 'dustem_create_continuum_2'] ;Intensity of NIR continuum
  321 + 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum
326 322 iv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001]
327 323 Npar=n_elements(pd)
328 324 ulimed=replicate(0,Npar)
... ... @@ -338,7 +334,7 @@ CASE use_model OF
338 334 '(*!dustem_params).grains(2).mdust_o_mh', $ ;Gra
339 335 '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra
340 336 '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil
341   - 'dustem_create_continuum_2'] ;Intensity of NIR continuum
  337 + 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum
342 338 iv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001]
343 339 Npar=n_elements(pd)
344 340 ulimed=replicate(0,Npar)
... ...
src/idl/dustem_plugin_modify_dust_polar_angle.pro
  Diff comments   View file @7a1d28b
... ... @@ -58,9 +58,6 @@ IF keyword_set(key) THEN BEGIN
58 58 a=where(key EQ 1,count1)
59 59 IF count1 NE 0 THEN psi=val[a[0]] ; setting psi from pd
60 60  
61   - print, 'psi in modify_dust_polar_angle is:'
62   - print, psi
63   -
64 61 ENDIF
65 62  
66 63  
... ... @@ -77,13 +74,6 @@ frac(tes)=0.
77 74  
78 75 polar_ippsi2iqu,I,Q,U,frac,replicate(psi,Nwaves)
79 76  
80   -print, 'U in modify_dust_polar_angle is:'
81   -print, U
82   -print, 'Q in modify_dust_polar_angle is:'
83   -print, Q
84   -
85   -
86   -
87 77 out=fltarr(Nwaves,2)
88 78  
89 79 out[*,0]=Q
... ...
src/idl/dustem_plugin_modify_spinning_sed.pro 0 → 100644
  Diff comments   View file @7a1d28b
... ... @@ -0,0 +1,86 @@
  1 +FUNCTION dustem_plugin_modify_spinning_sed, key=key, val=val, scope=scope, help=help
  2 +
  3 +;+
  4 +; NAME:
  5 +; dustem_create_stext
  6 +; PURPOSE:
  7 +; MODIFIES THE SPINNING DUST SED BY INTRODUCING A POLARIZATION FRACTION p AND A POLARIZATION ANGLE psi
  8 +; CATEGORY:
  9 +; DUSTEM Wrapper
  10 +; CALLING SEQUENCE:
  11 +; dustem_plugin_modify_spinning_sed(st,key=key,val=val)
  12 +; INPUTS:
  13 +; st (st = dustem_run(p_dim))
  14 +; OPTIONAL INPUT PARAMETERS:
  15 +; key = input parameter number
  16 +; val = input parameter value
  17 +; OUTPUTS:
  18 +; out = array containing the extinction stokes parameters
  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 +; None
  27 +; RESTRICTIONS:
  28 +; The dustem fortran code must be installed
  29 +; The dustem idl wrapper must be installed
  30 +; PROCEDURE:
  31 +; This is a dustem plugin
  32 +;-
  33 +
  34 +IF keyword_set(help) THEN BEGIN
  35 + doc_library,'dustem_plugin_modify_spinning_sed'
  36 + goto,the_end
  37 +ENDIF
  38 +
  39 +
  40 +scope='REPLACE_POLSED+REPLACE_QSED+REPLACE_USED'
  41 +
  42 +IF ptr_valid((*!dustem_fit).CURRENT_PARAM_VALUES) THEN BEGIN
  43 + p_dii=(*(*!dustem_fit).CURRENT_PARAM_VALUES)
  44 + st=dustem_run(p_dii)
  45 +
  46 +ENDIF ELSE BEGIN
  47 + p_dii=(*(*!dustem_fit).PARAM_INIT_VALUES)
  48 + st=dustem_run(p_di)
  49 +
  50 + ENDELSE
  51 +
  52 +
  53 +psi=0.
  54 +smallp=0. ;supposing the spinning dust is
  55 +IF keyword_set(key) THEN BEGIN
  56 + ind1=where(key EQ 1,count1)
  57 + ind2=where(key EQ 2,count2)
  58 + IF count1 NE 0 THEN smallp=val[ind1[0]] ; setting psi from pd
  59 + IF count2 NE 0 THEN psi=val[ind2[0]] ; setting psi from pd
  60 +ENDIF
  61 +
  62 +
  63 +fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/(st.polsed).wav)*1.e20/1.e7
  64 +
  65 +
  66 +SpinningI=(st.sed.em_grain_1)+(st.sed.em_grain_2)*fact
  67 +SpinningP=SpinningI*smallp
  68 +
  69 +Nwaves=(size(SpinningI))[1]
  70 +
  71 +polar_ippsi2iqu,SpinningI,SpinningQ,SpinningU,replicate(smallp,Nwaves),replicate(psi,Nwaves)
  72 +
  73 +out=fltarr(Nwaves,4)
  74 +
  75 +out[*,0]=SpinningI
  76 +out[*,1]=SpinningQ
  77 +out[*,2]=SpinningU
  78 +out[*,3]=SpinningP
  79 +
  80 +return, out
  81 +
  82 +the_end:
  83 +
  84 +
  85 +end
  86 +
... ...
src/idl/dustem_sed_plot.pro
  Diff comments   View file @7a1d28b
... ... @@ -16,7 +16,7 @@ IF not keyword_set(function_name) THEN BEGIN
16 16  
17 17 ;I do not understand the reason begind the stop below - either way polarization is handled by other procedures
18 18 IF keyword_set(pol) THEN BEGIN
19   - stop
  19 + ;stop
20 20 ;dustem_polsed=dustem_compute_polsed(p_min,st=pst,cont=cont,freefree=freefree,synchrotron=synchrotron,out_st=out_st)
21 21 dustem_plot_polar,st,ps=ps,_Extra=extra,help=help,UV=UV,SED=SED,Pfrac=Pfrac,print_ratio=print_ratio,win=win,cont=cont,xr=xr,yr=yr,donotclose=donotclose,aligned=aligned,noerrbars=noerrbars,multi=multi,almabands=almabands,nsmooth=nsmooth
22 22 ENDIF
... ...
src/idl/dustem_write_gas.pro
  Diff comments   View file @7a1d28b
... ... @@ -49,38 +49,38 @@ free_lun,unit
49 49 END
50 50  
51 51  
52   -;==== LEFT OVERS
53   -
54   -PRO dustem_write_gas,dir,st
55   -
56   -
57   -;Writes ALL the GAS_xxx.DAT files
58   -;stop
59   -
60   -comments=st.comments
61   -Ncomments=n_elements(comments)
62   -
63   -frmt='(A)'
64   -
65   -;Get filename
66   -ffile=st.file
67   -fv=str_sep(ffile,'/')
68   -file=dir+fv(n_elements(fv)-1)
69   -;Open file
70   -message,'Will write GAS file '+file,/continue
71   -openw,unit,file,/get_lun
72   -;Print comments
73   -IF Ncomments NE 0 THEN BEGIN
74   - FOR ii=0,Ncomments-1 DO BEGIN
75   - printf,unit,comments[ii]
76   - ENDFOR
77   -ENDIF
78   -;print data
79   -FOR k=0L,n_elements(st.fgas)-1 DO BEGIN
80   - printf,unit,(st.fgas)[k],format=frmt
81   -ENDFOR
82   -close,unit
83   -free_lun,unit
84   -
85   -END
86   -
  52 +; ;==== LEFT OVERS
  53 +;
  54 +; PRO dustem_write_gas,dir,st
  55 +;
  56 +;
  57 +; ;Writes ALL the GAS_xxx.DAT files
  58 +; ;stop
  59 +;
  60 +; comments=st.comments
  61 +; Ncomments=n_elements(comments)
  62 +;
  63 +; frmt='(A)'
  64 +;
  65 +; ;Get filename
  66 +; ffile=st.file
  67 +; fv=str_sep(ffile,'/')
  68 +; file=dir+fv(n_elements(fv)-1)
  69 +; ;Open file
  70 +; message,'Will write GAS file '+file,/continue
  71 +; openw,unit,file,/get_lun
  72 +; ;Print comments
  73 +; IF Ncomments NE 0 THEN BEGIN
  74 +; FOR ii=0,Ncomments-1 DO BEGIN
  75 +; printf,unit,comments[ii]
  76 +; ENDFOR
  77 +; ENDIF
  78 +; ;print data
  79 +; FOR k=0L,n_elements(st.fgas)-1 DO BEGIN
  80 +; printf,unit,(st.fgas)[k],format=frmt
  81 +; ENDFOR
  82 +; close,unit
  83 +; free_lun,unit
  84 +;
  85 +; END
  86 +;
... ...