OV-GSO-DC / dustem-wrapper_idl

Commit 0b5fce89a201ed72950a686f8f0b4829ae792ee3

Authored by Jean-Philippe Bernard
1 parent 46980f16
Exists in master

added provision for a test on Yvsg random variations

Showing 3 changed files with 245 additions and 265 deletions   Show diff stats
Inline Side-by-side
LabTools/IRAP/JPB/phangs_brute_force_fit_with_isrf_grid.pro
  Diff comments   View file @0b5fce8
... ... @@ -24,6 +24,7 @@ PRO phangs_brute_force_fit_with_isrf_grid,source_name=source_name $
24 24 ,marginalize_method=marginalize_method $
25 25 ,do_postscripts=do_postscripts $
26 26 ,postscripts_dir=postscripts_dir $
  27 + ,perturbe_yvsg=perturbe_yvsg $
27 28 ,silent=silent
28 29  
29 30 ;+
... ... @@ -93,6 +94,11 @@ IF keyword_set(set_yvsg) and not keyword_set(fit_G0) THEN BEGIN
93 94 stop
94 95 ENDIF
95 96  
  97 +perturbed_yvsg_str=''
  98 +IF keyword_set(perturbe_yvsg) THEN BEGIN
  99 + perturbed_yvsg_str='perturbed_yvsg'
  100 +ENDIF
  101 +
96 102 dustem_define_la_common
97 103 dustem_init
98 104  
... ... @@ -641,6 +647,13 @@ ENDIF
641 647 IF use_yvsg_option EQ 'MAP' THEN BEGIN
642 648 message,'Making map of Yvsg fixed values',/continue
643 649 new_map=project2(YVSG_fits_header,YVSG_fits_map,href,/round)
  650 + ;add a random map to Yvsg map, from 0 to perturbe_yvsg*Yvsg
  651 + IF keyword_set(perturbe_yvsg) THEN BEGIN
  652 + pert=perturbe_yvsg*randomu(seed,(size(new_map))[1],(size(new_map))[2]) ;map with min max from 0 to perturbe_yvsg
  653 + ;pert=la_mul(Yg3_map,pert)
  654 + ;Yg3_map=la_add(Yg3_map,pert)
  655 + new_map=la_add(new_map,la_mul(new_map,pert))
  656 + ENDIF
644 657 yvsgs_from_map=fltarr(Nvor)
645 658 FOR vid=0L,Nvor -1 DO BEGIN
646 659 yvsgs_from_map[vid]=la_mean(new_map[*all_seds_indices[vid]])
... ... @@ -967,12 +980,12 @@ ENDFOR
967 980 stats_st=st
968 981 IF keyword_set(save) THEN BEGIN
969 982 dustem_grid_saved=!dustem_grid
970   - 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'
  983 + 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+perturbed_yvsg_str+'.sav'
971 984 save,all_filters,use_fit_filters,dustem_grid_saved,stellar_normalizations,stats_st,hit_st,hit_st_perc,EBVs,EBVs_uv,tausV,taus_UV,F0s,F0s_predicted,G0s,G0s_predicted, $
972 985 Ygrain1,Ygrain2,Ygrain3,facts,chi2s,rchi2s,dF0s,dYgrain1,dYgrain2,dYgrain3, $
973 986 F0_hit,Ygrain1_hit,Ygrain2_hit,Ygrain3_hit,NH_values,best_fit_seds,observed_seds,file=file_save
974 987 message,'Wrote '+file_save,/continue
975   - file_save=data_dir+'seds_'+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'
  988 + file_save=data_dir+'seds_'+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+perturbed_yvsg_str+'.sav'
976 989 save,all_wavs,fit_wavs,saved_seds_st,file=file_save
977 990 message,'Wrote '+file_save,/continue
978 991 ;stop
... ... @@ -1058,7 +1071,7 @@ FOR vid=0LL,Nvor-1 DO BEGIN
1058 1071 ENDFOR
1059 1072  
1060 1073 IF keyword_set(save) THEN BEGIN
1061   - file_save=data_dir+use_source_name+use_model+'_logn_JWST_G0_YPAH_YVSG_maps_with-ISRFclasses'+norm_str+mathis_str+reso_str+fit_G0_str+set_yvsg_str+'.sav'
  1074 + file_save=data_dir+use_source_name+use_model+'_logn_JWST_G0_YPAH_YVSG_maps_with-ISRFclasses'+norm_str+mathis_str+reso_str+fit_G0_str+set_yvsg_str+perturbed_yvsg_str+'.sav'
1062 1075 save,NHCO,NHI,NHmap,norm_stars_map,EBV_map,EBV_uv_map,tauV_map,tau_uv_map,N_ebv_uv_filters,wwav_star,wwav_uv,F0_map,G0_map,F0_predicted_map,G0_predicted_map,Yg1_map,Yg2_map,Yg3_map,fact_map,chi2_map,rchi2_map,dF0_map,dYg1_map,dYg2_map,dYg3_map,F0_hit_map,Yg1_hit_map,Yg2_hit_map,Yg3_hit_map,NH_values_map,residual_maps,residual_maps_perc,file=file_save
1063 1076 message,'Wrote '+file_save,/continue
1064 1077 ENDIF
... ...
LabTools/IRAP/JPB/phangs_isrf_pipeline.pro
  Diff comments   View file @0b5fce8
... ... @@ -623,9 +623,40 @@ phangs_make_figures_from_restore,model='DL07',source_name='ngc0628',/normalize,f
623 623  
624 624  
625 625  
  626 +;This is an attempt at explaining increased PAH abundances by a factor of 2 when going from PACS3 to 1arcsec resolution
  627 +;by seeing the effect of a change of Yvsg from the low to the high resolution, by a random factor of 2
  628 +;ngc0628_1arcsec_B1_NEAREST (1arcsec_B1_NEAREST)
  629 +define_la_common
  630 +table_name_root='_logn_MuseISRF_JWST_adaptedG0_Ypah1added_Yvsgadded_4Phangs_noionis_F0upext'
  631 +class_min_max=0
  632 +data_dir=!phangs_data_dir+'/ISRF/WORK/'
  633 +source_name='ngc0628'
  634 +model='DL07'
  635 +norm_str='_normalized'
  636 +mathis_str=''
  637 +reso_str='_PACS3'
  638 +fit_G0_str=''
  639 +set_yvsg_str=''
  640 +file_save=data_dir+source_name+model+'_JWST_G0_YPAH_YVSG_maps_with-ISRFclasses'+norm_str+mathis_str+reso_str+fit_G0_str+set_yvsg_str+'.sav'
  641 +restore,file_save,/verb
  642 +;=== here modifying YG3 map
  643 +perturbe_yvsg=-0.5 ;This will decrease Yvsg map by a random fraction from 0% to 50%
  644 +restore,data_dir+source_name+'_ref_header'+reso_str+'.sav',/verb
  645 +nostop=1
  646 +show_sed=0
  647 +marginalize_method='NEAREST'
  648 +phangs_brute_force_fit_with_isrf_grid,model='DL07',source_name='ngc0628',/include_herschel,/normalize,fit_G0=1,set_yvsg=1,/save,resolution_value=1./60.^2 $
  649 + ,/subtract_star_light,show_sed=show_sed,/nostop,table_name_root=table_name_root,class_min_max=class_min_max $
  650 + ,yvsg_option='MAP',yvsg_fits_map=Yg3_map,yvsg_fits_header=href,marginalize_method=marginalize_method,/silent,perturbe_yvsg=perturbe_yvsg
