OV-GSO-DC / dustem-wrapper_idl

Commit 27bec12ae15b1b21e9fffc4d44eb8e7b7763f986

Authored by Annie Hughes
1 parent 8929b443
Exists in master

tidied up help and cleaned code to make it clearer

Showing 1 changed file with 75 additions and 179 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_stellarpopisrf_example.pro
  Diff comments   View file @27bec12
... ... @@ -4,6 +4,7 @@ PRO dustem_stellarpopisrf_example,model=model $
4 4 ,postscript=postscript $
5 5 ,fits_save_and_restore=fits_save_and_restore $
6 6 ,wait=wait $
  7 + ,noobj=noobj $
7 8 ,verbose=verbose $
8 9 ,help=help
9 10  
... ... @@ -13,18 +14,17 @@ PRO dustem_stellarpopisrf_example,model=model $
13 14 ;
14 15 ; PURPOSE:
15 16 ; This routine is an example of how to fit an observational SED
16   -; (StokesI only) with DustEM and DustEMWrap. The objective is to
17   -; illustrate how to use the stellar population plugin and not to do science -- the fit
18   -; obtained by running this example is likely to be poor.
19   -;
20   -; For this example, the code uses the SED in the file example_SED_2.xcat,
21   -; which is distributed in the Data/EXAMPLE_OBSDATA/ directory
  17 +; (StokesI only) with DustEM and DustEMWrap, and an ISRF that is due
  18 +; to a user-defined population of nearby main sequence stars.
22 19 ;
23   -; The example SED has Stokes I photometric data points from
24   -; IRAC, MIPS and IRAS. Examples illustrating running DustEMWrap to
25   -; fit spectral data, polarisation data and extinction data
26   -; are provided in other _example routines in the src/idl/
27   -; directory. See the DustEMWrap User Guide for more information.
  20 +; DustEMWrap reads information about the stellar spectral types (effective
  21 +; temperature, radius) from the EEM_dwarf_UBVIJHK_colors_Teff.txt file
  22 +; that is located in the Data/STELLARPOPS/ directory. This data file was
  23 +; authored by Prof. Erik Mamajek (see the file for more details).
  24 +;
  25 +; For this example, we generate an SED using an input model and then
  26 +; launch the fit with a starting guess that has been shifted away from
  27 +; the true parameter values.
28 28 ;
29 29 ; CATEGORY:
30 30 ; DustEMWrap, Distributed, High-Level, User Example
... ... @@ -66,6 +66,7 @@ PRO dustem_stellarpopisrf_example,model=model $
66 66 ; wait = if set, wait this many seconds between each step of
67 67 ; the code (for illustration purposes)
68 68 ; verbose = if set, subroutines will run in verbose mode
  69 +; noobj = if set, runs with no object graphics
69 70 ;
70 71 ; COMMON BLOCKS:
71 72 ; None
... ... @@ -102,10 +103,16 @@ ENDIF ELSE BEGIN
102 103 use_model='J13' ;Default is last dustem model
103 104 ENDELSE
104 105  
  106 +known_mdls=['MC10','DBP90','DL01','WD01_RV5P5B','DL07','J13','G17_MODELA','G17_MODELB','G17_MODELC','G17_MODELD']
  107 +test_model = where(known_mdls eq use_model,ct)
  108 +if ct eq 0 then begin
  109 + message,'ISM dust model '+use_model+' unknown',/continue
  110 + message,'Known models are MC10,DBP90,DL01,WD01_RV5P5B,DL07,J13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue
  111 + stop
  112 +end
105 113  
106 114  
107   -
108   -
  115 +use_polarization=0 ; initialize Dustemwrap in no polarization mode
