Commit cd2aaf3e1ffb2d35eededc837c3de94180518f3e

Authored by Jean-Philippe Bernard
1 parent 8b306a32
Exists in master

improved

LabTools/IRAP/JPB/dustem_merge_grids.pro
... ... @@ -4,13 +4,14 @@ dustem_read_grid_table2arrays,input_table_name1 $
4 4 ,Ngrains1,model1,Ncomb1,NFparams1,INCSP1,polar1,double1 $
5 5 ,parameter_description1,parameter_pmin1,parameter_pmax1,parameter_pnval1,parameter_plog1,fparameter_description1,fparameter_value1 $
6 6 ,param_values1,Nparams1,seds1,filters1,Nfilters1,wavs1,Nwavs1 $
7   - ,I_tot_vec1,ext_tot_vec1,I_grains_vec1,ext_grains_vec1
  7 + ,I_tot_vec1,ext_tot_vec1,I_grains_vec1,ext_grains_vec1,grain_keywords=grain_keywords1
  8 +
8 9  
9 10 dustem_read_grid_table2arrays,input_table_name2 $
10 11 ,Ngrains2,model2,Ncomb2,NFparams2,INCSP2,polar2,double2 $
11 12 ,parameter_description2,parameter_pmin2,parameter_pmax2,parameter_pnval2,parameter_plog2,fparameter_description2,fparameter_value2 $
12 13 ,param_values2,Nparams2,seds2,filters2,Nfilters2,wavs2,Nwavs2 $
13   - ,I_tot_vec2,ext_tot_vec2,I_grains_vec2,ext_grains_vec2
  14 + ,I_tot_vec2,ext_tot_vec2,I_grains_vec2,ext_grains_vec2,grain_keywords=grain_keywords2
14 15  
15 16 help,Ngrains1,model1,Ncomb1,NFparams1,INCSP1,polar1,double1
16 17 help,Ngrains2,model2,Ncomb2,NFparams2,INCSP2,polar2,double2
... ... @@ -30,7 +31,7 @@ Nfparams1=n_elements(fparameter_description1)
30 31 Nfparams2=n_elements(fparameter_description2)
31 32  
32 33 ;stop
33   -;==== This below is to merge tables with complmentray values of one free-parameter
  34 +;==== This below is to merge tables with complementary values of one free-parameter
34 35 IF model1 NE model2 THEN BEGIN
35 36 message,'Models in the two grids are different',/continue
36 37 stop
... ... @@ -59,6 +60,12 @@ IF double1 NE double2 THEN BEGIN
59 60 message,'data type in the two grids are different',/continue
60 61 stop
61 62 ENDIF
  63 +FOR i=0L,Ngrains1-1 DO BEGIN
  64 + IF grain_keywords1[i] NE grain_keywords2[i] THEN BEGIN
  65 + message,'grain keywords are different in the two grids',/continue
  66 + stop
  67 + ENDIF
  68 +ENDFOR
62 69  
63 70 Nparams=Nparams1
64 71 Nfparams=Nfparams1
... ... @@ -72,6 +79,13 @@ FOR i=0L,Nparams-1 DO BEGIN
72 79 ENDIF
73 80 ENDFOR
74 81  
  82 +FOR i=0L,Nparams-1 DO BEGIN
  83 + IF parameter_plog1[i] NE parameter_plog2[i] THEN BEGIN
  84 + message,'Found different parameter log flags in the two grids',/continue
  85 + stop
  86 + ENDIF
  87 +ENDFOR
  88 +
75 89 FOR i=0L,Nfparams-1 DO BEGIN
76 90 IF fparameter_description1[i] NE fparameter_description2[i] THEN BEGIN
77 91 message,'Found different fixed parameter description in the two grids',/continue
... ... @@ -102,6 +116,7 @@ FOR i=0L,Nparams-1 DO BEGIN
102 116 parameter_pmin[i]=min([parameter_pmin1[i],parameter_pmin2[i]])
103 117 parameter_pmax[i]=max([parameter_pmax1[i],parameter_pmax2[i]])
104 118 parameter_pnval[i]=parameter_pnval1[i]+parameter_pnval2[i]
  119 + parameter_plog[i]=parameter_plog1[i] ;which has to be equal to parameter_plog2[i]
