OV-GSO-DC / dustem-wrapper_idl

Commit 78db85de6e95a68dac4ec6b5545b1a2b576f4841

Authored by Jean-Philippe Bernard
1 parent 4ddd208e
Exists in master

modified for different flavors of DL07 dealing with APH1 and PAH0

Showing 4 changed files with 200 additions and 24 deletions   Show diff stats
Inline Side-by-side
LabTools/IRAP/JPB/make_sed_phangs_tables.pro
  Diff comments   View file @78db85d
... ... @@ -405,6 +405,146 @@ CASE use_grid_type OF
405 405 dustem_add_linear_params2grid,table_name,pd,iv_min,iv_max,iv_Nvalues,out_filename=output_table_name,plog=plog,show_seds=show_seds
406 406 goto,the_end
407 407 END
  408 + 222: BEGIN ;This is same as 221, but using the standard DL07 GRAIN.DAT, ie, Ypah+ set by Fortran and PAH keywords set to mix-logn
  409 + ;Note: This does not really work, because PAH+ contribution does not get multiplied linearly by abundance ...
  410 + !quiet=1
  411 + model='DL07' ;This is to use the DBP90 dust model
  412 + ;======= This is to produce a grid with the above model for G0 and PAH abundance and ISRF from the Muse data
  413 + ;======= Note that g0 is fixed parameters is actually not used, replaced by dustem_plugin_phangs_class_isrf_2
  414 + pd = [ $
  415 + 'dustem_plugin_phangs_class_isrf_2' $ ;G0 factor for Phangs ISRF classes
  416 + ]
  417 + iv_min = [0.1]
  418 + iv_max = [10.]
  419 + plog=[1]
  420 + IF keyword_set(test) THEN BEGIN
  421 + iv_min = [.1]
  422 + iv_max = [10.]
  423 + plog=[1]
  424 + ENDIF
  425 + fpd=['(*!dustem_params).gas.g0','(*!dustem_params).g0','dustem_plugin_phangs_class_isrf_1','dustem_plugin_phangs_class_isrf_4'] ;ISRF class to be used
  426 + dustem_init,model=model
  427 + fortran_user=dustem_set_up_fortran(/random_name) ;use a random fortran number
  428 + !dustem_verbose=0
  429 + (*!dustem_params).KEYWORDS='quiet '+(*!dustem_params).KEYWORDS ;This makes Fortran be quiet too
  430 + ;stop
  431 + (*!dustem_params).grains[0].TYPE_KEYWORDS='logn-mix'
  432 + (*!dustem_params).grains[1].TYPE_KEYWORDS='logn-mix'
  433 + ;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn'
  434 + ;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn'
  435 + grain_keywords=[(*!dustem_params).grains.TYPE_KEYWORDS] ;This is to inform dustem_make_sed_table, which actually runs the models
  436 + ;=== select filters to be used for the grid (HOOPS, we need MUSE also, and maybe all filters ...)
  437 + ;filters=dustem_get_all_filter_names()
  438 + filters=get_phangs_filter_names(/order_by_wavelengths)
  439 + IF keyword_set(isrf_class) THEN BEGIN ;do just one table
  440 + use_isrf_class=isrf_class
  441 + ENDIF ELSE BEGIN
  442 + use_isrf_class=0 ;This is for Mathis field
  443 + ENDELSE
  444 + fiv=[-1.,1.,use_isrf_class,0.] ;This sets the ISRF class to the requested value. also removes ionising photons.
  445 + isrf_class_str='_isrfclass'+strtrim(use_isrf_class,2)
  446 + ;define the number of free parameters in the grid and the grid fits table name
  447 + IF keyword_set(test) THEN BEGIN
  448 + iv_Nvalues=[4]
  449 + ;table_name=dir+test_str+model+'_MuseISRF_JWST_G0_4Phangs_noionis'+isrf_class_str+'.fits'
  450 + ENDIF ELSE BEGIN
  451 + ;iv_Nvalues=[19]
  452 + iv_Nvalues=[12]
  453 + ENDELSE
  454 + table_name=dir+test_str+model+'_orig_MuseISRF_JWST_adaptedG0_4Phangs_noionis'+isrf_class_str+'.fits'
  455 + ;stop
  456 + 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, $
  457 + grain_keywords=grain_keywords
  458 + ;adding linear parameters (small grains abundances)
  459 + message,'Adding Dust abundances to grid',/continue
  460 + ;pd=['(*!dustem_params).grains[0].MDUST_O_MH','(*!dustem_params).grains[1].MDUST_O_MH','(*!dustem_params).grains[2].MDUST_O_MH']
  461 + pd=['(*!dustem_params).grains[0].MDUST_O_MH','(*!dustem_params).grains[2].MDUST_O_MH']
  462 + default_value=(dustem_get_param_values_default(pd))
  463 + iv_min=[default_value/10.]
  464 + ;iv_max=[default_value[0:1]*10.,default_value[2]*100.]
  465 + iv_max=[default_value[0]*10.,default_value[1]*100.]
  466 + ;iv_Nvalues=[11,11,11]
  467 + ;iv_Nvalues=[13,13,13]
  468 + iv_Nvalues=[13,13]
  469 + ;plog=[1,1,1]
  470 + plog=[1,1]
  471 + IF keyword_set(test) THEN BEGIN
  472 + iv_min=[default_value/2.]
  473 + iv_max=[default_value*2.]
  474 + ;iv_Nvalues=[3,3,3]
  475 + ;plog=[1,1,1]
  476 + iv_Nvalues=[3,3]
  477 + plog=[1,1]
  478 + ENDIF
  479 + output_table_name=dir+test_str+model+'_orig_MuseISRF_JWST_adaptedG0_Ypah1added_Yvsgadded_4Phangs_noionis'+isrf_class_str+'.fits'
  480 + dustem_add_linear_params2grid,table_name,pd,iv_min,iv_max,iv_Nvalues,out_filename=output_table_name,plog=plog,show_seds=show_seds
  481 + goto,the_end
  482 + END
  483 + 223: BEGIN ;This is same as 221, but using the standard DL07 GRAIN.DAT, ie, YPAH0 and Ypah1 keywords set to logn
  484 + ;Note: This assumes a variable Xion which is different from DL07 and no size depednence of Xion
  485 + !quiet=1
  486 + model='DL07' ;This is to use the DBP90 dust model
  487 + ;======= This is to produce a grid with the above model for G0 and PAH abundance and ISRF from the Muse data
  488 + ;======= Note that g0 is fixed parameters is actually not used, replaced by dustem_plugin_phangs_class_isrf_2
  489 + pd = [ $
  490 + 'dustem_plugin_phangs_class_isrf_2' $ ;G0 factor for Phangs ISRF classes
  491 + ]
  492 + iv_min = [0.1]
  493 + iv_max = [10.]
  494 + plog=[1]
  495 + IF keyword_set(test) THEN BEGIN
  496 + iv_min = [.1]
  497 + iv_max = [10.]
  498 + plog=[1]
  499 + ENDIF
  500 + fpd=['(*!dustem_params).gas.g0','(*!dustem_params).g0','dustem_plugin_phangs_class_isrf_1','dustem_plugin_phangs_class_isrf_4'] ;ISRF class to be used
  501 + dustem_init,model=model
  502 + fortran_user=dustem_set_up_fortran(/random_name) ;use a random fortran number
  503 + !dustem_verbose=0
  504 + (*!dustem_params).KEYWORDS='quiet '+(*!dustem_params).KEYWORDS ;This makes Fortran be quiet too
  505 + ;stop
  506 + (*!dustem_params).grains[0].TYPE_KEYWORDS='logn'
  507 + (*!dustem_params).grains[1].TYPE_KEYWORDS='logn'
  508 + grain_keywords=[(*!dustem_params).grains.TYPE_KEYWORDS] ;This is to inform dustem_make_sed_table, which actually runs the models
  509 + ;=== select filters to be used for the grid (HOOPS, we need MUSE also, and maybe all filters ...)
  510 + ;filters=dustem_get_all_filter_names()
  511 + filters=get_phangs_filter_names(/order_by_wavelengths)
  512 + IF keyword_set(isrf_class) THEN BEGIN ;do just one table
  513 + use_isrf_class=isrf_class
  514 + ENDIF ELSE BEGIN
  515 + use_isrf_class=0 ;This is for Mathis field
  516 + ENDELSE
  517 + fiv=[-1.,1.,use_isrf_class,0.] ;This sets the ISRF class to the requested value. also removes ionising photons.
  518 + isrf_class_str='_isrfclass'+strtrim(use_isrf_class,2)
  519 + ;define the number of free parameters in the grid and the grid fits table name
  520 + IF keyword_set(test) THEN BEGIN
  521 + iv_Nvalues=[4]
  522 + ENDIF ELSE BEGIN
  523 + iv_Nvalues=[12]
  524 + ENDELSE
  525 + table_name=dir+test_str+model+'_logn_MuseISRF_JWST_adaptedG0_4Phangs_noionis'+isrf_class_str+'.fits'
  526 + ;stop
  527 + 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, $
  528 + grain_keywords=grain_keywords
  529 + ;adding linear parameters (small grains abundances)
  530 + message,'Adding Dust abundances to grid',/continue
  531 + pd=['(*!dustem_params).grains[0].MDUST_O_MH','(*!dustem_params).grains[1].MDUST_O_MH','(*!dustem_params).grains[2].MDUST_O_MH']
  532 + ;pd=['(*!dustem_params).grains[0].MDUST_O_MH','(*!dustem_params).grains[2].MDUST_O_MH']
  533 + default_value=(dustem_get_param_values_default(pd))
  534 + iv_min=[default_value/10.]
  535 + iv_max=[default_value[0:1]*10.,default_value[2]*100.]
  536 + iv_Nvalues=[13,13,13]
  537 + plog=[1,1,1]
  538 + IF keyword_set(test) THEN BEGIN
  539 + iv_min=[default_value/2.]
  540 + iv_max=[default_value*2.]
  541 + iv_Nvalues=[3,3,3]
  542 + plog=[1,1,1]
  543 + ENDIF
  544 + output_table_name=dir+test_str+model+'_logn_MuseISRF_JWST_adaptedG0_Ypah1added_Yvsgadded_4Phangs_noionis'+isrf_class_str+'.fits'
  545 + dustem_add_linear_params2grid,table_name,pd,iv_min,iv_max,iv_Nvalues,out_filename=output_table_name,plog=plog,show_seds=show_seds
  546 + goto,the_end
  547 + END