626 651  
627 652  
628   -
  653 +;doing postscripts
  654 +table_name_root='_logn_MuseISRF_JWST_adaptedG0_Ypah1added_Yvsgadded_4Phangs_noionis_F0upext'
  655 +postscripts_dir=!phangs_data_dir+'/ISRF/WORK/PS/ngc0628_1arcsec_B1_NEAREST'
  656 +nostop=1
  657 +phangs_make_figures_from_restore,model='DL07',source_name='ngc0628',/normalize,fit_G0=1,set_yvsg=1,resolution_value=1./60.^2 $
  658 + ,nostop=nostop,table_name_root=table_name_root $
  659 + ,/do_postscripts,postscripts_dir=postscripts_dir,which_contours='ENVNH',/zoom_maps,/perturbe_yvsg
629 660  
630 661  
631 662  
... ...
LabTools/IRAP/JPB/phangs_make_figures_from_restore.pro
  Diff comments   View file @0b5fce8
... ... @@ -16,10 +16,10 @@ PRO phangs_make_figures_from_restore,source_name=source_name $
16 16 ,plot_these_seds=plot_these_seds $
17 17 ,which_contours=which_contours $
18 18 ,zoom_maps=zoom_maps $
  19 + ,perturbe_yvsg=perturbe_yvsg $
19 20 ,silent=silent $
20 21 ,help=help
21 22  
22   -
23 23 ;+
24 24 ; NAME:
25 25 ; phangs_make_figures_from_restore
... ... @@ -79,11 +79,14 @@ IF keyword_set(help) THEN BEGIN
79 79 goto,the_end
80 80 ENDIF
81 81  
82   -obp=[1.1,0.,1.15,1]
83 82 win=0L
84 83 wxs=800
85 84 wys=900
86 85  
  86 +obp=[1.1,0.,1.15,1]
  87 +cleanplot
  88 +image_color_table='jpbloadct,/whitebackground'
  89 +
87 90 ;weightning='1/chi2'
88 91 ;weightning='likelyhood'
89 92 weightning='rellikelyhood' ;weightning strategy used for weighted parameters
... ... @@ -105,6 +108,11 @@ IF keyword_set(set_yvsg) and not keyword_set(fit_G0) THEN BEGIN
105 108 stop
106 109 ENDIF
107 110  
  111 +perturbed_yvsg_str=''
  112 +IF keyword_set(perturbe_yvsg) THEN BEGIN
  113 + perturbed_yvsg_str='perturbed_yvsg'
  114 +ENDIF
  115 +
108 116 dustem_define_la_common
109 117 dustem_init
110 118  
... ... @@ -157,7 +165,7 @@ ENDELSE
157 165  
158 166 ;stop
159 167  
160   -file=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'
  168 +file=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+perturbed_yvsg_str+'.sav'
161 169 st_info=file_info(file)
162 170 IF st_info.exists NE 1 THEN BEGIN
163 171 message,'Could not find '+file,/continue
... ... @@ -195,7 +203,7 @@ restore,file,/verb
195 203 ; % RESTORE: Restored variable: NH_VALUES.
196 204 ; % RESTORE: Restored variable: BEST_FIT_SEDS.
197 205 ; % RESTORE: Restored variable: OBSERVED_SEDS.
198   -file=data_dir+use_source_name+use_model+'_logn_JWST_G0_YPAH_YVSG_maps_with-ISRFclasses'+norm_str+mathis_str+reso_str+fit_G0_str+set_yvsg_str+'.sav'
  206 +file=data_dir+use_source_name+use_model+'_logn_JWST_G0_YPAH_YVSG_maps_with-ISRFclasses'+norm_str+mathis_str+reso_str+fit_G0_str+set_yvsg_str+perturbed_yvsg_str+'.sav'
199 207 st_info=file_info(file)
200 208 IF st_info.exists NE 1 THEN BEGIN
201 209 message,'Could not find '+file,/continue
... ... @@ -354,8 +362,17 @@ IF keyword_set(zoom_maps) THEN BEGIN
354 362 Yg2_hit_map=project2(href,Yg2_hit_map,hout,/silent,/previous_lut)
355 363 Yg3_hit_map=project2(href,Yg3_hit_map,hout,/silent,/previous_lut)
356 364 NH_values_map=project2(href,NH_values_map,hout,/silent,/previous_lut)
  365 + Nresiduals=(size(residual_maps))[3]
  366 + new_residual_maps=fltarr(sxpar(hout,'NAXIS1'),sxpar(hout,'NAXIS2'),Nresiduals)
  367 + new_residual_maps_perc=fltarr(sxpar(hout,'NAXIS1'),sxpar(hout,'NAXIS2'),Nresiduals)
  368 + FOR i=0L,Nresiduals-1 DO BEGIN
  369 + new_residual_maps[*,*,i]=project2(href,residual_maps[*,*,i],hout,/silent,/previous_lut)
  370 + new_residual_maps_perc[*,*,i]=project2(href,residual_maps_perc[*,*,i],hout,/silent,/previous_lut)
  371 + ENDFOR
  372 + residual_maps=new_residual_maps
  373 + residual_maps_perc=new_residual_maps_perc
357 374 href=hout
358   - stop
  375 + ;stop
359 376 ENDIF
360 377  
361 378 Xions=la_div(Ygrain2,la_add(Ygrain1,Ygrain2))
... ... @@ -500,8 +517,10 @@ IF not keyword_set(no_MC) THEN BEGIN
500 517 dir=!phangs_data_dir+'/phangs_drive/Archive/Products/cloud_catalogs/team_release/by_galaxy/'+use_source_name+'/'
501 518 file=dir+use_source_name+'_12m+7m+tp_co21_native_binmask2D.fits.bz2'
502 519 d=mrdfits(file,0,hmcs)
  520 + ind=where(d EQ 0)
  521 + d[ind]=-1
