From 78db85de6e95a68dac4ec6b5545b1a2b576f4841 Mon Sep 17 00:00:00 2001 From: Jean-Philippe Bernard Date: Tue, 17 Sep 2024 14:16:16 +0200 Subject: [PATCH] modified for different flavors of DL07 dealing with APH1 and PAH0 --- LabTools/IRAP/JPB/make_sed_phangs_tables.pro | 140 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ LabTools/IRAP/JPB/phangs_brute_force_fit_with_isrf_grid.pro | 57 +++++++++++++++++++++++++++++++++++++-------------------- LabTools/IRAP/JPB/phangs_isrf_pipeline.pro | 23 +++++++++++++++++++++-- LabTools/IRAP/JPB/srun/make_phangs_grids.pro | 4 ++-- 4 files changed, 200 insertions(+), 24 deletions(-) diff --git a/LabTools/IRAP/JPB/make_sed_phangs_tables.pro b/LabTools/IRAP/JPB/make_sed_phangs_tables.pro index ceee59b..2a14ddf 100644 --- a/LabTools/IRAP/JPB/make_sed_phangs_tables.pro +++ b/LabTools/IRAP/JPB/make_sed_phangs_tables.pro @@ -405,6 +405,146 @@ CASE use_grid_type OF dustem_add_linear_params2grid,table_name,pd,iv_min,iv_max,iv_Nvalues,out_filename=output_table_name,plog=plog,show_seds=show_seds goto,the_end END + 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 + ;Note: This does not really work, because PAH+ contribution does not get multiplied linearly by abundance ... + !quiet=1 + model='DL07' ;This is to use the DBP90 dust model + ;======= This is to produce a grid with the above model for G0 and PAH abundance and ISRF from the Muse data + ;======= Note that g0 is fixed parameters is actually not used, replaced by dustem_plugin_phangs_class_isrf_2 + pd = [ $ + 'dustem_plugin_phangs_class_isrf_2' $ ;G0 factor for Phangs ISRF classes + ] + iv_min = [0.1] + iv_max = [10.] + plog=[1] + IF keyword_set(test) THEN BEGIN + iv_min = [.1] + iv_max = [10.] + plog=[1] + ENDIF + 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 + dustem_init,model=model + fortran_user=dustem_set_up_fortran(/random_name) ;use a random fortran number + !dustem_verbose=0 + (*!dustem_params).KEYWORDS='quiet '+(*!dustem_params).KEYWORDS ;This makes Fortran be quiet too + ;stop + (*!dustem_params).grains[0].TYPE_KEYWORDS='logn-mix' + (*!dustem_params).grains[1].TYPE_KEYWORDS='logn-mix' + ;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn' + ;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn' + grain_keywords=[(*!dustem_params).grains.TYPE_KEYWORDS] ;This is to inform dustem_make_sed_table, which actually runs the models + ;=== select filters to be used for the grid (HOOPS, we need MUSE also, and maybe all filters ...) + ;filters=dustem_get_all_filter_names() + filters=get_phangs_filter_names(/order_by_wavelengths) + IF keyword_set(isrf_class) THEN BEGIN ;do just one table + use_isrf_class=isrf_class + ENDIF ELSE BEGIN + use_isrf_class=0 ;This is for Mathis field + ENDELSE + fiv=[-1.,1.,use_isrf_class,0.] ;This sets the ISRF class to the requested value. also removes ionising photons. + isrf_class_str='_isrfclass'+strtrim(use_isrf_class,2) + ;define the number of free parameters in the grid and the grid fits table name + IF keyword_set(test) THEN BEGIN + iv_Nvalues=[4] + ;table_name=dir+test_str+model+'_MuseISRF_JWST_G0_4Phangs_noionis'+isrf_class_str+'.fits' + ENDIF ELSE BEGIN + ;iv_Nvalues=[19] + iv_Nvalues=[12] + ENDELSE + table_name=dir+test_str+model+'_orig_MuseISRF_JWST_adaptedG0_4Phangs_noionis'+isrf_class_str+'.fits' + ;stop + 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, $ + grain_keywords=grain_keywords + ;adding linear parameters (small grains abundances) + message,'Adding Dust abundances to grid',/continue + ;pd=['(*!dustem_params).grains[0].MDUST_O_MH','(*!dustem_params).grains[1].MDUST_O_MH','(*!dustem_params).grains[2].MDUST_O_MH'] + pd=['(*!dustem_params).grains[0].MDUST_O_MH','(*!dustem_params).grains[2].MDUST_O_MH'] + default_value=(dustem_get_param_values_default(pd)) + iv_min=[default_value/10.] + ;iv_max=[default_value[0:1]*10.,default_value[2]*100.] + iv_max=[default_value[0]*10.,default_value[1]*100.] + ;iv_Nvalues=[11,11,11] + ;iv_Nvalues=[13,13,13] + iv_Nvalues=[13,13] + ;plog=[1,1,1] + plog=[1,1] + IF keyword_set(test) THEN BEGIN + iv_min=[default_value/2.] + iv_max=[default_value*2.] + ;iv_Nvalues=[3,3,3] + ;plog=[1,1,1] + iv_Nvalues=[3,3] + plog=[1,1] + ENDIF + output_table_name=dir+test_str+model+'_orig_MuseISRF_JWST_adaptedG0_Ypah1added_Yvsgadded_4Phangs_noionis'+isrf_class_str+'.fits' + dustem_add_linear_params2grid,table_name,pd,iv_min,iv_max,iv_Nvalues,out_filename=output_table_name,plog=plog,show_seds=show_seds + goto,the_end + END + 223: BEGIN ;This is same as 221, but using the standard DL07 GRAIN.DAT, ie, YPAH0 and Ypah1 keywords set to logn + ;Note: This assumes a variable Xion which is different from DL07 and no size depednence of Xion + !quiet=1 + model='DL07' ;This is to use the DBP90 dust model + ;======= This is to produce a grid with the above model for G0 and PAH abundance and ISRF from the Muse data + ;======= Note that g0 is fixed parameters is actually not used, replaced by dustem_plugin_phangs_class_isrf_2 + pd = [ $ + 'dustem_plugin_phangs_class_isrf_2' $ ;G0 factor for Phangs ISRF classes + ] + iv_min = [0.1] + iv_max = [10.] + plog=[1] + IF keyword_set(test) THEN BEGIN + iv_min = [.1] + iv_max = [10.] + plog=[1] + ENDIF + 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 + dustem_init,model=model + fortran_user=dustem_set_up_fortran(/random_name) ;use a random fortran number + !dustem_verbose=0 + (*!dustem_params).KEYWORDS='quiet '+(*!dustem_params).KEYWORDS ;This makes Fortran be quiet too + ;stop + (*!dustem_params).grains[0].TYPE_KEYWORDS='logn' + (*!dustem_params).grains[1].TYPE_KEYWORDS='logn' + grain_keywords=[(*!dustem_params).grains.TYPE_KEYWORDS] ;This is to inform dustem_make_sed_table, which actually runs the models + ;=== select filters to be used for the grid (HOOPS, we need MUSE also, and maybe all filters ...) + ;filters=dustem_get_all_filter_names() + filters=get_phangs_filter_names(/order_by_wavelengths) + IF keyword_set(isrf_class) THEN BEGIN ;do just one table + use_isrf_class=isrf_class + ENDIF ELSE BEGIN + use_isrf_class=0 ;This is for Mathis field + ENDELSE + fiv=[-1.,1.,use_isrf_class,0.] ;This sets the ISRF class to the requested value. also removes ionising photons. + isrf_class_str='_isrfclass'+strtrim(use_isrf_class,2) + ;define the number of free parameters in the grid and the grid fits table name + IF keyword_set(test) THEN BEGIN + iv_Nvalues=[4] + ENDIF ELSE BEGIN + iv_Nvalues=[12] + ENDELSE + table_name=dir+test_str+model+'_logn_MuseISRF_JWST_adaptedG0_4Phangs_noionis'+isrf_class_str+'.fits' + ;stop + 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, $ + grain_keywords=grain_keywords + ;adding linear parameters (small grains abundances) + message,'Adding Dust abundances to grid',/continue + pd=['(*!dustem_params).grains[0].MDUST_O_MH','(*!dustem_params).grains[1].MDUST_O_MH','(*!dustem_params).grains[2].MDUST_O_MH'] + ;pd=['(*!dustem_params).grains[0].MDUST_O_MH','(*!dustem_params).grains[2].MDUST_O_MH'] + default_value=(dustem_get_param_values_default(pd)) + iv_min=[default_value/10.] + iv_max=[default_value[0:1]*10.,default_value[2]*100.] + iv_Nvalues=[13,13,13] + plog=[1,1,1] + IF keyword_set(test) THEN BEGIN + iv_min=[default_value/2.] + iv_max=[default_value*2.] + iv_Nvalues=[3,3,3] + plog=[1,1,1] + ENDIF + output_table_name=dir+test_str+model+'_logn_MuseISRF_JWST_adaptedG0_Ypah1added_Yvsgadded_4Phangs_noionis'+isrf_class_str+'.fits' + dustem_add_linear_params2grid,table_name,pd,iv_min,iv_max,iv_Nvalues,out_filename=output_table_name,plog=plog,show_seds=show_seds + goto,the_end + END ENDCASE diff --git a/LabTools/IRAP/JPB/phangs_brute_force_fit_with_isrf_grid.pro b/LabTools/IRAP/JPB/phangs_brute_force_fit_with_isrf_grid.pro index 83332a3..d4b6494 100644 --- a/LabTools/IRAP/JPB/phangs_brute_force_fit_with_isrf_grid.pro +++ b/LabTools/IRAP/JPB/phangs_brute_force_fit_with_isrf_grid.pro @@ -582,6 +582,8 @@ IF use_yvsg_option EQ 'MAP' THEN BEGIN ENDIF FOR isrf_class=class_min,class_max DO BEGIN ;isrf_class runs from 1 + ;print,isrf_class + ;stop ;==== get G0prime of the template IF isrf_class NE 0 THEN BEGIN this_G0prime=G0prime[isrf_class] @@ -785,19 +787,28 @@ FOR isrf_class=class_min,class_max DO BEGIN ;isrf_class runs from 1 ;diffs=la_sub(*best_fit_seds[vid],*observed_seds[vid]) ;stop Ygrain1[vid]=weighted_params[1];/fact ;Caution, extenssive quantities must be divided by normalization factor. We don't here - Ygrain2[vid]=weighted_params[2];/fact - Ygrain3[vid]=weighted_params[3];/fact + ;=== below is if PAH+ are in the table + ;Ygrain2[vid]=weighted_params[2];/fact + ;Ygrain3[vid]=weighted_params[3];/fact + ;=== below is if PAH+ are NOT in the table + Ygrain3[vid]=weighted_params[2];/fact facts[vid]=fact chi2s[vid]=chi2 rchi2s[vid]=rchi2 F0_hit[vid]=params_hit[0] Ygrain1_hit[vid]=params_hit[1] - Ygrain2_hit[vid]=params_hit[2] - Ygrain3_hit[vid]=params_hit[3] + ;=== below is if PAH+ are in the table + ;Ygrain2_hit[vid]=params_hit[2] + ;Ygrain3_hit[vid]=params_hit[3] + ;=== below is if PAH+ are NOT in the table + Ygrain3_hit[vid]=params_hit[2] dF0s[vid]=0. ;should use difference between fixed param and nearest table neighbor dYgrain1[vid]=params_uncertainties[0]/fact ;extenssive quantities must be devided by normalization factor - dYgrain2[vid]=params_uncertainties[1]/fact - dYgrain3[vid]=params_uncertainties[2]/fact + ;=== below is if PAH+ are in the table + ;dYgrain2[vid]=params_uncertainties[1]/fact + ;dYgrain3[vid]=params_uncertainties[2]/fact + ;=== below is if PAH+ are NOT in the table + dYgrain3[vid]=params_uncertainties[1]/fact ;print,params[0],params[1],params[2],fact,chi2 ENDIF ELSE BEGIN 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 cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue' cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red' -ind=where(Ygrain2 NE la_undef()) -default_value=dustem_get_param_values_default('(*!dustem_params).grains[1].MDUST_O_MH') -parval_min=((*!dustem_grid).PMIN_VALUES)[2] -parval_max=((*!dustem_grid).PMAX_VALUES)[2] -;xrange=alog10([default_value/range_fact,default_value*range_fact]) -xrange=enlarge_range(alog10([parval_min,parval_max]),0.1) -yrange=[1,max(res)] -res=histogram(alog10(Ygrain2[ind]),locations=xv,Nbins=100) -cgplot,xv,res,psym=10,title='Y Grain2 histogram',xtit='log10(Y)',ytit='Number',/ylog,/ysty,yrange=yrange,xrange=xrange,/xsty -cgoplot,replicate(alog10(default_value),2),10^(!y.crange),linestyle=2,color='black' -cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue' -cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red' +ind=where(Ygrain2 NE la_undef(),count) +IF count NE 0 THEN BEGIN + default_value=dustem_get_param_values_default('(*!dustem_params).grains[1].MDUST_O_MH') + parval_min=((*!dustem_grid).PMIN_VALUES)[2] + parval_max=((*!dustem_grid).PMAX_VALUES)[2] + ;xrange=alog10([default_value/range_fact,default_value*range_fact]) + xrange=enlarge_range(alog10([parval_min,parval_max]),0.1) + yrange=[1,max(res)] + res=histogram(alog10(Ygrain2[ind]),locations=xv,Nbins=100) + cgplot,xv,res,psym=10,title='Y Grain2 histogram',xtit='log10(Y)',ytit='Number',/ylog,/ysty,yrange=yrange,xrange=xrange,/xsty + cgoplot,replicate(alog10(default_value),2),10^(!y.crange),linestyle=2,color='black' + cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue' + cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red' +ENDIF ind=where(Ygrain3 NE la_undef(),count3) IF count3 NE 0 THEN BEGIN default_value=dustem_get_param_values_default('(*!dustem_params).grains[2].MDUST_O_MH') - parval_min=((*!dustem_grid).PMIN_VALUES)[3] - parval_max=((*!dustem_grid).PMAX_VALUES)[3] + ;=== IF PAH+ are in table + ;parval_min=((*!dustem_grid).PMIN_VALUES)[3] + ;parval_max=((*!dustem_grid).PMAX_VALUES)[3] + ;=== IF PAH+ are not in table + parval_min=((*!dustem_grid).PMIN_VALUES)[2] + parval_max=((*!dustem_grid).PMAX_VALUES)[2] ;xrange=alog10([default_value/range_fact,default_value*range_fact]) xrange=enlarge_range(alog10([parval_min,parval_max]),0.1) yrange=[1,max(res)] diff --git a/LabTools/IRAP/JPB/phangs_isrf_pipeline.pro b/LabTools/IRAP/JPB/phangs_isrf_pipeline.pro index 9f4d72f..6c38e53 100644 --- a/LabTools/IRAP/JPB/phangs_isrf_pipeline.pro +++ b/LabTools/IRAP/JPB/phangs_isrf_pipeline.pro @@ -46,9 +46,14 @@ make_phangs_isrf_classes,bidon,source_name='ngc0628',/do_postscripts ;make_sed_phangs_tables,/test,/show_seds ;test of plugin dustem_plugin_phangs_class_isrf.pro ;Generally very long. Do it using nohup (see make_phangs_grid.pro) -make_sed_phangs_tables,grid_type=2,isrf_class=15 -make_sed_phangs_tables,grid_type=2,isrf_class=0 +make_sed_phangs_tables,grid_type=223,isrf_class=15,/test +make_sed_phangs_tables,grid_type=222,isrf_class=14 +make_sed_phangs_tables,grid_type=222,isrf_class=13 +make_sed_phangs_tables,grid_type=222,isrf_class=12 ;... +FOR i=4L,26 DO make_sed_phangs_tables,grid_type=223,isrf_class=i +FOR i=16L,26 DO make_sed_phangs_tables,grid_type=222,isrf_class=i +make_sed_phangs_tables,grid_type=222,isrf_class=0 ;=== Fit the SEDs with the ISRF Grids 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 make_phangs_isrf_classes,bidon,source_name='ngc0628',/save,resolution_filter='PACS3' +;===test new tables +dustem_init +table_name_root='_orig_MuseISRF_JWST_adaptedG0_Ypah1added_Yvsgadded_4Phangs_noionis' +class_min_max=[12,15] +phangs_brute_force_fit_with_isrf_grid,model='DL07',source_name='ngc0628',/include_herschel,/normalize,/fit_G0,/save,resolution_filter='PACS3' $ + ,/subtract_star_light,/show,/nostop,table_name_root=table_name_root,class_min_max=class_min_max + + table_name_root='_MuseISRF_JWST_adaptedG0_Ypah1added_Ypah2added_adaptedYvsgadded_4Phangs_noionis' class_min_max=0 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 make_phangs_isrf_classes,bidon,source_name='ngc0628',/save,resolution_filter='PACS1' +dustem_init +table_name_root='_orig_MuseISRF_JWST_adaptedG0_Ypah1added_Yvsgadded_4Phangs_noionis' +;class_min_max=[12,15] +phangs_brute_force_fit_with_isrf_grid,model='DL07',source_name='ngc0628',/include_herschel,/normalize,/fit_G0,/save,resolution_filter='PACS1' $ + ,/subtract_star_light,/show,/nostop,table_name_root=table_name_root,class_min_max=class_min_max + table_name_root='_MuseISRF_JWST_adaptedG0_Ypah1added_Ypah2added_adaptedYvsgadded_4Phangs_noionis' class_min_max=0 phangs_brute_force_fit_with_isrf_grid,model='DL07',source_name='ngc0628',/include_herschel,/normalize,/fit_G0,/save,resolution_filter='PACS1' $ diff --git a/LabTools/IRAP/JPB/srun/make_phangs_grids.pro b/LabTools/IRAP/JPB/srun/make_phangs_grids.pro index 5100500..304b41e 100644 --- a/LabTools/IRAP/JPB/srun/make_phangs_grids.pro +++ b/LabTools/IRAP/JPB/srun/make_phangs_grids.pro @@ -19,8 +19,8 @@ t0=systime(/sec) bidon=1 ;make_sed_phangs_tables,grid_type=221,isrf_class=0 ;alma1 (For Mathis field) -FOR i=5L,26 DO make_sed_phangs_tables,grid_type=221,isrf_class=i ;alma1 (For grid ISRF field) - +;FOR i=5L,26 DO make_sed_phangs_tables,grid_type=221,isrf_class=i ;alma1 (For grid ISRF field) +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 ;make grid ;=== This is DL07 model single G0s -- libgit2 0.21.2