408 548 ENDCASE
409 549  
410 550  
... ...
LabTools/IRAP/JPB/phangs_brute_force_fit_with_isrf_grid.pro
  Diff comments   View file @78db85d
... ... @@ -582,6 +582,8 @@ IF use_yvsg_option EQ 'MAP' THEN BEGIN
582 582 ENDIF
583 583  
584 584 FOR isrf_class=class_min,class_max DO BEGIN ;isrf_class runs from 1
  585 + ;print,isrf_class
  586 + ;stop
585 587 ;==== get G0prime of the template
586 588 IF isrf_class NE 0 THEN BEGIN
587 589 this_G0prime=G0prime[isrf_class]
... ... @@ -785,19 +787,28 @@ FOR isrf_class=class_min,class_max DO BEGIN ;isrf_class runs from 1
785 787 ;diffs=la_sub(*best_fit_seds[vid],*observed_seds[vid])
786 788 ;stop
787 789 Ygrain1[vid]=weighted_params[1];/fact ;Caution, extenssive quantities must be divided by normalization factor. We don't here
788   - Ygrain2[vid]=weighted_params[2];/fact
789   - Ygrain3[vid]=weighted_params[3];/fact
  790 + ;=== below is if PAH+ are in the table
  791 + ;Ygrain2[vid]=weighted_params[2];/fact
  792 + ;Ygrain3[vid]=weighted_params[3];/fact
  793 + ;=== below is if PAH+ are NOT in the table
  794 + Ygrain3[vid]=weighted_params[2];/fact
