OV-GSO-DC / dustem-wrapper_idl

Commit 9a95b8f382a9ec94201485ac46767e8613916e0f

Authored by Jean-Philippe Bernard
2 parents cd0dc454 f6389d8f
Exists in master

Merge branch 'master' of https://gitlab.irap.omp.eu/OV-GSO-DC/dustem-wrapper_idl

Showing 4 changed files with 1743 additions and 0 deletions   Show diff stats
Inline Side-by-side
LabTools/IRAP/RM/phangs_fits_seds_sizedist.pro 0 → 100644
  Diff comments   View file @9a95b8f
... ... @@ -0,0 +1,143 @@
  1 +PRO phangs_fits_seds_sizedist,show_sed=show_sed $
  2 + ,normalize=normalize $
  3 + ,weighted=weighted
  4 +
  5 +;RUN EXAMPLE : phangs_fits_seds_sizedist, /show_sed, /norma,/weighted
  6 +
  7 +dustem_define_la_common
  8 +
  9 +dir=!phangs_data_dir+'/ISRF/GRIDS/'
  10 +dir_to_fit = '/home/rmaris/Desktop/Soft_Libraries/dustem-wrapper_idl/LabTools/IRAP/RM/tmp/sed_size.fits'
  11 +
  12 +isrf_class = 0
  13 +isrf_class_str='_isrfclass'+strtrim(isrf_class,2)
  14 +model='DL07'
  15 +;table_name=dir+model+'_MuseISRF_JWST_G0_YPAH_YVSG_4Phangs'+isrf_class_str+'.fits'
  16 +table_name=dir+'DL07_logn_MuseISRF_JWST_adaptedG0_Ypah1added_Yvsgadded_4Phangs_noionis'+isrf_class_str+'.fits'
  17 +table_to_fit = dir_to_fit
  18 +; table_name = dir + 'DBP90_MuseISRF_JWST_G0_YPAH_YVSG_4Phangs_noionis'+isrf_class_str+'.fits'
  19 +stp=mrdfits(table_to_fit ,1,hp)
  20 +st=mrdfits(table_to_fit ,2,h)
  21 +Nst=n_elements(st)
  22 +Nseds = n_elements(st.(0))
  23 +Nparams=4
  24 +Ntrueparams=2
  25 +predicted_params=fltarr(Nseds,Nparams)
  26 +true_params=fltarr(Nseds,Ntrueparams)
  27 +rchi2_res = fltarr(Nseds)
  28 +
  29 +use_filters_names=['F300M','F335M','F360M','F0770W','F1000W','F1130W','F2100W','PACS_BLUE','PACS_GREEN','PACS_RED','PSW','PMW','PLW']
  30 +filters=dustem_filter_names2filters(use_filters_names)
  31 +Nfilters=n_elements(filters)
  32 +
  33 +array_seds = fltarr(Nseds,Nfilters)
  34 +
  35 +FOR j=0L,Nseds-1 DO BEGIN
  36 +
  37 + sed=dustem_initialize_sed(Nfilters)
  38 + ;sed=fltarr(Nfilters)
  39 + tags=tag_names(st)
  40 + FOR i=0L,Nfilters-1 DO BEGIN
  41 + ind=where(tags EQ 'I'+filters[i],count)
  42 + sed[i].filter=filters[i]
  43 + sed[i].stokesI=st[j].(ind[0])
  44 + ENDFOR
  45 + sed.sigmaii=(sed.stokesI*0.1)^2
  46 + FOR i=0L,Ntrueparams-1 DO BEGIN
  47 + true_params[j,i]=stp[j].(i)
  48 + ENDFOR
  49 + vsg_id=2L
  50 + vsg_id_str=strtrim(vsg_id,2)
  51 + fixed_param_desc= ['dustem_plugin_phangs_class_isrf_2'];,'(*!dustem_params).grains['+vsg_id_str+'].MDUST_O_MH']
  52 + fixed_param_val= [1.];,0.0016600000]
  53 +
  54 + array_seds[j,*]=sed.stokesI
  55 +
  56 + params=dustem_brute_force_fit(sed,table_name,filters,fact=fact,chi2=chi2,normalize=normalize,rchi2=rchi2,show_sed=show_sed $
  57 + ,params_hit=params_hit,params_uncertainties=params_uncertainties,params_min=params_min,params_max=params_max $
  58 + ,fixed_parameters_description=fixed_param_desc,fixed_parameters_values=fixed_param_val $
  59 + ,weighted_params=weighted_params,reset=reset,best_sed=best_sed,best_grid_sed=best_grid_sed,status=status,weightning=weightning $
  60 + ,weighted_sed=weighted_sed,weighted_chi2=weighted_chi2,rweighted_chi2=rweighted_chi2,marginalize_method=marginalize_method)
  61 +
  62 + ;print,params,fact
  63 + ;IF fact GT 2 THEN BEGIN
  64 + ; print,ii,fact,chi2
  65 + ; stop
  66 + ;ENDIF
  67 +
  68 + print,' true_params = ', true_params[j,*]
  69 + print,'fitted_params= ', params
  70 + print,'rweighted_chi2= ', rchi2
  71 + ;print,'best_fit_seds = ', best_sed
  72 + print,'============================================================================='
  73 +
  74 +
  75 +
  76 + FOR i=0L,Nparams-1 DO BEGIN
  77 + predicted_params[j,*]=weighted_params
  78 + ENDFOR
  79 +
  80 + rchi2_res[j] = rchi2
  81 +
  82 +
  83 +ENDFOR
  84 +
  85 +Ypah0 = predicted_params[*,1]
  86 +Ypah1 = predicted_params[*,2]
  87 +Ypah = Ypah0 + Ypah1
  88 +Xion =Ypah1/Ypah
  89 +amax_Ypah0 = true_params[*,0]
  90 +amax_Ypah1 = true_params[*,1]
  91 +
  92 +
  93 +
  94 +amax = amax_Ypah0 + amax_Ypah1
  95 +
  96 +ind = where(amax_Ypah0 EQ amax_Ypah1,count)
  97 +
  98 +
  99 +xsize=500
  100 +ysize=500
  101 +!P.Position = [1., 1., 1., 1.]
  102 +
  103 +window, 0 ,xsize=500,ysize=500
  104 +cgplot,amax,Xion,psym='Filled Circle',color='blue',/ysty,xtit='amax',ytit='Xion',/xlog,yrange=[la_min(Xion),la_max(Xion)]
  105 +cgplot,amax[ind],Xion[ind],psym='Filled Circle',color='red',/ysty,xtit='amax',ytit='Xion',/xlog,yrange=[la_min(Xion),la_max(Xion)],/overplot
  106 +
  107 +ind=where(amax_Ypah1 EQ amax_Ypah1[0] ,count)
  108 +window, 1,xsize=500,ysize=500
  109 +cgplot,amax_Ypah0[ind],Xion,psym='Filled Circle',color='blue',/ysty,xtit='amax_Ypah0',ytit='Xion',/xlog,yrange=[la_min(Xion),la_max(Xion)]
  110 +
  111 +ind=where(amax_Ypah0 EQ amax_Ypah0[0] ,count)
  112 +window, 2,xsize=500,ysize=500
  113 +cgplot,amax_Ypah1[ind],Xion,psym='Filled Circle',color='blue',/ysty,xtit='amax_Ypah1',ytit='Xion',/xlog,yrange=[la_min(Xion),la_max(Xion)]
  114 +
  115 +window, 3,xsize=500,ysize=500
  116 +cgplot,amax,Ypah1,psym='Filled Circle',color='red',/ysty,xtit='amax',ytit='Ypah1',/xlog ,yrange=[la_min(Ypah1),la_max(Ypah1)]
  117 +
  118 +window, 4,xsize=500,ysize=500
  119 +cgplot,amax,Ypah0,psym='Filled Circle',color='red',/ysty,xtit='amax',ytit='Ypah0',/xlog ,yrange=[la_min(Ypah0),la_max(Ypah0)]
  120 +
  121 +window, 5,xsize=500,ysize=500
  122 +cgplot,amax,rchi2_res,psym='Filled Circle',color='black',xtit='amax',ytit='rchi2',/xlog,yrange=[la_min(rchi2_res),la_max(rchi2_res)]
  123 +cgplot,amax[ind],rchi2_res[ind],psym='Filled Circle',color='red',/ysty,xtit='amax',ytit='rchi2',/xlog,yrange=[la_min(rchi2_res),la_max(rchi2_res)],/overplot
  124 +
  125 +window, 6,xsize=500,ysize=500
  126 +cgplot,amax,array_seds[*,1],psym='Filled Circle',color='black',xtit='amax',ytit='F335M',/xlog,yrange=[la_min(array_seds[*,1]),la_max(array_seds[*,1])]
  127 +
  128 +window, 7,xsize=500,ysize=500
  129 +cgplot,amax,array_seds[*,3],psym='Filled Circle',color='black',xtit='amax',ytit='F770W',/xlog,yrange=[la_min(array_seds[*,3]),la_max(array_seds[*,3])]
  130 +
  131 +
  132 +window, 8,xsize=500,ysize=500
  133 +cgplot,amax,array_seds[*,5],psym='Filled Circle',color='black',xtit='amax',ytit='F1130W',/xlog,yrange=[la_min(array_seds[*,5]),la_max(array_seds[*,5])]
  134 +
  135 +ratio = array_seds[*,5]/array_seds[*,3]
  136 +
  137 +window, 9,xsize=500,ysize=500
  138 +cgplot,amax,ratio,psym='Filled Circle',color='black',xtit='amax',ytit='ratio',/xlog,yrange=[la_min(ratio),la_max(ratio)]
  139 +
  140 +
  141 +END
  142 +
  143 +
... ...
LabTools/IRAP/RM/phangs_get_observed_seds.pro 0 → 100644
  Diff comments   View file @9a95b8f