109 116 use_verbose=0
110 117 if keyword_set(verbose) then use_verbose=1
111 118 use_Nitermax=5 ; maximum number of iterations for the fit
... ... @@ -117,183 +124,72 @@ dustem_define_la_common
117 124 ;=== Refer to the DustEM and DustEMWrap userguides for an explanation
118 125 ; of the different grain types
119 126  
120   -CASE use_model OF
121   - 'DBP90':BEGIN
122   - pd = [ $
123   - '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
124   - '(*!dustem_params).grains(1).mdust_o_mh',$ ;VSG mass fraction
125   - '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx
126   - 'dustem_plugin_stellar_population_O7V3',$ ;distance to stellar population
127   - 'dustem_plugin_stellar_population_O5V4'] ;number of stars
128   - rv = [4.3e-4, 4.7e-4,6.4e-3,1.5,100.] ;rv because we're generating fake data and try to retrieve the original values
129   -
130   - iv = rv+[5.00E-4,5.00E-4,5.00E-4,-1.,-30] ;###setting the initial parameter vector
131   -
132   - Npar=n_elements(pd)
133   - ulimed=replicate(0,Npar)
134   - llimed=replicate(1,Npar)
135   - llims=replicate(0.,Npar)
136   -
137   -
138   - fpd=[ $
139   - 'dustem_plugin_stellar_population_O7V4',$ ;number of stars of stellar population' , $ ;multiplicative factor to total ISRF
140   - 'dustem_plugin_stellar_population_O5V3'$ ; stellar distance
141   - ]
142   - ;initial parameter values for fixed parameters
143   - fiv=[ $
144   - 2. , $ ;number of stars of pop O7
145   - 45$ ;dustance of pop O5
146   - ]
147   -
148   - END
149   - 'DL01':BEGIN
150   - pd = [ $
151   - '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
152   - '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
153   - '(*!dustem_params).grains(2).mdust_o_mh', $ ;Gra
154   - '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra
155   - '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil
156   - 'dustem_plugin_stellar_population_O5V4'] ;number of stars of pop O5
157   -
158   - rv = [5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,50]
159   - iv = rv+[5.e-4, 5.e-4,5.e-5,5.e-3,5.e-3,100]
160   -
161   - fpd=[ $
162   - 'dustem_plugin_stellar_population_O5V3'$ ; stellar distance
163   - ]
164   - ;initial parameter values for fixed parameters
165   - fiv=[ $
166   - 45$ ;fixed stellar distance
167   - ]
  127 +;; ; example parameter initialisation for DBP90
  128 +;; pd = [ $
  129 +;; '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
  130 +;; '(*!dustem_params).grains(1).mdust_o_mh',$ ;VSG mass fraction
  131 +;; '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx
  132 +;; 'dustem_plugin_stellar_population_O7V3',$ ;distance to O7V stars
  133 +;; 'dustem_plugin_stellar_population_O5V4'] ;number of O5V stars
168 134  
  135 +;; rv = [4.3e-4, 4.7e-4,6.4e-3,1.5,100.] ;rv because we're generating fake data and try to retrieve the original values
  136 +;; iv = rv+[5.00E-4,5.00E-4,5.00E-4,-1.,-30] ;###setting the initial parameter vector
169 137  
170   - Npar=n_elements(pd)
171   - ulimed=replicate(0,Npar)
172   - llimed=replicate(1,Npar)
173   - llims=replicate(0.,Npar)
174   - END
175   - 'WD01_RV5p5B':BEGIN
176   -;; ***COMMENT AH***
177   -;; we need to implement this, or remove
178   - message, 'WD01 model not yet implemented in DustEMWrap',/info
179   -;; ***END COMMENT AH***
180   - END
181   - 'DL07':BEGIN
182   - pd = [ $
183   - '(*!dustem_params).grains(0).mdust_o_mh', $ ;PAH0 mass fraction
184   - '(*!dustem_params).grains(1).mdust_o_mh', $ ;PAH1 mass fraction
185   - '(*!dustem_params).grains(2).mdust_o_mh', $ ;Gra
186   - '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra
187   - '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil
188   - 'dustem_plugin_stellar_population_O5V3' $ ;distance to stellar pop O5V
189   - ]
190   -
191   - rv = [5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,5.];,10,1.,10.,1.]
192   - iv = rv+ [5.e-4, 5.e-4,5.e-4,5.e-3,5.e-3,25.]
193   -
194   - Npar=n_elements(pd)
195   - ulimed=replicate(0,Npar)
196   - llimed=replicate(1,Npar)
197   - llims=replicate(0.,Npar)
198   -
199   - fpd=[ $
200   - 'dustem_plugin_stellar_population_O5V4'$ ; stellar number of pop O5V
201   - ]
202   - ;initial parameter values for fixed parameters
203   - fiv=[ $
204   - 30$ ;fixed stellar population number
205   - ]
206   -
207   - END
208   - 'MC10':BEGIN
209   -
210   - pd = [ $
211   - 'dustem_plugin_stellar_population_O5V4', $ ;stellar number of pop O5V
212   - '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
213   - '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
214   - '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx
215   - '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx
216   - '(*!dustem_params).grains(4).mdust_o_mh' $ ;aSilx
217   - ]
218   - ;initial parameter values for parameters to be fitted
219   - rv = [ $
220   - 50, $
221   - 7.8e-4, $ ;mass fraction of PAH0
222   - 7.8e-4, $ ;mass fraction of PAH1
223   - 1.65e-4, $ ;mass fraction of amCBEx
224   - 1.45e-3, $ ;mass fraction of amCBEx
225   - 7.8e-3 $ ;mass fraction of aSilx
226   - ]
227   -
228   - iv = rv + [-30,5.E-4,5.E-4,5.E-4,5.E-3,5.E-3]
229   -
230   - fpd = ['dustem_plugin_stellar_population_O5V3'] ;distance to stellar population O5V
231   - fiv = [25]
  138 +;; Npar=n_elements(pd)
  139 +;; ulimed=replicate(0,Npar)
  140 +;; llimed=replicate(1,Npar)
  141 +;; llims=replicate(0.,Npar)