503 522 mcs=project2(hmcs,d,href)
504   - ;image_cont20,mcs,href,/square,/silent,image_color_table=image_color_table,axis_color_table=1,off_bar=obp,levels=use_level_st
  523 + image_cont20,mcs,href,/square,/silent,image_color_table=image_color_table,axis_color_table=1,off_bar=obp,levels=use_level_st
505 524 ;stop
506 525 is_mc=mcs
507 526 ind=where(mcs NE 0 and mcs NE la_undef(),count)
... ... @@ -530,6 +549,10 @@ IF not keyword_set(no_HII) THEN BEGIN
530 549 sxaddpar,hhii,'EQUINOX',2000.
531 550 ;image_cont20,hii,hhii,/square,/silent
532 551 hii=project2(hhii,hii,href,/silent,/round)
  552 + is_hii=hii
  553 + ind=where(hii NE -1 and hii NE la_undef(),count)
  554 + is_hii[ind]=1
  555 + is_hii[ind_bad]=la_undef()
533 556 images=fltarr(Nx,Ny,2)
534 557 images[*,*,0]=hii
535 558 images[*,*,1]=hii
... ... @@ -622,263 +645,156 @@ IF keyword_set(plot_these_seds) THEN BEGIN
622 645 ENDFOR
623 646 ENDIF
624 647  
625   -;=== compute statistics on HII regions
626   -IF not keyword_set(no_HII) THEN BEGIN
627   - ;window,win,xsize=wxs,ysize=wys & win=win+1
628   - ;image_cont20,hii,href,/square,/silent,image_color_table='jpbloadct,/white',axis_color_table=1,off_bar=obp,levels=use_level_st
629   - ;=== do a simple HII region mask (1 if an HII, 0 otherwise) and a contour structure
630   - is_hii=hii
631   - ind=where(hii NE -1 and hii NE la_undef(),count)
632   - is_hii[ind]=1
633   - is_hii[ind_bad]=la_undef()
634   - ;stop
635   - hii_m=fltarr(Nx,Ny,2)
636   - hii_m[*,*,0]=is_hii
637   - hii_m[*,*,1]=is_hii
638   - levs_HII=[0.5,1]
639   - hii_st=create_contour_st(hii_m,lev1=levs_HII,lev2=levs_HII)
640   - hii_st.(0).color=0
641   - hii_st.(0).thick=2.
642   - hii_st.(1).color=100
643   - hii_st.(1).thick=1.0
644   - ;hii[*]=la_undef()
645   - ;image_cont20,hii,href,/square,/silent,levels=hii_st,image_color_table='jpbloadct,/no',axis_color_table=1,off_bar=obp
646   - ;image_cont20,is_hii,href,/square,/silent,levels=hii_st,image_color_table='jpbloadct',axis_color_table=1
647   - ;stop
648   - ;=== do a mask per HII region and a surounding mask
649   - N_hii=long(max(hii))
650   - hii_indices=ptrarr(N_hii)
651   - hii_count=lonarr(N_hii)
652   - hii_numbers=lonarr(N_hii)
653   - FOR i=0L,N_hii-1 DO BEGIN
654   - ind_hii=where(hii EQ i,count)
655   - IF count NE 0 THEN BEGIN
656   - hii_indices[i]=ptr_new(ind_hii)
657   - ENDIF
658   - hii_count[i]=count
659   - hii_numbers[i]=i+1
660   - ENDFOR
661   - ;=== remove unfound HII regions
662   - ind=where(hii_count NE 0,count)
663   - IF count NE 0 THEN BEGIN
664   - hii_count=hii_count[ind]
665   - hii_indices=hii_indices[ind]
666   - hii_numbers=hii_numbers[ind] ;Original HII region number
667   - N_hii=n_elements(hii_indices) ;new number of HII regions
668   - ENDIF
669   - ;=== order HII regions by size
670   - order=sort(hii_count)
671   - hii_count=hii_count[order]
672   - hii_indices=hii_indices[order]
673   - hii_numbers=hii_numbers[order]
674   - ;=== do a surounding mask per HII region
675   - around_hii_indices=ptrarr(N_hii)
676   - Npix=5 ;number of pixels by which to extend
677   - count_env=lonarr(N_hii)
678   - maskk=intarr(Nx,Ny)
679   - ind=where(is_hii EQ 1)
680   - maskk[ind]=1
681   - FOR i=0L,N_hii-1 DO BEGIN
682   - indd=enlarge(*hii_indices[i],Npix,Nx,Ny,/remove_seed,remove_these=maskk,count=count) ;This enlarges but does remove coliding hii regions.
683   - IF count NE 0 THEN around_hii_indices[i]=ptr_new(indd)
684   - count_env[i]=count
685   - ENDFOR
686   - ;stop
  648 +;== This is for local statistics around regions (HII, MCs, nebula)
  649 +one_stat={region:'',med_NH:la_undef(4),std_NH:la_undef(4),med_G0:la_undef(4),std_G0:la_undef(4),med_YG1:la_undef(4),med_YG2:la_undef(4),med_YG3:la_undef(4),med_Xion:la_undef(4),med_rchi2:la_undef(4)}
687 650  
688   - ;=== compute global statistics of fitted parameters in various regions
689   - one_stat={region:'',count:0L,med_G0:la_undef(4),med_YG1:la_undef(4),med_YG2:la_undef(4),med_YG3:la_undef(4),med_Xion:la_undef(4),med_rchi2:la_undef(4)}
690   - regions=['ALL','HII','not-HII']
691   - Nregions=n_elements(regions)
692   - regional_stats=replicate(one_stat,Nregions)
693   - regional_stats.region=regions
694   - FOR i=0L,Nregions-1 DO BEGIN
695   - CASE regions[i] OF
696   - 'ALL': BEGIN
697   - ind=where(G0_map NE la_undef(),count)
698   - END
699   - 'HII': BEGIN
700   - ind=where(G0_map NE la_undef() AND is_hii EQ 1 AND is_hii NE la_undef(),count)
701   - END
702   - 'not-HII': BEGIN
703   - ind=where(G0_map NE la_undef() AND is_hii EQ 0 AND is_hii NE la_undef(),count)
704   - END
705   - ENDCASE
706   - regional_stats[i].med_G0=la_median(G0_map[ind])
707   - regional_stats[i].count=count
708   - regional_stats[i].med_YG1=la_median(Yg1_map[ind])
709   - regional_stats[i].med_YG2=la_median(Yg2_map[ind])
710   - regional_stats[i].med_YG3=la_median(Yg3_map[ind])
711   - regional_stats[i].med_Xion=la_median(Xion_map[ind])
712   - regional_stats[i].med_rchi2=la_median(rchi2_map[ind])
713   - ENDFOR
  651 +images=fltarr(Nx,Ny,7)
  652 +images[*,*,0]=NH_values_map
  653 +images[*,*,1]=G0_map
  654 +images[*,*,2]=Yg1_map
  655 +images[*,*,3]=Yg2_map
  656 +images[*,*,4]=Yg3_map
  657 +images[*,*,5]=Xion_map
  658 +images[*,*,6]=rchi2_map
