check_phangs_ssps_against_uv_saved.pro
7.54 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
181
182
183
184
185
186
187
188
189
190
191
192
PRO check_phangs_ssps_against_uv,show_map=show_map
;check_phangs_ssps_against_uv
;===== This computes a map of prediction by dustem_plugin_phangs_stellar_continuum
;===== constrained on JWST NIR data to astrosat data in the UV for NGC0628.
;===== uses the plugin in a mode where the fortran is not run, but yet, the output files of the fortran are read at each pixel,
;===== which is not efficient (supposed to run in 46hrs .... See check_phangs_ssps_against_uv_fast for an alternative use.
window,0,xsize=900,ysize=1000
pdp_define_la_common
obp=[1.1,0,1.15,1]
source_name='ngc0628'
data_dir='/Volumes/PILOT_FLIGHT1/PHANGS-JWST/DR1/'
;dustem_init,show_plots=show_plots
;needed only for NHCO
restore,data_dir+'ngc0628_jwst_images.sav',/verb
;needed for stellar parameters and voronoi bin info
restore,data_dir+'ngc0628_muse_images.sav',/verb
;restore,data_dir+'ngc0628_muse_NH.sav',/verb
;restore,data_dir+'ngc0628_astrosat_data.sav',/verb
;stop
;=== restore SEDs in voronoi bins
restore,data_dir+'ngc0628_jwst_seds_muse_pixels.sav',/verb
jwst_seds=all_seds
jwst_filters=(*jwst_seds[0]).filter
restore,data_dir+'ngc0628_muse_seds_muse_pixels.sav',/verb
muse_seds=all_seds
muse_filters=(*muse_seds[0]).filter
restore,data_dir+'ngc0628_astrosat_seds_muse_pixels.sav',/verb
astrosat_seds=all_seds
astrosat_filters=(*astrosat_seds[0]).filter
all_filters=[astrosat_filters,muse_filters,jwst_filters]
Nvor=max(voronoi_id)
Nfilters=n_elements(jwst_filters)+n_elements(astrosat_filters)+n_elements(muse_filters)
use_NHmap=NHCO ;used NH map in 1.e21 (from CO)
;galaxy_distance=10. ;MPc
use_model='DBP90' ;Example with default keywords uses the DBP90 model
;use_model='MC10' ;Example with default keywords uses the DBP90 model
;use_model='DL07' ;Example with default keywords uses the DBP90 model
use_polarization=0 ; initialize Dustemwrap in no polarization mode
use_window=2 ; default graphics window number to use for plotting the results
use_verbose=0
use_Nitermax=5 ; maximum number of iterations for the fit
dustem_define_la_common
;parameters to fit
pd = [ $
'(*!dustem_params).G0', $ ;G0
'(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
'(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
;'(*!dustem_params).grains(2).mdust_o_mh', $ ;VSGs
;'(*!dustem_params).grains(3).mdust_o_mh', $ ;VSGs
;'(*!dustem_params).grains(4).mdust_o_mh', $ ;BGs
'dustem_plugin_phangs_stellar_continuum_1'] ;stellar continuum amplitude
;iv = [1.6, [2.2e-3, 2.2e-3, 5.7e-3, 5.7e-3],2.e-3]
iv = [1.6, [2.2e-3, 2.2e-3],5.e-5]
;iv = [1.6, [2.2e-3, 2.2e-3,2e-3, 2e-3],5.e-5]
Npar=n_elements(pd)
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar)
;llims=[1.e-2,[1.e-4,1.e-4,1.e-4,1.e-4],1.e-10]
llims=[1.e-2,[1.e-4,1.e-4],1.e-10]
;llims=[1.e-2,[1.e-4,1.e-4,1.e-4,1.e-4],1.e-10]
pd_to_update=['dustem_plugin_phangs_stellar_continuum_1','(*!dustem_params).grains(0).mdust_o_mh']
pd_filter_names=['NIRCAM11','MIRI2']
;== INITIALISE DUSTEM
dustem_init,model=use_model,polarization=use_polarization,show_plots=show_plots
;stop
!quiet=1
!dustem_nocatch=1
!dustem_verbose=use_verbose
IF keyword_set(noobj) THEN !dustem_noobj=1
;!EXCEPT=2 ; for debugging (not sure hat this is)
;=== INFORMATION TO MAKE THE PLOT
;yr=[1.00e-2,10.0] ; y-axis limits
yr=[1.00e-3,100.0] ; y-axis limits
;xr=[0.5,1.00e2] ; x-axis limits
xr=[0.1,1.00e2] ; x-axis limits
tit='FIT INTENSITY EXAMPLE' ; plot title
ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') ; y-axis title
xtit=textoidl('\lambda (\mum)') ; x-axis title
;voronoi bins for test. Last series are randomly picked
one_sed=(*jwst_seds[0])[0]
Nfilters=n_elements(astrosat_filters)+n_elements(muse_filters)+n_elements(filters)
seds=replicate(one_sed,Nfilters)
diff_seds=fltarr(Nvor,Nfilters)
diff_maps=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),Nfilters) ;difference maps
diff_maps[*]=la_undef()
;imrange=[-1.,1.] ;MJy/sr
imrange=[-0.3,0.2] ;MJy/sr
t1=systime(0,/sec)
first_vid=0L
first_vid=5960L
FOR vid=first_vid,Nvor-1 DO BEGIN
IF vid mod 10 EQ 0 THEN BEGIN
t2=systime(0,/sec)
delta_t_hr=(t2-t1)/60.^2 ;elapsed time [hr]
frac_perc=vid/Nvor*100
delta_t_remain=delta_t_hr/frac_perc*100.-delta_t_hr ;remaining time [hr]
message,'done '+strtrim(frac_perc)+' % in '+strtrim(delta_t_hr,2)+' hr remaining '+strtrim(delta_t_remain,2)+' hr',/continue
ind=where(diff_maps[*,*,0] NE la_undef(),count)
print,la_min(diff_maps[*,*,0]),la_max(diff_maps[*,*,0]),count
image_cont20,diff_maps[*,*,0],href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,off_bar_pos=obp,imrange=imrange
;stop
ENDIF
;index=vid
;===== get the SSP weights for the given voronoi bins
weights=phangs_binid2weights(st_muse_weights,vid,st_templates,age_values,metalicity_values,reddening=reddening)
;===== extract SED
Nf=0L & Nt=Nf+n_elements(astrosat_filters)-1
seds[Nf:Nt]=*astrosat_seds[vid]
Nf=Nt+1 & Nt=Nf+n_elements(muse_filters)-1
seds[Nf:Nt]=*muse_seds[vid]
Nf=Nt+1 & Nt=Nf+n_elements(jwst_filters)-1
seds[Nf:Nt]=*jwst_seds[vid]
sed=seds
;needed to add this because astrosat had 0 variances
sed.sigmaii=(5./100.*sed.stokesI)^2
;===== set the fixed parameters and initial values
fpd=['dustem_plugin_phangs_stellar_continuum_2'] ;This is E(B-V)
fiv=[reddening] ;This is E(B-V) value
Nage=n_elements(age_values)
NZ=n_elements(metalicity_values)
;NH_value=la_mul(la_mean(use_NHmap[*all_seds_indices[vid]]),10.) ;NH value in 1.e20 H/cm2
;sed.stokesI=sed.stokesI/NH_value ;normalize SED to 1.e20 H/cm2
;sed.sigmaII=sed.sigmaII/NH_value^2 ;normalize SED to 1.e20 H/cm2
;
FOR i=0L,Nage-1 DO BEGIN
FOR j=0L,NZ-1 DO BEGIN
IF weights[i,j] NE 0 AND bins[i,j] NE la_undef() THEN BEGIN
ij=lonarr(1,2) & ij[0,0]=i & ij[0,1]=j
ind=ij2index(ij,[Nage,Nz])+1
fpd=[fpd,'dustem_plugin_phangs_stellar_continuum_'+strtrim(ind[0]+2,2)] ;+2 is because amplitude,E(B-V) is parameter 0 of pluggin and bins start at zero
fiv=[fiv,weights[i,j]]
ENDIF
ENDFOR
ENDFOR
;all_params=[pd,fpd]
;all_values=[iv,fiv]
;try=dustem_parameter_description2type(all_params[0],string_name=string_name)
;stop
;key=lindgen(n_elements(fpd))
;val=fiv
;toto=dustem_plugin_phangs_stellar_continuum(key=key,val=val)
;== SET THE OBSERVATIONAL STRUCTURE
dustem_set_data,m_fit=sed,m_show=sed
stop
IF vid EQ first_vid THEN BEGIN
input_dustem_st=0
use_previous_fortran = 0
new_iv=dustem_first_guess(pd,iv,sed,pd_to_update,pd_filter_names,fpd=fpd,fiv=fiv,pol=pol, $
use_previous_fortran=use_previous_fortran,output_dustem_st=dustem_st,input_dustem_st=input_dustem_st,dustem_predicted_sed=dustem_predicted_sed)
;pb after this !dustem_plugin)[i].spec not set
ENDIF ELSE BEGIN
use_previous_fortran = 1
new_iv=dustem_first_guess(pd,iv,sed,pd_to_update,pd_filter_names,fpd=fpd,fiv=fiv,pol=pol, $
use_previous_fortran=use_previous_fortran,input_dustem_st=dustem_st,dustem_predicted_sed=dustem_predicted_sed)
ENDELSE
;stop
;=== get a first guess at parameters (must be done after dustem_set_data)
;== compute predicted sed
;dustem_predicted_sed=dustem_compute_sed(new_iv,st=st)
;== compute difference between predicted and observed SED
diff_sed=la_sub(sed.stokesI,dustem_predicted_sed)
diff_seds[vid,*]=diff_sed
;=== update difference maps
FOR j=0L,Nfilters-1 DO BEGIN
map=diff_maps[*,*,j]
map[*all_seds_indices[vid]]=diff_seds[vid,j]
diff_maps[*,*,j]=map
ENDFOR
ENDFOR
message,'========== Finished computing difference SEDs',/continue
stop
END