232 142  
233   -
234   - Npar=n_elements(pd)
  143 +;; ;fixed parameters
  144 +;; fpd=[ $
  145 +;; 'dustem_plugin_stellar_population_O7V4',$ ;number of O7V stars
  146 +;; 'dustem_plugin_stellar_population_O5V3'$ ; distance to O5V stars
  147 +;; ]
  148 +;; ;initial parameter values for fixed parameters
  149 +;; fiv=[2., 45.]
  150 +
  151 +;'J13':BEGIN
  152 +pd = [ $
  153 + '(*!dustem_params).grains(1).mdust_o_mh',$ ;CM20 mass fraction
  154 + '(*!dustem_params).grains(0).mdust_o_mh',$ ;CM20 mass fraction
  155 + '(*!dustem_params).grains(2).mdust_o_mh', $ ;aPyM5
  156 + '(*!dustem_params).grains(3).mdust_o_mh', $ ;aPyM5
  157 + 'dustem_plugin_stellar_population_O5V4'] ; number of pop O5V stars
  158 +
  159 +rv = [7.8e-4, 7.8e-4,1.65e-4,1.45e-3,35] ; real values
  160 +iv = rv + [2.e-4,-1.e-4,3.e-5,5.e-4,-5] ; starting guess values
235 161  
236   - ulimed=replicate(0,Npar)
237   - llimed=replicate(1,Npar)
238   - llims=replicate(0.,Npar)
239   -
240   - END
241   - 'J13':BEGIN
242   - pd = [ $
243   - '(*!dustem_params).grains(1).mdust_o_mh',$ ;CM20 mass fraction
244   - '(*!dustem_params).grains(0).mdust_o_mh',$ ;CM20 mass fraction
245   - '(*!dustem_params).grains(2).mdust_o_mh', $ ;aPyM5
246   - '(*!dustem_params).grains(3).mdust_o_mh', $ ;aPyM5
247   - 'dustem_plugin_stellar_population_O5V4'] ;stellar number of pop O5V
248   -
249   - rv = [7.8e-4, 7.8e-4,1.65e-4,1.45e-3,35]
250   - iv = rv+ [5.e-4,5.e-4,5.e-4,5.e-3,35]
251   -
252   - Npar=n_elements(pd)
253   - ulimed=replicate(0,Npar)
254   - llimed=replicate(1,Npar)
255   - llims=replicate(0.,Npar)
256   -
257   - fpd = ['dustem_plugin_stellar_population_O5V3'] ;distance to stellar population O5V
258   - fiv = [5]
  162 +Npar=n_elements(pd)
  163 +ulimed=replicate(0,Npar)
  164 +llimed=replicate(1,Npar)
  165 +llims=replicate(0.,Npar)
  166 +
  167 +;fixed parameters and their values
  168 +fpd = ['dustem_plugin_stellar_population_O5V3' $ ;distance to stellar population O5V
  169 + ,'(*!dustem_params).G0'] ; Mathis field
  170 +fiv = [5.,1.e-12]