714 659  
715   - write_xcat,regional_stats,'/tmp/bidon',/nowrite
716   - ;stop
  660 +struct_tags=['NH','G0','YG1','YG2','YG3','Xion','rchi2'] ;must correspond to images above
717 661  
718   - ;=== compute local statistics of fitted parameters locally around HII regions
719   - one_stat={region:'',med_NH:la_undef(4),std_NH:la_undef(4),med_G0:la_undef(4),std_G0:la_undef(4),med_YG1:la_undef(4),med_YG2:la_undef(4),med_YG3:la_undef(4),med_Xion:la_undef(4),med_rchi2:la_undef(4)}
720   - regions=['HII']
721   - Nregions=n_elements(regions)
722   - hii_stats=replicate(one_stat,Nregions)
723   - hii_stats.region=regions
724   - med_NH_in=fltarr(N_hii) & med_NH_out=fltarr(N_hii)
725   - med_G0_in=fltarr(N_hii) & med_G0_out=fltarr(N_hii)
726   - med_YG1_in=fltarr(N_hii) & med_YG1_out=fltarr(N_hii)
727   - med_YG2_in=fltarr(N_hii) & med_YG2_out=fltarr(N_hii)
728   - med_YG3_in=fltarr(N_hii) & med_YG3_out=fltarr(N_hii)
729   - med_Xion_in=fltarr(N_hii) & med_Xion_out=fltarr(N_hii)
730   - med_rchi2_in=fltarr(N_hii) & med_rchi2_out=fltarr(N_hii)
731   - FOR i=0L,N_hii-1 DO BEGIN
732   - med_NH_in[i]=la_median(NH_values_map[*hii_indices[i]])
733   - med_NH_out[i]=la_median(NH_values_map[*around_hii_indices[i]])
734   - med_G0_in[i]=la_median(G0_map[*hii_indices[i]])
735   - med_G0_out[i]=la_median(G0_map[*around_hii_indices[i]])
736   - med_YG1_in[i]=la_median(YG1_map[*hii_indices[i]])
737   - med_YG1_out[i]=la_median(YG1_map[*around_hii_indices[i]])
738   - med_YG2_in[i]=la_median(YG2_map[*hii_indices[i]])
739   - med_YG2_out[i]=la_median(YG2_map[*around_hii_indices[i]])
740   - med_YG3_in[i]=la_median(YG3_map[*hii_indices[i]])
741   - med_YG3_out[i]=la_median(YG3_map[*around_hii_indices[i]])
742   - med_Xion_in[i]=la_median(Xion_map[*hii_indices[i]])
743   - med_Xion_out[i]=la_median(Xion_map[*around_hii_indices[i]])
744   - med_rchi2_in[i]=la_median(rchi2_map[*hii_indices[i]])
745   - med_rchi2_out[i]=la_median(rchi2_map[*around_hii_indices[i]])
746   - ENDFOR
747   - ;=== compute average HII statistics per HII region size
748   - ind=where(med_NH_in NE la_undef() and med_NH_out NE la_undef(),count)
749   - hii_count=hii_count[ind]
750   - med_NH_in=med_NH_in[ind] & med_NH_out=med_NH_out[ind]
751   - med_G0_in=med_G0_in[ind] & med_G0_out=med_G0_out[ind]
752   - med_YG1_in=med_YG1_in[ind] & med_YG1_out=med_YG1_out[ind]
753   - med_YG2_in=med_YG2_in[ind] & med_YG2_out=med_YG2_out[ind]
754   - med_YG3_in=med_YG3_in[ind] & med_YG3_out=med_YG3_out[ind]
755   - med_Xion_in=med_Xion_in[ind] & med_Xion_out=med_Xion_out[ind]
756   - med_rchi2_in=med_rchi2_in[ind] & med_rchi2_out=med_rchi2_out[ind]
757   -
758   - order=sort(hii_count)
759   - toto=hii_count[order]
760   - un=uniq(toto)
761   - un_sizes=toto[un]
762   - Nun_sizes=n_elements(un_sizes)
763   - one_stat={size:0.,NHII:0L, $
764   - med_d_NH:la_undef(4),std_d_NH:la_undef(4), $
765   - med_d_G0:la_undef(4),std_d_G0:la_undef(4), $
766   - med_d_YG1:la_undef(4),std_d_YG1:la_undef(4), $
767   - med_d_YG2:la_undef(4),std_d_YG2:la_undef(4), $
768   - med_d_YG3:la_undef(4),std_d_YG3:la_undef(4), $
769   - med_d_Xion:la_undef(4),std_d_Xion:la_undef(4), $
770   - med_d_rchi2:la_undef(4),std_d_rchi2:la_undef(4)}
771   - size_hii_stat=replicate(one_stat,Nun_sizes)
772   - FOR i=0L,Nun_sizes-1 DO BEGIN
773   - ind=where(hii_count EQ un_sizes[i],count)
774   - size_hii_stat[i].size=la_median(hii_count[ind])
775   - size_hii_stat[i].NHII=count
776   - vv=(med_NH_in[ind]-med_NH_out[ind])/med_NH_out[ind]*100.
777   - size_hii_stat[i].med_d_NH=la_median(vv) & size_hii_stat[i].std_d_NH=la_sigma(vv)
778   - vv=(med_G0_in[ind]-med_G0_out[ind])/med_G0_out[ind]*100.
779   - size_hii_stat[i].med_d_G0=la_median(vv) & size_hii_stat[i].std_d_G0=la_sigma(vv)
780   - vv=(med_YG1_in[ind]-med_YG1_out[ind])/med_YG1_out[ind]*100.
781   - size_hii_stat[i].med_d_YG1=la_median(vv) & size_hii_stat[i].std_d_YG1=la_sigma(vv)
782   - vv=(med_YG2_in[ind]-med_YG2_out[ind])/med_YG2_out[ind]*100.
783   - size_hii_stat[i].med_d_YG2=la_median(vv) & size_hii_stat[i].std_d_YG2=la_sigma(vv)
784   - vv=(med_YG3_in[ind]-med_YG3_out[ind])/med_YG3_out[ind]*100.
785   - size_hii_stat[i].med_d_YG3=la_median(vv) & size_hii_stat[i].std_d_YG3=la_sigma(vv)
786   - vv=(med_Xion_in[ind]-med_Xion_out[ind])/med_Xion_out[ind]*100.
787   - size_hii_stat[i].med_d_Xion=la_median(vv) & size_hii_stat[i].std_d_Xion=la_sigma(vv)
788   - vv=(med_rchi2_in[ind]-med_rchi2_out[ind])/med_rchi2_out[ind]*100.
789   - size_hii_stat[i].med_d_rchi2=la_median(vv) & size_hii_stat[i].std_d_rchi2=la_sigma(vv)
790   - ENDFOR
  662 +;======================================================================================================
  663 +;=== compute local statistics of fitted parameters locally around HII regions
  664 +;======================================================================================================
  665 +
  666 +hii_local_stats=regional_stats(hii,images,/median,/do_print $
  667 + ,struct_tags=struct_tags,stdev_obj_st=hii_median_stdev_obj_st,stdev_back_st=hii_median_stdev_back_st $
  668 + ,order_by='OBJECT_COUNT',local_obj_st=hii_median_local_stats,per_count_bin_local_st=hii_per_count_bin_local_st)