... ... @@ -0,0 +1,1536 @@
  1 +PRO phangs_get_observed_seds,source_name=source_name $
  2 + ,model=model $
  3 + ,normalize=normalize $
  4 + ,from_restore=from_restore $
  5 + ,save=save $
  6 + ,force_mathis=force_mathis $
  7 + ,fit_G0=fit_G0 $
  8 + ,G0_option=G0_option $
  9 + ,G0_fits_map=G0_fits_map $
  10 + ,G0_fits_header=G0_fits_header $
  11 + ,set_yvsg=set_yvsg $
  12 + ,yvsg_option=yvsg_option $
  13 + ,yvsg_fits_map=yvsg_fits_map $
  14 + ,yvsg_fits_header=yvsg_fits_header $
  15 + ,help=help $
  16 + ,include_herschel=include_herschel $
  17 + ,show_seds=show_seds $
  18 + ,resolution_filter=resolution_filter $
  19 + ,resolution_value=resolution_value $
  20 + ,subtract_star_light=subtract_star_light $
  21 + ,table_name_root=table_name_root $
  22 + ,class_min_max=class_min_max $
  23 + ,nostop=nostop $
  24 + ,marginalize_method=marginalize_method $
  25 + ,do_postscripts=do_postscripts $
  26 + ,postscripts_dir=postscripts_dir $
  27 + ,silent=silent
  28 +
  29 +;+
  30 +; NAME:
  31 +; phangs_brute_force_fit_with_isrf_grid
  32 +; CALLING SEQUENCE:
  33 +; phangs_brute_force_fit_with_isrf_grid[,source_name=][,/normalize][,/save][,/help][,/from_restore][,/force_mathis]
  34 +; PURPOSE:
  35 +; Brute-Force fit SEDs of a given phangs Galaxy using ISRF grid
  36 +; INPUTS:
  37 +; None
  38 +; OPTIONAL KEYWORDS:
  39 +; source_name = source name (default='ngc0628')
  40 +;subtract_star_light = if set, subtract star light contribution before the fit
  41 +; normalize = if set, SEDs are normalized to NH=1.e20 H/cm2 (from HI+CO) and are then brute-force fitted against the table without automatic rescaling
  42 +; from_restore = redo final plots from previoulsy saved file
  43 +; save = if set, save the classes
  44 +; fit_G0 = if set, G0 is fitted. Otherwise assumed equal to predicted (or taken from a ower resolution map
  45 +; set_yvsg = if set, sets Ysvg. assumed equal to predicted (or taken from a ower resolution map
  46 +; help = if set, print this help
  47 +; OUTPUTS:
  48 +; None
  49 +; OPTIONAL INPUT:
  50 +; source_name = source name (default='ngc0628')
  51 +; resolution_filter = dustemwrap filter name at the resolution of which the fit should be done
  52 +; table_name_root = string giving the grid fits file name root (default = see code)
  53 +; OPTIONAL OUTPUT:
  54 +; None
  55 +; PROCEDURE AND SUBROUTINE USED
  56 +;
  57 +; SIDE EFFECTS:
  58 +; produces file *_DBP90_JWST_G0_YPAH_YVSG_on-voronoi_with-ISRFclasses*.sav
  59 +; EXAMPLE:
  60 +; phangs_brute_force_fit_with_isrf_grid,source_name='ngc0628',/normalize,/save
  61 +; phangs_brute_force_fit_with_isrf_grid,source_name='ngc0628',/normalize,/save,/include_herschel,/fit_G0,/show_seds
  62 +; phangs_brute_force_fit_with_isrf_grid,source_name='ngc0628',/normalize,/save,/include_herschel,/fit_G0,/force_Mathis
  63 +; MODIFICATION HISTORY:
  64 +; written by Jean-Philippe Bernard
  65 +;-
  66 +
  67 +IF keyword_set(help) THEN BEGIN
  68 + doc_library,'phangs_brute_force_fit_with_isrf_grid'
  69 + goto,the_end
  70 +ENDIF
  71 +
  72 +obp=[1.1,0.,1.15,1]
  73 +win=0L
  74 +
  75 +;weightning='1/chi2'
  76 +;weightning='likelyhood'
  77 +weightning='rellikelyhood' ;weightning strategy used for weighted parameters
  78 +
  79 +use_source_name='ngc0628'
  80 +IF keyword_set(source_name) THEN use_source_name=source_name
  81 +
  82 +norm_str=''
  83 +IF keyword_set(normalize) THEN BEGIN
  84 + norm_str='_normalized'
  85 +ENDIF
  86 +mathis_str=''
  87 +IF keyword_set(force_mathis) THEN BEGIN
  88 + mathis_str='_mathis'
  89 +ENDIF
  90 +
  91 +IF keyword_set(set_yvsg) and not keyword_set(fit_G0) THEN BEGIN
  92 + message,'/fit_G0 should be set with /set_yvsg.',/continue
  93 + stop
  94 +ENDIF
  95 +
  96 +dustem_define_la_common
  97 +
  98 +reso_str=''
  99 +IF keyword_set(resolution_filter) OR keyword_set(resolution_value) THEN BEGIN
  100 + IF keyword_set(resolution_filter) THEN BEGIN
  101 + reso_str='_'+resolution_filter
  102 + reso=dustem_filter2reso(resolution_filter)
  103 + ENDIF ELSE BEGIN
  104 + reso_str='_'+strtrim(resolution_value,2)
  105 + reso=resolution_value
  106 + ENDELSE
  107 +ENDIF
  108 +
  109 +use_model='DBP90'
  110 +IF keyword_set(model) THEN BEGIN
  111 + use_model=model
  112 +ENDIF
  113 +;stop
  114 +;=== define what we call vsg, for /set_yvsg option
  115 +CASE use_model OF
  116 + 'DL07': vsg_id=2L
  117 + 'DBP90': vsg_id=1L
  118 +ENDCASE
  119 +vsg_id_str=strtrim(vsg_id,2)
  120 +;stop
  121 +
  122 +data_dir=!phangs_data_dir+'/ISRF/WORKJP/'
  123 +grids_data_dir=!phangs_data_dir+'/ISRF/GRIDS/'
  124 +
  125 +IF keyword_set(do_postscripts) THEN BEGIN
  126 + ps_dir=!phangs_data_dir+'/ISRF/WORKJP/PS/'
  127 + IF keyword_set(postscripts_dir) THEN BEGIN
  128 + ps_dir=postscripts_dir+'/'
  129 + ENDIF
  130 + force_mkdir,ps_dir
  131 +ENDIF
  132 +
  133 +IF keyword_set(from_restore) THEN GOTO,from_restore
  134 +
  135 +dustem_init,show_plots=show_plots
  136 +
  137 +;==== needed for NHCO,NHHI
  138 +file=data_dir+use_source_name+'_CO_images'+reso_str+'.sav'
  139 +st_info=file_info(file)
  140 +IF st_info.exists NE 1 THEN BEGIN
  141 + message,'Could not find '+file,/continue
  142 + stop
  143 +ENDIF
  144 +restore,file,/verb
  145 +;% RESTORE: Restored variable: NHCO.
  146 +;% RESTORE: Restored variable: HREF.
  147 +file=data_dir+use_source_name+'_HI_images'+reso_str+'.sav'
  148 +st_info=file_info(file)
  149 +IF st_info.exists NE 1 THEN BEGIN
  150 + message,'Could not find '+file,/continue
  151 + stop
  152 +ENDIF
  153 +restore,file,/verb
  154 +;% RESTORE: Restored variable: NHI.
  155 +;% RESTORE: Restored variable: HREF.
  156 +
  157 +;only using voronoi_id and only for its maximum value
  158 +file=data_dir+use_source_name+'_muse_data'+reso_str+'.sav'
  159 +st_info=file_info(file)
  160 +IF st_info.exists NE 1 THEN BEGIN
  161 + message,'Could not find '+file,/continue
  162 + stop
  163 +ENDIF
  164 +restore,file,/verb
  165 +;% RESTORE: Restored variable: ST_TEMPLATES.
  166 +;% RESTORE: Restored variable: ST_MUSE_WEIGHTS.
  167 +;% RESTORE: Restored variable: VORONOI_ID.
  168 +;% RESTORE: Restored variable: AGE_VALUES.
  169 +;% RESTORE: Restored variable: METALICITY_VALUES.
  170 +;% RESTORE: Restored variable: BINS.
  171 +;% RESTORE: Restored variable: HREF.
  172 +
  173 +;stop
  174 +
  175 +;==== apparently, classes is not used by this code
  176 +;file=data_dir+'isrf_classes_one_ratio_on_'+use_source_name+reso_str+'.fits'
  177 +;st_info=file_info(file)
  178 +;IF st_info.exists NE 1 THEN BEGIN
  179 +; message,'Could not find '+file,/continue
  180 +; stop
  181 +;ENDIF
  182 +;classes=mrdfits(file,2,h)
  183 +
  184 +;==== apparently, map_classes is not used by this code
  185 +;file=data_dir+use_source_name+'_isrf_classes_map'+reso_str+'.fits'
  186 +;st_info=file_info(file)
  187 +;IF st_info.exists NE 1 THEN BEGIN
  188 +; message,'Could not find '+file,/continue
  189 +; stop
  190 +;ENDIF
  191 +;map_classes=readfits(file,hc)
  192 +
  193 +file=data_dir+use_source_name+'_isrf_classes_voronoi'+reso_str+'.sav'
  194 +st_info=file_info(file)
  195 +IF st_info.exists NE 1 THEN BEGIN
  196 + message,'Could not find '+file,/continue
  197 + stop
  198 +ENDIF
  199 +restore,file,/verb
  200 +;% RESTORE: Restored variable: VOR_CLASS.
  201 +;==== Restore the ISRF predictions
  202 +file=data_dir+use_source_name+'_isrf_min_prediction'+reso_str+'.sav'
  203 +st_info=file_info(file)
  204 +IF st_info.exists NE 1 THEN BEGIN
  205 + message,'Could not find '+file,/continue
  206 + stop
  207 +ENDIF
  208 +restore,file,/verb
  209 +;% RESTORE: Restored variable: ISRFS.
  210 +;% RESTORE: Restored variable: F0S.
  211 +;% RESTORE: Restored variable: HREF.
  212 +;% RESTORE: Restored variable: USE_RESO_FILTER.
  213 +;stop
  214 +;==== somehow, the F0s in the save set have only one value. Not sure why ... So, gets recomputed below.
  215 +;F0s_predicted=F0s
  216 +;In case F0s are not in save set file:
  217 +NNvor=(size(ISRFs))[2]
  218 +F0s_predicted=fltarr(NNvor)
  219 +G0s_predicted=fltarr(NNvor)
  220 +;==== get Mathis field for F0 calculations
  221 +bidon=dustem_get_wavelengths(isrf_wavelengths=isrf_wavelengths)
  222 +file=!dustem_soft_dir+'/data/ISRF_MATHIS.DAT'
  223 +st_info=file_info(file)
  224 +IF st_info.exists NE 1 THEN BEGIN
  225 + message,'Could not find '+file,/continue
  226 + stop
  227 +ENDIF
  228 +readcol,file,Mathis_wavs,Mathis_ISRF
  229 +Mathis_ISRF=interpol(Mathis_ISRF,Mathis_wavs,isrf_wavelengths)
  230 +Mathis_1mic=interpol(Mathis_ISRF,isrf_wavelengths,1.0)
  231 +;==== compute F0s predicted by ISRF prediction
  232 +FOR i=0L,NNvor-1 DO BEGIN
  233 + F0s_predicted[i]=interpol(ISRFs[*,i],isrf_wavelengths,1.0)/Mathis_1mic
  234 +ENDFOR
  235 +;stop
  236 +;=== If using the Mathis field, F0s_predicted should be rescaled by G0prime of the class
  237 +IF keyword_set(force_mathis) THEN BEGIN
  238 + F0s_predicted=la_mul(F0s_predicted,G0prime[vor_class])
  239 +ENDIF
  240 +
  241 +ind=where(vor_class NE 0)
  242 +class_min=la_min(vor_class[ind])
  243 +class_max=la_max(vor_class[ind])
  244 +IF keyword_set(class_min_max) THEN BEGIN
  245 + class_min=class_min_max[0]
  246 + class_max=class_min_max[1]
  247 +ENDIF
  248 +;=== For test (vid=0 is in class 23 for NGC0628)
  249 +;class_min=23
  250 +;class_max=23
  251 +
  252 +;stop
  253 +help,nhco,nhi,voronoi_id,classes,map_classes,vor_class
  254 +;NHCO FLOAT = Array[3515, 4576]
  255 +;NHI FLOAT = Array[3515, 4576]
  256 +;VORONOI_ID LONG = Array[3515, 4576]
  257 +;CLASSES STRUCT = -> <Anonymous> Array[30]
  258 +;MAP_CLASSES LONG = Array[3515, 4576]
  259 +;VOR_CLASS LONG = Array[47069]
  260 +
  261 +help,classes,/str
  262 +;** Structure <50e057d8>, 9 tags, length=1648, data length=1640, refs=1:
  263 +; NUMBER LONG 0
  264 +; NAME STRING ' '
  265 +; ISRF_WAVELENGTHS
  266 +; FLOAT Array[200]
  267 +; ISRF_TEMPLATE FLOAT Array[200]
  268 +; FLUX_RATIO FLOAT 1.00000e-05
  269 +; FLUX_RATIO_WAV FLOAT 0.150000
  270 +; REFERENCE_WAVELENGTH
  271 +; FLOAT 1.00000
  272 +; FLUX_RATIO_MIN FLOAT 7.88046e-06
  273 +; FLUX_RATIO_MAX FLOAT 1.26896e-05
  274 +
  275 +;For test
  276 +;class_min=11 & class_max=26
  277 +;class_min=1 & class_max=30
  278 +;class_min=25 & class_max=25
  279 +
  280 +;==== If force_mathis is set, use only class 0 for the dust fit
  281 +IF keyword_set(force_mathis) THEN BEGIN
  282 + class_min=0L
  283 + class_max=0L
  284 +ENDIF
  285 +
  286 +;=== restore pre-computed SEDs
  287 +file=data_dir+use_source_name+'_jwst_seds_muse_pixels'+reso_str+'.sav'
  288 +st_info=file_info(file)
  289 +IF st_info.exists NE 1 THEN BEGIN
  290 + message,'Could not find '+file,/continue
  291 + stop
  292 +ENDIF
  293 +restore,file,/verb
  294 +;% RESTORE: Restored variable: ALL_SEDS.
  295 +jwst_seds=all_seds
  296 +jwst_filters=(*jwst_seds[0]).filter
  297 +file=data_dir+use_source_name+'_muse_seds_muse_pixels'+reso_str+'.sav'
  298 +st_info=file_info(file)
  299 +IF st_info.exists NE 1 THEN BEGIN
  300 + message,'Could not find '+file,/continue
  301 + stop
  302 +ENDIF
  303 +restore,file,/verb ;contains the Muse SEDs in voronoi bins
  304 +;% RESTORE: Restored variable: ALL_SEDS.
  305 +muse_seds=all_seds
  306 +muse_filters=(*muse_seds[0]).filter
  307 +file=data_dir+use_source_name+'_astrosat_seds_muse_pixels'+reso_str+'.sav'
  308 +st_info=file_info(file)
  309 +IF st_info.exists NE 1 THEN BEGIN
  310 + message,'Could not find '+file,/continue
  311 + stop
  312 +ENDIF
  313 +restore,file,/verb
  314 +;% RESTORE: Restored variable: ALL_SEDS.
  315 +astrosat_seds=all_seds
  316 +astrosat_filters=(*astrosat_seds[0]).filter
  317 +file=data_dir+use_source_name+'_herschel_seds_muse_pixels'+reso_str+'.sav'
  318 +st_info=file_info(file)
  319 +IF st_info.exists NE 1 THEN BEGIN
  320 + message,'Could not find '+file,/continue
  321 + stop
  322 +ENDIF
  323 +restore,file,/verb
  324 +;% RESTORE: Restored variable: ALL_SEDS.
  325 +herschel_seds=all_seds
  326 +herschel_filters=(*herschel_seds[0]).filter
  327 +file=data_dir+use_source_name+'_seds_indices'+reso_str+'.sav'
  328 +st_info=file_info(file)
  329 +IF st_info.exists NE 1 THEN BEGIN
  330 + message,'Could not find '+file,/continue
  331 + stop
  332 +ENDIF
  333 +restore,file,/verb
  334 +;% RESTORE: Restored variable: ALL_SEDS_INDICES.
  335 +
  336 +help,astrosat_seds,muse_seds,jwst_seds,herschel_seds
  337 +;stop
  338 +
  339 +Nvor=long(max(voronoi_id))+1 ;because Voronoi bins starts from 0
  340 +
  341 +Muse_factors_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  342 +Voronoi_surfaces_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  343 +
  344 +IF not keyword_set(resolution_filter) AND not keyword_set(resolution_value) THEN BEGIN ;because Muse factors are not saved at smoothed resolution ....
  345 + FOR vid=0LL,Nvor-1 DO BEGIN
  346 + IF ptr_valid(all_seds_indices[vid]) THEN BEGIN
  347 + Muse_factors_map[*all_seds_indices[vid]]=Muse_factors[vid]
  348 + Voronoi_surfaces_map[*all_seds_indices[vid]]=voronoi_surfaces[vid]
  349 + ENDIF
  350 + ENDFOR
  351 + window,win & win=win+1
  352 + tit='Muse Scale factors'
  353 + image_cont20,Muse_factors_map,href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,imrange=[0.5,2],tit=tit,off_bar_pos=obp
  354 + ;image_cont20,Muse_factors_map,href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1
  355 + ;stop
  356 + message,'Average Muse Scale factor ='+strtrim(la_mean(Muse_factors_map),2),/continue
  357 + window,win & win=win+1
  358 + tit='1/.Voronoi surface'
  359 + image_cont20,1./Voronoi_surfaces_map,href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,imrange=[-0.1,1]*1.e8,tit=tit,off_bar_pos=obp
  360 + ;stop
  361 +ENDIF
  362 +
  363 +;===== Compute seds pointer
  364 +one_sed=(*jwst_seds[0])
  365 +all_filters=[jwst_filters,astrosat_filters,muse_filters]
  366 +IF keyword_set(include_herschel) THEN BEGIN
  367 + all_filters=[all_filters,herschel_filters]
  368 +ENDIF
  369 +Nfilters=n_elements(all_filters)
  370 +seds_ptr=ptrarr(Nvor)
  371 +nperc=5./100.
  372 +
  373 +empty_sed=dustem_initialize_sed(Nfilters)
  374 +IF keyword_set(include_herschel) THEN BEGIN
  375 + FOR ii=0L,Nvor-1 DO BEGIN
  376 + IF ptr_valid(astrosat_seds[ii]) AND ptr_valid(muse_seds[ii]) AND ptr_valid(jwst_seds[ii]) AND ptr_valid(herschel_seds[ii]) THEN BEGIN
  377 + seds_ptr[ii]=ptr_new([*astrosat_seds[ii],*muse_seds[ii],*jwst_seds[ii],*herschel_seds[ii]])
  378 + ENDIF ELSE BEGIN
  379 + seds_ptr[ii]=ptr_new(empty_sed)
  380 + ENDELSE
  381 + ENDFOR
  382 +ENDIF ELSE BEGIN
  383 + FOR ii=0L,Nvor-1 DO BEGIN
  384 + IF ptr_valid(astrosat_seds[ii]) AND ptr_valid(muse_seds[ii]) AND ptr_valid(jwst_seds[ii]) THEN BEGIN
  385 + seds_ptr[ii]=ptr_new([*astrosat_seds[ii],*muse_seds[ii],*jwst_seds[ii]])
  386 + ENDIF ELSE BEGIN
  387 + seds_ptr[ii]=ptr_new(empty_sed)
  388 + ENDELSE
  389 + ENDFOR
  390 +ENDELSE
  391 +
  392 +;==== remove zero variances
  393 +FOR ii=0L,Nvor-1 DO BEGIN
  394 + ;=== check for 0 variances
  395 + ind=where((*seds_ptr[ii]).sigmaII EQ 0,count)
  396 + IF count NE 0 THEN BEGIN
  397 + this_sed=*seds_ptr[ii]
  398 + Ivalues=((*seds_ptr[ii]).stokesI)[ind]
  399 + sigmaii_values=(Ivalues*nperc)^2
  400 + this_sed[ind].sigmaii=sigmaii_values
  401 + seds_ptr[ii]=ptr_new(this_sed)
  402 + ENDIF
  403 +ENDFOR
  404 +
  405 +;==== Brute force fit using ISRF classes
  406 +
  407 +;Nvor=max(voronoi_id)
  408 +F0s=fltarr(Nvor)+la_undef()
  409 +G0s=fltarr(Nvor)+la_undef()
  410 +EBVs=fltarr(Nvor)+la_undef()
  411 +observed_seds=ptrarr(Nvor) ;All observed SEDs
  412 +best_fit_seds=ptrarr(Nvor) ;All best fit SEDs
  413 +Ygrain1=fltarr(Nvor)+la_undef()
  414 +Ygrain2=fltarr(Nvor)+la_undef()
  415 +Ygrain3=fltarr(Nvor)+la_undef()
  416 +facts=fltarr(Nvor)+la_undef()
  417 +chi2s=fltarr(Nvor)+la_undef()
  418 +rchi2s=fltarr(Nvor)+la_undef()
  419 +dF0s=fltarr(Nvor)+la_undef()
  420 +dYgrain1=fltarr(Nvor)+la_undef()
  421 +dYgrain2=fltarr(Nvor)+la_undef()
  422 +dYgrain3=fltarr(Nvor)+la_undef()
  423 +F0_hit=fltarr(Nvor)+la_undef()
  424 +Ygrain1_hit=fltarr(Nvor)+la_undef()
  425 +Ygrain2_hit=fltarr(Nvor)+la_undef()
  426 +Ygrain3_hit=fltarr(Nvor)+la_undef()
  427 +NH_values=fltarr(Nvor)+la_undef()
  428 +stellar_normalizations=fltarr(Nvor)+la_undef()
  429 +
  430 +NHmap=la_add(NHCO,NHI) ;used NH map in 1.e21 (from CO+HI)
  431 +
  432 +;==== This is all the filters in a given SED
  433 +all_filters=(*seds_ptr[0]).filter
  434 +all_instru=dustem_filter2instru(all_filters)
  435 +all_reso=dustem_filter2reso(all_filters) ;caution: resolution not good for astrosat and SDSS. Done by hand below
  436 +all_wavs=dustem_filter2wav(all_filters)
  437 +
  438 +;==== filters used for brute force fit
  439 +try=get_phangs_filter_names()
  440 +IF keyword_set(include_herschel) THEN BEGIN
  441 + use_fit_filters_names=['F300M','F335M','F360M','F0770W','F1000W','F1130W','F2100W','PACS_BLUE','PACS_GREEN','PACS_RED','PSW','PMW','PLW'] ;JWST+Herschel
  442 + ;use_fit_filters_names=['F200W','F300M','F335M','F360M','F0770W','F1000W','F1130W','F2100W','PACS_BLUE','PACS_GREEN','PACS_RED','PSW','PMW','PLW'] ;JWST+Herschel
  443 + ;use_fit_filters_names=['F360M','F0770W','F1000W','F1130W','F2100W','PACS_BLUE','PACS_GREEN','PACS_RED','PSW','PMW','PLW'] ;JWST+Herschel
  444 + ;use_fit_filters_names=['PACS_BLUE','PACS_GREEN','PACS_RED','PSW','PMW','PLW'] ;JWST+Herschel
  445 +ENDIF ELSE BEGIN
  446 + use_fit_filters_names=['F300M','F335M','F360M','F0770W','F1000W','F1130W','F2100W'] ;This is with JWST only
  447 + ;use_fit_filters_names=['F200W','F300M','F335M','F360M','F0770W','F1000W','F1130W','F2100W'] ;This is with JWST only
  448 + ;use_fit_filters_names=['F360M','F0770W','F1000W','F1130W','F2100W'] ;This is with JWST only
  449 +ENDELSE
  450 +use_fit_filters=dustem_filter_names2filters(use_fit_filters_names)
  451 +
  452 +;==== remove filters with improper resolution (wider than resolution_filter) from the SED to be fitted
  453 +IF keyword_set(resolution_filter) OR keyword_set(resolution_value) THEN BEGIN
  454 + use_fit_instru=dustem_filter2instru(use_fit_filters)
  455 + use_fit_reso=dustem_filter2reso(use_fit_filters) ;caution: resolution not good for astrosat and SDSS. Done by hand below
  456 + ind=where(use_fit_instru EQ 'ASTROSAT',count)
  457 + ;IF count NE 0 THEN use_fit_reso[ind]=4./60.^2 ;CAUTION: This is made up
  458 + IF count NE 0 THEN use_fit_reso[ind]=1./60.^2 ;CAUTION: This is made up to keep astrosat all the way down to 1 arcsec
  459 + ind=where(use_fit_instru EQ 'SDSS',count)
  460 + IF count NE 0 THEN use_fit_reso[ind]=1./60.^2
  461 + ind=where(use_fit_reso LE reso,count)
  462 + use_fit_filters=use_fit_filters[ind]
  463 + ;stop
  464 +ENDIF
  465 +
  466 +fit_wavs=dustem_filter2wav(use_fit_filters)
  467 +N_fit_filters=n_elements(use_fit_filters)
  468 +;=== compute filter index used to restrict SEDs to the filters above.
  469 +ind_fit_filters=[-1]
  470 +FOR i=0L,N_fit_filters-1 DO BEGIN
  471 + indd=where(all_filters EQ use_fit_filters[i],count)
  472 + IF count NE 0 THEN ind_fit_filters=[ind_fit_filters,indd]
  473 +ENDFOR
  474 +ind_fit_filters=ind_fit_filters[1:*]
  475 +
  476 +show_sed=0
  477 +
  478 +use_f0_option='DEFAULT'
  479 +IF keyword_set(G0_option) THEN use_f0_option=G0_option
  480 +use_yvsg_option='DEFAULT'
  481 +IF keyword_set(yvsg_option) THEN use_yvsg_option=yvsg_option
  482 +
  483 +IF not keyword_set(fit_G0) THEN BEGIN
  484 + fit_G0_str='_GOfixed'
  485 +ENDIF ELSE BEGIN
  486 + fit_G0_str=''
  487 +ENDELSE
  488 +
  489 +IF keyword_set(set_yvsg) THEN BEGIN
  490 + set_yvsg_str='_YVSGfixed'
  491 +ENDIF ELSE BEGIN
  492 + set_yvsg_str=''
  493 +ENDELSE
  494 +
  495 +;IF not keyword_set(force_mathis) AND not keyword_set(fit_G0) THEN BEGIN
  496 +IF not keyword_set(fit_G0) THEN BEGIN
  497 + fixed_parameters_description=['dustem_plugin_phangs_class_isrf_2'] ;This is the ISRF F0 factor, fixed for grid
  498 +ENDIF ELSE BEGIN
  499 + IF keyword_set(set_yvsg) THEN BEGIN
  500 + fixed_parameters_description=['(*!dustem_params).grains['+vsg_id_str+'].MDUST_O_MH']
  501 + ENDIF ELSE BEGIN
  502 + ;left blanck on purpose. GO free if fit_G0=1
  503 + ENDELSE
  504 +ENDELSE
  505 +
  506 +;file=!dustem_soft_dir+'/data/ISRF_MATHIS.DAT'
  507 +;readcol,file,Mathis_wavs,Mathis_ISRF
  508 +;Mathis_ISRF=interpol(Mathis_ISRF,Mathis_wavs,isrf_wavelengths)
  509 +;Mathis_1mic=interpol(Mathis_ISRF,isrf_wavelengths,1.0)
  510 +
  511 +ind=where(NHmap NE la_undef(),count)
  512 +NH_min=la_abs(la_min(NHmap[ind]))
  513 +
  514 +;stop
  515 +
  516 +;=== stuff for stellar subtraction
  517 +;starlight_filter_names=['F200W'] ;These are the filters that will be used to normalize stellar light (ie assumed to be dominated by stellar emission)
  518 +starlight_filter_names=['F300M']
  519 +N_starlight_filters=n_elements(starlight_filter_names)
  520 +use_starlight_filters=dustem_filter_names2filters(starlight_filter_names)
  521 +wav_star=dustem_filter2wav(use_starlight_filters)
  522 +
  523 +;starlight_ebv_names=['F148W'] ;Astrosat1, used to compute extinction
  524 +starlight_ebv_names=['sdss_g'] ;MUSE filter shortest wavelength, used to compute extinction
  525 +N_ebv_filters=n_elements(starlight_ebv_names)
  526 +use_ebv_filters=dustem_filter_names2filters(starlight_ebv_names)
  527 +wav_ebv_star=dustem_filter2wav(use_ebv_filters)
  528 +
  529 +indstar=[-1L]
  530 +indebv=[-1L]
  531 +FOR i=0L,N_starlight_filters-1 DO BEGIN
  532 + ind=where(all_filters EQ use_starlight_filters[i],count)
  533 + IF count NE 0 THEN BEGIN
  534 + indstar=[indstar,ind[0]]
  535 + ENDIF
  536 + ind=where(all_filters EQ use_ebv_filters[i],count)
  537 + IF count NE 0 THEN BEGIN
  538 + indebv=[indebv,ind[0]]
  539 + ENDIF
  540 +ENDFOR
  541 +indstar=indstar[1:*]
  542 +indebv=indebv[1:*]
  543 +
  544 +one_st={isrf_class:0L, $
  545 + Nvor:0L, $
  546 + min_G0:0.,med_G0:0.,max_G0:0., $
  547 + min_Yg1:0.,med_Yg1:0.,max_Yg1:0., $
  548 + min_Yg2:0., med_Yg2:0., max_Yg2:0., $
  549 + min_Yg3:0.,med_Yg3:0.,max_Yg3:0., $
  550 + min_rchi2:0.,med_rchi2:0.,max_rchi2:0.}
  551 +st=replicate(one_st,30)
  552 +
  553 +;indcalz=where(all_wavs GE 0.912 AND all_wavs LT 2.2,count_calz) ;Where Calzeti's correction is valid
  554 +indcalz=where(all_wavs LT 2.2,count_calz) ;Where Calzeti's correction is valid
  555 +
  556 +;This is for the DBP90 tables
  557 +;table_name_root='_MuseISRF_JWST_G0_YPAH_YVSG_4Phangs_noionis'
  558 +;This is for the new DL07 tables
  559 +use_table_name_root='_MuseISRF_JWST_G0_Ypah1added_Ypah2added_Yvsgadded_4Phangs_noionis'
  560 +IF keyword_set(table_name_root) THEN BEGIN
  561 + use_table_name_root=table_name_root
  562 +ENDIF
  563 +
  564 +;=== get fixed parameter values from map if needed
  565 +IF use_f0_option EQ 'MAP' THEN BEGIN
  566 + message,'Making map of F0 fixed values',/continue
  567 + new_map=project2(G0_fits_header,G0_fits_map,href)
  568 + F0s_from_map=fltarr(Nvor)
  569 + FOR vid=0L,Nvor -1 DO BEGIN
  570 + F0s_from_map[vid]=la_mean(new_map[*all_seds_indices[vid]])
  571 + ENDFOR
  572 +ENDIF
  573 +IF use_yvsg_option EQ 'MAP' THEN BEGIN
  574 + message,'Making map of Yvsg fixed values',/continue
  575 + new_map=project2(YVSG_fits_header,YVSG_fits_map,href)
  576 + yvsgs_from_map=fltarr(Nvor)
  577 + FOR vid=0L,Nvor -1 DO BEGIN
  578 + yvsgs_from_map[vid]=la_mean(new_map[*all_seds_indices[vid]])
  579 + ENDFOR
  580 +ENDIF
  581 +
  582 +FOR isrf_class=class_min,class_max DO BEGIN ;isrf_class runs from 1
  583 + ;print,isrf_class
  584 + ;stop
  585 + ;==== get G0prime of the template
  586 + IF isrf_class NE 0 THEN BEGIN
  587 + this_G0prime=G0prime[isrf_class]
  588 + ENDIF ELSE BEGIN
  589 + this_G0prime=1.
  590 + ENDELSE
  591 + ;stop
  592 + ;==== This is to use the grid for this ISRF class
  593 + isrf_class_str='_isrfclass'+strtrim(isrf_class,2)
  594 + table_name=grids_data_dir+use_model+use_table_name_root+isrf_class_str+'.fits'
  595 + ;select Voronoi bins with that ISRF class
  596 + IF not keyword_set(force_mathis) THEN BEGIN
  597 + ind_class=where(vor_class EQ isrf_class,Nvor_class)
  598 + ENDIF ELSE BEGIN
  599 + ;stop
  600 + ind_class=lindgen(Nvor)
  601 + Nvor_class=Nvor
  602 + ENDELSE
  603 + reset=1
  604 + IF Nvor_class NE 0 THEN BEGIN
  605 + FOR vvid=0LL,Nvor_class-1 DO BEGIN
  606 + IF vvid mod 100 EQ 0 THEN BEGIN
  607 + message,'isrf#'+strtrim(isrf_class,2)+' Doing sed '+strtrim(vvid,2)+' '+strtrim(1.*vvid/Nvor_class*100.,2)+' %',/continue
  608 + ENDIF
  609 + vid=ind_class[vvid]
  610 + ;IF vid EQ 37255 THEN stop ;caution, this is a test only
  611 + ;=== Get the predicted F0 value for that vid
  612 + F0=F0s_predicted[vid]
  613 + ;print,F0s_predicted[vid],predicted_F0s[vid]
  614 + ;stop
  615 + IF not keyword_set(fit_G0) THEN BEGIN
  616 + ;fixed_parameters_values=[F0]
  617 + ;REM: This value could ultimately be taken from a map at low resolution
  618 + IF use_f0_option EQ 'DEFAULT' THEN BEGIN
  619 + use_F0=F0*10. ;plots for NG0628 show that F0 is about 10*F0_predicted
  620 + ;use_F0=F0 ;plots for NG0628 show that F0 is about 10*F0_predicted
  621 + ENDIF ELSE BEGIN ;get value from lower resolution map
  622 + use_F0=F0s_from_map[vid]
  623 + ENDELSE
  624 + fixed_parameters_values=[use_F0]
  625 + ENDIF
  626 + IF keyword_set(set_yvsg) THEN BEGIN
  627 + ;REM: This value could ultimately be taken from a map at low resolution
  628 + IF use_yvsg_option EQ 'DEFAULT' THEN BEGIN
  629 + use_yvsg=dustem_get_param_values_default('(*!dustem_params).grains['+vsg_id_str+'].MDUST_O_MH')
  630 + ENDIF ELSE BEGIN ;get value from lower resolution map
  631 + use_yvsg=yvsgs_from_map[vid]
  632 + ENDELSE
  633 + fixed_parameters_values=[use_yvsg] ;plots for NG0628 show that F0 is about 10*F0_predicted
  634 + ENDIF
  635 +;stop
  636 + sed=*seds_ptr[vid] ;full sed
  637 + sed_norm=*seds_ptr[vid] ;full sed normalized by NH (below)
  638 + total_predicted_sed=*seds_ptr[vid]
  639 + ;=== restrict sed to requested filters
  640 + sed_fit=sed[ind_fit_filters]
  641 + ;=== normalize SED to NH=1.e20 H/cm2, if requested
  642 + IF keyword_set(normalize) THEN BEGIN
  643 + IF ptr_valid(all_seds_indices[vid]) THEN BEGIN
  644 + NH_value=la_mul(la_mean(NHmap[*all_seds_indices[vid]]),10.) ;NH value in 1.e20 H/cm2
  645 + IF NH_value LE 0 THEN NH_value=NH_value+NH_min
  646 + ENDIF ELSE BEGIN
  647 + NH_value=la_undef()
  648 + ENDELSE
  649 + sed_fit.stokesI=la_div(sed_fit.stokesI,NH_value)
  650 + ;sed_fit.sigmaII=la_div(sed_fit.sigmaII,NH_value) ;This would keep the original variances
  651 + sed_fit.sigmaII=(sed_fit.stokesI*0.1)^2 ;This instead assumes assumes 10% uncertainties
  652 + ;This is to force smaller uncertainties where pic emission of large grains is
  653 + ;stop
  654 + indwav=where(sed_fit.wave GT 70.,countwav)
  655 + IF countwav NE 0 THEN BEGIN
  656 + sed_fit[indwav].sigmaII=(sed_fit[indwav].stokesI*0.05)^2
  657 + ENDIF
  658 + sed_norm.stokesI=la_div(sed_norm.stokesI,NH_value)
  659 + sed_fit_normalize=0
  660 + NH_values[vid]=NH_value
  661 + ENDIF ELSE BEGIN
  662 + sed_fit_normalize=1
  663 + ENDELSE
  664 + ;=== Subtract star light from the SED, if needed
  665 + ;=== caution this is done by simple interpolation normalized on one filter value.
  666 + total_sed_fit=sed_fit
  667 + IF keyword_set(subtract_star_light) THEN BEGIN
  668 + toto=interpol(ISRFs[*,vid],isrf_wavelengths,wav_star) ;ISRF interpolated at stellar wavelengths in [ergs/s/cm2/Hz] scaled so that Inustar matches the muse data
  669 + fact_mjy2ergs=1.e-20*1.e7*1.e-4
  670 + toto=toto/fact_mjy2ergs ;Now in MJy/sr
  671 + toto=toto/NH_value
  672 + norm_stars=la_mean((sed_norm.stokesI)[indstar]/toto[indstar]) ;un-redenned stellar normalisation factor
  673 + ;print,norm_stars,NH_value,norm_stars*NH_value
  674 + ;message,'Stellar normalization factor='+strtrim(norm_stars,2),/continue
  675 + ;stop
  676 + ;sed_stars=interpol(ISRFs[*,vid],isrf_wavelengths,all_wavs)/fact_mjy2ergs
  677 + sed_stars=interpol(ISRFs[*,vid],isrf_wavelengths,all_wavs)/fact_mjy2ergs/NH_value
  678 + sed_stars=la_mul(sed_stars,norm_stars) ;This is the stellar Full SED
  679 +
  680 + wwav=all_wavs[indebv[0]]*1.e4 ;AA
  681 + ;flux_ratio=(sed_norm.stokesI/sed_stars)[indebv[0]] ;flux ratio between unredened and redenned (observed) flux
  682 + ;flux_ratio=sed_norm[indebv[0]].stokesI/sed_stars[indebv[0]] ;flux ratio between unredened and redenned (observed) flux
  683 + flux_ratio=sed_stars[indebv[0]]/sed_norm[indebv[0]].stokesI ;flux ratio between unredened and redenned (observed) flux
  684 + EBV=calz_guess_ebv(wwav,flux_ratio,R_v=R_v)
  685 + ;stop
  686 + sed_stars_reddened=sed_stars
  687 + ;We want to reden, so, use -1.*EBV
  688 + calz_unred, all_wavs[indcalz]*1.e4, sed_stars[indcalz], -1.*EBV, result, R_V = R_V
  689 + sed_stars_reddened[indcalz]=result
  690 + ;calz_unred, all_wavs*1.e4, sed_stars, -1.*EBV, result, R_V = R_V
  691 + norm_stars_redenned=la_mean((sed_norm.stokesI)[indstar]/sed_stars_reddened[indstar]) ;redenned stellar normalisation factor
  692 +
  693 + ;stop
  694 + sed_fit_stars=interpol(ISRFs[*,vid],isrf_wavelengths,fit_wavs)/fact_mjy2ergs/NH_value
  695 + sed_fit_stars=la_mul(sed_fit_stars,norm_stars) ;stellar SED at fit wavelengths, normalized to match
  696 + sed_fit.stokesI=la_sub(sed_fit.stokesI,sed_fit_stars) ;This is the dust-only SED to be fit
  697 + ;stop
  698 + stellar_normalizations[vid]=norm_stars
  699 + ENDIF
  700 + ind=where(sed_fit.stokesI EQ la_undef(),count)
  701 + ;=== brute force fit the SED
  702 + IF count EQ 0 THEN BEGIN
  703 + ;stop
  704 +
  705 + seds_stars[vid]=ptr_new(sed_norm.stokesI)
  706 + best_fit_seds[vid]=ptr_new(total_predicted_sed.stokesI)
  707 +
  708 + ;print,params[0],params[1],params[2],fact,chi2
  709 +
  710 + ENDIF
  711 + ENDFOR
  712 + ;ind=where(vor_class EQ isrf_class,Nvor_class)
  713 + st[isrf_class].Nvor=Nvor_class
  714 + st[isrf_class].isrf_class=isrf_class
  715 + st[isrf_class].min_G0=la_min(G0s[ind_class])
  716 + st[isrf_class].med_G0=la_median(G0s[ind_class])
  717 + st[isrf_class].max_G0=la_max(G0s[ind_class])
  718 + st[isrf_class].min_Yg1=la_min(Ygrain1[ind_class])
  719 + st[isrf_class].med_Yg1=la_median(Ygrain1[ind_class])
  720 + st[isrf_class].max_Yg1=la_max(Ygrain1[ind_class])
  721 + st[isrf_class].min_Yg2=la_min(Ygrain2[ind_class])
  722 + st[isrf_class].med_Yg2=la_median(Ygrain2[ind_class])
  723 + st[isrf_class].max_Yg2=la_max(Ygrain2[ind_class])
  724 + st[isrf_class].min_Yg3=la_min(Ygrain3[ind_class])
  725 + st[isrf_class].med_Yg3=la_median(Ygrain3[ind_class])
  726 + st[isrf_class].max_Yg3=la_max(Ygrain3[ind_class])
  727 + st[isrf_class].min_rchi2=la_min(rchi2s[ind_class])
  728 + st[isrf_class].med_rchi2=la_median(rchi2s[ind_class])
  729 + st[isrf_class].max_rchi2=la_max(rchi2s[ind_class])
  730 +; write_xcat,st,'/tmp/'+getenv('LOGNAME')+'/essai.xcat'
  731 + write_xcat,st,data_dir+use_source_name+use_model+'_'+norm_str+mathis_str+reso_str+fit_G0_str+set_yvsg_str+'.xcat'
  732 + ENDIF ELSE BEGIN
  733 + message,'No Voronoi bins found with ISRF class '+strtrim(isrf_class,2),/continue
  734 + ENDELSE
  735 +ENDFOR
  736 +
  737 +stats_st=st
  738 +IF keyword_set(save) THEN BEGIN
  739 + dustem_grid_saved=!dustem_grid
  740 + 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'
  741 + save,all_filters,use_fit_filters,dustem_grid_saved,stellar_normalizations,stats_st,EBVs,F0s,F0s_predicted,G0s,G0s_predicted, $
  742 + Ygrain1,Ygrain2,Ygrain3,facts,chi2s,rchi2s,dF0s,dYgrain1,dYgrain2,dYgrain3, $
  743 + F0_hit,Ygrain1_hit,Ygrain2_hit,Ygrain3_hit,NH_values,best_fit_seds,observed_seds,file=file_save
  744 + message,'Wrote '+file_save,/continue
  745 +ENDIF
  746 +
  747 +;stop
  748 +;=== make maps of parameters
  749 +EBV_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  750 +chi2_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  751 +rchi2_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  752 +F0_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  753 +G0_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  754 +F0_predicted_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  755 +G0_predicted_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  756 +Yg1_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  757 +Yg2_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  758 +Yg3_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  759 +dF0_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  760 +dYg1_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  761 +dYg2_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  762 +dYg3_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  763 +fact_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  764 +F0_hit_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  765 +Yg1_hit_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  766 +Yg2_hit_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  767 +Yg3_hit_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  768 +NH_values_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  769 +norm_stars_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
  770 +;residual maps
  771 +residual_maps=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),Nfilters)+la_undef()
  772 +residual_maps_perc=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),Nfilters)+la_undef()
  773 +
  774 +FOR vid=0LL,Nvor-1 DO BEGIN
  775 + IF vid mod 100 EQ 0 THEN message,'Done maps :'+strtrim((1.*vid)/Nvor*100.,2)+' %',/continue
  776 + ;ind=where(voronoi_id EQ vid)
  777 + IF ptr_valid(all_seds_indices[vid]) THEN BEGIN
  778 + IF ptr_valid(best_fit_seds[vid]) THEN BEGIN
  779 + diffs=la_sub(*best_fit_seds[vid],*observed_seds[vid])
  780 + diffs_perc=la_mul(la_div(diffs,*observed_seds[vid]),100.) ;difference in percent
  781 + ;stop
  782 + FOR i=0L,Nfilters-1 DO BEGIN
  783 + residual_map=residual_maps[*,*,i]
  784 + residual_map[*all_seds_indices[vid]]=diffs[i]
  785 + residual_maps[*,*,i]=residual_map
  786 + residual_map=residual_maps_perc[*,*,i]
  787 + residual_map[*all_seds_indices[vid]]=diffs_perc[i]
  788 + residual_maps_perc[*,*,i]=residual_map
  789 + ENDFOR
  790 + ENDIF
  791 + EBV_map[*all_seds_indices[vid]]=EBVs[vid]
  792 + F0_map[*all_seds_indices[vid]]=F0s[vid]
  793 + G0_map[*all_seds_indices[vid]]=G0s[vid]
  794 + F0_predicted_map[*all_seds_indices[vid]]=F0s_predicted[vid]
  795 + G0_predicted_map[*all_seds_indices[vid]]=G0s_predicted[vid]
  796 + Yg1_map[*all_seds_indices[vid]]=Ygrain1[vid]
  797 + Yg2_map[*all_seds_indices[vid]]=Ygrain2[vid]
  798 + Yg3_map[*all_seds_indices[vid]]=Ygrain3[vid]
  799 + dF0_map[*all_seds_indices[vid]]=dF0s[vid]
  800 + dYg1_map[*all_seds_indices[vid]]=dYgrain1[vid]
  801 + dYg2_map[*all_seds_indices[vid]]=dYgrain2[vid]
  802 + dYg3_map[*all_seds_indices[vid]]=dYgrain3[vid]
  803 + chi2_map[*all_seds_indices[vid]]=chi2s[vid]
  804 + rchi2_map[*all_seds_indices[vid]]=rchi2s[vid]
  805 + fact_map[*all_seds_indices[vid]]=facts[vid]
  806 + F0_hit_map[*all_seds_indices[vid]]=F0_hit[vid]
  807 + Yg1_hit_map[*all_seds_indices[vid]]=Ygrain1_hit[vid]
  808 + Yg2_hit_map[*all_seds_indices[vid]]=Ygrain2_hit[vid]
  809 + Yg3_hit_map[*all_seds_indices[vid]]=Ygrain3_hit[vid]
  810 + NH_values_map[*all_seds_indices[vid]]=NH_values[vid]
  811 + norm_stars_map[*all_seds_indices[vid]]=stellar_normalizations[vid]
  812 + IF not keyword_set(resolution_filter) AND not keyword_set(resolution_value)THEN BEGIN
  813 + Muse_factors_map[*all_seds_indices[vid]]=Muse_factors[vid]
  814 + ENDIF
  815 + ENDIF
  816 +ENDFOR
  817 +
  818 +IF keyword_set(save) THEN BEGIN
  819 +; dustem_grid_saved=!dustem_grid
  820 +; file_save=data_dir+use_source_name+use_model+'_JWST_G0_YPAH_YVSG_on-voronoi_with-ISRFclasses'+norm_str+mathis_str+reso_str+fit_G0_str+set_yvsg_str+'.sav'
  821 +; save,all_filters,use_fit_filters,dustem_grid_saved,stellar_normalizations,stats_st,EBVs,F0s,F0s_predicted,G0s,G0s_predicted, $
  822 +; Ygrain1,Ygrain2,Ygrain3,facts,chi2s,rchi2s,dF0s,dYgrain1,dYgrain2,dYgrain3, $
  823 +; F0_hit,Ygrain1_hit,Ygrain2_hit,Ygrain3_hit,Yvsg_hit,NH_values,best_fit_seds,observed_seds,file=file_save
  824 +; message,'Wrote '+file_save,/continue
  825 + 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'
  826 + save,NHCO,NHI,NHmap,norm_stars_map,EBV_map,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
  827 + message,'Wrote '+file_save,/continue
  828 +ENDIF
  829 +
  830 +from_restore:
  831 +
  832 +IF not keyword_set(fit_G0) THEN BEGIN
  833 + fit_G0_str='_GOfixed'
  834 +ENDIF ELSE BEGIN
  835 + fit_G0_str=''
  836 +ENDELSE
  837 +
  838 +IF keyword_set(set_yvsg) THEN BEGIN
  839 + set_yvsg_str='_YVSGfixed'
  840 +ENDIF ELSE BEGIN
  841 + set_yvsg_str=''
  842 +ENDELSE
  843 +
  844 +;stop
  845 +
  846 +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'
  847 +st_info=file_info(file)
  848 +IF st_info.exists NE 1 THEN BEGIN
  849 + message,'Could not find '+file,/continue
  850 + stop
  851 +ENDIF
  852 +restore,file,/verb
  853 +; % RESTORE: Restored variable: STATS_ST.
  854 +; % RESTORE: Restored variable: F0S.
  855 +; % RESTORE: Restored variable: F0S_PREDICTED.
  856 +; % RESTORE: Restored variable: YPAHS.
  857 +; % RESTORE: Restored variable: YVSGS.
  858 +; % RESTORE: Restored variable: FACTS.
  859 +; % RESTORE: Restored variable: CHI2S.
  860 +; % RESTORE: Restored variable: RCHI2S.
  861 +; % RESTORE: Restored variable: DF0S.
  862 +; % RESTORE: Restored variable: DYPAHS.
  863 +; % RESTORE: Restored variable: DYVSGS.
  864 +; % RESTORE: Restored variable: F0_HIT.
  865 +; % RESTORE: Restored variable: YPAH_HIT.
  866 +; % RESTORE: Restored variable: YVSG_HIT.
  867 +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'
  868 +st_info=file_info(file)
  869 +IF st_info.exists NE 1 THEN BEGIN
  870 + message,'Could not find '+file,/continue
  871 + stop
  872 +ENDIF
  873 +restore,file,/verb
  874 +; % RESTORE: Restored variable: NHCO.
  875 +; % RESTORE: Restored variable: NHI.
  876 +; % RESTORE: Restored variable: NHMAP.
  877 +; % RESTORE: Restored variable: NORM_STARS_MAP.
  878 +; % RESTORE: Restored variable: EBV_MAP.
  879 +; % RESTORE: Restored variable: F0_MAP.
  880 +; % RESTORE: Restored variable: G0_MAP.
  881 +; % RESTORE: Restored variable: F0_PREDICTED_MAP.
  882 +; % RESTORE: Restored variable: G0_PREDICTED_MAP.
  883 +; % RESTORE: Restored variable: YG1_MAP.
  884 +; % RESTORE: Restored variable: YG2_MAP.
  885 +; % RESTORE: Restored variable: YG3_MAP.
  886 +; % RESTORE: Restored variable: FACT_MAP.
  887 +; % RESTORE: Restored variable: CHI2_MAP.
  888 +; % RESTORE: Restored variable: RCHI2_MAP.
  889 +; % RESTORE: Restored variable: DF0_MAP.
  890 +; % RESTORE: Restored variable: DYG1_MAP.
  891 +; % RESTORE: Restored variable: DYG2_MAP.
  892 +; % RESTORE: Restored variable: DYG3_MAP.
  893 +; % RESTORE: Restored variable: F0_HIT_MAP.
  894 +; % RESTORE: Restored variable: YG1_HIT_MAP.
  895 +; % RESTORE: Restored variable: YG2_HIT_MAP.
  896 +; % RESTORE: Restored variable: YG3_HIT_MAP.
  897 +; % RESTORE: Restored variable: NH_VALUES_MAP.
  898 +; % RESTORE: Restored variable: RESIDUAL_MAPS.
  899 +; % RESTORE: Restored variable: RESIDUAL_MAPS_PERC.
  900 +file=data_dir+use_source_name+'_isrf_classes_voronoi'+reso_str+'.sav'
  901 +st_info=file_info(file)
  902 +IF st_info.exists NE 1 THEN BEGIN
  903 + message,'Could not find '+file,/continue
  904 + stop
  905 +ENDIF
  906 +restore,file,/verb
  907 +;% RESTORE: Restored variable: VOR_CLASS.
  908 +file=data_dir+use_source_name+'_ref_header'+reso_str+'.sav'
  909 +st_info=file_info(file)
  910 +IF st_info.exists NE 1 THEN BEGIN
  911 + message,'Could not find '+file,/continue
  912 + stop
  913 +ENDIF
  914 +restore,file,/verb
  915 +;% RESTORE: Restored variable: HREF.
  916 +
  917 +IF defined(dustem_grid_saved) THEN defsysv,'!dustem_grid',dustem_grid_saved
  918 +;=== get the stellar mask
  919 +file=data_dir+use_source_name+'_star_mask'+reso_str+'.sav'
  920 +restore,file,/verb
  921 +;% RESTORE: Restored variable: STAR_MASK.
  922 +;% RESTORE: Restored variable: HSTARMASK.
  923 +;% RESTORE: Restored variable: STAR_MASK_REF.
  924 +
  925 +;stop
  926 +
  927 +;==== set masked values in all maps
  928 +ind_masked=where(star_mask_ref NE 0 and star_mask_ref NE la_undef(),count_masked)
  929 +IF count_masked NE 0 THEN BEGIN
  930 + FOR i=0L,n_elements(all_filters)-1 DO BEGIN
  931 + map=residual_maps_perc[*,*,i]
  932 + map[ind_masked]=la_undef()
  933 + residual_maps_perc[*,*,i]=map
  934 + map=residual_maps[*,*,i]
  935 + map[ind_masked]=la_undef()
  936 + residual_maps[*,*,i]=map
  937 + ENDFOR
  938 + NHmap[ind_masked]=la_undef()
  939 + norm_stars_map[ind_masked]=la_undef()
  940 + EBV_map[ind_masked]=la_undef()
  941 + F0_map[ind_masked]=la_undef()
  942 + G0_map[ind_masked]=la_undef()
  943 + F0_predicted_map[ind_masked]=la_undef()
  944 + G0_predicted_map[ind_masked]=la_undef()
  945 + Yg1_map[ind_masked]=la_undef()
  946 + Yg2_map[ind_masked]=la_undef()
  947 + Yg3_map[ind_masked]=la_undef()
  948 + fact_map[ind_masked]=la_undef()
  949 + chi2_map[ind_masked]=la_undef()
  950 + rchi2_map[ind_masked]=la_undef()
  951 + dF0_map[ind_masked]=la_undef()
  952 + dYg1_map[ind_masked]=la_undef()
  953 + dYg2_map[ind_masked]=la_undef()
  954 + dYg3_map[ind_masked]=la_undef()
  955 + F0_hit_map[ind_masked]=la_undef()
  956 + Yg1_hit_map[ind_masked]=la_undef()
  957 + Yg2_hit_map[ind_masked]=la_undef()
  958 + Yg3_hit_map[ind_masked]=la_undef()
  959 + NH_values_map[ind_masked]=la_undef()
  960 +ENDIF
  961 +
  962 +cleanplot
  963 +win=0L
  964 +
  965 +;This is to allow dustem to return default parameter values
  966 +dustem_init,model=use_model
  967 +
  968 +Nx=sxpar(href,'NAXIS1')
  969 +Ny=sxpar(href,'NAXIS2')
  970 +;=== This is to make a common mask for all maps
  971 +mask=fltarr(Nx,Ny)
  972 +mask[*]=1
  973 +ind=where(G0_map EQ la_undef(),count)
  974 +IF count NE 0 THEN BEGIN
  975 + mask[ind]=0
  976 +ENDIF
  977 +ind_bad=where(mask NE 1,count_bad)
  978 +;=== apply mask to NH maps
  979 +IF count_bad NE 0 THEN BEGIN
  980 + NHmap[ind_bad]=la_undef()
  981 + F0_map[ind_bad]=la_undef()
  982 + F0_predicted_map[ind_bad]=la_undef()
  983 +ENDIF
  984 +
  985 +;=== make an overlay with contours of the total NH
  986 +;=== We smooth the NH map to 5 arcesec, or else, too many contours are shown.
  987 +IF reso LE 5./60.^2 THEN BEGIN
  988 + use_NHmap=degrade_res(NHmap,href,reso,5./60.^2,hout)
  989 + use_NHmap=project2(hout,use_NHmap,href,/silent)
  990 +ENDIF ELSE BEGIN
  991 + use_NHmap=NHmap
  992 +ENDELSE
  993 +levs=[la_undef(),-0.5,-0.25,0.,0.25,0.5,0.75,1.]
  994 +;levs=[-0.5,-0.25,0.,0.25,0.5,0.75,1.]
  995 +NNHmap=fltarr(Nx,Ny,2)
  996 +NNHmap[*,*,0]=la_log10(use_NHmap)
  997 +NNHmap[*,*,1]=la_log10(use_NHmap)
  998 +cont_st=create_contour_st(NNHmap,lev1=levs)
  999 +cont_st.(0).color=0
  1000 +cont_st.(0).thick=2.
  1001 +cont_st.(1).color=100
  1002 +cont_st.(1).thick=1.0
  1003 +
  1004 +
  1005 +;===== plot nhits vs isrf_class
  1006 +Ncl=n_elements(stats_st)
  1007 +F0_hitp_per_class=lonarr(Ncl)
  1008 +F0_hitm_per_class=lonarr(Ncl)
  1009 +Ygrain1_hitp_per_class=lonarr(Ncl)
  1010 +Ygrain1_hitm_per_class=lonarr(Ncl)
  1011 +Ygrain2_hitp_per_class=lonarr(Ncl)
  1012 +Ygrain2_hitm_per_class=lonarr(Ncl)
  1013 +Ygrain3_hitp_per_class=lonarr(Ncl)
  1014 +Ygrain3_hitm_per_class=lonarr(Ncl)
  1015 +FOR i=0L,Ncl-1 DO BEGIN
  1016 + ind=where(vor_class EQ stats_st[i].isrf_class AND F0_hit EQ 1,count)
  1017 + F0_hitp_per_class[i]=count
  1018 + ind=where(vor_class EQ stats_st[i].isrf_class AND F0_hit EQ -1,count)
  1019 + F0_hitm_per_class[i]=count
  1020 +
  1021 + ind=where(vor_class EQ stats_st[i].isrf_class AND Ygrain1_hit EQ 1,count)
  1022 + Ygrain1_hitp_per_class[i]=count
  1023 + ind=where(vor_class EQ stats_st[i].isrf_class AND Ygrain1_hit EQ -1,count)
  1024 + Ygrain1_hitm_per_class[i]=count
  1025 +
  1026 + ind=where(vor_class EQ stats_st[i].isrf_class AND Ygrain2_hit EQ 1,count)
  1027 + Ygrain2_hitp_per_class[i]=count
  1028 + ind=where(vor_class EQ stats_st[i].isrf_class AND Ygrain2_hit EQ -1,count)
  1029 + Ygrain2_hitm_per_class[i]=count
  1030 +
  1031 + ind=where(vor_class EQ stats_st[i].isrf_class AND Ygrain3_hit EQ 1,count)
  1032 + Ygrain3_hitp_per_class[i]=count
  1033 + ind=where(vor_class EQ stats_st[i].isrf_class AND Ygrain3_hit EQ -1,count)
  1034 + Ygrain3_hitm_per_class[i]=count
  1035 +ENDFOR
  1036 +window,win & win=win+1
  1037 +!p.multi=[0,2,2]
  1038 +xtit='ISRF class'
  1039 +ytit='N bins'
  1040 +tit='F0_hit'
  1041 +cgplot,stats_st.isrf_class,F0_hitp_per_class,color='red',title=tit,xtit=xtit,ytit=ytit
  1042 +cgoplot,stats_st.isrf_class,F0_hitm_per_class,color='blue'
  1043 +tit='Ygrain1_hit'
  1044 +cgplot,stats_st.isrf_class,Ygrain1_hitp_per_class,color='red',title=tit,xtit=xtit,ytit=ytit
  1045 +cgoplot,stats_st.isrf_class,Ygrain1_hitm_per_class,color='blue'
  1046 +tit='Ygrain2_hit'
  1047 +cgplot,stats_st.isrf_class,Ygrain2_hitp_per_class,color='red',title=tit,xtit=xtit,ytit=ytit
  1048 +cgoplot,stats_st.isrf_class,Ygrain2_hitm_per_class,color='blue'
  1049 +tit='Ygrain3_hit'
  1050 +cgplot,stats_st.isrf_class,Ygrain3_hitp_per_class,color='red',title=tit,xtit=xtit,ytit=ytit
  1051 +cgoplot,stats_st.isrf_class,Ygrain3_hitm_per_class,color='blue'
  1052 +
  1053 +window,win & win=win+1
  1054 +;==== plot hits vs F0hit
  1055 +!p.multi=[0,2,2]
  1056 +xtit='F0_hit'
  1057 +ytit='grain hit'
  1058 +tit='Ygrain1_hit'
  1059 +xr=[-1.1,1.1]
  1060 +cgplot,F0_hitp_per_class,Ygrain1_hitp_per_class,color='red',title=tit,xtit=xtit,ytit=ytit,psym='Open Circle',xrange=xrange,/xsty
  1061 +cgoplot,F0_hitp_per_class,Ygrain1_hitm_per_class,color='blue',psym='Open Circle'
  1062 +tit='Ygrain2_hit'
  1063 +cgplot,F0_hitp_per_class,Ygrain2_hitp_per_class,color='red',title=tit,xtit=xtit,ytit=ytit,psym='Open Circle',xrange=xrange,/xsty
  1064 +cgoplot,F0_hitp_per_class,Ygrain2_hitm_per_class,color='blue',psym='Open Circle'
  1065 +tit='Ygrain3_hit'
  1066 +cgplot,F0_hitp_per_class,Ygrain3_hitp_per_class,color='red',title=tit,xtit=xtit,ytit=ytit,psym='Open Circle',xrange=xrange,/xsty
  1067 +cgoplot,F0_hitp_per_class,Ygrain3_hitm_per_class,color='blue',psym='Open Circle'
  1068 +
  1069 +;stop
  1070 +
  1071 +;===== histograms
  1072 +window,win & win=win+1
  1073 +ind=where(vor_class NE la_undef() and vor_class NE 0)
  1074 +res=histogram(vor_class[ind],locations=xv)
  1075 +cgplot,xv,res,psym=10,title='ISRF class histogram',xtit='class',ytit='Number',/ylog,yrange=[1,max(res)],/ysty
  1076 +
  1077 +;===== histograms 4 pannels
  1078 +range_fact=15.
  1079 +F0_range_fact=10. ;%
  1080 +
  1081 +window,win,xsize=700,ysize=900 & win=win+1
  1082 +!p.multi=[0,2,2]
  1083 +
  1084 +ind=where(F0s NE la_undef() AND F0s GT 0,count)
  1085 +res=histogram(alog10(F0s[ind]),locations=xv,Nbins=100)
  1086 +yrange=[1,max(res)]
  1087 +F0_cent=1.
  1088 +parval_min=((*!dustem_grid).PMIN_VALUES)[0]
  1089 +parval_max=((*!dustem_grid).PMAX_VALUES)[0]
  1090 +;xrange=alog10([F0_min/F0_range_fact,F0_max*F0_range_fact])
  1091 +xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
  1092 +cgplot,xv,res,psym=10,title='F0s histogram',xtit='log10(F0)',ytit='Number',/ylog,yrange=yrange,/ysty,xrange=xrange,/xsty
  1093 +cgoplot,replicate(alog10(F0_cent),2),10^(!y.crange),linestyle=2,color='black'
  1094 +cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
  1095 +cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  1096 +
  1097 +ind=where(Ygrain1 NE la_undef())
  1098 +default_value=dustem_get_param_values_default('(*!dustem_params).grains[0].MDUST_O_MH')
  1099 +parval_min=((*!dustem_grid).PMIN_VALUES)[1]
  1100 +parval_max=((*!dustem_grid).PMAX_VALUES)[1]
  1101 +;xrange=alog10([default_value/range_fact,default_value*range_fact])
  1102 +xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
  1103 +yrange=[1,max(res)]
  1104 +res=histogram(alog10(Ygrain1[ind]),locations=xv,Nbins=100)
  1105 +cgplot,xv,res,psym=10,title='Y Grain1 histogram',xtit='log10(Y)',ytit='Number',/ylog,/ysty,yrange=yrange,xrange=xrange,/xsty
  1106 +cgoplot,replicate(alog10(default_value),2),10^(!y.crange),linestyle=2,color='black'
  1107 +cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
  1108 +cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  1109 +
  1110 +ind=where(Ygrain2 NE la_undef(),count)
  1111 +IF count NE 0 THEN BEGIN
  1112 + default_value=dustem_get_param_values_default('(*!dustem_params).grains[1].MDUST_O_MH')
  1113 + parval_min=((*!dustem_grid).PMIN_VALUES)[2]
  1114 + parval_max=((*!dustem_grid).PMAX_VALUES)[2]
  1115 + ;xrange=alog10([default_value/range_fact,default_value*range_fact])
  1116 + xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
  1117 + yrange=[1,max(res)]
  1118 + res=histogram(alog10(Ygrain2[ind]),locations=xv,Nbins=100)
  1119 + cgplot,xv,res,psym=10,title='Y Grain2 histogram',xtit='log10(Y)',ytit='Number',/ylog,/ysty,yrange=yrange,xrange=xrange,/xsty
  1120 + cgoplot,replicate(alog10(default_value),2),10^(!y.crange),linestyle=2,color='black'
  1121 + cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
  1122 + cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  1123 +ENDIF
  1124 +
  1125 +ind=where(Ygrain3 NE la_undef(),count3)
  1126 +IF count3 NE 0 THEN BEGIN
  1127 + default_value=dustem_get_param_values_default('(*!dustem_params).grains[2].MDUST_O_MH')
  1128 + ;=== IF PAH+ are in table
  1129 + parval_min=((*!dustem_grid).PMIN_VALUES)[3]
  1130 + parval_max=((*!dustem_grid).PMAX_VALUES)[3]
  1131 + ;=== IF PAH+ are not in table
  1132 + ;parval_min=((*!dustem_grid).PMIN_VALUES)[2]
  1133 + ;parval_max=((*!dustem_grid).PMAX_VALUES)[2]
  1134 + ;xrange=alog10([default_value/range_fact,default_value*range_fact])
  1135 + xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
  1136 + yrange=[1,max(res)]
  1137 + res=histogram(alog10(Ygrain3[ind]),locations=xv,Nbins=100)
  1138 + cgplot,xv,res,psym=10,title='Y Grain3 histogram',xtit='log10(Y)',ytit='Number',/ylog,/ysty,yrange=yrange,xrange=xrange,/xsty
  1139 + cgoplot,replicate(alog10(default_value),2),10^(!y.crange),linestyle=2,color='black'
  1140 + cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
  1141 + cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  1142 +ENDIF
  1143 +
  1144 +IF keyword_set(do_postscripts) THEN BEGIN
  1145 + cleanplot
  1146 +
  1147 + ps_file=ps_dir+use_source_name+'_histo_F0'+reso_str+fit_G0_str+mathis_str+'.ps'
  1148 + set_plot,'PS'
  1149 + device,filename=ps_file
  1150 +
  1151 + ind=where(F0s NE la_undef() AND F0s GT 0,count)
  1152 + res=histogram(alog10(F0s[ind]),locations=xv,Nbins=100)
  1153 + yrange=[1,max(res)]
  1154 + F0_cent=1.
  1155 + parval_min=((*!dustem_grid).PMIN_VALUES)[0]
  1156 + parval_max=((*!dustem_grid).PMAX_VALUES)[0]
  1157 + ;xrange=alog10([F0_min/F0_range_fact,F0_max*F0_range_fact])
  1158 + xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
  1159 + cgplot,xv,res,psym=10,title='F0s histogram',xtit='log10(F0)',ytit='Number',/ylog,yrange=yrange,/ysty,xrange=xrange,/xsty
  1160 + cgoplot,replicate(alog10(F0_cent),2),10^(!y.crange),linestyle=2,color='black'
  1161 + cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
  1162 + cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  1163 +
  1164 + device,/close
  1165 + set_plot,'X'
  1166 + message,'Wrote '+ps_file,/continue
  1167 + ;stop
  1168 +
  1169 + ps_file=ps_dir+use_source_name+'_histo_YG1'+reso_str+fit_G0_str+mathis_str+'.ps'
  1170 + set_plot,'PS'
  1171 + device, filename=ps_file
  1172 +
  1173 + ind=where(Ygrain1 NE la_undef())
  1174 + default_value=dustem_get_param_values_default('(*!dustem_params).grains[0].MDUST_O_MH')
  1175 + parval_min=((*!dustem_grid).PMIN_VALUES)[1]
  1176 + parval_max=((*!dustem_grid).PMAX_VALUES)[1]
  1177 + ;xrange=alog10([default_value/range_fact,default_value*range_fact])
  1178 + xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
  1179 + yrange=[1,max(res)]
  1180 + res=histogram(alog10(Ygrain1[ind]),locations=xv,Nbins=100)
  1181 + cgplot,xv,res,psym=10,title='Y Grain1 histogram',xtit='log10(Y)',ytit='Number',/ylog,/ysty,yrange=yrange,xrange=xrange,/xsty
  1182 + cgoplot,replicate(alog10(default_value),2),10^(!y.crange),linestyle=2,color='black'
  1183 + cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
  1184 + cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  1185 +
  1186 + device,/close
  1187 + set_plot,'X'
  1188 +
  1189 + ps_file=ps_dir+use_source_name+'_histo_YG2'+reso_str+fit_G0_str+mathis_str+'.ps'
  1190 + set_plot,'PS'
  1191 + device, filename=ps_file
  1192 +
  1193 + ind=where(Ygrain2 NE la_undef())
  1194 + default_value=dustem_get_param_values_default('(*!dustem_params).grains[1].MDUST_O_MH')
  1195 + parval_min=((*!dustem_grid).PMIN_VALUES)[2]
  1196 + parval_max=((*!dustem_grid).PMAX_VALUES)[2]
  1197 + ;xrange=alog10([default_value/range_fact,default_value*range_fact])
  1198 + xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
  1199 + yrange=[1,max(res)]
  1200 + res=histogram(alog10(Ygrain2[ind]),locations=xv,Nbins=100)
  1201 + cgplot,xv,res,psym=10,title='Y Grain2 histogram',xtit='log10(Y)',ytit='Number',/ylog,/ysty,yrange=yrange,xrange=xrange,/xsty
  1202 + cgoplot,replicate(alog10(default_value),2),10^(!y.crange),linestyle=2,color='black'
  1203 + cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
  1204 + cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  1205 +
  1206 + device,/close
  1207 + set_plot,'X'
  1208 +
  1209 + ps_file=ps_dir+use_source_name+'_histo_YG3'+reso_str+fit_G0_str+mathis_str+'.ps'
  1210 + set_plot,'PS'
  1211 + device, filename=ps_file
  1212 +
  1213 + ind=where(Ygrain3 NE la_undef(),count3)
  1214 + IF count3 NE 0 THEN BEGIN
  1215 + default_value=dustem_get_param_values_default('(*!dustem_params).grains[2].MDUST_O_MH')
  1216 + parval_min=((*!dustem_grid).PMIN_VALUES)[3]
  1217 + parval_max=((*!dustem_grid).PMAX_VALUES)[3]
  1218 + ;xrange=alog10([default_value/range_fact,default_value*range_fact])
  1219 + xrange=enlarge_range(alog10([parval_min,parval_max]),0.1)
  1220 + yrange=[1,max(res)]
  1221 + res=histogram(alog10(Ygrain3[ind]),locations=xv,Nbins=100)
  1222 + cgplot,xv,res,psym=10,title='Y Grain3 histogram',xtit='log10(Y)',ytit='Number',/ylog,/ysty,yrange=yrange,xrange=xrange,/xsty
  1223 + cgoplot,replicate(alog10(default_value),2),10^(!y.crange),linestyle=2,color='black'
  1224 + cgoplot,replicate(alog10(parval_min),2),10^(!y.crange),linestyle=2,color='blue'
  1225 + cgoplot,replicate(alog10(parval_max),2),10^(!y.crange),linestyle=2,color='red'
  1226 + ENDIF
  1227 +
  1228 + device,/close
  1229 + set_plot,'X'
  1230 +
  1231 +ENDIF
  1232 +;stop
  1233 +
  1234 +;==== show maps
  1235 +obp=[1.1,0.,1.15,1]
  1236 +cleanplot
  1237 +image_color_table='jpbloadct,/whitebackground'
  1238 +
  1239 +window,win,xsize=800,ysize=900 & win=win+1
  1240 +imrange=[-1,1.5]
  1241 +title='NH map'
  1242 +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
  1243 +IF keyword_set(do_postscripts) THEN BEGIN
  1244 + ps_file=ps_dir+use_source_name+'_NHmap'+reso_str+fit_G0_str+mathis_str+'.ps'
  1245 + 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, $
  1246 + postscript=ps_file,orientation='portrait'
  1247 + 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, $
  1248 + postscript=ps_file,orientation='portrait'
  1249 + message,'Wrote '+ps_file,/continue
  1250 +ENDIF
  1251 +;stop
  1252 +
  1253 +window,win,xsize=800,ysize=900 & win=win+1
  1254 +imrange=1.+[-0.5,+0.5]
  1255 +title='stellar normalization'
  1256 +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
  1257 +IF keyword_set(do_postscripts) THEN BEGIN
  1258 + ps_file=ps_dir+use_source_name+'_norm_stars_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  1259 + 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, $
  1260 + postscript=ps_file,orientation='portrait',/nologo
  1261 + message,'Wrote '+ps_file,/continue
  1262 +ENDIF
  1263 +
  1264 +window,win,xsize=800,ysize=900 & win=win+1
  1265 +imrange=[0.1,0.25]
  1266 +title='E(B-V)'
  1267 +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
  1268 +IF keyword_set(do_postscripts) THEN BEGIN
  1269 + ps_file=ps_dir+use_source_name+'_EBV_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  1270 + 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, $
  1271 + postscript=ps_file,orientation='portrait',/nologo
  1272 + message,'Wrote '+ps_file,/continue
  1273 +ENDIF
  1274 +
  1275 +window,win,xsize=800,ysize=900 & win=win+1
  1276 +imrange=[-1.5,1.2]
  1277 +title='log10(F0_predicted)'
  1278 +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
  1279 +IF keyword_set(do_postscripts) THEN BEGIN
  1280 + ps_file=ps_dir+use_source_name+'_F0_predicted_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  1281 + 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, $
  1282 + postscript=ps_file,orientation='portrait',/nologo
  1283 + message,'Wrote '+ps_file,/continue
  1284 +ENDIF
  1285 +
  1286 +window,win,xsize=800,ysize=900 & win=win+1
  1287 +imrange=[-1.,+1.0] ;with Herschel data
  1288 +title='log10(F0)'
  1289 +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
  1290 +IF keyword_set(do_postscripts) THEN BEGIN
  1291 + ps_file=ps_dir+use_source_name+'_F0_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  1292 + 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, $
  1293 + postscript=ps_file,orientation='portrait',/nologo
  1294 + message,'Wrote '+ps_file,/continue
  1295 +ENDIF
  1296 +
  1297 +window,win,xsize=800,ysize=900 & win=win+1
  1298 +imrange=[-1.,+1.0]
  1299 +title='log10(G0)'
  1300 +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
  1301 +IF keyword_set(do_postscripts) THEN BEGIN
  1302 + ps_file=ps_dir+use_source_name+'_G0_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  1303 + 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, $
  1304 + postscript=ps_file,orientation='portrait',/nologo
  1305 + message,'Wrote '+ps_file,/continue
  1306 +ENDIF
  1307 +
  1308 +window,win,xsize=800,ysize=900 & win=win+1
  1309 +imrange=[-1.,+1.0]
  1310 +title='log10(G0_predicted)'
  1311 +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
  1312 +IF keyword_set(do_postscripts) THEN BEGIN
  1313 + ps_file=ps_dir+use_source_name+'_G0_predicted_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  1314 + 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, $
  1315 + postscript=ps_file,orientation='portrait',/nologo
  1316 + message,'Wrote '+ps_file,/continue
  1317 +ENDIF
  1318 +
  1319 +window,win,xsize=800,ysize=900 & win=win+1
  1320 +;imrange=[-3.5,-2.]
  1321 +imrange=[-3.8,-2.9]
  1322 +title='log10(YG1)'
  1323 +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
  1324 +IF keyword_set(do_postscripts) THEN BEGIN
  1325 + ps_file=ps_dir+use_source_name+'_YG1_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  1326 + 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, $
  1327 + postscript=ps_file,orientation='portrait',/nologo
  1328 + message,'Wrote '+ps_file,/continue
  1329 +ENDIF
  1330 +;stop
  1331 +
  1332 +window,win,xsize=800,ysize=900 & win=win+1
  1333 +;imrange=[-4,-2.]
  1334 +imrange=[-3.8,-2.9]
  1335 +title='log10(YG2)'
  1336 +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
  1337 +IF keyword_set(do_postscripts) THEN BEGIN
  1338 + ps_file=ps_dir+use_source_name+'_YG2_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  1339 + 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, $
  1340 + postscript=ps_file,orientation='portrait',/nologo
  1341 + message,'Wrote '+ps_file,/continue
  1342 +ENDIF
  1343 +;stop
  1344 +
  1345 +window,win,xsize=800,ysize=900 & win=win+1
  1346 +;imrange=[-4,-2.]
  1347 +imrange=[-3.8,-2.5]
  1348 +title='log10(YG3)'
  1349 +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
  1350 +IF keyword_set(do_postscripts) THEN BEGIN
  1351 + ps_file=ps_dir+use_source_name+'_YG3_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  1352 + 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, $
  1353 + postscript=ps_file,orientation='portrait',/nologo
  1354 + message,'Wrote '+ps_file,/continue
  1355 +ENDIF
  1356 +
  1357 +;window,win,xsize=800,ysize=900 & win=win+1
  1358 +;imrange=[1.5,2.2]
  1359 +;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
  1360 +
  1361 +window,win,xsize=800,ysize=900 & win=win+1
  1362 +;imrange=[1.2,1.8]
  1363 +;imrange=[-0.1,2.]
  1364 +imrange=[-1.,2.]
  1365 +title='log10(rChi2/med(rchi2))'
  1366 +map=rchi2_map
  1367 +;ind=where(rchi2_map EQ la_min(rchi2_map))
  1368 +;map=la_div(map,rchi2_map[ind[0]])
  1369 +map=la_div(map,la_median(rchi2_map))
  1370 +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
  1371 +IF keyword_set(do_postscripts) THEN BEGIN
  1372 + ps_file=ps_dir+use_source_name+'_rchi2_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  1373 + 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, $
  1374 + postscript=ps_file,orientation='portrait',/nologo
  1375 + message,'Wrote '+ps_file,/continue
  1376 +ENDIF
  1377 +;stop
  1378 +
  1379 +;window,win,xsize=800,ysize=900 & win=win+1
  1380 +;imrange=[-3,2] ;normalize=0
  1381 +;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
  1382 +
  1383 +;stop
  1384 +window,win,xsize=800,ysize=900 & win=win+1
  1385 +imrange=[-2,2]
  1386 +title='F0 hit'
  1387 +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
  1388 +IF keyword_set(do_postscripts) THEN BEGIN
  1389 + ps_file=ps_dir+use_source_name+'_F0_hit_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  1390 + 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, $
  1391 + postscript=ps_file,orientation='portrait',/nologo
  1392 + message,'Wrote '+ps_file,/continue
  1393 +ENDIF
  1394 +
  1395 +window,win,xsize=800,ysize=900 & win=win+1
  1396 +imrange=[-2,2]
  1397 +title='YG1 hit'
  1398 +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
  1399 +IF keyword_set(do_postscripts) THEN BEGIN
  1400 + ps_file=ps_dir+use_source_name+'_YG1_hit_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  1401 + 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, $
  1402 + postscript=ps_file,orientation='portrait',/nologo
  1403 + message,'Wrote '+ps_file,/continue
  1404 +ENDIF
  1405 +
  1406 +window,win,xsize=800,ysize=900 & win=win+1
  1407 +imrange=[-2,2]
  1408 +title='YG2 hit'
  1409 +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
  1410 +IF keyword_set(do_postscripts) THEN BEGIN
  1411 + ps_file=ps_dir+use_source_name+'_YG2_hit_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  1412 + 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, $
  1413 + postscript=ps_file,orientation='portrait',/nologo
  1414 + message,'Wrote '+ps_file,/continue
  1415 +ENDIF
  1416 +
  1417 +window,win,xsize=800,ysize=900 & win=win+1
  1418 +imrange=[-2,2]
  1419 +title='YG3 hit'
  1420 +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
  1421 +IF keyword_set(do_postscripts) THEN BEGIN
  1422 + ps_file=ps_dir+use_source_name+'_YG3_hit_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  1423 + 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, $
  1424 + postscript=ps_file,orientation='portrait',/nologo
  1425 + message,'Wrote '+ps_file,/continue
  1426 +ENDIF
  1427 +
  1428 +;window,win,xsize=800,ysize=900 & win=win+1
  1429 +;imrange=0
  1430 +;imrange=[-2.0,-1.45]
  1431 +;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
  1432 +
  1433 +window,win,xsize=800,ysize=900 & win=win+1
  1434 +rap_F0=la_div(F0_map,F0_predicted_map)
  1435 +imrange=[-0.1,1.5]
  1436 +tit='log(G0/G0_predicted)'
  1437 +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
  1438 +IF keyword_set(do_postscripts) THEN BEGIN
  1439 + ps_file=ps_dir+use_source_name+'_G0_ratio_map'+reso_str+fit_G0_str+mathis_str+'.ps'
  1440 + 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, $
  1441 + postscript=ps_file,orientation='portrait',/nologo
  1442 + message,'Wrote '+ps_file,/continue
  1443 +ENDIF
  1444 +
  1445 +;==== Show residual maps
  1446 +window,win,xsize=800,ysize=900 & win=win+1
  1447 +residual_range=[-100,100.]/2.
  1448 +FOR i=0L,n_elements(all_filters)-1 DO BEGIN
  1449 + title='Residual percent Filter '+all_filters[i]+' '+strtrim(dustem_filter2wav(all_filters[i]),2)
  1450 + filter_str='_'+all_filters[i]
  1451 + map=residual_maps_perc[*,*,i]
  1452 + ind=where(finite(map) NE 1,count)
  1453 + IF count NE 0 THEN map[ind]=la_undef()
  1454 + ;print,la_min(map),la_max(map)
  1455 + IF la_min(map) NE la_max(map) THEN use_obp=obp ELSE use_obp=0
  1456 + 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
  1457 + IF keyword_set(do_postscripts) THEN BEGIN
  1458 + ps_file=ps_dir+use_source_name+'_residual_map'+reso_str+fit_G0_str+mathis_str+filter_str+'.ps'
  1459 + 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, $
  1460 + postscript=ps_file,orientation='portrait',/nologo
  1461 + message,'Wrote '+ps_file,/continue
  1462 + ENDIF
  1463 + IF not keyword_set(nostop) THEN stop
  1464 +ENDFOR
  1465 +
  1466 +;=== other plots
  1467 +NH_range=[0.1,1000.]
  1468 +F0_range=[0.1,10]
  1469 +F0_range_pred=[0.001,10]
  1470 +Ypah_range=[8.e-5,2.e-1]
  1471 +Yvsg_range=[8.e-5,2.e-1]
  1472 +chi2_range=[1.e3,1.e7]/1.e2
  1473 +
  1474 +window,win,xsize=800,ysize=900 & win=win+1
  1475 +cgplot,F0s_predicted,F0s,/xlog,/ylog,xrange=F0_range_pred,/xsty,yrange=F0_range,/ysty,psym=3,xtit='F0 predicted',ytit='F0'
  1476 +cgoplot,[1.e-3,1.e3],[1.e-3,1.e3],linesyle=2,color='red'
  1477 +cgoplot,[1.e-3,1.e3],[1.e-3,1.e3]*10.,linesyle=2,color='blue'
  1478 +cgoplot,[1.e-3,1.e3],[1.e-3,1.e3]*100.,linesyle=2,color='blue'
  1479 +
  1480 +window,win,xsize=800,ysize=900 & win=win+1
  1481 +!p.multi=[0,1,4]
  1482 +cgplot,F0s,Ygrain1,/xlog,/ylog,xr=F0_range,/xsty,psym=3,xtit='F0',ytit='Y',title='Grain 1'
  1483 +cgplot,F0s,Ygrain2,/xlog,/ylog,xr=F0_range,/xsty,psym=3,xtit='F0',ytit='Y',title='Grain 2'
  1484 +cgplot,F0s,Ygrain3,/xlog,/ylog,xr=F0_range,/xsty,psym=3,xtit='F0',ytit='Y',title='Grain 3'
  1485 +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'
  1486 +
  1487 +;stop
  1488 +
  1489 +window,win,xsize=800,ysize=900 & win=win+1
  1490 +!p.multi=[0,1,4]
  1491 +NH_range=[0.1,1000.]
  1492 +F0_range=[0.1,10]
  1493 +Ypah_range=[8.e-5,2.e-1]
  1494 +Yvsg_range=[8.e-5,2.e-1]
  1495 +cgplot,NH_values,F0s,/xlog,/ylog,xr=NH_range,yr=F0_range,/xsty,psym=3,xtit='NH',ytit='F0'
  1496 +cgplot,NH_values,Ygrain1,/xlog,/ylog,xr=NH_range,yr=Ypah_range,/xsty,psym=3,xtit='NH',ytit='Y',title='Grain 1'
  1497 +cgplot,NH_values,Ygrain2,/xlog,/ylog,xr=NH_range,yr=Ypah_range,/xsty,psym=3,xtit='NH',ytit='Y',title='Grain 2'
  1498 +cgplot,NH_values,Ygrain3,/xlog,/ylog,xr=NH_range,yr=Yvsg_range,/xsty,psym=3,xtit='NH',ytit='Y',title='Grain 3'
  1499 +
  1500 +window,win,xsize=800,ysize=900 & win=win+1
  1501 +!p.multi=[0,1,4]
  1502 +cgplot,chi2s,F0s,/xlog,/ylog,xr=chi2_range,yr=F0_range,/xsty,psym=3,xtit='chi2',ytit='F0'
  1503 +cgplot,chi2s,Ygrain1,/xlog,/ylog,xr=chi2_range,yr=Ypah_range,/xsty,psym=3,xtit='chi2',ytit='Y',title='Grain 1'
  1504 +cgplot,chi2s,Ygrain2,/xlog,/ylog,xr=chi2_range,yr=Ypah_range,/xsty,psym=3,xtit='chi2',ytit='Y',title='Grain 2'
  1505 +cgplot,chi2s,Ygrain3,/xlog,/ylog,xr=chi2_range,yr=Yvsg_range,/xsty,psym=3,xtit='chi2',ytit='Y',title='Grain 3'
  1506 +
  1507 +goto,the_end
  1508 +stop
  1509 +
  1510 +jwst_4_coutours=fltarr((size(jwst_images))[1],(size(jwst_images))[2],2)
  1511 +jwst_4_coutours[*,*,0]=jwst_images[*,*,0,5]
  1512 +jwst_4_coutours[*,*,1]=jwst_images[*,*,0,5]
  1513 +levs=create_contour_st(jwst_4_coutours,lev1=[1,5])
  1514 +levs.(0).color=100
  1515 +levs.(1).color=200
  1516 +levs.(1).thick=2
  1517 +
  1518 +obp=[1.1,0.,1.15,1]
  1519 +image_cont20,jwst_4_coutours[*,*,0],href,/square,imrange=[-1.,10],image_color_table='jpbloadct',/silent,title='',off_bar=obp,axis_color_table=1,levels=levs
  1520 +
  1521 +window,0,xsize=800,ysize=1000 & win=win+1
  1522 +aa=alog10(F0_map)
  1523 +wset,0 & image_cont20,F0_map,href,/square,imrange=[-2.,100],image_color_table='jpbloadct',/silent,title='F0',off_bar=obp,axis_color_table=1
  1524 +wset,0 & image_cont20,aa,href,/square,imrange=[0.5,2.5],image_color_table='jpbloadct',/silent,title='F0',off_bar=obp,axis_color_table=1,levels=levs
  1525 +window,1,xsize=800,ysize=1000
  1526 +wset,1 & image_cont20,Ypah_map,href,/square,imrange=[-0.005,0.5],image_color_table='jpbloadct',/silent,title='Ypah',off_bar=obp,axis_color_table=1,levels=levs
  1527 +window,2,xsize=800,ysize=1000
  1528 +wset,2 & image_cont20,Yvsg_map,href,/square,imrange=[-0.005,0.5],image_color_table='jpbloadct',/silent,title='Yvsg',off_bar=obp,axis_color_table=1,levels=levs
  1529 +window,3,xsize=800,ysize=1000
  1530 +wset,3 & image_cont20,alog10(chi2_map),href,/square,imrange=[2,6],image_color_table='jpbloadct',/silent,title='chi2',off_bar=obp,axis_color_table=1,levels=levs
  1531 +window,4,xsize=800,ysize=1000
  1532 +image_cont20,alog10(fact_map),href,/square,imrange=[-1,5],image_color_table='jpbloadct',/silent,title='fact',off_bar=obp,axis_color_table=1,levels=levs
  1533 +
  1534 +the_end:
  1535 +
  1536 +END