259 171  
260   - END
261   - 'ELSE':BEGIN
262   - message,'model '+model+' unknown',/continue
263   - message,'Known models are MC10,DBP90,DL01,DL07,J13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue
264   - stop
265   - END
266   -ENDCASE
267   -
268   -
269   -if keyword_set(wait) then begin
270   - message,'Finished setting dust model and plug-in parameters: '+use_model,/info
271   - wait,wait
272   -end
273   -
274   -
275   -dustem_init,model=use_model
  172 +dustem_init,model=use_model,polarization=use_polarization
276 173  
277 174 !dustem_nocatch=1
278 175 !dustem_verbose=1
279 176 ;!dustem_show_plot=1
280   -!EXCEPT=2 ;so I can locate the plotting error...
281   -;!dustem_dim=1 ; this option is to dim the stellar population ISRF with the current Dustem extinction
282   -;We're fitting total optical depths so the dimmed-ISRF scenario is a reasonable assumption.
283   -;Because of this the final parameter values are slightly lower than the initial (real) values.
284   -;But this is because the first run was not extinct with the current dustem extinction
285   -;since it comes prior to the dustem run. Meaning the first ISRF has to be extinct to see if the paramters are retrieved.
286   -
  177 +;; !EXCEPT=2 ;so I can locate the plotting error...
  178 +;; ;!dustem_dim=1 ; this option is to dim the stellar population ISRF with the current Dustem extinction
  179 +;; ;We're fitting total optical depths so the dimmed-ISRF scenario is a reasonable assumption.
  180 +;; ;Because of this the final parameter values are slightly lower than the initial (real) values.
  181 +;; ;But this is because the first run was not extinct with the current dustem extinction
  182 +;; ;since it comes prior to the dustem run. Meaning the first ISRF has to be extinct to see if the paramters are retrieved.
  183 +IF keyword_set(noobj) THEN !dustem_noobj=1
287 184  
288 185 !dustem_isrf_file=ptr_new(!dustem_soft_dir+'data/ISRF_MATHIS.DAT')
289 186  
  187 +;=== READ EXAMPLE DATA:
290 188  
291   -;=== READ EXAMPLE DATA: EXTINCTION
292   -
293   -;NB:HERE WE ARE NOT FITTING THE DATA. WE'RE JUST USING THE EXISTING FILTERS.
294   -
  189 +;NB: HERE WE ARE READING AN SED FILE JUST TO SET-UP THE SED STRUCTURE
  190 +;AND FILTERS. WE REPLACE THE STOKES I and STOKES I UNCERTAINTIES LATER
295 191 dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
296   -file=dir+'example_SED_2.xcat'
  192 +file=dir+'example_SED_1.xcat'
297 193 IF keyword_set(sed_file) THEN file=sed_file
298 194 sed=read_xcat(file,/silent)
299 195  
... ... @@ -303,7 +199,7 @@ if keyword_set(wait) then begin
303 199 end
304 200  
305 201  
306   -;=== GENERATE EXAMPLE DATA: EXTINCTION
  202 +;=== GENERATE EXAMPLE DATA:
307 203  
308 204 ;=== initializing IQU and associated errors to avoid problems when checking SED in dustem_set_data.pro
309 205 ;sed.StokesI=1.
... ... @@ -331,7 +227,7 @@ dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,u
331 227 sed.StokesI = dustem_compute_sed(rv,st=st)
332 228  
333 229 ;== SET THE OBSERVATIONAL STRUCTURE
334   -dustem_set_data, sed,sed
  230 +dustem_set_data, sed, sed
335 231  
336 232  
337 233 if keyword_set(wait) then begin
... ... @@ -412,7 +308,7 @@ IF keyword_set(fits_save_and_restore) THEN BEGIN
412 308 ENDIF
413 309 ENDIF
414 310  
415   -message,'Finished dustem_fit_intensity_example',/info
  311 +message,'Finished dustem_stellarpopisrf_example',/info
416 312  
417 313  
418 314 the_end:
... ...