790 795 facts[vid]=fact
791 796 chi2s[vid]=chi2
792 797 rchi2s[vid]=rchi2
793 798 F0_hit[vid]=params_hit[0]
794 799 Ygrain1_hit[vid]=params_hit[1]
795   - Ygrain2_hit[vid]=params_hit[2]
796   - Ygrain3_hit[vid]=params_hit[3]
  800 + ;=== below is if PAH+ are in the table
  801 + ;Ygrain2_hit[vid]=params_hit[2]
  802 + ;Ygrain3_hit[vid]=params_hit[3]
  803 + ;=== below is if PAH+ are NOT in the table
  804 + Ygrain3_hit[vid]=params_hit[2]
797 805 dF0s[vid]=0. ;should use difference between fixed param and nearest table neighbor
798 806 dYgrain1[vid]=params_uncertainties[0]/fact ;extenssive quantities must be devided by normalization factor
799   - dYgrain2[vid]=params_uncertainties[1]/fact
800   - dYgrain3[vid]=params_uncertainties[2]/fact
  807 + ;=== below is if PAH+ are in the table
  808 + ;dYgrain2[vid]=params_uncertainties[1]/fact
  809 + ;dYgrain3[vid]=params_uncertainties[2]/fact
  810 + ;=== below is if PAH+ are NOT in the table
  811 + dYgrain3[vid]=params_uncertainties[1]/fact
