PRO check_phangs_ssps_isrf_prediction,show_map=show_map,from_restore=from_restore,from_classes_restore=from_classes_restore,save=save ;check_phangs_ssps_isrf_prediction ;check_phangs_ssps_isrf_prediction,/from_restore ;check_phangs_ssps_isrf_prediction,/from_classes_restore,/save ;must have run make_phangs_ssps_isrf_prediction.pro before from_restore: ;make_phangs_isrf_classes,/save make_phangs_isrf_classes stop from_classes_restore: data_dir=!phangs_data_dir+'/ISRF/WORK/' ;file=data_dir+'isrf_classes_one_ratio.sav' restore,data_dir+'isrf_classes_one_ratio.sav',/verbose restore,data_dir+'ngc0628_isrf_min_prediction.sav',/verb ;% RESTORE: Restored variable: ISRFS. ;% RESTORE: Restored variable: OBJECT_DISTANCE. ;% RESTORE: Restored variable: OBJECT_THICKNESS. ;% RESTORE: Restored variable: SOURCE_NAME. restore,data_dir+'ngc0628_astrosat_voronoi_prediction_fast.sav',/verb ;% RESTORE: Restored variable: ACTUAL_SEDS. ;% RESTORE: Restored variable: PREDICTED_SEDS. ;% RESTORE: Restored variable: HREF. ;% RESTORE: Restored variable: ALL_SEDS_INDICES. ;% RESTORE: Restored variable: DO_NORMALIZE. ;% RESTORE: Restored variable: ALL_FILTERS. ;% RESTORE: Restored variable: NHVOR. Nclasses=n_elements(classes) refs_wavs=[classes[0].FLUX_RATIO_WAV] reference_wavelength=classes[0].reference_wavelength ratios=classes.flux_ratio lambir=dustem_get_wavelengths() Nlambir=(size(ISRFS))[1] Nvor=(size(ISRFS))[2] n_ISRFS=fltarr(Nlambir,Nvor) n_wavs=fltarr(Nlambir,Nvor) ind=where(lambir GT 1.) indd=ind[0] FOR vid=0L,Nvor-1 DO BEGIN n_ISRFS[*,vid]=ISRFS[*,vid]/ISRFS[indd,vid] n_wavs[*,vid]=lambir ENDFOR ;compute flux ratio for each voronoi bin vor_ratios=fltarr(Nvor) FOR i=0L,Nvor-1 DO BEGIN vor_ratios[i]=interpol(n_ISRFs[*,i],lambir,refs_wavs[0]) ENDFOR ;fill in template ISRF of each class ;method='one_spectrum' method='average_spectrum' empty_classes=intarr(Nclasses) CASE method OF 'one_spectrum':BEGIN ;This is selection based on a single ISRF FOR i=0L,Nclasses-1 DO BEGIN ind=where(vor_ratios GT ratios[i],count) classes[i].ISRF_template=ptr_new(n_ISRFS[*,ind[0]]) ENDFOR END 'average_spectrum':BEGIN ;This is based on an average delta=1.0 ;This is a bit arbitrary. Chosen so that there is no empty classes FOR i=0L,Nclasses-1 DO BEGIN ;ind=where(vor_ratios GE ratios[i]-ratios[i]*delta and vor_ratios LE ratios[i]+ratios[i]*delta,count) ratio_min=classes[i].flux_ratio_min ratio_max=classes[i].flux_ratio_max ;ind=where(vor_ratios GT ratios_min[i] and vor_ratios LE ratios_max[i],count) ind=where(vor_ratios GT ratio_min and vor_ratios LE ratio_max,count) IF count NE 0 THEN BEGIN classes[i].ISRF_template=ptr_new(la_mean(n_ISRFS[*,ind],dim=-1)) ENDIF ELSE BEGIN message,'Found no voronoi bin in class '+strtrim(i+1,2),/continue ;stop empty_classes[i]=1 ENDELSE ENDFOR END ENDCASE ;===== remove empty classes ind=where(empty_classes EQ 0,count) classes=classes[ind] Nclasses=count ;ratios_min=ratios_min[ind] ;ratios_max=ratios_max[ind] ;===== Fill in class number and name FOR i=0L,Nclasses-1 DO BEGIN classes[i].number=i+1 ;because class 0 will be for Mathis ISRF str='I('+strtrim(refs_wavs[0],2)+'/I('+strtrim(reference_wavelength,2)+')='+strtrim(ratios[i],2) classes[i].name=str ENDFOR stop ;classify Voronoi bin according to flux ratio ;This could also be done on chi2, instead of just ratio ... vor_class=lonarr(Nvor) FOR i=0L,Nvor-1 DO BEGIN ;ratio=interpol(n_ISRFs[*,i],lambir,refs_wavs[0]) this_ratio=vor_ratios[i] ind=where(classes.flux_ratio_min LT this_ratio and classes.flux_ratio_max GE this_ratio,count) ;ind=where(ratios_min LT this_ratio and ratios_max GE this_ratio,count) xv=ind[0] ;xv=round(interpol(xratios,ratios,this_ratio)) IF ptr_valid(classes[xv].voronoi_members) THEN BEGIN classes[xv].voronoi_members=ptr_new([*classes[xv].voronoi_members,i]) ENDIF ELSE BEGIN classes[xv].voronoi_members=ptr_new([i]) ENDELSE classes[xv].Nvor=classes[xv].Nvor+1 this_chi2=total((n_ISRFS[*,i])^2-(*classes[xv].ISRF_template)^2) IF ptr_valid(classes[xv].voronoi_chi2) THEN BEGIN classes[xv].voronoi_chi2=ptr_new([*classes[xv].voronoi_chi2,this_chi2]) ENDIF ELSE BEGIN classes[xv].voronoi_chi2=ptr_new([this_chi2]) ENDELSE ;This is a voronoi vector of the belonging class vor_class[i]=xv ENDFOR ;classes[xv].voronoi_chi2=ptr_new((*classes[xv].voronoi_chi2)[1:*]) ;classes[xv].voronoi_members=ptr_new((*classes[xv].voronoi_members)[1:*]) ;==== plot histogram of voronoi classes res=histogram(vor_class,min=0,max=Nclasses-1,Nbin=Nclasses,locations=x) window,win & win=win+1 !p.multi=[0,1,2] cgplot,x,res,psym=10,xrange=[-1,Nclasses+1],/xstyle,xtit='Class',ytit='Npix' cgplot,ratios,res,psym=10,xrange=[1.e-5,10],/xstyle,xtit='Ratio',ytit='Npix',/xlog ;==== make a sky map of the classes window,win & win=win+1 map_classes=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef() FOR i=0L,Nvor-1 DO BEGIN map_classes[*all_seds_indices[i]]=vor_class[i] ENDFOR obp=[1.1,0.,1.15,1] image_cont20,map_classes,href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,imrange=[-1,Nclasses+1],off_bar=obp,tit='Class map' ;=== plot ISRF templates of all classes window,win & win=win+1 xtit='Wavelength [mic]' ytit='Normalized ISRF' tit='ISRF classes' loadct,13 !p.multi=0 med_isrf=la_median(n_ISRFS,dim=-1) cgplot,lambir,med_isrf,color='black',/xlog,/ylog,xrange=[1.e-2,1.e2],/xsty,yrange=[1.e-5,10.],/nodata,xtit=xtit,ytit=ytit,tit=tit colors=fix(range_gen(Nclasses,[10,250])) FOR i=0L,Nclasses-1 DO cgoplot,lambir,*classes[i].ISRF_template,color=colors[i] cgoplot,lambir,med_isrf,color='black',thick=2. cgoplot,lambir,*class0.isrf_template,color='black',thick=2. cgoplot,replicate(refs_wavs,2),10^!y.crange,linestyle=2 cgoplot,replicate(reference_wavelength,2),10^!y.crange,linestyle=2 ;stop IF keyword_set(save) THEN BEGIN file_save=data_dir+'ngc0682_isrf_classes_one_ratio.sav' message,'About to write '+file_save,/continue message,'Make sure you are allowed to do that and type .c to proceed.',/continue save,classes,vor_class,map_classes,file=file_save,/verb message,'Saved '+file_save,/continue ENDIF message,'========== Finished computing ISRF classes',/continue stop END