OV-GSO-DC / dustem-wrapper_idl

Commit f3236ab4ebe9bab351572d7370c8cbbecf7e9cbd

Authored by Jean-Philippe Bernard
1 parent 9a95b8f3
Exists in master

first commit

Showing 2 changed files with 1225 additions and 0 deletions   Show diff stats
Inline Side-by-side
LabTools/IRAP/JPB/phangs_make_figures_from_restore.pro 0 → 100644
  Diff comments   View file @f3236ab
... ... @@ -0,0 +1,1143 @@
  1 +PRO phangs_make_figures_from_restore,source_name=source_name $
  2 + ,model=model $
  3 + ,normalize=normalize $
  4 + ,force_mathis=force_mathis $
  5 + ,fit_G0=fit_G0 $
  6 + ,set_yvsg=set_yvsg $
  7 + ,help=help $
  8 + ,resolution_filter=resolution_filter $
  9 + ,resolution_value=resolution_value $
  10 + ,table_name_root=table_name_root $
  11 + ,nostop=nostop $
  12 + ,marginalize_method=marginalize_method $
  13 + ,do_postscripts=do_postscripts $
  14 + ,postscripts_dir=postscripts_dir $
  15 + ,silent=silent
  16 +
  17 +obp=[1.1,0.,1.15,1]
  18 +win=0L
  19 +wxs=800
  20 +wys=900
  21 +
  22 +;weightning='1/chi2'
  23 +;weightning='likelyhood'
  24 +weightning='rellikelyhood' ;weightning strategy used for weighted parameters
  25 +
  26 +use_source_name='ngc0628'
  27 +IF keyword_set(source_name) THEN use_source_name=source_name
  28 +
  29 +norm_str=''
  30 +IF keyword_set(normalize) THEN BEGIN
  31 + norm_str='_normalized'
  32 +ENDIF
  33 +mathis_str=''
  34 +IF keyword_set(force_mathis) THEN BEGIN
  35 + mathis_str='_mathis'
  36 +ENDIF
  37 +
  38 +IF keyword_set(set_yvsg) and not keyword_set(fit_G0) THEN BEGIN
  39 + message,'/fit_G0 should be set with /set_yvsg.',/continue
  40 + stop
  41 +ENDIF
  42 +
  43 +pdp_define_la_common
  44 +dustem_init
  45 +
  46 +reso_str=''
  47 +IF keyword_set(resolution_filter) OR keyword_set(resolution_value) THEN BEGIN
  48 + IF keyword_set(resolution_filter) THEN BEGIN
  49 + reso_str='_'+resolution_filter
  50 + reso=dustem_filter2reso(resolution_filter)
  51 + ENDIF ELSE BEGIN
  52 + reso_str='_'+strtrim(resolution_value,2)
  53 + reso=resolution_value
  54 + ENDELSE
  55 +ENDIF
  56 +
  57 +use_model='DBP90'
  58 +IF keyword_set(model) THEN BEGIN
  59 + use_model=model
  60 +ENDIF
  61 +;stop
  62 +;=== define what we call vsg, for /set_yvsg option
  63 +CASE use_model OF
  64 + 'DL07': vsg_id=2L
  65 + 'DBP90': vsg_id=1L
  66 +ENDCASE
  67 +vsg_id_str=strtrim(vsg_id,2)
  68 +;stop
  69 +
  70 +data_dir=!phangs_data_dir+'/ISRF/WORK/'
  71 +grids_data_dir=!phangs_data_dir+'/ISRF/GRIDS/'
  72 +
  73 +IF keyword_set(do_postscripts) THEN BEGIN
  74 + ps_dir=!phangs_data_dir+'/ISRF/WORK/PS/'
  75 + IF keyword_set(postscripts_dir) THEN BEGIN
  76 + ps_dir=postscripts_dir+'/'
  77 + ENDIF
  78 + force_mkdir,ps_dir
  79 +ENDIF
  80 +
  81 +IF not keyword_set(fit_G0) THEN BEGIN
  82 + fit_G0_str='_GOfixed'
  83 +ENDIF ELSE BEGIN
  84 + fit_G0_str=''
  85 +ENDELSE
  86 +
  87 +IF keyword_set(set_yvsg) THEN BEGIN
  88 + set_yvsg_str='_YVSGfixed'
  89 +ENDIF ELSE BEGIN
  90 + set_yvsg_str=''
  91 +ENDELSE
  92 +
  93 +;stop
  94 +
  95 +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'
  96 +st_info=file_info(file)
  97 +IF st_info.exists NE 1 THEN BEGIN
  98 + message,'Could not find '+file,/continue
  99 + stop
  100 +ENDIF
  101 +restore,file,/verb
  102 +; % RESTORE: Restored variable: STATS_ST.
  103 +; % RESTORE: Restored variable: F0S.
  104 +; % RESTORE: Restored variable: F0S_PREDICTED.
  105 +; % RESTORE: Restored variable: YPAHS.
  106 +; % RESTORE: Restored variable: YVSGS.
  107 +; % RESTORE: Restored variable: FACTS.
  108 +; % RESTORE: Restored variable: CHI2S.
  109 +; % RESTORE: Restored variable: RCHI2S.
  110 +; % RESTORE: Restored variable: DF0S.
  111 +; % RESTORE: Restored variable: DYPAHS.
  112 +; % RESTORE: Restored variable: DYVSGS.
  113 +; % RESTORE: Restored variable: F0_HIT.
  114 +; % RESTORE: Restored variable: YPAH_HIT.
  115 +; % RESTORE: Restored variable: YVSG_HIT.
  116 +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'
  117 +st_info=file_info(file)
  118 +IF st_info.exists NE 1 THEN BEGIN
  119 + message,'Could not find '+file,/continue
  120 + stop
  121 +ENDIF
  122 +restore,file,/verb
  123 +; % RESTORE: Restored variable: NHCO.
  124 +; % RESTORE: Restored variable: NHI.
  125 +; % RESTORE: Restored variable: NHMAP.
  126 +; % RESTORE: Restored variable: NORM_STARS_MAP.
  127 +; % RESTORE: Restored variable: EBV_MAP.
  128 +; % RESTORE: Restored variable: F0_MAP.
  129 +; % RESTORE: Restored variable: G0_MAP.
  130 +; % RESTORE: Restored variable: F0_PREDICTED_MAP.
  131 +; % RESTORE: Restored variable: G0_PREDICTED_MAP.
  132 +; % RESTORE: Restored variable: YG1_MAP.
  133 +; % RESTORE: Restored variable: YG2_MAP.
  134 +; % RESTORE: Restored variable: YG3_MAP.
  135 +; % RESTORE: Restored variable: FACT_MAP.
  136 +; % RESTORE: Restored variable: CHI2_MAP.
  137 +; % RESTORE: Restored variable: RCHI2_MAP.
  138 +; % RESTORE: Restored variable: DF0_MAP.
  139 +; % RESTORE: Restored variable: DYG1_MAP.
  140 +; % RESTORE: Restored variable: DYG2_MAP.
  141 +; % RESTORE: Restored variable: DYG3_MAP.
  142 +; % RESTORE: Restored variable: F0_HIT_MAP.
  143 +; % RESTORE: Restored variable: YG1_HIT_MAP.
  144 +; % RESTORE: Restored variable: YG2_HIT_MAP.
  145 +; % RESTORE: Restored variable: YG3_HIT_MAP.
  146 +; % RESTORE: Restored variable: NH_VALUES_MAP.
  147 +; % RESTORE: Restored variable: RESIDUAL_MAPS.
  148 +; % RESTORE: Restored variable: RESIDUAL_MAPS_PERC.
  149 +file=data_dir+use_source_name+'_isrf_classes_voronoi'+reso_str+'.sav'
  150 +st_info=file_info(file)
  151 +IF st_info.exists NE 1 THEN BEGIN
  152 + message,'Could not find '+file,/continue
  153 + stop
  154 +ENDIF
  155 +restore,file,/verb
  156 +;% RESTORE: Restored variable: VOR_CLASS.
  157 +file=data_dir+use_source_name+'_ref_header'+reso_str+'.sav'
  158 +st_info=file_info(file)
  159 +IF st_info.exists NE 1 THEN BEGIN
  160 + message,'Could not find '+file,/continue
  161 + stop
  162 +ENDIF
  163 +restore,file,/verb
  164 +;% RESTORE: Restored variable: HREF.
  165 +
  166 +IF defined(dustem_grid_saved) THEN defsysv,'!dustem_grid',dustem_grid_saved
  167 +;=== get the stellar mask
  168 +file=data_dir+use_source_name+'_star_mask'+reso_str+'.sav'
  169 +restore,file,/verb
  170 +;% RESTORE: Restored variable: STAR_MASK.
  171 +;% RESTORE: Restored variable: HSTARMASK.
  172 +;% RESTORE: Restored variable: STAR_MASK_REF.
  173 +
  174 +;stop
  175 +
  176 +;==== restore mask of low brightness regions
  177 +file=data_dir+use_source_name+'_brightness_mask'+reso_str+'.sav'
  178 +restore,file,/verb
  179 +;% RESTORE: Restored variable: BRIGHTNESS_MASK.
  180 +;stop
  181 +
  182 +;==== set masked values in all maps
  183 +;ind_masked=where(star_mask_ref NE 0 and star_mask_ref NE la_undef(),count_masked)
  184 +ind_masked=where((star_mask_ref NE 0 and star_mask_ref NE la_undef()) OR $
  185 + (brightness_mask NE 0 and brightness_mask NE la_undef()),count_masked)
  186 +;stop
  187 +IF count_masked NE 0 THEN BEGIN
  188 + FOR i=0L,n_elements(all_filters)-1 DO BEGIN
  189 + map=residual_maps_perc[*,*,i]
  190 + map[ind_masked]=la_undef()
  191 + residual_maps_perc[*,*,i]=map
  192 + map=residual_maps[*,*,i]
  193 + map[ind_masked]=la_undef()
  194 + residual_maps[*,*,i]=map
  195 + ENDFOR
  196 + NHmap[ind_masked]=la_undef()
  197 + norm_stars_map[ind_masked]=la_undef()
  198 + EBV_map[ind_masked]=la_undef()
  199 + F0_map[ind_masked]=la_undef()
  200 + G0_map[ind_masked]=la_undef()
  201 + F0_predicted_map[ind_masked]=la_undef()
  202 + G0_predicted_map[ind_masked]=la_undef()
  203 + Yg1_map[ind_masked]=la_undef()
  204 + Yg2_map[ind_masked]=la_undef()
  205 + Yg3_map[ind_masked]=la_undef()
  206 + fact_map[ind_masked]=la_undef()
  207 + chi2_map[ind_masked]=la_undef()
  208 + rchi2_map[ind_masked]=la_undef()
  209 + dF0_map[ind_masked]=la_undef()
  210 + dYg1_map[ind_masked]=la_undef()
  211 + dYg2_map[ind_masked]=la_undef()
  212 + dYg3_map[ind_masked]=la_undef()
  213 + F0_hit_map[ind_masked]=la_undef()
  214 + Yg1_hit_map[ind_masked]=la_undef()
  215 + Yg2_hit_map[ind_masked]=la_undef()
  216 + Yg3_hit_map[ind_masked]=la_undef()
  217 + NH_values_map[ind_masked]=la_undef()
  218 +ENDIF
  219 +
  220 +Xion_map=la_div(Yg2_map,la_add(Yg1_map,Yg2_map))
  221 +
  222 +cleanplot
  223 +win=0L
  224 +
  225 +;This is to allow dustem to return default parameter values
  226 +dustem_init,model=use_model
  227 +
  228 +Nfilters=n_elements(all_filters)
  229 +wavs=dustem_filter2wav(all_filters)
  230 +NH_thr=2. ;2e21 H/cm2
  231 +
  232 +;=== compute the average and median residuals
  233 +residuals_avg=fltarr(Nfilters)+la_undef()
  234 +residuals_median=fltarr(Nfilters)+la_undef()
  235 +residuals_avg_highNH=fltarr(Nfilters)+la_undef()
  236 +residuals_median_highNH=fltarr(Nfilters)+la_undef()
  237 +FOR i=0L,Nfilters-1 DO BEGIN
  238 + vv=residual_maps_perc[*,*,i]
  239 + ind=where(vv ne la_undef() AND finite(vv))
  240 + residuals_avg[i]=la_mean(vv[ind])
  241 + residuals_median[i]=la_median(vv[ind])
  242 + ind=where(vv ne la_undef() AND finite(vv) AND NHmap GT NH_thr,count)
  243 + residuals_avg_highNH[i]=la_mean(vv[ind])
  244 + residuals_median_highNH[i]=la_median(vv[ind])
  245 +ENDFOR
  246 +
  247 +;window,win & win=win+1
  248 +;ind=where(residuals_avg NE la_undef() and finite(residuals_avg))
  249 +;cgplot,wavs[ind],abs(residuals_avg[ind]),/xlog,/ylog
  250 +;plot_negative_log,wavs[ind],residuals_avg[ind]
  251 +
  252 +all_resos=dustem_filter2reso(all_filters)
  253 +ind=where(residuals_median NE la_undef() and finite(residuals_median) and all_resos LE reso,count)
  254 +xtit='Wavelength [mic]'
  255 +ytit='(model - data)/data [%]'
  256 +yrange=[-40.,40.]
  257 +yrange=[-100.,100.]
  258 +cgplot,wavs[ind],residuals_median[ind],/xlog,xrange=[0.1,500],yrange=yrange,/xsty,/ysty,/nodata,xtit=xtit,ytit=ytit
  259 +cgoplot,wavs[ind],residuals_median[ind],psym='Filled Circle',color='red'
  260 +cgoplot,wavs[ind],residuals_median[ind],linestyle=2,color='red'
  261 +cgoplot,wavs[ind],residuals_avg[ind],psym='Filled Circle',color='blue'
  262 +cgoplot,wavs[ind],residuals_avg[ind],linestyle=2,color='blue'
  263 +cgoplot,wavs[ind],residuals_median_highNH[ind],psym='Filled Square',color='red'
  264 +cgoplot,wavs[ind],residuals_median_highNH[ind],linestyle=2,color='red'
  265 +cgoplot,wavs[ind],residuals_avg_highNH[ind],psym='Filled Square',color='blue'
  266 +cgoplot,wavs[ind],residuals_avg_highNH[ind],linestyle=2,color='blue'
  267 +cgoplot,10^!x.crange,[0.,0.],linestyle=2,color='black'
  268 +;plot_negative_log,wavs[ind],residuals_avg[ind]
  269 +print,residuals_median[ind]
  270 +print,residuals_avg[ind]
  271 +
  272 +;stop
  273 +
  274 +IF keyword_set(do_postscripts) THEN BEGIN
  275 + cleanplot
  276 +
  277 + ps_file=ps_dir+use_source_name+'_residual_sed_'+reso_str+fit_G0_str+mathis_str+'.ps'
  278 + set_plot,'PS'
  279 + device,filename=ps_file
  280 +
  281 + cgplot,wavs[ind],residuals_median[ind],/xlog,xrange=[0.1,500],yrange=[-40.,40.],/xsty,/ysty,/nodata,xtit=xtit,ytit=ytit
  282 + cgoplot,wavs[ind],residuals_median[ind],psym='Filled Circle',color='red'
  283 + cgoplot,wavs[ind],residuals_median[ind],linestyle=2,color='red'
  284 + cgoplot,wavs[ind],residuals_avg[ind],psym='Filled Circle',color='blue'
  285 + cgoplot,wavs[ind],residuals_avg[ind],linestyle=2,color='blue'
  286 + cgoplot,10^!x.crange,[0.,0.],linestyle=2,color='black'
  287 +
  288 + device,/close
  289 + set_plot,'X'
  290 + message,'Wrote '+ps_file,/continue
  291 + ;stop
  292 +ENDIF
  293 +
  294 +Nx=sxpar(href,'NAXIS1')
  295 +Ny=sxpar(href,'NAXIS2')
  296 +;=== This is to make a common mask of undefined values
  297 +mask=fltarr(Nx,Ny)
  298 +mask[*]=1
  299 +ind=where(G0_map EQ la_undef(),count)
  300 +IF count NE 0 THEN BEGIN
  301 + mask[ind]=0
  302 +ENDIF
  303 +ind_bad=where(mask NE 1,count_bad)
  304 +;=== apply mask to NH maps
  305 +IF count_bad NE 0 THEN BEGIN
  306 + NHmap[ind_bad]=la_undef()
  307 + F0_map[ind_bad]=la_undef()
  308 + F0_predicted_map[ind_bad]=la_undef()
  309 +ENDIF
  310 +
  311 +;=== make an overlay with contours of the total NH
  312 +;=== We smooth the NH map to 5 arcesec, or else, too many contours are shown.
  313 +IF reso LE 5./60.^2 THEN BEGIN
  314 + use_NHmap=degrade_res(NHmap,href,reso,5./60.^2,hout) ;in 1.e21 H/cm2
  315 + use_NHmap=project2(hout,use_NHmap,href,/silent)
  316 +ENDIF ELSE BEGIN
  317 + use_NHmap=NHmap
  318 +ENDELSE
  319 +;Nlevs=7
  320 +;levs=range_gen(Nlevs,[-0.5,1])
  321 +Nlevs=9
  322 +levs=range_gen(Nlevs,[-0.5,2])
  323 +
  324 +;levs=[la_undef(),-0.5,-0.25,0.,0.25,0.5,0.75,1.]
  325 +;levs=[-0.5,-0.25,0.,0.25,0.5,0.75,1.]
  326 +NNHmap=fltarr(Nx,Ny,2)
  327 +NNHmap[*,*,0]=la_log10(use_NHmap)
  328 +NNHmap[*,*,1]=la_log10(use_NHmap)
  329 +cont_st=create_contour_st(NNHmap,lev1=levs)
  330 +cont_st.(0).color=0
  331 +cont_st.(0).thick=2.
  332 +cont_st.(1).color=100
  333 +cont_st.(1).thick=1.0
  334 +
  335 +;===== get HII region mask
  336 +obj_name=strupcase(source_name)
  337 +IF obj_name EQ 'NGC0628' THEN obj_name='NGC628'
  338 +dir=!phangs_data_dir+'/phangs_drive/Archive/Products/Nebulae_catalogs/Nebulae_catalogue_v1/spatial_masks/'
  339 +file=dir+obj_name+'_HIIreg_mask.fits'
  340 +hii=readfits(file,hhii)
  341 +ind=where(finite(hii) NE 1,count)
  342 +IF count NE 0 THEN hii[ind]=la_undef()
  343 +hhii=cd2astro_header(hhii)
  344 +sxaddpar,hhii,'EQUINOX',2000.
  345 +image_cont20,hii,hhii,/square,/silent
  346 +hii=project2(hhii,hii,href,/silent,/round)
  347 +window,win,xsize=wxs,ysize=wys & win=win+1
  348 +image_cont20,hii,href,/square,/silent,image_color_table='jpbloadct,/white',axis_color_table=1,off_bar=obp
  349 +;=== do a simple HII region mask (1 if an HII, 0 otherwise) and a contour structure
  350 +is_hii=hii
  351 +ind=where(hii NE -1 and hii NE la_undef(),count)
  352 +is_hii[ind]=1
  353 +is_hii[ind_bad]=la_undef()
  354 +;stop
  355 +hii_m=fltarr(Nx,Ny,2)
  356 +hii_m[*,*,0]=is_hii
  357 +hii_m[*,*,1]=is_hii
  358 +levs_HII=[0.5,1]
  359 +hii_st=create_contour_st(hii_m,lev1=levs_HII,lev2=levs_HII)
  360 +hii_st.(0).color=0
  361 +hii_st.(0).thick=2.
  362 +hii_st.(1).color=100
  363 +hii_st.(1).thick=1.0
  364 +;hii[*]=la_undef()
  365 +image_cont20,hii,href,/square,/silent,levels=hii_st,image_color_table='jpbloadct,/no',axis_color_table=1,off_bar=obp
  366 +;image_cont20,is_hii,href,/square,/silent,levels=hii_st,image_color_table='jpbloadct',axis_color_table=1
  367 +;stop
  368 +;=== do a mask per HII region and a surounding mask
  369 +N_hii=long(max(hii))
  370 +hii_indices=ptrarr(N_hii)
  371 +hii_count=lonarr(N_hii)
  372 +hii_numbers=lonarr(N_hii)
  373 +FOR i=0L,N_hii-1 DO BEGIN
  374 + ind_hii=where(hii EQ i,count)
  375 + IF count NE 0 THEN BEGIN
  376 + hii_indices[i]=ptr_new(ind_hii)
  377 + ENDIF
  378 + hii_count[i]=count
  379 + hii_numbers[i]=i+1
  380 +ENDFOR
  381 +;=== remove unfound HII regions
  382 +ind=where(hii_count NE 0,count)
  383 +IF count NE 0 THEN BEGIN
  384 + hii_count=hii_count[ind]
  385 + hii_indices=hii_indices[ind]
  386 + hii_numbers=hii_numbers[ind] ;Original HII region number
  387 + N_hii=n_elements(hii_indices) ;new number of HII regions
  388 +ENDIF
  389 +;=== order HII regions by size
  390 +order=sort(hii_count)
  391 +hii_count=hii_count[order]
  392 +hii_indices=hii_indices[order]
  393 +hii_numbers=hii_numbers[order]
  394 +;=== do a surounding mask per HII region
  395 +around_hii_indices=ptrarr(N_hii)
  396 +Npix=5 ;number of pixels by which to extend
  397 +count_env=lonarr(N_hii)
  398 +maskk=intarr(Nx,Ny)
  399 +ind=where(is_hii EQ 1)
  400 +maskk[ind]=1
  401 +FOR i=0L,N_hii-1 DO BEGIN
  402 + indd=enlarge(*hii_indices[i],Npix,Nx,Ny,/remove_seed,remove_these=maskk,count=count) ;This enlarges but does remove coliding hii regions.
  403 + IF count NE 0 THEN around_hii_indices[i]=ptr_new(indd)
  404 + count_env[i]=count
  405 +ENDFOR
  406 +stop
  407 +
  408 +;=== compute statistics of fitted parameters in various regions
  409 +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)}
  410 +regions=['ALL','HII','not-HII']
  411 +Nregions=n_elements(regions)
  412 +regional_stats=replicate(one_stat,Nregions)
  413 +regional_stats.region=regions
  414 +FOR i=0L,Nregions-1 DO BEGIN
  415 + CASE regions[i] OF
  416 + 'ALL': BEGIN
  417 + ind=where(G0_map NE la_undef(),count)
  418 + END
  419 + 'HII': BEGIN
  420 + ind=where(G0_map NE la_undef() AND is_hii NE 0 AND is_hii NE la_undef(),count)
  421 + END
  422 + 'not-HII': BEGIN
  423 + ind=where(G0_map NE la_undef() AND is_hii EQ 1 AND is_hii NE la_undef(),count)
  424 + END
  425 + ENDCASE
  426 + regional_stats[i].med_G0=la_median(G0_map[ind])
  427 + regional_stats[i].count=count
  428 + regional_stats[i].med_YG1=la_median(Yg1_map[ind])
  429 + regional_stats[i].med_YG2=la_median(Yg2_map[ind])
  430 + regional_stats[i].med_YG3=la_median(Yg3_map[ind])
  431 + regional_stats[i].med_Xion=la_median(Xion_map[ind])
  432 + regional_stats[i].med_rchi2=la_median(rchi2_map[ind])
  433 +ENDFOR
  434 +
  435 +write_xcat,regional_stats,'/tmp/bidon',/nowrite
  436 +
  437 +;=== compute statistics of fitted parameters locally around HII regions
  438 +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)}
  439 +regions=['HII']
  440 +Nregions=n_elements(regions)
  441 +hii_stats=replicate(one_stat,Nregions)
  442 +hii_stats.region=regions
  443 +med_NH_in=fltarr(N_hii) & med_NH_out=fltarr(N_hii)
  444 +med_G0_in=fltarr(N_hii) & med_G0_out=fltarr(N_hii)
  445 +med_YG1_in=fltarr(N_hii) & med_YG1_out=fltarr(N_hii)
  446 +med_YG2_in=fltarr(N_hii) & med_YG2_out=fltarr(N_hii)
  447 +med_YG3_in=fltarr(N_hii) & med_YG3_out=fltarr(N_hii)
  448 +med_Xion_in=fltarr(N_hii) & med_Xion_out=fltarr(N_hii)
  449 +med_rchi2_in=fltarr(N_hii) & med_rchi2_out=fltarr(N_hii)
  450 +FOR i=0L,N_hii-1 DO BEGIN
  451 + med_NH_in[i]=la_median(NH_values_map[*hii_indices[i]])
  452 + med_NH_out[i]=la_median(NH_values_map[*around_hii_indices[i]])
  453 + med_G0_in[i]=la_median(G0_map[*hii_indices[i]])
  454 + med_G0_out[i]=la_median(G0_map[*around_hii_indices[i]])
  455 + med_YG1_in[i]=la_median(YG1_map[*hii_indices[i]])
  456 + med_YG1_out[i]=la_median(YG1_map[*around_hii_indices[i]])
  457 + med_YG2_in[i]=la_median(YG2_map[*hii_indices[i]])
  458 + med_YG2_out[i]=la_median(YG2_map[*around_hii_indices[i]])
  459 + med_YG3_in[i]=la_median(YG3_map[*hii_indices[i]])
  460 + med_YG3_out[i]=la_median(YG3_map[*around_hii_indices[i]])
  461 + med_Xion_in[i]=la_median(Xion_map[*hii_indices[i]])
  462 + med_Xion_out[i]=la_median(Xion_map[*around_hii_indices[i]])
  463 + med_rchi2_in[i]=la_median(rchi2_map[*hii_indices[i]])
  464 + med_rchi2_out[i]=la_median(rchi2_map[*around_hii_indices[i]])
  465 +ENDFOR
  466 +window,win,xsize=wxs,ysize=wys & win=win+1
  467 +!p.multi=[0,3,3]
  468 +ind=where(med_NH_in NE la_undef() and med_NH_out NE la_undef(),count)
  469 +vv=(med_NH_in[ind]-med_NH_out[ind])/med_NH_out[ind]*100.
  470 +xrange=[1,max(hii_count)]
  471 +cgplot,hii_count[ind],vv,tit='NH',xrange=xrange,/xlog & cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
  472 +hii_stats.med_NH=la_median(vv) & hii_stats.std_NH=hii_stats.med_NH/la_sigma(vv)
  473 +ind=where(med_G0_in NE la_undef() and med_G0_out NE la_undef(),count)
  474 +vv=(med_G0_in[ind]-med_G0_out[ind])/med_G0_out[ind]*100.
  475 +cgplot,hii_count[ind],vv,tit='G0',xrange=xrange,/xlog & cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
  476 +hii_stats.med_G0=la_median(vv) & hii_stats.std_G0=hii_stats.med_G0/la_sigma(vv)
  477 +ind=where(med_YG1_in NE la_undef() and med_YG1_out NE la_undef(),count)
  478 +vv=(med_YG1_in[ind]-med_YG1_out[ind])/med_YG1_out[ind]*100.
  479 +cgplot,hii_count[ind],vv,tit='YG1',xrange=xrange,/xlog & cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
  480 +hii_stats.med_YG1=la_median(vv)
  481 +ind=where(med_YG2_in NE la_undef() and med_YG2_out NE la_undef(),count)
  482 +vv=(med_YG2_in[ind]-med_YG2_out[ind])/med_YG2_out[ind]*100.
  483 +cgplot,hii_count[ind],vv,tit='YG2',xrange=xrange,/xlog & cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
  484 +hii_stats.med_YG2=la_median(vv)
  485 +ind=where(med_YG3_in NE la_undef() and med_YG3_out NE la_undef(),count)
  486 +vv=(med_YG3_in[ind]-med_YG3_out[ind])/med_YG3_out[ind]*100.
  487 +cgplot,hii_count[ind],vv,tit='YG3',xrange=xrange,/xlog & cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
  488 +hii_stats.med_YG3=la_median(vv)
  489 +ind=where(med_Xion_in NE la_undef() and med_Xion_out NE la_undef(),count)
  490 +vv=(med_Xion_in[ind]-med_Xion_out[ind])/med_Xion_out[ind]*100.
  491 +cgplot,hii_count[ind],vv,tit='Xion',xrange=xrange,/xlog & cgoplot,10^!x.crange,[0,0.],linestyle=2,color='red'
  492 +hii_stats.med_Xion=la_median(vv)
  493 +ind=where(med_rchi2_in NE la_undef() and med_rchi2_out NE la_undef(),count)
  494 +vv=(med_rchi2_in[ind]-med_rchi2_out[ind])/med_rchi2_out[ind]*100.
  495 +cgplot,hii_count[ind],vv,tit='rchi2',xrange=xrange,/xlog,/ylog & cgoplot,10^!x.crange,[1,1.],linestyle=2,color='red'
  496 +hii_stats.med_rchi2=la_median(vv)
  497 +
  498 +write_xcat,hii_stats,'/tmp/bidon',/nowrite
  499 +
  500 +cleanplot
  501 +
  502 +stop
  503 +
  504 +;===== plot nhits vs isrf_class
  505 +Ncl=n_elements(stats_st)
  506 +F0_hitp_per_class=lonarr(Ncl)
  507 +F0_hitm_per_class=lonarr(Ncl)
  508 +Ygrain1_hitp_per_class=lonarr(Ncl)
  509 +Ygrain1_hitm_per_class=lonarr(Ncl)
  510 +Ygrain2_hitp_per_class=lonarr(Ncl)
  511 +Ygrain2_hitm_per_class=lonarr(Ncl)
  512 +Ygrain3_hitp_per_class=lonarr(Ncl)
  513 +Ygrain3_hitm_per_class=lonarr(Ncl)
  514 +FOR i=0L,Ncl-1 DO BEGIN
  515 + ind=where(vor_class EQ stats_st[i].isrf_class AND F0_hit EQ 1,count)
  516 + F0_hitp_per_class[i]=count
  517 + ind=where(vor_class EQ stats_st[i].isrf_class AND F0_hit EQ -1,count)
  518 + F0_hitm_per_class[i]=count
  519 +
  520 + ind=where(vor_class EQ stats_st[i].isrf_class AND Ygrain1_hit EQ 1,count)
  521 + Ygrain1_hitp_per_class[i]=count
  522 + ind=where(vor_class EQ stats_st[i].isrf_class AND Ygrain1_hit EQ -1,count)
  523 + Ygrain1_hitm_per_class[i]=count
  524 +
  525 + ind=where(vor_class EQ stats_st[i].isrf_class AND Ygrain2_hit EQ 1,count)
  526 + Ygrain2_hitp_per_class[i]=count
  527 + ind=where(vor_class EQ stats_st[i].isrf_class AND Ygrain2_hit EQ -1,count)
  528 + Ygrain2_hitm_per_class[i]=count
  529 +
  530 + ind=where(vor_class EQ stats_st[i].isrf_class AND Ygrain3_hit EQ 1,count)
  531 + Ygrain3_hitp_per_class[i]=count
  532 + ind=where(vor_class EQ stats_st[i].isrf_class AND Ygrain3_hit EQ -1,count)
  533 + Ygrain3_hitm_per_class[i]=count
  534 +ENDFOR
  535 +window,win & win=win+1
  536 +!p.multi=[0,2,2]
  537 +xtit='ISRF class'
  538 +ytit='N bins'
  539 +tit='F0_hit'
  540 +cgplot,stats_st.isrf_class,F0_hitp_per_class,color='red',title=tit,xtit=xtit,ytit=ytit
  541 +cgoplot,stats_st.isrf_class,F0_hitm_per_class,color='blue'
  542 +tit='Ygrain1_hit'
  543 +cgplot,stats_st.isrf_class,Ygrain1_hitp_per_class,color='red',title=tit,xtit=xtit,ytit=ytit
  544 +cgoplot,stats_st.isrf_class,Ygrain1_hitm_per_class,color='blue'
  545 +tit='Ygrain2_hit'
  546 +cgplot,stats_st.isrf_class,Ygrain2_hitp_per_class,color='red',title=tit,xtit=xtit,ytit=ytit
  547 +cgoplot,stats_st.isrf_class,Ygrain2_hitm_per_class,color='blue'
  548 +tit='Ygrain3_hit'
  549 +cgplot,stats_st.isrf_class,Ygrain3_hitp_per_class,color='red',title=tit,xtit=xtit,ytit=ytit
  550 +cgoplot,stats_st.isrf_class,Ygrain3_hitm_per_class,color='blue'
  551 +
  552 +window,win & win=win+1
  553 +;==== plot hits vs F0hit
  554 +!p.multi=[0,2,2]
  555 +xtit='F0_hit'
  556 +ytit='grain hit'
  557 +tit='Ygrain1_hit'
  558 +xr=[-1.1,1.1]
  559 +cgplot,F0_hitp_per_class,Ygrain1_hitp_per_class,color='red',title=tit,xtit=xtit,ytit=ytit,psym='Open Circle',xrange=xrange,/xsty
  560 +cgoplot,F0_hitp_per_class,Ygrain1_hitm_per_class,color='blue',psym='Open Circle'
  561 +tit='Ygrain2_hit'
  562 +cgplot,F0_hitp_per_class,Ygrain2_hitp_per_class,color='red',title=tit,xtit=xtit,ytit=ytit,psym='Open Circle',xrange=xrange,/xsty
  563 +cgoplot,F0_hitp_per_class,Ygrain2_hitm_per_class,color='blue',psym='Open Circle'
  564 +tit='Ygrain3_hit'
  565 +cgplot,F0_hitp_per_class,Ygrain3_hitp_per_class,color='red',title=tit,xtit=xtit,ytit=ytit,psym='Open Circle',xrange=xrange,/xsty
  566 +cgoplot,F0_hitp_per_class,Ygrain3_hitm_per_class,color='blue',psym='Open Circle'
  567 +
  568 +;stop
  569 +
  570 +;===== histograms
  571 +window,win & win=win+1
  572 +ind=where(vor_class NE la_undef() and vor_class NE 0)
  573 +res=histogram(vor_class[ind],locations=xv)
  574 +cgplot,xv,res,psym=10,title='ISRF class histogram',xtit='class',ytit='Number',/ylog,yrange=[1,max(res)],/ysty
  575 +
  576 +;===== histograms 4 pannels
  577 +range_fact=15.
  578 +F0_range_fact=10. ;%
  579 +
  580 +window,win,xsize=wxs,ysize=wys & win=win+1
  581 +!p.multi=[0,2,2]
  582 +
  583 +ind=where(F0s NE la_undef() AND F0s GT 0,count)
  584 +res=histogram(alog10(F0s[ind]),locations=xv,Nbins=100)
  585 +yrange=[1,max(res)]
  586 +F0_cent=1.
  587 +parval_min=((*!dustem_grid).PMIN_VALUES)[0]
  588 +parval_max=((*!dustem_grid).PMAX_VALUES)[0]
  589 +;xrange=alog10([F0_min/F0_range_fact,F0_max*F0_range_fact])
  590 +xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
  591 +cgplot,xv,res,psym=10,title='F0s histogram',xtit='log10(F0)',ytit='Number',/ylog,yrange=yrange,/ysty,xrange=xrange,/xsty
  592 +cgoplot,replicate(alog10(F0_cent),2),10^(!y.crange),linestyle=2,color='black'
  593 +cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
  594 +cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  595 +
  596 +ind=where(Ygrain1 NE la_undef())
  597 +default_value=dustem_get_param_values_default('(*!dustem_params).grains[0].MDUST_O_MH')
  598 +parval_min=((*!dustem_grid).PMIN_VALUES)[1]
  599 +parval_max=((*!dustem_grid).PMAX_VALUES)[1]
  600 +;xrange=alog10([default_value/range_fact,default_value*range_fact])
  601 +xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
  602 +yrange=[1,max(res)]
  603 +res=histogram(alog10(Ygrain1[ind]),locations=xv,Nbins=100)
  604 +cgplot,xv,res,psym=10,title='Y Grain1 histogram',xtit='log10(Y)',ytit='Number',/ylog,/ysty,yrange=yrange,xrange=xrange,/xsty
  605 +cgoplot,replicate(alog10(default_value),2),10^(!y.crange),linestyle=2,color='black'
  606 +cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
  607 +cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  608 +
  609 +ind=where(Ygrain2 NE la_undef(),count)
  610 +IF count NE 0 THEN BEGIN
  611 + default_value=dustem_get_param_values_default('(*!dustem_params).grains[1].MDUST_O_MH')
  612 + parval_min=((*!dustem_grid).PMIN_VALUES)[2]
  613 + parval_max=((*!dustem_grid).PMAX_VALUES)[2]
  614 + ;xrange=alog10([default_value/range_fact,default_value*range_fact])
  615 + xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
  616 + yrange=[1,max(res)]
  617 + res=histogram(alog10(Ygrain2[ind]),locations=xv,Nbins=100)
  618 + cgplot,xv,res,psym=10,title='Y Grain2 histogram',xtit='log10(Y)',ytit='Number',/ylog,/ysty,yrange=yrange,xrange=xrange,/xsty
  619 + cgoplot,replicate(alog10(default_value),2),10^(!y.crange),linestyle=2,color='black'
  620 + cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
  621 + cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  622 +ENDIF
  623 +
  624 +ind=where(Ygrain3 NE la_undef(),count3)
  625 +IF count3 NE 0 THEN BEGIN
  626 + default_value=dustem_get_param_values_default('(*!dustem_params).grains[2].MDUST_O_MH')
  627 + ;=== IF PAH+ are in table
  628 + parval_min=((*!dustem_grid).PMIN_VALUES)[3]
  629 + parval_max=((*!dustem_grid).PMAX_VALUES)[3]
  630 + ;=== IF PAH+ are not in table
  631 + ;parval_min=((*!dustem_grid).PMIN_VALUES)[2]
  632 + ;parval_max=((*!dustem_grid).PMAX_VALUES)[2]
  633 + ;xrange=alog10([default_value/range_fact,default_value*range_fact])
  634 + xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
  635 + yrange=[1,max(res)]
  636 + res=histogram(alog10(Ygrain3[ind]),locations=xv,Nbins=100)
  637 + cgplot,xv,res,psym=10,title='Y Grain3 histogram',xtit='log10(Y)',ytit='Number',/ylog,/ysty,yrange=yrange,xrange=xrange,/xsty
  638 + cgoplot,replicate(alog10(default_value),2),10^(!y.crange),linestyle=2,color='black'
  639 + cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
  640 + cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  641 +ENDIF
  642 +
  643 +IF keyword_set(do_postscripts) THEN BEGIN
  644 + cleanplot
  645 +
  646 + ps_file=ps_dir+use_source_name+'_histo_F0'+reso_str+fit_G0_str+mathis_str+'.ps'
  647 + set_plot,'PS'
  648 + device,filename=ps_file
  649 +
  650 + ind=where(F0s NE la_undef() AND F0s GT 0,count)
  651 + res=histogram(alog10(F0s[ind]),locations=xv,Nbins=100)
  652 + yrange=[1,max(res)]
  653 + F0_cent=1.
  654 + parval_min=((*!dustem_grid).PMIN_VALUES)[0]
  655 + parval_max=((*!dustem_grid).PMAX_VALUES)[0]
  656 + ;xrange=alog10([F0_min/F0_range_fact,F0_max*F0_range_fact])
  657 + xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
  658 + cgplot,xv,res,psym=10,title='F0s histogram',xtit='log10(F0)',ytit='Number',/ylog,yrange=yrange,/ysty,xrange=xrange,/xsty
  659 + cgoplot,replicate(alog10(F0_cent),2),10^(!y.crange),linestyle=2,color='black'
  660 + cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
  661 + cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  662 +
  663 + device,/close
  664 + set_plot,'X'
  665 + message,'Wrote '+ps_file,/continue
  666 + ;stop
  667 +
  668 + ps_file=ps_dir+use_source_name+'_histo_YG1'+reso_str+fit_G0_str+mathis_str+'.ps'
  669 + set_plot,'PS'
  670 + device, filename=ps_file
  671 +
  672 + ind=where(Ygrain1 NE la_undef())
  673 + default_value=dustem_get_param_values_default('(*!dustem_params).grains[0].MDUST_O_MH')
  674 + parval_min=((*!dustem_grid).PMIN_VALUES)[1]
  675 + parval_max=((*!dustem_grid).PMAX_VALUES)[1]
  676 + ;xrange=alog10([default_value/range_fact,default_value*range_fact])
  677 + xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
  678 + yrange=[1,max(res)]
  679 + res=histogram(alog10(Ygrain1[ind]),locations=xv,Nbins=100)
  680 + cgplot,xv,res,psym=10,title='Y Grain1 histogram',xtit='log10(Y)',ytit='Number',/ylog,/ysty,yrange=yrange,xrange=xrange,/xsty
  681 + cgoplot,replicate(alog10(default_value),2),10^(!y.crange),linestyle=2,color='black'
  682 + cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
  683 + cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  684 +
  685 + device,/close
  686 + set_plot,'X'
  687 +
  688 + ps_file=ps_dir+use_source_name+'_histo_YG2'+reso_str+fit_G0_str+mathis_str+'.ps'
  689 + set_plot,'PS'
  690 + device, filename=ps_file
  691 +
  692 + ind=where(Ygrain2 NE la_undef())
  693 + default_value=dustem_get_param_values_default('(*!dustem_params).grains[1].MDUST_O_MH')
  694 + parval_min=((*!dustem_grid).PMIN_VALUES)[2]
  695 + parval_max=((*!dustem_grid).PMAX_VALUES)[2]
  696 + ;xrange=alog10([default_value/range_fact,default_value*range_fact])
  697 + xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
  698 + yrange=[1,max(res)]
  699 + res=histogram(alog10(Ygrain2[ind]),locations=xv,Nbins=100)
  700 + cgplot,xv,res,psym=10,title='Y Grain2 histogram',xtit='log10(Y)',ytit='Number',/ylog,/ysty,yrange=yrange,xrange=xrange,/xsty
  701 + cgoplot,replicate(alog10(default_value),2),10^(!y.crange),linestyle=2,color='black'
  702 + cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
  703 + cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  704 +
  705 + device,/close
  706 + set_plot,'X'
  707 +
  708 + ps_file=ps_dir+use_source_name+'_histo_YG3'+reso_str+fit_G0_str+mathis_str+'.ps'
  709 + set_plot,'PS'
  710 + device, filename=ps_file
  711 +
  712 + ind=where(Ygrain3 NE la_undef(),count3)
  713 + IF count3 NE 0 THEN BEGIN
  714 + default_value=dustem_get_param_values_default('(*!dustem_params).grains[2].MDUST_O_MH')
  715 + parval_min=((*!dustem_grid).PMIN_VALUES)[3]
  716 + parval_max=((*!dustem_grid).PMAX_VALUES)[3]
  717 + ;xrange=alog10([default_value/range_fact,default_value*range_fact])
  718 + xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
  719 + yrange=[1,max(res)]
  720 + res=histogram(alog10(Ygrain3[ind]),locations=xv,Nbins=100)
  721 + cgplot,xv,res,psym=10,title='Y Grain3 histogram',xtit='log10(Y)',ytit='Number',/ylog,/ysty,yrange=yrange,xrange=xrange,/xsty
  722 + cgoplot,replicate(alog10(default_value),2),10^(!y.crange),linestyle=2,color='black'
  723 + cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
  724 + cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  725 + ENDIF
  726 +
  727 + device,/close
  728 + set_plot,'X'
  729 +
  730 +ENDIF
  731 +;stop
  732 +
  733 +;==== show maps
  734 +obp=[1.1,0.,1.15,1]
  735 +cleanplot
  736 +image_color_table='jpbloadct,/whitebackground'
  737 +
  738 +window,win,xsize=wxs,ysize=wys & win=win+1
  739 +CASE use_source_name OF
  740 + 'ngc0628':imrange=[-1,1.5]
  741 + 'ngc4321':imrange=[-1,2.]
  742 + 'ngc3351':imrange=[-1,1.5]
  743 + 'ngc1566':imrange=[-1,1.5]
  744 +ENDCASE
  745 +title='NH map'
  746 +image_cont20,la_log10(NHmap),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st
  747 +;stop
  748 +IF keyword_set(do_postscripts) THEN BEGIN
  749 + ps_file=ps_dir+use_source_name+'_NHmap'+reso_str+fit_G0_str+mathis_str+'.ps'
  750 + image_cont20,la_log10(NHmap),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st, $
  751 + postscript=ps_file,orientation='portrait'
  752 + image_cont20,la_log10(NHmap),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st, $
  753 + postscript=ps_file,orientation='portrait'
  754 + message,'Wrote '+ps_file,/continue
  755 +ENDIF
  756 +;stop
  757 +
  758 +window,win,xsize=wxs,ysize=wys & win=win+1
  759 +CASE use_source_name OF
  760 + 'ngc0628':imrange=1.+[-0.5,+0.5]
  761 + 'ngc4321':imrange=1.+[-0.5,+0.5]
  762 + 'ngc3351':imrange=2.+[-0.5,+0.5]
  763 + 'ngc1566':imrange=1.+[-0.5,+0.5]
  764 +ENDCASE
  765 +;imrange=1.+[-0.5,+0.5]
  766 +title='stellar normalization'
  767 +image_cont20,norm_stars_map,Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st
  768 +IF keyword_set(do_postscripts) THEN BEGIN
  769 + ps_file=ps_dir+use_source_name+'_norm_stars_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  770 + image_cont20,norm_stars_map,Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st, $
  771 + postscript=ps_file,orientation='portrait',/nologo
  772 + message,'Wrote '+ps_file,/continue
  773 +ENDIF
  774 +
  775 +window,win,xsize=wxs,ysize=wys & win=win+1
  776 +CASE use_source_name OF
  777 + 'ngc0628':imrange=[0.1,0.25]
  778 + 'ngc4321':imrange=[0.1,0.25]
  779 + 'ngc3351':imrange=[0.1,0.6]
  780 + 'ngc1566':imrange=[0.1,0.25]
  781 +ENDCASE
  782 +title='E(B-V)'
  783 +image_cont20,EBV_map,Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st
  784 +IF keyword_set(do_postscripts) THEN BEGIN
  785 + ps_file=ps_dir+use_source_name+'_EBV_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  786 + image_cont20,EBV_map,Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st, $
  787 + postscript=ps_file,orientation='portrait',/nologo
  788 + message,'Wrote '+ps_file,/continue
  789 +ENDIF
  790 +;stop
  791 +
  792 +window,win,xsize=wxs,ysize=wys & win=win+1
  793 +imrange=[-1.5,1.2]
  794 +title='log10(F0_predicted)'
  795 +image_cont20,la_log10(F0_predicted_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st
  796 +IF keyword_set(do_postscripts) THEN BEGIN
  797 + ps_file=ps_dir+use_source_name+'_F0_predicted_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  798 + image_cont20,la_log10(F0_predicted_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st, $
  799 + postscript=ps_file,orientation='portrait',/nologo
  800 + message,'Wrote '+ps_file,/continue
  801 +ENDIF
  802 +
  803 +window,win,xsize=wxs,ysize=wys & win=win+1
  804 +imrange=[-1.,+1.0] ;with Herschel data
  805 +imrange=[-1.,+1.5] ;with Herschel data
  806 +title='log10(F0)'
  807 +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
  808 +;stop
  809 +IF keyword_set(do_postscripts) THEN BEGIN
  810 + ps_file=ps_dir+use_source_name+'_F0_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  811 + 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, $
  812 + postscript=ps_file,orientation='portrait',/nologo
  813 + message,'Wrote '+ps_file,/continue
  814 +ENDIF
  815 +
  816 +window,win,xsize=wxs,ysize=wys & win=win+1
  817 +imrange=[-1.,+1.0]
  818 +title='log10(G0)'
  819 +image_cont20,la_log10(G0_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st
  820 +IF keyword_set(do_postscripts) THEN BEGIN
  821 + ps_file=ps_dir+use_source_name+'_G0_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  822 + image_cont20,la_log10(G0_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st, $
  823 + postscript=ps_file,orientation='portrait',/nologo
  824 + message,'Wrote '+ps_file,/continue
  825 +ENDIF
  826 +;stop
  827 +
  828 +window,win,xsize=wxs,ysize=wys & win=win+1
  829 +imrange=[-1.,+1.0]
  830 +title='log10(G0)'
  831 +image_cont20,la_log10(G0_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=hii_st
  832 +;stop
  833 +IF keyword_set(do_postscripts) THEN BEGIN
  834 + ps_file=ps_dir+use_source_name+'_G0_map_hiicont'+reso_str+fit_G0_str+mathis_str+'.ps'
  835 + image_cont20,la_log10(G0_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=hii_st, $
  836 + postscript=ps_file,orientation='portrait',/nologo
  837 + message,'Wrote '+ps_file,/continue
  838 +ENDIF
  839 +;image_cont20,hii,Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=hii_st
  840 +
  841 +window,win,xsize=wxs,ysize=wys & win=win+1
  842 +imrange=[-1.,+1.0]
  843 +title='log10(G0_predicted)'
  844 +image_cont20,la_log10(G0_predicted_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st
  845 +IF keyword_set(do_postscripts) THEN BEGIN
  846 + ps_file=ps_dir+use_source_name+'_G0_predicted_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  847 + image_cont20,la_log10(G0_predicted_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st, $
  848 + postscript=ps_file,orientation='portrait',/nologo
  849 + message,'Wrote '+ps_file,/continue
  850 +ENDIF
  851 +
  852 +window,win,xsize=wxs,ysize=wys & win=win+1
  853 +CASE use_source_name OF
  854 + 'ngc0628':imrange=[-3.8,-2.9]
  855 + 'ngc4321':imrange=[-3.8,-2.9]
  856 + 'ngc3351':imrange=[-3.8,-2.9]-0.3
  857 + 'ngc1566':imrange=[-3.8,-2.9]
  858 +ENDCASE
  859 +title='log10(YG1)'
  860 +image_cont20,la_log10(Yg1_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st
  861 +IF keyword_set(do_postscripts) THEN BEGIN
  862 + ps_file=ps_dir+use_source_name+'_YG1_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  863 + image_cont20,la_log10(Yg1_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st, $
  864 + postscript=ps_file,orientation='portrait',/nologo
  865 + message,'Wrote '+ps_file,/continue
  866 +ENDIF
  867 +title='log10(YG1)'
  868 +image_cont20,la_log10(Yg1_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=hii_st
  869 +IF keyword_set(do_postscripts) THEN BEGIN
  870 + ps_file=ps_dir+use_source_name+'_YG1_map_hiicont'+reso_str+fit_G0_str+mathis_str+'.ps'
  871 + image_cont20,la_log10(Yg1_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=hii_st, $
  872 + postscript=ps_file,orientation='portrait',/nologo
  873 + message,'Wrote '+ps_file,/continue
  874 +ENDIF
  875 +stop
  876 +
  877 +window,win,xsize=wxs,ysize=wys & win=win+1
  878 +CASE use_source_name OF
  879 + 'ngc0628':imrange=[-3.8,-2.9]
  880 + 'ngc4321':imrange=[-3.8,-2.9]
  881 + 'ngc3351':imrange=[-3.8,-2.9]-0.3
  882 + 'ngc1566':imrange=[-3.8,-2.9]
  883 +ENDCASE
  884 +title='log10(YG2)'
  885 +image_cont20,la_log10(Yg2_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st
  886 +IF keyword_set(do_postscripts) THEN BEGIN
  887 + ps_file=ps_dir+use_source_name+'_YG2_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  888 + image_cont20,la_log10(Yg2_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st, $
  889 + postscript=ps_file,orientation='portrait',/nologo
  890 + message,'Wrote '+ps_file,/continue
  891 +ENDIF
  892 +;stop
  893 +
  894 +window,win,xsize=wxs,ysize=wys & win=win+1
  895 +;imrange=[-4,-2.]
  896 +imrange=[-3.8,-2.5]
  897 +title='log10(YG3)'
  898 +image_cont20,la_log10(Yg3_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st
  899 +IF keyword_set(do_postscripts) THEN BEGIN
  900 + ps_file=ps_dir+use_source_name+'_YG3_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  901 + image_cont20,la_log10(Yg3_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st, $
  902 + postscript=ps_file,orientation='portrait',/nologo
  903 + message,'Wrote '+ps_file,/continue
  904 +ENDIF
  905 +
  906 +;window,win,xsize=800,ysize=900 & win=win+1
  907 +;imrange=[1.5,2.2]
  908 +;image_cont20,la_log10(chi2_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table='jpbloadct',/silent,title='log10(Chi2)',off_bar=obp
  909 +
  910 +window,win,xsize=wxs,ysize=wys & win=win+1
  911 +;imrange=[1.2,1.8]
  912 +;imrange=[-0.1,2.]
  913 +imrange=[-1.,2.]
  914 +title='log10(rChi2/med(rchi2))'
  915 +map=rchi2_map
  916 +;ind=where(rchi2_map EQ la_min(rchi2_map))
  917 +;map=la_div(map,rchi2_map[ind[0]])
  918 +map=la_div(map,la_median(rchi2_map))
  919 +image_cont20,la_log10(map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st
  920 +IF keyword_set(do_postscripts) THEN BEGIN
  921 + ps_file=ps_dir+use_source_name+'_rchi2_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  922 + image_cont20,la_log10(map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,off_bar=obp,title=title,levels=cont_st, $
  923 + postscript=ps_file,orientation='portrait',/nologo
  924 + message,'Wrote '+ps_file,/continue
  925 +ENDIF
  926 +;stop
  927 +
  928 +;window,win,xsize=800,ysize=900 & win=win+1
  929 +;imrange=[-3,2] ;normalize=0
  930 +;image_cont20,la_log10(fact_map),Href,/square,imrange=imrange,axis_color_table=1,image_color_table='jpbloadct',/silent,title='log10(fact)',off_bar=obp
  931 +
  932 +;stop
  933 +window,win,xsize=wxs,ysize=wys & win=win+1
  934 +imrange=[-2,2]
  935 +title='F0 hit'
  936 +image_cont20,F0_hit_map,Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,title=title,off_bar=obp,levels=cont_st
  937 +IF keyword_set(do_postscripts) THEN BEGIN
  938 + ps_file=ps_dir+use_source_name+'_F0_hit_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  939 + image_cont20,F0_hit_map,Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,title=title,off_bar=obp,levels=cont_st, $
  940 + postscript=ps_file,orientation='portrait',/nologo
  941 + message,'Wrote '+ps_file,/continue
  942 +ENDIF
  943 +
  944 +window,win,xsize=wxs,ysize=wys & win=win+1
  945 +imrange=[-2,2]
  946 +title='YG1 hit'
  947 +image_cont20,Yg1_hit_map,Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,title=title,off_bar=obp,levels=cont_st
  948 +IF keyword_set(do_postscripts) THEN BEGIN
  949 + ps_file=ps_dir+use_source_name+'_YG1_hit_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  950 + image_cont20,Yg1_hit_map,Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,title=title,off_bar=obp,levels=cont_st, $
  951 + postscript=ps_file,orientation='portrait',/nologo
  952 + message,'Wrote '+ps_file,/continue
  953 +ENDIF
  954 +
  955 +window,win,xsize=wxs,ysize=wys & win=win+1
  956 +imrange=[-2,2]
  957 +title='YG2 hit'
  958 +image_cont20,Yg2_hit_map,Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,title=title,off_bar=obp,levels=cont_st
  959 +IF keyword_set(do_postscripts) THEN BEGIN
  960 + ps_file=ps_dir+use_source_name+'_YG2_hit_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  961 + image_cont20,Yg2_hit_map,Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,title=title,off_bar=obp,levels=cont_st, $
  962 + postscript=ps_file,orientation='portrait',/nologo
  963 + message,'Wrote '+ps_file,/continue
  964 +ENDIF
  965 +
  966 +window,win,xsize=wxs,ysize=wys & win=win+1
  967 +imrange=[-2,2]
  968 +title='YG3 hit'
  969 +image_cont20,Yg3_hit_map,Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,title=title,off_bar=obp,levels=cont_st
  970 +IF keyword_set(do_postscripts) THEN BEGIN
  971 + ps_file=ps_dir+use_source_name+'_YG3_hit_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  972 + image_cont20,Yg3_hit_map,Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,title=title,off_bar=obp,levels=cont_st, $
  973 + postscript=ps_file,orientation='portrait',/nologo
  974 + message,'Wrote '+ps_file,/continue
  975 +ENDIF
  976 +
  977 +;window,win,xsize=800,ysize=900 & win=win+1
  978 +;imrange=0
  979 +;imrange=[-2.0,-1.45]
  980 +;image_cont20,la_log10(la_mul(Ypah_map,fact_map)),Href,/square,imrange=imrange,axis_color_table=1,image_color_table='jpbloadct',/silent,title='Ypah *fact',off_bar=obp
  981 +
  982 +;stop
  983 +window,win,xsize=wxs,ysize=wys & win=win+1
  984 +CASE use_source_name OF
  985 + 'ngc0628':imrange=[0.5-0.2,0.5+0.2]
  986 + 'ngc4321':imrange=[0.5-0.2,0.5+0.2]
  987 + 'ngc3351':imrange=[0.5-0.05,0.5+0.05]
  988 + 'ngc1566':imrange=[0.5-0.2,0.5+0.2]
  989 +ENDCASE
  990 +tit='Xion'
  991 +image_cont20,Xion_map,Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,title=tit,off_bar=obp,levels=cont_st
  992 +IF keyword_set(do_postscripts) THEN BEGIN
  993 + ps_file=ps_dir+use_source_name+'_Xion_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  994 + image_cont20,Xion_map,Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,title=tit,off_bar=obp,levels=cont_st, $
  995 + postscript=ps_file,orientation='portrait',/nologo
  996 + message,'Wrote '+ps_file,/continue
  997 +ENDIF
  998 +image_cont20,Xion_map,Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,title=tit,off_bar=obp,levels=hii_st
  999 +IF keyword_set(do_postscripts) THEN BEGIN
  1000 + ps_file=ps_dir+use_source_name+'_Xion_map_hiicont'+reso_str+fit_G0_str+mathis_str+'.ps'
  1001 + image_cont20,Xion_map,Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,title=tit,off_bar=obp,levels=hii_st, $
  1002 + postscript=ps_file,orientation='portrait',/nologo
  1003 + message,'Wrote '+ps_file,/continue
  1004 +ENDIF
  1005 +stop
  1006 +
  1007 +window,win,xsize=wxs,ysize=wys & win=win+1
  1008 +CASE use_source_name OF
  1009 + 'ngc0628':imrange=[-0.1,1.5]
  1010 + 'ngc4321':imrange=[-1.0,1.5]
  1011 + 'ngc3351':imrange=[-0.1,1.5]
  1012 + 'ngc1566':imrange=[-0.1,1.5]
  1013 +ENDCASE
  1014 +rap_F0=la_div(F0_map,F0_predicted_map)
  1015 +tit='log(G0/G0_predicted)'
  1016 +image_cont20,la_log10(rap_F0),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,title=tit,off_bar=obp,levels=cont_st
  1017 +IF keyword_set(do_postscripts) THEN BEGIN
  1018 + ps_file=ps_dir+use_source_name+'_G0_ratio_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  1019 + image_cont20,la_log10(rap_F0),Href,/square,imrange=imrange,axis_color_table=1,image_color_table=image_color_table,/silent,title=tit,off_bar=obp,levels=cont_st, $
  1020 + postscript=ps_file,orientation='portrait',/nologo
  1021 + message,'Wrote '+ps_file,/continue
  1022 +ENDIF
  1023 +;stop
  1024 +
  1025 +;==== Show residual maps
  1026 +window,win,xsize=wxs,ysize=wys & win=win+1
  1027 +residual_range=[-100,100.]/2.
  1028 +FOR i=0L,n_elements(all_filters)-1 DO BEGIN
  1029 + title='Residual percent Filter '+all_filters[i]+' '+strtrim(dustem_filter2wav(all_filters[i]),2)
  1030 + filter_str='_'+all_filters[i]
  1031 + map=residual_maps_perc[*,*,i]
  1032 + ind=where(finite(map) NE 1,count)
  1033 + IF count NE 0 THEN map[ind]=la_undef()
  1034 + ;print,la_min(map),la_max(map)
  1035 + IF la_min(map) NE la_max(map) THEN use_obp=obp ELSE use_obp=0
  1036 + image_cont20,residual_maps_perc[*,*,i],href,/square,/silent,off_bar_pos=use_obp,image_color_table=image_color_table,axis_color_table=1,title=title,levels=cont_st,imrange=residual_range
  1037 + IF keyword_set(do_postscripts) THEN BEGIN
  1038 + ps_file=ps_dir+use_source_name+'_residual_map'+reso_str+fit_G0_str+mathis_str+filter_str+'.ps'
  1039 + image_cont20,residual_maps_perc[*,*,i],href,/square,/silent,off_bar_pos=use_obp,image_color_table=image_color_table,axis_color_table=1,title=title,levels=cont_st,imrange=residual_range, $
  1040 + postscript=ps_file,orientation='portrait',/nologo
  1041 + message,'Wrote '+ps_file,/continue
  1042 + ENDIF
  1043 + IF not keyword_set(nostop) THEN stop
  1044 +ENDFOR
  1045 +
  1046 +;=== other plots
  1047 +NH_range=[0.1,1000.]
  1048 +F0_range=[0.1,10]
  1049 +F0_range_pred=[0.001,10]
  1050 +Ypah_range=[8.e-5,2.e-1]
  1051 +Yvsg_range=[8.e-5,2.e-1]
  1052 +chi2_range=[1.e3,1.e7]/1.e2
  1053 +
  1054 +window,win,xsize=wxs,ysize=wys & win=win+1
  1055 +cgplot,F0s_predicted,F0s,/xlog,/ylog,xrange=F0_range_pred,/xsty,yrange=F0_range,/ysty,psym=3,xtit='F0 predicted',ytit='F0'
  1056 +cgoplot,[1.e-3,1.e3],[1.e-3,1.e3],linesyle=2,color='red'
  1057 +cgoplot,[1.e-3,1.e3],[1.e-3,1.e3]*10.,linesyle=2,color='blue'
  1058 +cgoplot,[1.e-3,1.e3],[1.e-3,1.e3]*100.,linesyle=2,color='blue'
  1059 +
  1060 +window,win,xsize=wxs,ysize=wys & win=win+1
  1061 +!p.multi=[0,1,4]
  1062 +cgplot,F0s,Ygrain1,/xlog,/ylog,xr=F0_range,/xsty,psym=3,xtit='F0',ytit='Y',title='Grain 1'
  1063 +cgplot,F0s,Ygrain2,/xlog,/ylog,xr=F0_range,/xsty,psym=3,xtit='F0',ytit='Y',title='Grain 2'
  1064 +cgplot,F0s,Ygrain3,/xlog,/ylog,xr=F0_range,/xsty,psym=3,xtit='F0',ytit='Y',title='Grain 3'
  1065 +cgplot,Ygrain3,Ygrain1,/xlog,/ylog,xr=[1.e-4,1.e-1],yr=[1.e-4,1.e-1],/xsty,psym=3,xtit='Yvsg',ytit='Ypah'
  1066 +
  1067 +;stop
  1068 +
  1069 +window,win,xsize=wxs,ysize=wys & win=win+1
  1070 +!p.multi=[0,1,4]
  1071 +NH_range=[0.1,1000.]
  1072 +F0_range=[0.1,10]
  1073 +Ypah_range=[8.e-5,2.e-1]
  1074 +Yvsg_range=[8.e-5,2.e-1]
  1075 +cgplot,NH_values,F0s,/xlog,/ylog,xr=NH_range,yr=F0_range,/xsty,psym=3,xtit='NH',ytit='F0'
  1076 +cgplot,NH_values,Ygrain1,/xlog,/ylog,xr=NH_range,yr=Ypah_range,/xsty,psym=3,xtit='NH',ytit='Y',title='Grain 1'
  1077 +cgplot,NH_values,Ygrain2,/xlog,/ylog,xr=NH_range,yr=Ypah_range,/xsty,psym=3,xtit='NH',ytit='Y',title='Grain 2'
  1078 +cgplot,NH_values,Ygrain3,/xlog,/ylog,xr=NH_range,yr=Yvsg_range,/xsty,psym=3,xtit='NH',ytit='Y',title='Grain 3'
  1079 +
  1080 +window,win,xsize=wxs,ysize=wys & win=win+1
  1081 +!p.multi=[0,1,4]
  1082 +cgplot,chi2s,F0s,/xlog,/ylog,xr=chi2_range,yr=F0_range,/xsty,psym=3,xtit='chi2',ytit='F0'
  1083 +cgplot,chi2s,Ygrain1,/xlog,/ylog,xr=chi2_range,yr=Ypah_range,/xsty,psym=3,xtit='chi2',ytit='Y',title='Grain 1'
  1084 +cgplot,chi2s,Ygrain2,/xlog,/ylog,xr=chi2_range,yr=Ypah_range,/xsty,psym=3,xtit='chi2',ytit='Y',title='Grain 2'
  1085 +cgplot,chi2s,Ygrain3,/xlog,/ylog,xr=chi2_range,yr=Yvsg_range,/xsty,psym=3,xtit='chi2',ytit='Y',title='Grain 3'
  1086 +
  1087 +;stop
  1088 +Xions=la_div(Ygrain2,la_add(Ygrain1,Ygrain2))
  1089 +
  1090 +window,win,xsize=wxs,ysize=wys & win=win+1
  1091 +!p.multi=[0,1,4]
  1092 +G0_range=[0.1,500]
  1093 +cgplot,G0s,Ygrain1,/xlog,/ylog,xr=G0_range,/xsty,psym=3,xtit='G0',ytit='Y',title='Grain 1'
  1094 +cgplot,G0s,Ygrain2,/xlog,/ylog,xr=G0_range,/xsty,psym=3,xtit='G0',ytit='Y',title='Grain 2'
  1095 +cgplot,G0s,Ygrain3,/xlog,/ylog,xr=G0_range,/xsty,psym=3,xtit='G0',ytit='Y',title='Grain 3'
  1096 +cgplot,G0s,Xions,/xlog,/ylog,xr=G0_range,/xsty,psym=3,xtit='G0',ytit='Xion',title='Xion'
  1097 +
  1098 +window,win,xsize=wxs,ysize=wys & win=win+1
  1099 +Nx=100 & Ny=100
  1100 +xmin=min(NH_range)
  1101 +xmax=max(NH_range)
  1102 +ymin=min(G0_range)
  1103 +ymax=max(G0_range)
  1104 +logxbox=1
  1105 +logybox=1
  1106 +plot_bin_correl,G0s,NH_values,Nx,Ny,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,Nlevels=Nlevels, $
  1107 + fact=fact,logxbox=logxbox,logybox=logybox,range=range, $
  1108 + binned_map=binned_map,binned_xaxis=binned_xaxis,binned_yaxis=binned_yaxis,no_plot=no_plot,xtit='NH',ytit='G0', $
  1109 + /median_plot,upper_perc_plot=95.,lower_perc_plot=95.
  1110 +
  1111 +window,win,xsize=wxs,ysize=wys & win=win+1
  1112 +Xion_range=[0.1,0.9]
  1113 +;Nx=50 & Ny=50
  1114 +xmin=min(G0_range)
  1115 +xmax=max(G0_range)
  1116 +ymin=min(Xion_range)
  1117 +ymax=max(Xion_range)
  1118 +logxbox=1
  1119 +logybox=0
  1120 +plot_bin_correl,G0s,Xions,Nx,Ny,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,Nlevels=Nlevels, $
  1121 + fact=fact,logxbox=logxbox,logybox=logybox,range=range, $
  1122 + binned_map=binned_map,binned_xaxis=binned_xaxis,binned_yaxis=binned_yaxis,no_plot=no_plot,xtit='G0',ytit='Xion', $
  1123 + /median_plot,upper_perc_plot=95.,lower_perc_plot=95.
  1124 +
  1125 +window,win,xsize=wxs,ysize=wys & win=win+1
  1126 +Xion_range=[0.1,0.9]
  1127 +;Nx=50 & Ny=50
  1128 +xmin=min(NH_range)
  1129 +xmax=max(NH_range)
  1130 +ymin=min(Xion_range)
  1131 +ymax=max(Xion_range)
  1132 +logxbox=1
  1133 +logybox=0
  1134 +plot_bin_correl,NH_values,Xions,Nx,Ny,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,Nlevels=Nlevels, $
  1135 + fact=fact,logxbox=logxbox,logybox=logybox,range=range, $
  1136 + binned_map=binned_map,binned_xaxis=binned_xaxis,binned_yaxis=binned_yaxis,no_plot=no_plot,xtit='NH',ytit='Xion', $
  1137 + /median_plot,upper_perc_plot=95.,lower_perc_plot=95.
  1138 +
  1139 +goto,the_end
  1140 +
  1141 +the_end:
  1142 +
  1143 +END
0 1144 \ No newline at end of file
... ...
LabTools/IRAP/JPB/phangs_make_intensity_mask.pro 0 → 100644
  Diff comments   View file @f3236ab
... ... @@ -0,0 +1,82 @@
  1 +PRO phangs_make_intensity_mask,source_name=source_name,resolution_filter=resolution_filter,resolution_value=resolution_value,save=save,show_map=show_map
  2 +
  3 +;https://iopscience.iop.org/article/10.3847/1538-4357/ad54bd/pdf
  4 +;Les cutoff ont été appliqués sur trois filtres je crois : 0.09M Jy/sr pour F770W, 0.13 Mjy/sr pour F1130W et 0.29 MJy/sr pour F2100W.
  5 +
  6 +dustem_define_la_common
  7 +
  8 +brightness_thresholds_min=[0.09,0.13,0.29]
  9 +brightness_thresholds_filters=['F770W','F1130W','F2100W']
  10 +Nthrs=n_elements(brightness_thresholds_min)
  11 +
  12 +data_dir=!phangs_data_dir+'/ISRF/WORK/'
  13 +
  14 +use_source_name='ngc0628'
  15 +IF keyword_set(source_name) THEN use_source_name=source_name
  16 +
  17 +reso_str=''
  18 +IF keyword_set(resolution_filter) OR keyword_set(resolution_value) THEN BEGIN
  19 + IF keyword_set(resolution_filter) THEN BEGIN
  20 + reso_str='_'+resolution_filter
  21 + ENDIF ELSE BEGIN
  22 + reso_str='_'+strtrim(resolution_value,2)
  23 + ENDELSE
  24 +ENDIF
  25 +
  26 +file=data_dir+use_source_name+'_muse_data'+reso_str+'.sav'
  27 +st_info=file_info(file)
  28 +IF st_info.exists NE 1 THEN BEGIN
  29 + message,'Could not find '+file,/continue
  30 + stop
  31 +ENDIF
  32 +restore,file,/verb
  33 +;% RESTORE: Restored variable: ST_TEMPLATES.
  34 +;% RESTORE: Restored variable: ST_MUSE_WEIGHTS.
  35 +;% RESTORE: Restored variable: VORONOI_ID.
  36 +;% RESTORE: Restored variable: AGE_VALUES.
  37 +;% RESTORE: Restored variable: METALICITY_VALUES.
  38 +;% RESTORE: Restored variable: BINS.
  39 +;% RESTORE: Restored variable: HREF.
  40 +file=data_dir+use_source_name+'_jwst_images'+reso_str+'.sav'
  41 +st_info=file_info(file)
  42 +IF st_info.exists NE 1 THEN BEGIN
  43 + message,'Could not find '+file,/continue
  44 + stop
  45 +ENDIF
  46 +restore,file,/verb
  47 +;% RESTORE: Restored variable: JWST_IMAGES.
  48 +;% RESTORE: Restored variable: FILTERS.
  49 +;% RESTORE: Restored variable: HREF.
  50 +
  51 +;restore,file_save_indices,/verb ;just to get the all_seds_indices variable
  52 +
  53 +Nx=sxpar(href,'NAXIS1')
  54 +Ny=sxpar(href,'NAXIS2')
  55 +brightness_mask=intarr(Nx,Ny)
  56 +FOR i=0L,Nthrs-1 DO BEGIN
  57 + ind_filt=where(filters EQ brightness_thresholds_filters[i],count_filt)
  58 + im=jwst_images(*,*,0,ind_filt[0])
  59 + ind_remove=where(im LE brightness_thresholds_min[i],count)
  60 + IF count NE 0 THEN BEGIN
  61 + brightness_mask[ind_remove]=1
  62 + message,'Flagged '+strtrim(count,2)+' pixels on brightness in '+brightness_thresholds_filters[i],/continue
  63 + ENDIF
  64 + ind_indef=where(im EQ la_undef(),count)
  65 + IF count NE 0 THEN BEGIN
  66 + brightness_mask[ind_indef]=la_undef()
  67 + ENDIF
  68 +ENDFOR
  69 +
  70 +IF keyword_set(show_map) THEN BEGIN
  71 + image_cont20,brightness_mask,href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,off_bar_pos=[1.1,0,1.15,1],title=use_source_name
  72 +ENDIF
  73 +
  74 +IF keyword_set(save) THEN BEGIN
  75 + file_save=data_dir+use_source_name+'_brightness_mask'+reso_str+'.sav'
  76 + save,brightness_mask,file=file_save
  77 + message,'Wrote '+file_save,/continue
  78 +ENDIF
  79 +
  80 +stop
  81 +
  82 +END
0 83 \ No newline at end of file
... ...