801 812 ;print,params[0],params[1],params[2],fact,chi2
802 813 ENDIF ELSE BEGIN
803 814 message,'dustem_brute_force_fit status ='+strtrim(status,2),/continue
... ... @@ -1210,24 +1221,30 @@ cgoplot,replicate(alog10(default_value),2),10^(!y.crange),linestyle=2,color='bla
1210 1221 cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
1211 1222 cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
1212 1223  
1213   -ind=where(Ygrain2 NE la_undef())
1214   -default_value=dustem_get_param_values_default('(*!dustem_params).grains[1].MDUST_O_MH')
1215   -parval_min=((*!dustem_grid).PMIN_VALUES)[2]
1216   -parval_max=((*!dustem_grid).PMAX_VALUES)[2]
1217   -;xrange=alog10([default_value/range_fact,default_value*range_fact])
1218   -xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
1219   -yrange=[1,max(res)]
1220   -res=histogram(alog10(Ygrain2[ind]),locations=xv,Nbins=100)
1221   -cgplot,xv,res,psym=10,title='Y Grain2 histogram',xtit='log10(Y)',ytit='Number',/ylog,/ysty,yrange=yrange,xrange=xrange,/xsty
1222   -cgoplot,replicate(alog10(default_value),2),10^(!y.crange),linestyle=2,color='black'
1223   -cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
1224   -cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  1224 +ind=where(Ygrain2 NE la_undef(),count)
  1225 +IF count NE 0 THEN BEGIN
  1226 + default_value=dustem_get_param_values_default('(*!dustem_params).grains[1].MDUST_O_MH')
  1227 + parval_min=((*!dustem_grid).PMIN_VALUES)[2]
  1228 + parval_max=((*!dustem_grid).PMAX_VALUES)[2]
  1229 + ;xrange=alog10([default_value/range_fact,default_value*range_fact])
  1230 + xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
  1231 + yrange=[1,max(res)]
  1232 + res=histogram(alog10(Ygrain2[ind]),locations=xv,Nbins=100)
  1233 + cgplot,xv,res,psym=10,title='Y Grain2 histogram',xtit='log10(Y)',ytit='Number',/ylog,/ysty,yrange=yrange,xrange=xrange,/xsty
  1234 + cgoplot,replicate(alog10(default_value),2),10^(!y.crange),linestyle=2,color='black'
  1235 + cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
  1236 + cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  1237 +ENDIF
1225 1238  
1226 1239 ind=where(Ygrain3 NE la_undef(),count3)
1227 1240 IF count3 NE 0 THEN BEGIN
1228 1241 default_value=dustem_get_param_values_default('(*!dustem_params).grains[2].MDUST_O_MH')
1229   - parval_min=((*!dustem_grid).PMIN_VALUES)[3]
1230   - parval_max=((*!dustem_grid).PMAX_VALUES)[3]
  1242 + ;=== IF PAH+ are in table
  1243 + ;parval_min=((*!dustem_grid).PMIN_VALUES)[3]
  1244 + ;parval_max=((*!dustem_grid).PMAX_VALUES)[3]
  1245 + ;=== IF PAH+ are not in table
  1246 + parval_min=((*!dustem_grid).PMIN_VALUES)[2]
  1247 + parval_max=((*!dustem_grid).PMAX_VALUES)[2]
1231 1248 ;xrange=alog10([default_value/range_fact,default_value*range_fact])
1232 1249 xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
1233 1250 yrange=[1,max(res)]
... ...
LabTools/IRAP/JPB/phangs_isrf_pipeline.pro
  Diff comments   View file @78db85d
... ... @@ -46,9 +46,14 @@ make_phangs_isrf_classes,bidon,source_name='ngc0628',/do_postscripts
46 46  
47 47 ;make_sed_phangs_tables,/test,/show_seds ;test of plugin dustem_plugin_phangs_class_isrf.pro
48 48 ;Generally very long. Do it using nohup (see make_phangs_grid.pro)
49   -make_sed_phangs_tables,grid_type=2,isrf_class=15
50   -make_sed_phangs_tables,grid_type=2,isrf_class=0
  49 +make_sed_phangs_tables,grid_type=223,isrf_class=15,/test
  50 +make_sed_phangs_tables,grid_type=222,isrf_class=14
  51 +make_sed_phangs_tables,grid_type=222,isrf_class=13
  52 +make_sed_phangs_tables,grid_type=222,isrf_class=12