... ...
LabTools/IRAP/RM/test_size_distribution.pro 0 → 100644
  Diff comments   View file @9a95b8f
... ... @@ -0,0 +1,32 @@
  1 +PRO test_size_distribution
  2 +
  3 +fpd=['(*!dustem_params).gas.g0','(*!dustem_params).g0','dustem_plugin_phangs_class_isrf_1','dustem_plugin_phangs_class_isrf_4'] ;ISRF class to be used
  4 +model='DL07'
  5 +dustem_init,model=model
  6 +fortran_user=dustem_set_up_fortran(/random_name) ;use a random fortran number
  7 +!dustem_verbose=0
  8 +(*!dustem_params).KEYWORDS='quiet '+(*!dustem_params).KEYWORDS ;This makes Fortran be quiet too
  9 +;stop
  10 +(*!dustem_params).grains[0].TYPE_KEYWORDS='logn'
  11 +(*!dustem_params).grains[1].TYPE_KEYWORDS='logn'
  12 +grain_keywords=[(*!dustem_params).grains.TYPE_KEYWORDS] ;This is to inform dustem_make_sed_table, which actually runs the models
  13 +
  14 +
  15 +params = ['(*!dustem_params).grains(0).amax','(*!dustem_params).grains(1).amax']
  16 +fiv=[-1., 1,0.,0] ;This sets the ISRF class to the requested value. also removes ionising photons.
  17 +; use_fit_filters_names=[dustem_instru2filters('NIRCAM'),dustem_instru2filters('MIRI'),dustem_instru2filters('PACS'),dustem_instru2filters('SPIRE')]
  18 +
  19 +use_fit_filters_names = dustem_filter_names2filters(['F300M','F335M','F360M','F0770W','F1000W','F1130W','F2100W','PACS_BLUE','PACS_GREEN','PACS_RED','PSW','PMW','PLW'])
  20 +path_save = 'tmp/sed_size.fits'
  21 +
  22 +
  23 +dustem_make_sed_table,model,params, [1.2000E-07,1.2000E-07],[3.300E-08,3.300E-08],[5,5],fpd=fpd,fiv=fiv,plog=[0,0],/show_seds,filters=use_fit_filters_names,filename=path_save, $
  24 + grain_keywords=grain_keywords
  25 +
  26 +
  27 +
  28 +
  29 +the_end:
  30 +message,'test_size_distribution',/info$
  31 +
  32 +END