791 669  
792   - write_xcat,size_hii_stat,/nowrite
793   - ;stop
794   -ENDIF
795 670 IF not keyword_set(nostop) THEN stop
796 671  
  672 +;stop
  673 +
797 674 ;==========================================================================================
798 675 ;======= plot statistics on HII regions
799 676 ;==========================================================================================
800 677 IF not keyword_set(no_HII) THEN BEGIN
801 678 window,win,xsize=wxs*2,ysize=wys & win=win+1
802 679 !p.multi=[0,3,2]
803   - xrange=[1,max(hii_count)]
  680 + xrange=[1,max(hii_median_local_stats.count)]
  681 + xtit='Surface [pix]'
804 682  
805   - ind=where(med_NH_in NE la_undef() and med_NH_out NE la_undef(),count)
806   - vv=(med_NH_in[ind]-med_NH_out[ind])/med_NH_out[ind]*100.
807   - cgplot,hii_count,vv,tit='NH',xrange=xrange,/xlog
  683 + ;ind=where(med_NH_in NE la_undef() and med_NH_out NE la_undef(),count)
  684 + ;vv=(med_NH_in[ind]-med_NH_out[ind])/med_NH_out[ind]*100.
  685 + cgplot,hii_median_local_stats.count,hii_median_local_stats.NH,tit='Delta-NH [%]',xrange=xrange,/xlog,xtit=xtit
808 686 cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
809   - cgoplot,size_hii_stat.size,size_hii_stat.med_d_NH,color='blue',psym='filled circle'
810   - hii_stats.med_NH=la_median(size_hii_stat.med_d_NH)
811   - cgoplot,10^!x.crange,replicate(hii_stats.med_NH,2),linestyle=2,color='blue'
812   - hii_stats.std_NH=hii_stats.med_NH/la_sigma(vv)
813   -
814   - ind=where(med_G0_in NE la_undef() and med_G0_out NE la_undef(),count)
815   - vv=(med_G0_in[ind]-med_G0_out[ind])/med_G0_out[ind]*100.
816   - cgplot,hii_count[ind],vv,tit='G0',xrange=xrange,/xlog
  687 + cgoplot,hii_per_count_bin_local_st.count,hii_per_count_bin_local_st.NH,color='blue',psym='filled circle'
  688 + ;hii_stats.med_NH=la_median(size_hii_stat.med_d_NH)
  689 + cgoplot,10^!x.crange,replicate(la_median(hii_median_local_stats.NH),2),linestyle=2,color='blue'
  690 +
  691 + ind=where(hii_median_local_stats.G0 NE la_undef())
  692 + cgplot,hii_median_local_stats[ind].count,hii_median_local_stats[ind].G0,tit='Delta-G0 [%]',xrange=xrange,/xlog,xtit=xtit
817 693 cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
818   - cgoplot,size_hii_stat.size,size_hii_stat.med_d_G0,color='blue',psym='filled circle'
819   - hii_stats.med_G0=la_median(size_hii_stat.med_d_G0)
820   - cgoplot,10^!x.crange,replicate(la_median(size_hii_stat.med_d_G0),2),linestyle=2,color='blue'
821   - hii_stats.std_G0=hii_stats.med_G0/la_sigma(vv)
822   -
823   - ind=where(med_YG1_in NE la_undef() and med_YG1_out NE la_undef(),count)
824   - vv=(med_YG1_in[ind]-med_YG1_out[ind])/med_YG1_out[ind]*100.
825   - cgplot,hii_count[ind],vv,tit='YG1',xrange=xrange,/xlog
  694 + cgoplot,hii_per_count_bin_local_st.count,hii_per_count_bin_local_st.G0,color='blue',psym='filled circle'
  695 + ;hii_stats.med_NH=la_median(size_hii_stat.med_d_NH)
  696 + cgoplot,10^!x.crange,replicate(la_median(hii_median_local_stats.G0),2),linestyle=2,color='blue'
  697 +
  698 + ind=where(hii_median_local_stats.YG1 NE la_undef())
  699 + cgplot,hii_median_local_stats[ind].count,hii_median_local_stats[ind].YG1,tit='Delta-YG1 [%]',xrange=xrange,/xlog,xtit=xtit
826 700 cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
827   - cgoplot,size_hii_stat.size,size_hii_stat.med_d_YG1,color='blue',psym='filled circle'
828   - hii_stats.med_YG1=la_median(size_hii_stat.med_d_YG1)
829   - cgoplot,10^!x.crange,replicate(hii_stats.med_YG1,2),linestyle=2,color='blue'
  701 + cgoplot,hii_per_count_bin_local_st.count,hii_per_count_bin_local_st.YG1,color='blue',psym='filled circle'
  702 + ;hii_stats.med_NH=la_median(size_hii_stat.med_d_NH)
  703 + cgoplot,10^!x.crange,replicate(la_median(hii_median_local_stats.YG1),2),linestyle=2,color='blue'
830 704  
831   - ind=where(med_YG2_in NE la_undef() and med_YG2_out NE la_undef(),count)
832   - vv=(med_YG2_in[ind]-med_YG2_out[ind])/med_YG2_out[ind]*100.
833   - cgplot,hii_count[ind],vv,tit='YG2',xrange=xrange,/xlog
  705 + ind=where(hii_median_local_stats.YG2 NE la_undef())
  706 + cgplot,hii_median_local_stats[ind].count,hii_median_local_stats[ind].YG2,tit='Delta-YG2 [%]',xrange=xrange,/xlog,xtit=xtit
