make_phangs_isrf_classes.pro
5.38 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
PRO make_phangs_isrf_classes,bidon,save=save,help=help
;+
; NAME:
; make_phangs_isrf_classes
; CALLING SEQUENCE:
; make_phangs_isrf_classes[,/save]
; PURPOSE:
; generates the ISRF classes used in fiting PHANGs data
; INPUTS:
; None
; OPTIONAL KEYWORDS:
; save = if set, save the classes
; help = if set, print this help
; OUTPUTS:
; None
; OPTIONAL INPUT:
; None
; OPTIONAL OUTPUT:
; None
; PROCEDURE AND SUBROUTINE USED
; None
; SIDE EFFECTS:
; None
; EXAMPLE:
; make_phangs_isrf_classes,/save
; MODIFICATION HISTORY:
; written by Jean-Philippe Bernard
;-
IF keyword_set(help) THEN BEGIN
doc_library,'make_phangs_isrf_classes'
goto,the_end
ENDIF
win=0L
;data_dir='/Volumes/PILOT_FLIGHT1/PHANGS-JWST/DR1/'
data_dir=!phangs_data_dir+'/ISRF/WORK/'
use_model='DBP90' ;Example with default keywords uses the DBP90 model
use_polarization=0 ; initialize Dustemwrap in no polarization mode
;== INITIALISE DUSTEM
dustem_init,model=use_model,polarization=use_polarization,show_plots=show_plots
;==== restore needed data
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.
;stop
;=== attempt to classify ISRFs
;normalized ISRFs
lambir=dustem_get_wavelengths(isrf_wavelengths=isrf_wavelengths)
Nlambir=(size(ISRFS))[1]
Nvor=(size(ISRFS))[2]
n_ISRFS=fltarr(Nlambir,Nvor)
n_wavs=fltarr(Nlambir,Nvor)
ind=where(isrf_wavelengths GT 1.)
indd=ind[0]
FOR vid=0L,Nvor-1 DO BEGIN
n_ISRFS[*,vid]=ISRFS[*,vid]/ISRFS[indd,vid]
n_wavs[*,vid]=isrf_wavelengths
ENDFOR
;=== 2D histogram of normalized ISRFs in NGC0628
x=reform(n_wavs,Nlambir*Nvor)
y=reform(n_ISRFS,Nlambir*Nvor)
;jpbloadct
loadct,13
Nx=200 & Ny=200
Nx=100 & Ny=100
xmin=-1.5 & xmax=2.
ymin=-5 & ymax=1.
xtit='log10(wav [mic])'
ytit='log10(ISRF/ISRF(1 mic))'
range=[-1,5]
window,win & win=win+1
plot_bin_correl,alog10(x),alog10(y),Nx,Ny,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,xtit=xtit,ytit=ytit,range=range;,/xlog ;,/ylog
;compute the median ISRFS
med_isrf=la_median(n_ISRFS,dim=-1)
;refs_wavs=[0.11] ;This is set to be just before the Lymann break
reference_wavelength=1. ;all ISRFs are normalized to 1 mic
refs_wavs=[0.15] ;This is set to be just before the Lymann break
ref_ratios=interpol(med_isrf,isrf_wavelengths,refs_wavs)
;These are the ratios which will be used to classify ISRFs
one_class={number:0L $;class number
,name:'' $ ;class name
,ISRF_template:ptr_new() $ ;isrf template
,flux_ratio:0. $ ;flux ratio
,flux_ratio_wav:0. $ ;flux ratio wavelength
,reference_wavelength:reference_wavelength $ ;ISRF normalisation wavelength
,flux_ratio_min:0. $ ;minimum flux ratio to belong to this class
,flux_ratio_max:0. $ ;maximum flux ratio to belong to this class
,flux_ref:0. $ ;
,Nvor:0L $ ;number of Voronoi bins in the class
,voronoi_members:ptr_new() $ ;voronoi pixel members of that class
,voronoi_chi2:ptr_new()} ;chi2 between ref ISRF and voronoi member's ISRF
;class_ratio_range=[8.e-5,3.]
;class_ratio_range=[5.e-5,3.]
;Nclasses=20L
class_ratio_range=[1.e-5,10.]
Nclasses=30L
;Nclasses=200L ;for test
ratios=range_gen(Nclasses,class_ratio_range,/log,/get_minmax,min_values=ratios_min,max_values=ratios_max)
xratios=indgen(Nclasses)
classes=replicate(one_class,Nclasses)
;fill in class values
classes.number=lindgen(Nclasses)
classes.flux_ratio=ratios
classes.flux_ratio_wav=refs_wavs[0]
classes.flux_ratio_min=ratios_min
classes.flux_ratio_max=ratios_max
FOR i=0L,Nclasses-1 DO BEGIN
classes[i].ISRF_template=ptr_new(n_ISRFs[*,i])
ENDFOR
;=== Do a class (class 0) with a Mathis field
class0=one_class
class0.number=0L
class0.name='Mathis'
class0.reference_wavelength=1.
class0.flux_ratio_wav=refs_wavs[0]
file=!dustem_soft_dir+'/data/ISRF_MATHIS.DAT'
;file='/Users/jpb/Soft_Libraries/dustem_fortran/data/ISRF.DAT'
readcol,file,Mathis_wavs,Mathis_ISRF
Mathis_ISRF=interpol(Mathis_ISRF,Mathis_wavs,isrf_wavelengths)
n_Mathis_ISRF=Mathis_ISRF/interpol(Mathis_ISRF,isrf_wavelengths,class0.reference_wavelength)
class0.ISRF_template=ptr_new(n_Mathis_ISRF)
ratio=interpol(n_Mathis_ISRF,isrf_wavelengths,class0.flux_ratio_wav)
class0.flux_ratio=ratio
class0.flux_ratio_min=ratio*(1.-0.1) ;arbitrary +-10% around ratio
class0.flux_ratio_max=ratio*(1.+0.1) ;arbitrary +-10% around ratio
Mathis_1mic=interpol(Mathis_ISRF,isrf_wavelengths,1.0)
message,'Mathis (G0=1) has 4*pi*Inu='+strtrim(Mathis_1mic,2)+' ergs/s/cm2/Hz',/continue
message,'Mathis (G0=1) has Inu(0.15 mic)/Inu(1 mic)='+strtrim(ratio,2)+' ',/continue
;stop
;save empty classes defined above
IF keyword_set(save) THEN BEGIN
file_save=data_dir+'isrf_classes_one_ratio.sav'
message,'Will save '+file_save,/continue
message,'.c to proceed ',/continue
stop
save,classes,class0,file=file_save,/verb
message,'Saved '+file_save,/continue
ENDIF
the_end:
END