0 33 \ No newline at end of file
... ...
LabTools/IRAP/RM/test_vsg.pro 0 → 100644
  Diff comments   View file @9a95b8f
... ... @@ -0,0 +1,32 @@
  1 +PRO test_vsg
  2 +
  3 +fpd=['(*!dustem_params).gas.g0','(*!dustem_params).g0','dustem_plugin_phangs_class_isrf_1','dustem_plugin_phangs_class_isrf_4'] ;ISRF class to be used
  4 +model='DL07'
  5 +dustem_init,model=model
  6 +fortran_user=dustem_set_up_fortran(/random_name) ;use a random fortran number
  7 +!dustem_verbose=0
  8 +(*!dustem_params).KEYWORDS='quiet '+(*!dustem_params).KEYWORDS ;This makes Fortran be quiet too
  9 +;stop
  10 +(*!dustem_params).grains[0].TYPE_KEYWORDS='logn'
  11 +(*!dustem_params).grains[1].TYPE_KEYWORDS='logn'
  12 +grain_keywords=[(*!dustem_params).grains.TYPE_KEYWORDS] ;This is to inform dustem_make_sed_table, which actually runs the models
  13 +
  14 +
  15 +params = ['(*!dustem_params).grains[4].MDUST_O_MH']
  16 +fiv=[-1., 1,0.,0] ;This sets the ISRF class to the requested value. also removes ionising photons.
  17 +; use_fit_filters_names=[dustem_instru2filters('NIRCAM'),dustem_instru2filters('MIRI'),dustem_instru2filters('PACS'),dustem_instru2filters('SPIRE')]
  18 +
  19 +use_fit_filters_names = dustem_filter_names2filters(['F300M','F335M','F360M','F0770W','F1000W','F1130W','F2100W','PACS_BLUE']);,'PACS_GREEN','PACS_RED','PSW','PMW','PLW'])
  20 +path_save = 'tmp/sed_vsg.fits'
  21 +
  22 +
  23 +dustem_make_sed_table,model,params, [1.6600E-04],[1.0E-3],[10],fpd=fpd,fiv=fiv,plog=[0],/show_seds,filters=use_fit_filters_names,filename=path_save, $
  24 + grain_keywords=grain_keywords
  25 +
  26 +
  27 +
  28 +
  29 +the_end:
  30 +message,'test_size_distribution',/info$
  31 +
  32 +END
0 33 \ No newline at end of file
... ...