105 120 ENDFOR
106 121  
107 122 fparameter_description=fparameter_description1
... ...
LabTools/IRAP/JPB/phangs_brute_force_fit_with_isrf_grid.pro
... ... @@ -554,6 +554,8 @@ st=replicate(one_st,30)
554 554  
555 555 one_st={isrf_class:0L,Nvor:0L,nhitF0up:0,nhitF0down:0,nhitYG1up:0,nhitYG1down:0,nhitYG2up:0,nhitYG2down:0,nhitYG3up:0,nhitYG3down:0}
556 556 hit_st=replicate(one_st,30)
  557 +one_st={isrf_class:0L,Nvor:0L,nhitF0up:0.,nhitF0down:0.,nhitYG1up:0.,nhitYG1down:0.,nhitYG2up:0.,nhitYG2down:0.,nhitYG3up:0.,nhitYG3down:0.}
  558 +hit_st_perc=replicate(one_st,30)
557 559  
558 560 ;indcalz=where(all_wavs GE 0.912 AND all_wavs LT 2.2,count_calz) ;Where Calzeti's correction is valid
559 561 indcalz=where(all_wavs LT 2.2,count_calz) ;Where Calzeti's correction is valid
... ... @@ -846,24 +848,35 @@ FOR isrf_class=class_min,class_max DO BEGIN ;isrf_class runs from 1
846 848 st[isrf_class].max_rchi2=la_max(rchi2s[ind_class])
847 849 hit_st[isrf_class].isrf_class=isrf_class
848 850 hit_st[isrf_class].Nvor=Nvor_class
  851 + hit_st_perc[isrf_class].isrf_class=isrf_class
  852 + hit_st_perc[isrf_class].Nvor=Nvor_class
849 853 ind=where(F0_hit[ind_class] EQ 1,count)
850 854 hit_st[isrf_class].nhitF0up=count
  855 + hit_st_perc[isrf_class].nhitF0up=(1.*count)/Nvor_class*100.
851 856 ind=where(F0_hit[ind_class] EQ -1,count)
852 857 hit_st[isrf_class].nhitF0down=count
  858 + hit_st_perc[isrf_class].nhitF0down=(1.*count)/Nvor_class*100.
853 859 ind=where(Ygrain1_hit[ind_class] EQ 1,count)
854 860 hit_st[isrf_class].nhitYG1up=count
  861 + hit_st_perc[isrf_class].nhitYG1up=(1.*count)/Nvor_class*100.
855 862 ind=where(Ygrain1_hit[ind_class] EQ -1,count)
856 863 hit_st[isrf_class].nhitYG1down=count
  864 + hit_st_perc[isrf_class].nhitYG1down=(1.*count)/Nvor_class*100.
857 865 ind=where(Ygrain2_hit[ind_class] EQ 1,count)
858 866 hit_st[isrf_class].nhitYG2up=count
  867 + hit_st_perc[isrf_class].nhitYG2up=(1.*count)/Nvor_class*100.
859 868 ind=where(Ygrain2_hit[ind_class] EQ -1,count)
860 869 hit_st[isrf_class].nhitYG2down=count
  870 + hit_st_perc[isrf_class].nhitYG2down=(1.*count)/Nvor_class*100.
861 871 ind=where(Ygrain3_hit[ind_class] EQ 1,count)
862 872 hit_st[isrf_class].nhitYG3up=count
  873 + hit_st_perc[isrf_class].nhitYG3up=(1.*count)/Nvor_class*100.
863 874 ind=where(Ygrain3_hit[ind_class] EQ -1,count)
864 875 hit_st[isrf_class].nhitYG3down=count
  876 + hit_st_perc[isrf_class].nhitYG3down=(1.*count)/Nvor_class*100.
865 877 write_xcat,st,'/tmp/'+getenv('LOGNAME')+'/essai.xcat'
866   - write_xcat,hit_st,'/tmp/'+getenv('LOGNAME')+'/essai.xcat'
  878 + ;write_xcat,hit_st,'/tmp/'+getenv('LOGNAME')+'/essai.xcat'
  879 + write_xcat,hit_st_perc,'/tmp/'+getenv('LOGNAME')+'/essai.xcat'