51 53 ;...
  54 +FOR i=4L,26 DO make_sed_phangs_tables,grid_type=223,isrf_class=i
  55 +FOR i=16L,26 DO make_sed_phangs_tables,grid_type=222,isrf_class=i
  56 +make_sed_phangs_tables,grid_type=222,isrf_class=0
52 57  
53 58 ;=== Fit the SEDs with the ISRF Grids
54 59 table_name_root='_MuseISRF_JWST_adaptedG0_Ypah1added_Ypah2added_adaptedYvsgadded_4Phangs_noionis'
... ... @@ -117,6 +122,14 @@ phangs_smooth_muse_isrf,'ngc0628',resolution_filter='PACS3',/save,/show,/nostop
117 122  
118 123 make_phangs_isrf_classes,bidon,source_name='ngc0628',/save,resolution_filter='PACS3'
119 124  
  125 +;===test new tables
  126 +dustem_init
  127 +table_name_root='_orig_MuseISRF_JWST_adaptedG0_Ypah1added_Yvsgadded_4Phangs_noionis'
  128 +class_min_max=[12,15]
  129 +phangs_brute_force_fit_with_isrf_grid,model='DL07',source_name='ngc0628',/include_herschel,/normalize,/fit_G0,/save,resolution_filter='PACS3' $
  130 + ,/subtract_star_light,/show,/nostop,table_name_root=table_name_root,class_min_max=class_min_max
  131 +
  132 +
120 133 table_name_root='_MuseISRF_JWST_adaptedG0_Ypah1added_Ypah2added_adaptedYvsgadded_4Phangs_noionis'
121 134 class_min_max=0
122 135 phangs_brute_force_fit_with_isrf_grid,model='DL07',source_name='ngc0628',/include_herschel,/normalize,/fit_G0,/save,resolution_filter='PACS3' $
... ... @@ -295,6 +308,12 @@ phangs_smooth_muse_isrf,'ngc0628',resolution_filter='PACS1',/save,/nostop
295 308  
296 309 make_phangs_isrf_classes,bidon,source_name='ngc0628',/save,resolution_filter='PACS1'
297 310  
  311 +dustem_init
  312 +table_name_root='_orig_MuseISRF_JWST_adaptedG0_Ypah1added_Yvsgadded_4Phangs_noionis'
  313 +;class_min_max=[12,15]
  314 +phangs_brute_force_fit_with_isrf_grid,model='DL07',source_name='ngc0628',/include_herschel,/normalize,/fit_G0,/save,resolution_filter='PACS1' $
  315 + ,/subtract_star_light,/show,/nostop,table_name_root=table_name_root,class_min_max=class_min_max
  316 +
298 317 table_name_root='_MuseISRF_JWST_adaptedG0_Ypah1added_Ypah2added_adaptedYvsgadded_4Phangs_noionis'
299 318 class_min_max=0
300 319 phangs_brute_force_fit_with_isrf_grid,model='DL07',source_name='ngc0628',/include_herschel,/normalize,/fit_G0,/save,resolution_filter='PACS1' $
... ...
LabTools/IRAP/JPB/srun/make_phangs_grids.pro
  Diff comments   View file @78db85d
... ... @@ -19,8 +19,8 @@ t0=systime(/sec)
19 19 bidon=1
20 20  
21 21 ;make_sed_phangs_tables,grid_type=221,isrf_class=0 ;alma1 (For Mathis field)
22   -FOR i=5L,26 DO make_sed_phangs_tables,grid_type=221,isrf_class=i ;alma1 (For grid ISRF field)
23   -
  22 +;FOR i=5L,26 DO make_sed_phangs_tables,grid_type=221,isrf_class=i ;alma1 (For grid ISRF field)
  23 +FOR i=5L,26 DO make_sed_phangs_tables,grid_type=223,isrf_class=i ;This is to compute DL07 with logn keyword for both PAH0 and PAH1
24 24  
25 25 ;make grid
26 26 ;=== This is DL07 model single G0s
... ...