834 707 cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
835   - cgoplot,size_hii_stat.size,size_hii_stat.med_d_YG2,color='blue',psym='filled circle'
836   - hii_stats.med_YG2=la_median(size_hii_stat.med_d_YG2)
837   - cgoplot,10^!x.crange,replicate(hii_stats.med_YG2,2),linestyle=2,color='blue'
838   -
839   - goto,jump_this_one
840   - ind=where(med_YG3_in NE la_undef() and med_YG3_out NE la_undef(),count)
841   - vv=(med_YG3_in[ind]-med_YG3_out[ind])/med_YG3_out[ind]*100.
842   - cgplot,hii_count[ind],vv,tit='YG3',xrange=xrange,/xlog
  708 + cgoplot,hii_per_count_bin_local_st.count,hii_per_count_bin_local_st.YG2,color='blue',psym='filled circle'
  709 + ;hii_stats.med_NH=la_median(size_hii_stat.med_d_NH)
  710 + cgoplot,10^!x.crange,replicate(la_median(hii_median_local_stats.YG2),2),linestyle=2,color='blue'
  711 +
  712 + ind=where(hii_median_local_stats.Xion NE la_undef())
  713 + cgplot,hii_median_local_stats[ind].count,hii_median_local_stats[ind].Xion,tit='Delta-Xion [%]',xrange=xrange,/xlog,xtit=xtit
843 714 cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
844   - cgoplot,size_hii_stat.size,size_hii_stat.med_d_YG3,color='blue',psym='filled circle'
845   - hii_stats.med_YG3=la_median(size_hii_stat.med_d_YG3)
846   - cgoplot,10^!x.crange,replicate(hii_stats.med_YG3,2),linestyle=2,color='blue'
847   - jump_this_one:
848   -
849   - ind=where(med_Xion_in NE la_undef() and med_Xion_out NE la_undef(),count)
850   - vv=(med_Xion_in[ind]-med_Xion_out[ind])/med_Xion_out[ind]*100.
851   - cgplot,hii_count[ind],vv,tit='Xion',xrange=xrange,/xlog
  715 + cgoplot,hii_per_count_bin_local_st.count,hii_per_count_bin_local_st.Xion,color='blue',psym='filled circle'
  716 + ;hii_stats.med_NH=la_median(size_hii_stat.med_d_NH)
  717 + cgoplot,10^!x.crange,replicate(la_median(hii_median_local_stats.Xion),2),linestyle=2,color='blue'
  718 +
  719 + ind=where(hii_median_local_stats.rchi2 NE la_undef())
  720 + cgplot,hii_median_local_stats[ind].count,hii_median_local_stats[ind].rchi2,tit='Delta-Rchi2 [%]',xrange=xrange,/xlog,xtit=xtit
852 721 cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
853   - cgoplot,size_hii_stat.size,size_hii_stat.med_d_Xion,color='blue',psym='filled circle'
854   - hii_stats.med_Xion=la_median(size_hii_stat.med_d_Xion)
855   - cgoplot,10^!x.crange,replicate(hii_stats.med_Xion,2),linestyle=2,color='blue'
  722 + cgoplot,hii_per_count_bin_local_st.count,hii_per_count_bin_local_st.rchi2,color='blue',psym='filled circle'
  723 + ;hii_stats.med_NH=la_median(size_hii_stat.med_d_NH)
  724 + cgoplot,10^!x.crange,replicate(la_median(hii_median_local_stats.rchi2),2),linestyle=2,color='blue'
856 725  
857   - ind=where(med_rchi2_in NE la_undef() and med_rchi2_out NE la_undef(),count)
858   - vv=(med_rchi2_in[ind]-med_rchi2_out[ind])/med_rchi2_out[ind]*100.
859   - cgplot,hii_count[ind],vv,tit='rchi2',xrange=xrange,/xlog,yrange=[-100,500],/ysty;,/ylog
860   - cgoplot,10^!x.crange,[1,1.],linestyle=2,color='red'
861   - cgoplot,size_hii_stat.size,size_hii_stat.med_d_rchi2,color='blue',psym='filled circle'
862   - hii_stats.med_rchi2=la_median(size_hii_stat.med_d_rchi2)
863   - cgoplot,10^!x.crange,replicate(hii_stats.med_rchi2,2),linestyle=2,color='blue'
  726 + cleanplot
  727 +ENDIF
  728 +
  729 +;stop
  730 +
  731 +mc_local_stats=regional_stats(mcs,images,backg_Npix=10,/median,/do_print $
  732 + ,struct_tags=struct_tags,stdev_obj_st=mc_median_stdev_obj_st,stdev_back_st=mc_median_stdev_back_st $
  733 + ,order_by='OBJECT_COUNT',local_obj_st=mc_median_local_stats,per_count_bin_local_st=mc_per_count_bin_local_st)
  734 +
  735 +
  736 +;==========================================================================================
  737 +;======= plot statistics on HII regions
  738 +;==========================================================================================
  739 +IF not keyword_set(no_mc) THEN BEGIN
  740 + window,win,xsize=wxs*2,ysize=wys & win=win+1
  741 + !p.multi=[0,3,2]
  742 + xrange=[1,max(mc_median_local_stats.count)]
  743 + xtit='Surface [pix]'
  744 + cgplot,mc_median_local_stats.count,mc_median_local_stats.NH,tit='Delta-NH [%]',xrange=xrange,/xlog,xtit=xtit
  745 + cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
  746 + cgoplot,mc_per_count_bin_local_st.count,mc_per_count_bin_local_st.NH,color='blue',psym='filled circle'
  747 + cgoplot,10^!x.crange,replicate(la_median(mc_median_local_stats.NH),2),linestyle=2,color='blue'
  748 +
  749 + ind=where(mc_median_local_stats.G0 NE la_undef())
  750 + cgplot,mc_median_local_stats[ind].count,mc_median_local_stats[ind].G0,tit='Delta-G0 [%]',xrange=xrange,/xlog,xtit=xtit
  751 + cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
  752 + cgoplot,mc_per_count_bin_local_st.count,mc_per_count_bin_local_st.G0,color='blue',psym='filled circle'
  753 + cgoplot,10^!x.crange,replicate(la_median(mc_median_local_stats.G0),2),linestyle=2,color='blue'
  754 +
  755 + ind=where(mc_median_local_stats.YG1 NE la_undef())
  756 + cgplot,mc_median_local_stats[ind].count,mc_median_local_stats[ind].YG1,tit='Delta-YG1 [%]',xrange=xrange,/xlog,xtit=xtit
  757 + cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
  758 + cgoplot,mc_per_count_bin_local_st.count,mc_per_count_bin_local_st.YG1,color='blue',psym='filled circle'
  759 + cgoplot,10^!x.crange,replicate(la_median(mc_median_local_stats.YG1),2),linestyle=2,color='blue'
  760 +
  761 + ind=where(mc_median_local_stats.YG2 NE la_undef())
  762 + cgplot,mc_median_local_stats[ind].count,mc_median_local_stats[ind].YG2,tit='Delta-YG2 [%]',xrange=xrange,/xlog,xtit=xtit
  763 + cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
  764 + cgoplot,mc_per_count_bin_local_st.count,mc_per_count_bin_local_st.YG2,color='blue',psym='filled circle'
  765 + cgoplot,10^!x.crange,replicate(la_median(mc_median_local_stats.YG2),2),linestyle=2,color='blue'
  766 +
  767 + ind=where(mc_median_local_stats.Xion NE la_undef())
  768 + cgplot,mc_median_local_stats[ind].count,mc_median_local_stats[ind].Xion,tit='Delta-Xion [%]',xrange=xrange,/xlog,xtit=xtit
  769 + cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
  770 + cgoplot,mc_per_count_bin_local_st.count,mc_per_count_bin_local_st.Xion,color='blue',psym='filled circle'
  771 + cgoplot,10^!x.crange,replicate(la_median(mc_median_local_stats.Xion),2),linestyle=2,color='blue'
