check_phangs_ssps_isrf_prediction.pro
6.13 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
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,bidon,/save
make_phangs_isrf_classes,bidon
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