867 880 ENDIF ELSE BEGIN
868 881 message,'No Voronoi bins found with ISRF class '+strtrim(isrf_class,2),/continue
869 882 ENDELSE
... ... @@ -873,7 +886,7 @@ stats_st=st
873 886 IF keyword_set(save) THEN BEGIN
874 887 dustem_grid_saved=!dustem_grid
875 888 file_save=data_dir+use_source_name+use_model+'_logn_JWST_G0_YPAH_YVSG_on-voronoi_with-ISRFclasses'+norm_str+mathis_str+reso_str+fit_G0_str+set_yvsg_str+'.sav'
876   - save,all_filters,use_fit_filters,dustem_grid_saved,stellar_normalizations,stats_st,EBVs,F0s,F0s_predicted,G0s,G0s_predicted, $
  889 + save,all_filters,use_fit_filters,dustem_grid_saved,stellar_normalizations,stats_st,hit_st,hit_st_perc,EBVs,F0s,F0s_predicted,G0s,G0s_predicted, $
877 890 Ygrain1,Ygrain2,Ygrain3,facts,chi2s,rchi2s,dF0s,dYgrain1,dYgrain2,dYgrain3, $
878 891 F0_hit,Ygrain1_hit,Ygrain2_hit,Ygrain3_hit,NH_values,best_fit_seds,observed_seds,file=file_save
879 892 message,'Wrote '+file_save,/continue
... ... @@ -1511,7 +1524,7 @@ imrange=[-1.,+1.0] ;with Herschel data
1511 1524 imrange=[-1.,+1.5] ;with Herschel data
1512 1525 title='log10(F0)'
1513 1526 image_cont20,la_log10(F0_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st
1514   -stop
  1527 +;stop
1515 1528 IF keyword_set(do_postscripts) THEN BEGIN
1516 1529 ps_file=ps_dir+use_source_name+'_F0_map'+reso_str+fit_G0_str+mathis_str+'.ps'
1517 1530 image_cont20,la_log10(F0_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st, $
... ... @@ -1761,8 +1774,60 @@ cgplot,chi2s,Ygrain1,/xlog,/ylog,xr=chi2_range,yr=Ypah_range,/xsty,psym=3,xtit='
1761 1774 cgplot,chi2s,Ygrain2,/xlog,/ylog,xr=chi2_range,yr=Ypah_range,/xsty,psym=3,xtit='chi2',ytit='Y',title='Grain 2'
1762 1775 cgplot,chi2s,Ygrain3,/xlog,/ylog,xr=chi2_range,yr=Yvsg_range,/xsty,psym=3,xtit='chi2',ytit='Y',title='Grain 3'
1763 1776  
  1777 +;stop
  1778 +Xions=la_div(Ygrain2,la_add(Ygrain1,Ygrain2))
  1779 +
  1780 +window,win,xsize=800,ysize=900 & win=win+1
  1781 +!p.multi=[0,1,4]
  1782 +G0_range=[0.1,500]
  1783 +cgplot,G0s,Ygrain1,/xlog,/ylog,xr=G0_range,/xsty,psym=3,xtit='G0',ytit='Y',title='Grain 1'
  1784 +cgplot,G0s,Ygrain2,/xlog,/ylog,xr=G0_range,/xsty,psym=3,xtit='G0',ytit='Y',title='Grain 2'
  1785 +cgplot,G0s,Ygrain3,/xlog,/ylog,xr=G0_range,/xsty,psym=3,xtit='G0',ytit='Y',title='Grain 3'
  1786 +cgplot,G0s,Xions,/xlog,/ylog,xr=G0_range,/xsty,psym=3,xtit='G0',ytit='Xion',title='Xion'
  1787 +
  1788 +window,win,xsize=800,ysize=900 & win=win+1
  1789 +Nx=100 & Ny=100
  1790 +xmin=min(NH_range)
  1791 +xmax=max(NH_range)
  1792 +ymin=min(G0_range)
  1793 +ymax=max(G0_range)
  1794 +logxbox=1
  1795 +logybox=1
  1796 +plot_bin_correl,G0s,NH_values,Nx,Ny,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,Nlevels=Nlevels, $
  1797 + fact=fact,logxbox=logxbox,logybox=logybox,range=range, $
  1798 + binned_map=binned_map,binned_xaxis=binned_xaxis,binned_yaxis=binned_yaxis,no_plot=no_plot,xtit='NH',ytit='G0', $
  1799 + /median_plot,upper_perc_plot=95.,lower_perc_plot=95.
  1800 +
  1801 +window,win,xsize=800,ysize=900 & win=win+1
  1802 +Xion_range=[0.1,0.9]
  1803 +;Nx=50 & Ny=50
  1804 +xmin=min(G0_range)
  1805 +xmax=max(G0_range)
  1806 +ymin=min(Xion_range)
  1807 +ymax=max(Xion_range)
  1808 +logxbox=1
  1809 +logybox=0
  1810 +plot_bin_correl,G0s,Xions,Nx,Ny,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,Nlevels=Nlevels, $
  1811 + fact=fact,logxbox=logxbox,logybox=logybox,range=range, $
  1812 + binned_map=binned_map,binned_xaxis=binned_xaxis,binned_yaxis=binned_yaxis,no_plot=no_plot,xtit='G0',ytit='Xion', $
  1813 + /median_plot,upper_perc_plot=95.,lower_perc_plot=95.
  1814 +
  1815 +window,win,xsize=800,ysize=900 & win=win+1
  1816 +Xion_range=[0.1,0.9]
  1817 +;Nx=50 & Ny=50
  1818 +xmin=min(NH_range)
  1819 +xmax=max(NH_range)
  1820 +ymin=min(Xion_range)
  1821 +ymax=max(Xion_range)
  1822 +logxbox=1
  1823 +logybox=0
  1824 +plot_bin_correl,NH_values,Xions,Nx,Ny,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,Nlevels=Nlevels, $
  1825 + fact=fact,logxbox=logxbox,logybox=logybox,range=range, $
  1826 + binned_map=binned_map,binned_xaxis=binned_xaxis,binned_yaxis=binned_yaxis,no_plot=no_plot,xtit='NH',ytit='Xion', $
  1827 + /median_plot,upper_perc_plot=95.,lower_perc_plot=95.
  1828 +
1764 1829 goto,the_end
1765   -stop
  1830 +;stop
1766 1831  
1767 1832 jwst_4_coutours=fltarr((size(jwst_images))[1],(size(jwst_images))[2],2)
1768 1833 jwst_4_coutours[*,*,0]=jwst_images[*,*,0,5]
... ...
LabTools/IRAP/JPB/phangs_isrf_pipeline.pro
... ... @@ -160,15 +160,7 @@ class_min_max=[15,20]
160 160 phangs_brute_force_fit_with_isrf_grid,model='DL07',source_name='ngc0628',/include_herschel,/normalize,/fit_G0,resolution_filter='PACS3' $
161 161 ,/subtract_star_light,/show,/nostop,table_name_root=table_name_root,class_min_max=class_min_max,/save
162 162  
163   -
164   -
165   -
166   -
167   -
168   -
169   -
170 163 ;ISRF classes
171   -;table_name_root='_logn_MuseISRF_JWST_adaptedG0_Ypah1added_Yvsgadded_4Phangs_noionis' ;rem: They also have YPAH0 added ....
172 164 table_name_root='_logn_MuseISRF_JWST_adaptedG0_Ypah1added_Yvsgadded_4Phangs_noionis_F0upext' ;Uses extended F0 range ....
173 165 class_min_max=0
174 166 show_sed=0
... ... @@ -582,7 +574,6 @@ make_phangs_isrf_classes,bidon,source_name='ngc0628',/save,resolution_value=1./6
582 574  
583 575 ;ngc0628_1arcsec_B1_NEAREST (1arcsec_B1_NEAREST)
584 576 table_name_root='_logn_MuseISRF_JWST_adaptedG0_Ypah1added_Yvsgadded_4Phangs_noionis_F0upext'
585   -;table_name_root='_logn_MuseISRF_JWST_adaptedG0_Ypah1added_Yvsgadded_4Phangs_noionis' ;rem: They also have YPAH0 added ....
586 577 class_min_max=0
587 578 data_dir=!phangs_data_dir+'/ISRF/WORK/'
588 579 source_name='ngc0628'
... ...
LabTools/IRAP/JPB/srun/make_phangs_grids.pro
... ... @@ -18,12 +18,14 @@ t0=systime(/sec)
18 18  
19 19 bidon=1
20 20  
  21 +;=== test
  22 +make_sed_phangs_tables,grid_type=223,isrf_class=4,/test
21 23 ;make_sed_phangs_tables,grid_type=221,isrf_class=0 ;alma1 (For Mathis field)
22 24 ;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
23 25 ;make_sed_phangs_tables,grid_type=223,isrf_class=0 ;Mathis
24 26  
25 27 ;make_sed_phangs_tables,grid_type='+2000',isrf_class=4
26   -FOR i=5L,26 DO make_sed_phangs_tables,grid_type='+2000',isrf_class=i ;This is to extend F0 by 2 twoards higher values DL07
  28 +;FOR i=5L,26 DO make_sed_phangs_tables,grid_type='+2000',isrf_class=i ;This is to extend F0 by 2 twoards higher values DL07
27 29 ;make_sed_phangs_tables,grid_type='+2000',isrf_class=0 ;alma1 (For Mathis field)
28 30  
29 31 ;make grid
... ...
src/idl/dustem_add_linear_params2grid.pro
... ... @@ -92,7 +92,7 @@ FOR ipd=0L,Npd-1 DO BEGIN
92 92 ,Ngrains,model,old_Ncomb,NFparams,INCSP,use_polar,use_double $
93 93 ,old_parameter_description,old_parameter_pmin,old_parameter_pmax,old_parameter_pnval,old_parameter_plog,old_fparameter_description,old_fparameter_value $
94 94 ,old_param_values,old_Nparams,old_seds,filters,Nfilters,wavs,Nwavs $
95   - ,old_I_tot_vec,old_ext_tot_vec,old_I_grains_vec,old_ext_grains_vec
  95 + ,old_I_tot_vec,old_ext_tot_vec,old_I_grains_vec,old_ext_grains_vec,grain_keywords=grain_keywords
96 96  
97 97 new_parameter_description=[old_parameter_description,pd[ipd]]
98 98 new_parameter_pmin=[old_parameter_pmin,iv_min[ipd]]
... ... @@ -307,7 +307,8 @@ FOR ipd=0L,Npd-1 DO BEGIN
307 307 fiv=old_fparameter_value, $
308 308 filename=use_output_table_name, $
309 309 use_polarisation=use_polarisation, $
310   - double=use_double
  310 + double=use_double, $
  311 + grain_keywords=grain_keywords
311 312 ;stop
312 313 ENDFOR
313 314  
... ...
src/idl/dustem_grids_extend_parameter_range.pro
... ... @@ -110,7 +110,6 @@ PRO dustem_grids_extend_parameter_range,input_table_name,pd_add,Nv_add $
110 110 add_iv_max[ind_add]=one_add_iv_max
111 111 add_iv_Nvalues=old_parameter_pnval
112 112 add_iv_Nvalues[ind_add]=Nv_add
113   - add_plog=old_add_parameter_plog
114 113 add_fpd=old_fparameter_description
115 114 add_fiv=old_fparameter_value
116 115 add_grain_keywords=old_grain_keywords
... ... @@ -128,44 +127,40 @@ PRO dustem_grids_extend_parameter_range,input_table_name,pd_add,Nv_add $
128 127 ;stop
129 128 add_table_name='/tmp/add_grid.fits'
130 129 ;=== decide which parameters are linear
131   - ;IF param_is_non_linear THEN BEGIN ;This is if the parameter to be added is non-linear
132   - ;== make table for non-linear parameters
133   - message,'========================== Making added grid for non-linear parameters ...',/continue
134   - ;In priciple, if non-linear parameters are already in the input table, there should be no reason to do that, but anyway ...
135   - dustem_make_sed_table,add_model $
136   - ,add_parameters_description[ind_non_linear] $
137   - ,add_iv_min[ind_non_linear] $
138   - ,add_iv_max[ind_non_linear] $
139   - ,add_iv_Nvalues[ind_non_linear] $
140   - ,fpd=add_fpd $
141   - ,fiv=add_fiv $
142   - ,grain_keywords=add_grain_keywords $
143   - ,filename=add_table_name $
144   - ,filters=add_filters $
145   - ,plog=add_plog[ind_non_linear] $
146   - ,show_seds=show_seds $
147   - ,print_params=print_params $
148   - ,help=help
149   - ;== extend table for linear parameter
150   - message,'========================== Making added grid for linear parameters ...',/continue
151   - dustem_add_linear_params2grid,add_table_name $
152   - ,add_parameters_description[ind_linear] $
153   - ,add_iv_min[ind_linear] $
154   - ,add_iv_max[ind_linear] $
155   - ,add_iv_Nvalues[ind_linear] $
156   - ,out_filename=add_table_name $
157   - ,plog=add_plog[ind_linear] $
158   - ,show_seds=show_seds
159   - ;ENDIF ELSE BEGIN
160   -
161   - ;ENDELSE
  130 + ;== make table for non-linear parameters
  131 + message,'========================== Making added grid for non-linear parameters ...',/continue
  132 + ;In priciple, if non-linear parameters are already in the input table, there should be no reason to do that, but anyway ...
  133 + dustem_make_sed_table,add_model $
  134 + ,add_parameters_description[ind_non_linear] $
  135 + ,add_iv_min[ind_non_linear] $
  136 + ,add_iv_max[ind_non_linear] $
  137 + ,add_iv_Nvalues[ind_non_linear] $
  138 + ,fpd=add_fpd $
  139 + ,fiv=add_fiv $
  140 + ,grain_keywords=add_grain_keywords $
  141 + ,filename=add_table_name $
  142 + ,filters=add_filters $
  143 + ,plog=add_plog[ind_non_linear] $
  144 + ,show_seds=show_seds $
  145 + ,print_params=print_params $
  146 + ,help=help
  147 + ;== extend table for linear parameters
  148 + message,'========================== Making added grid for linear parameters ...',/continue
  149 + dustem_add_linear_params2grid,add_table_name $
  150 + ,add_parameters_description[ind_linear] $
  151 + ,add_iv_min[ind_linear] $
  152 + ,add_iv_max[ind_linear] $
  153 + ,add_iv_Nvalues[ind_linear] $
  154 + ,out_filename=add_table_name $
  155 + ,plog=add_plog[ind_linear] $
  156 + ,show_seds=show_seds
162 157  
163 158 ;=== merge tables
164 159 use_output_table_name='/tmp/add_ext_grid.fits'
165 160 IF keyword_set(output_table_name) THEN use_output_table_name=output_table_name
166 161  
167 162 message,'========================== Merging grids ...',/continue
168   - dustem_merge_grids,input_table_name,add_table_name,output_table_name
  163 + dustem_merge_grids,input_table_name,add_table_name,use_output_table_name
169 164  
170 165 ;stop
171 166  
... ...
src/idl/dustem_make_grid_description_table.pro
1   -PRO dustem_make_grid_description_table,table_name
  1 +PRO dustem_make_grid_description_table,table_name,nolatex=nolatex
2 2  
3 3 ;dir=!phangs_data_dir+'/ISRF/GRIDS/'
4 4 ;isrf_class_str='_isrfclass'+strtrim(isrf_class,2)
... ... @@ -6,6 +6,9 @@ PRO dustem_make_grid_description_table,table_name
6 6 ;table_name=dir+model+'_MuseISRF_JWST_G0_YPAH_YVSG_4Phangs_noionis'+isrf_class_str+'.fits'
7 7 ;dustem_make_grid_description_table,
8 8  
  9 +;table_name='/Volumes/PILOT_FLIGHT1/PHANGS/ISRF/GRIDS/DL07_logn_MuseISRF_JWST_adaptedG0_Ypah1added_Yvsgadded_4Phangs_noionis_F0upext_isrfclass8.fits'
  10 +;dustem_make_grid_description_table,table_name
  11 +
9 12 dustem_define_grid,table_name,filters=filters,status=status
10 13  
11 14 Npar=(*!dustem_grid).NPARAMS
... ... @@ -20,8 +23,25 @@ st.plog=(*!dustem_grid).PLOG
20 23  
21 24 write_xcat,st,'/tmp/bidon.xcat'
22 25  
23   -struct2latex_table,st,'/tmp/bidon.tex'
24   -
  26 +IF not keyword_set(nolatex) THEN BEGIN
  27 + struct2latex_table,st,'/tmp/bidon.tex'
  28 +ENDIF
  29 +
  30 +;stop
  31 +;===== check parameter values
  32 +;stop
  33 +
  34 +;count parameters uniq values
  35 +n_diffs=lonarr(Npar)
  36 +FOR i=0L,Npar-1 DO BEGIN
  37 + vv=(*!dustem_grid).st_params.(3)
  38 + order=sort(vv)
  39 + vv=vv[order]
  40 + un=uniq(vv)
  41 + vv=vv[un]
  42 + n_diffs[i]=n_elements(vv)
  43 +ENDFOR
  44 +print,n_diffs
25 45 stop
26 46  
27 47 END
28 48 \ No newline at end of file
... ...
src/idl/dustem_marginalize_grid.pro
... ... @@ -80,6 +80,8 @@ FOR i=0L,Nfilters-1 DO BEGIN
80 80 ENDFOR
81 81 ;Nfilters=n_elements(filters)
82 82  
  83 +use_dist_min_perc=1. ;parameter values which are a this percent different will be considered equal in marginalization ...
  84 +
83 85 CASE use_method OF
84 86 'NEAREST': BEGIN
85 87 FOR i=0L,N_fixed_params-1 DO BEGIN
... ... @@ -87,11 +89,20 @@ CASE use_method OF
87 89 IF count NE 0 THEN BEGIN
88 90 ;fixed_parameter_value=fixed_parameters_values[ind[0]]
89 91 fixed_parameter_value=fixed_parameters_values[i]
90   - ;compute distance between table parameter and the requested fixed value
91   - distance=abs(fixed_parameter_value-params[*].(ind[0]))
  92 + ;compute absolute distance between table parameter and the requested fixed value
  93 + this_param_values=params[*].(ind[0])
  94 + ;=== This is a trick to make equal parameter values which are different below a percent level.
  95 + ;=== Needed to allow for slightly different values of a given parameter in merged tables (at numerical accuracy)
  96 + this_param_values=1.*round(this_param_values/this_param_values[0]*use_dist_min_perc*100.)/100.*this_param_values[0]
  97 + ;order=sort(this_param_values)
  98 + ;this_param_indiv_values=this_param_values[order]
  99 + ;un=uniq(this_param_indiv_values)
  100 + ;this_param_indiv_values=this_param_indiv_values[un]
  101 + ;print,this_param_indiv_values
  102 + ;compute distance between provided parameter and parameter values
  103 + distance=abs(fixed_parameter_value-this_param_values)
  104 + ;pick the smalest distance
92 105 indd=where(distance EQ min(distance),ccount)
93   - ;print,ccount
94   - ;2197
95 106 ;stop
96 107 IF ccount EQ 0 THEN BEGIN
97 108 message,'This should not happen',/continue
... ...
src/idl/dustem_write_sed_grid.pro
... ... @@ -280,7 +280,11 @@ IF keyword_set(grain_keywords) THEN BEGIN
280 280 FOR i=0L,Ngrains-1 DO BEGIN
281 281 sxaddpar,h0,'KEYWG'+strtrim(i+1,2),grain_keywords[i],'dustemwrap grain keywords for grain '+strtrim(i+1,2),before='END'
282 282 ENDFOR
283   -ENDIF
  283 +ENDIF ELSE BEGIN
  284 + FOR i=0L,Ngrains-1 DO BEGIN
  285 + sxaddpar,h0,'KEYWG'+strtrim(k+1,2),((*!dustem_params).grains.TYPE_KEYWORDS)[k],'Grain DustemWrap keywords Used',before='END'
  286 + ENDFOR
  287 +ENDELSE
284 288 sxaddpar,h0,'NHREF',*!dustem_HCD,'dustemwrap reference NH [H/cm2]',after='MODEL'
285 289 sxaddpar,h0,'POLAR',polar,'dustemwrap polarization ?',after='NHREF'
286 290 sxaddpar,h0,'DOUBLE',double4header,'double accuracy ?',after='POLAR'
... ... @@ -288,9 +292,6 @@ sxaddpar,h0,'COMMENT','**** Variable Dustemwrap PARAMETERS Used in the Grid ***'
288 292 sxaddpar,h0,'NPARAM',Nparams,'Number of free parameters',before='END'
289 293 sxaddpar,h0,'NCOMB',Nc,'Number of parameter values combinations',before='END'
290 294 sxaddpar,h0,'KEY',(*!dustem_params).keywords,'Global DustemWrap keywords Used',before='END'
291   -FOR k=0L,Ngrains-1 DO BEGIN
292   - sxaddpar,h0,'KEYG'+strtrim(k,2),((*!dustem_params).grains.TYPE_KEYWORDS)[k],'Grain DustemWrap keywords Used',before='END'
293   -ENDFOR
294 295 ;stop
295 296 FOR j=0L,Nparams-1 DO BEGIN
296 297 ;value=sxpar(*all_headers[0],'TTYPE'+strtrim(j+1,2))
... ...