864 772  
865   - write_xcat,hii_stats,'/tmp/bidon',/nowrite
  773 + ind=where(mc_median_local_stats.rchi2 NE la_undef())
  774 + cgplot,mc_median_local_stats[ind].count,mc_median_local_stats[ind].rchi2,tit='Delta-Rchi2 [%]',xrange=xrange,/xlog,xtit=xtit
  775 + cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
  776 + cgoplot,mc_per_count_bin_local_st.count,mc_per_count_bin_local_st.rchi2,color='blue',psym='filled circle'
  777 + cgoplot,10^!x.crange,replicate(la_median(mc_median_local_stats.rchi2),2),linestyle=2,color='blue'
866 778  
867 779 cleanplot
868 780 ENDIF
  781 +
  782 +;stop
  783 +
869 784 IF not keyword_set(nostop) THEN stop
870 785  
871   -;=== compute global statistics of fitted parameters in various regions
  786 +;=== compute global statistics of fitted parameters in various regions for Method paper
872 787 ;define NH Threshold values
873 788 NHthr_min=[0.01,0.1,1 ,10.] ;1.e21
874 789 NHthr_max=[0.1 ,1 ,10,100.] ;1.e21
875 790 one_stat={region:'',count:0L,med_NH:la_undef(4),med_G0:la_undef(4),med_YG1:la_undef(4),med_YG2:la_undef(4),med_YG3:la_undef(4),med_Xion:la_undef(4),med_rchi2:la_undef(4)}
876   -regions=['ALL','HII','MC','ARMS','INTERARMS','NHR0','NHR1','NHR2','NHR3']
  791 +regions=['ALL','HII','MC','ARMS','INTERARMS','BULGE','NHR0','NHR1','NHR2','NHR3','LOCAL_HII','LOCAL_MC']
877 792 Nregions=n_elements(regions)
878 793 regional_stats=replicate(one_stat,Nregions)
879 794 regional_stats_stdev=replicate(one_stat,Nregions)
880 795 regional_stats.region=regions
881 796 regional_stats_stdev.region=regions
  797 +scalings=[1.,1.,1.e4,1.e4,1.e4,1.,1.]
882 798 FOR i=0L,Nregions-1 DO BEGIN
883 799 CASE regions[i] OF
884 800 'ALL': BEGIN
... ... @@ -899,6 +815,10 @@ FOR i=0L,Nregions-1 DO BEGIN
899 815 'INTERARMS': BEGIN
900 816 ind=where(G0_map NE la_undef() AND arms_im EQ 0,count)
901 817 END
  818 + 'BULGE': BEGIN
  819 + ind=where(G0_map NE la_undef() AND bulge_im GE 1,count)
  820 + ;stop
  821 + END
902 822 'NHR0': BEGIN
903 823 ind=where(G0_map NE la_undef() AND NHmap GT NHthr_min[0] AND NHmap LE NHthr_max[0],count)
904 824 END
... ... @@ -911,28 +831,52 @@ FOR i=0L,Nregions-1 DO BEGIN
911 831 'NHR3': BEGIN
912 832 ind=where(G0_map NE la_undef() AND NHmap GT NHthr_min[3] AND NHmap LE NHthr_max[3],count)
913 833 END
  834 + 'LOCAL_HII': BEGIN
  835 + ind=where(G0_map NE la_undef(),count)
  836 + END
  837 + 'LOCAL_MC': BEGIN
  838 + ind=where(G0_map NE la_undef(),count)
  839 + END
914 840 ENDCASE
915   - regional_stats[i].med_NH=la_median(NHmap[ind])
916   - regional_stats_stdev[i].med_NH=la_sigma(NHmap[ind])
917   - regional_stats[i].med_G0=la_median(G0_map[ind])
918   - regional_stats_stdev[i].med_G0=la_sigma(G0_map[ind])
  841 + regional_stats[i].med_NH=la_mul(la_median(NHmap[ind]),scalings[0])
  842 + regional_stats_stdev[i].med_NH=la_mul(la_sigma(NHmap[ind]),scalings[0])
  843 + regional_stats[i].med_G0=la_mul(la_median(G0_map[ind]),scalings[1])
  844 + regional_stats_stdev[i].med_G0=la_mul(la_sigma(G0_map[ind]),scalings[1])
919 845 regional_stats[i].count=count
920 846 regional_stats_stdev[i].count=count
921   - regional_stats[i].med_YG1=la_median(Yg1_map[ind])
922   - regional_stats_stdev[i].med_YG1=la_sigma(Yg1_map[ind])
923   - regional_stats[i].med_YG2=la_median(Yg2_map[ind])
924   - regional_stats_stdev[i].med_YG2=la_sigma(Yg2_map[ind])
925   - regional_stats[i].med_YG3=la_median(Yg3_map[ind])
926   - regional_stats_stdev[i].med_YG3=la_sigma(Yg3_map[ind])
927   - regional_stats[i].med_Xion=la_median(Xion_map[ind])
928   - regional_stats_stdev[i].med_Xion=la_sigma(Xion_map[ind])
929   - regional_stats[i].med_rchi2=la_median(rchi2_map[ind])
930   - regional_stats_stdev[i].med_rchi2=la_sigma(rchi2_map[ind])
  847 + regional_stats[i].med_YG1=la_mul(la_median(Yg1_map[ind]),scalings[2])
  848 + regional_stats_stdev[i].med_YG1=la_mul(la_sigma(Yg1_map[ind]),scalings[2])
  849 + regional_stats[i].med_YG2=la_mul(la_median(Yg2_map[ind]),scalings[3])
  850 + regional_stats_stdev[i].med_YG2=la_mul(la_sigma(Yg2_map[ind]),scalings[3])
  851 + regional_stats[i].med_YG3=la_mul(la_median(Yg3_map[ind]),scalings[4])
  852 + regional_stats_stdev[i].med_YG3=la_mul(la_sigma(Yg3_map[ind]),scalings[4])
  853 + regional_stats[i].med_Xion=la_mul(la_median(Xion_map[ind]),scalings[5])
  854 + regional_stats_stdev[i].med_Xion=la_mul(la_sigma(Xion_map[ind]),scalings[5])
  855 + regional_stats[i].med_rchi2=la_mul(la_median(rchi2_map[ind]),scalings[6])
  856 + regional_stats_stdev[i].med_rchi2=la_mul(la_sigma(rchi2_map[ind]),scalings[6])
  857 + IF regions[i] EQ 'LOCAL_HII' THEN BEGIN
  858 + regional_stats[i].med_NH=la_median(hii_median_local_stats.NH)
  859 + regional_stats[i].med_G0=la_median(hii_median_local_stats.G0)
  860 + regional_stats[i].med_YG1=la_median(hii_median_local_stats.YG1)
  861 + regional_stats[i].med_YG2=la_median(hii_median_local_stats.YG2)
  862 + regional_stats[i].med_YG3=la_median(hii_median_local_stats.YG3)
  863 + regional_stats[i].med_XION=la_median(hii_median_local_stats.XION)
  864 + regional_stats[i].med_rchi2=la_median(hii_median_local_stats.rchi2)
  865 + ENDIF
  866 + IF regions[i] EQ 'LOCAL_MC' THEN BEGIN
  867 + regional_stats[i].med_NH=la_median(mc_median_local_stats.NH)
  868 + regional_stats[i].med_G0=la_median(mc_median_local_stats.G0)
  869 + regional_stats[i].med_YG1=la_median(mc_median_local_stats.YG1)
  870 + regional_stats[i].med_YG2=la_median(mc_median_local_stats.YG2)
  871 + regional_stats[i].med_YG3=la_median(mc_median_local_stats.YG3)
  872 + regional_stats[i].med_XION=la_median(mc_median_local_stats.XION)
  873 + regional_stats[i].med_rchi2=la_median(mc_median_local_stats.rchi2)
  874 + ENDIF
931 875 ENDFOR
932 876  
933 877 write_xcat,regional_stats,'/tmp/bidon',/nowrite
934 878  
935   -stop
  879 +;stop
936 880  
937 881 file=!phangs_data_dir+'/ISRF/WORK/'+use_source_name+'_regional_statistics'+reso_str+'.xcat'
938 882 write_xcat,regional_stats,file
... ... @@ -948,22 +892,15 @@ regional_stats=read_xcat(file)
948 892 file=!phangs_data_dir+'/ISRF/WORK/'+use_source_name+'_regional_statistics_stdev'+reso_str+'.xcat'
949 893 regional_stats_stdev=read_xcat(file)
950 894  
951   -use_regional_stats=regional_stats
952   -use_regional_stats_stdev=regional_stats_stdev
953   -use_regional_stats.med_G0=la_mul(use_regional_stats.med_G0,1.)
954   -use_regional_stats_stdev.med_G0=la_mul(use_regional_stats_stdev.med_G0,1.)
955   -use_regional_stats.med_YG1=la_mul(use_regional_stats.med_YG1,1.e4)
956   -use_regional_stats_stdev.med_YG1=la_mul(use_regional_stats_stdev.med_YG1,1.e4)
957   -use_regional_stats.med_YG2=la_mul(use_regional_stats.med_YG2,1.e4)
958   -use_regional_stats_stdev.med_YG2=la_mul(use_regional_stats_stdev.med_YG2,1.e4)
959   -use_regional_stats.med_YG3=la_mul(use_regional_stats.med_YG3,1.e4)
960   -use_regional_stats_stdev.med_YG3=la_mul(use_regional_stats_stdev.med_YG3,1.e4)
961   -use_regional_stats.med_Xion=la_mul(use_regional_stats.med_xion,1.)
962   -use_regional_stats_stdev.med_xion=la_mul(use_regional_stats_stdev.med_xion,1.)
963   -
964   -use_all_format='(A10, " & ",I6, " & ",E10.2, " & ",E10.2, " & ",E10.2, " & ",E10.2, " & ",E10.2, " & ",E10.2, " & ",E10.2," \\")'
  895 +;stop
  896 +use_regional_stats=structure_remove_column(structure_remove_column(regional_stats,'MED_YG3'),'MED_RCHI2')
  897 +use_regional_stats_stdev=structure_remove_column(structure_remove_column(regional_stats_stdev,'MED_YG3'),'MED_RCHI2')
  898 +
  899 +;use_all_format='(A10, " & ",I6, " & ",E10.2, " & ",E10.2, " & ",E10.2, " & ",E10.2, " & ",E10.2, " & ",E10.2, " & ",E10.2," \\")'
  900 +use_all_format='(A10, " & ",I6, " & ",E10.2, " & ",E10.2, " & ",E10.2, " & ",E10.2, " & ",E10.2, "\\")'
965 901 error_st=use_regional_stats_stdev
966   -erindex=[2,3,4,5,6,7]
  902 +;erindex=[2,3,4,5,6,7]
  903 +erindex=[2,3,4,5]
967 904 fileout=!phangs_data_dir+'/ISRF/WORK/'+use_source_name+'_regional_statistics_latex'+reso_str+'.txt'
968 905 struct2latex_table,use_regional_stats,fileout,exp_accur=exp_accur,silent=silent,formats=formats,indef_replace=indef_replace, $
969 906 error_st=error_st,use_all_format=use_all_format,erindex=erindex,parents=parents, $
... ... @@ -971,6 +908,8 @@ struct2latex_table,use_regional_stats,fileout,exp_accur=exp_accur,silent=silent,
971 908 replace_tagnames_only=replace_tagnames_only,small=small,tiny=tiny,landscape=landscape, $
972 909 long=long, units=units,forcepos=forcepos,label=label,two_columns=two_columns
973 910  
  911 +;stop
  912 +
974 913 IF not keyword_set(nostop) THEN stop
975 914  
976 915 goto,dontdothat1
... ... @@ -1369,9 +1308,6 @@ no_uv_plots:
1369 1308 IF not keyword_set(nostop) THEN stop
1370 1309  
1371 1310 ;==== show maps
1372   -obp=[1.1,0.,1.15,1]
1373   -cleanplot
1374   -image_color_table='jpbloadct,/whitebackground'
1375 1311  
1376 1312 ;==========================================================================================
1377 1313 ;======== map of NH
... ... @@ -1768,7 +1704,7 @@ IF not keyword_set(nostop) THEN stop
1768 1704 win=0L ;avoid "no more windows"
1769 1705  
1770 1706 ;=========================================================================================================
1771   -;====== maps of residuals in %
  1707 +;====== plot maps of residuals in %
1772 1708 ;=========================================================================================================
1773 1709 window,win,xsize=wxs,ysize=wys & win=win+1
1774 1710 residual_range=[-100,100.]/2.
... ...