dustem_plugin_phangs_stellar_continuum.pro
9.47 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
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
FUNCTION dustem_plugin_phangs_stellar_continuum,key=key $
,val=val $
,scope=scope $
,paramtag=paramtag $
,age_values=age_values $
,metalicity_values=metalicity_values $
,paramdefault=paramdefault $
; ,object_distance=object_distance $ ;not used
,Voronoi_Npix=Voronoi_Npix $
,help=help
;+
; NAME:
; dustem_plugin_phangs_stellar_continuum
; PURPOSE:
; DustEMWrap plugin to add EMILES stellar template emission to Dustemwrap SEDs
; CATEGORY:
; DustEM, Distributed, Mid-Level, Plugin
; CALLING SEQUENCE:
; stellar_continuum=dustem_plugin_phangs_stellar_continuum([,key=][,val=][,scope=][,paramtag=][paramdefault=][,/help])
; INPUTS:
; None
; OPTIONAL INPUT PARAMETERS:
; key = input parameter number
; First parameter: Amplitude factor for the emission
; Second parameter: E(B-V) for extinction to stars applied to the template, as used in Phangs
; Following 13*6 parameters are weights by which to multiply the MILES templates [Msun/pc^2]
; val = corresponding input parameter value
; OUTPUTS:
; stellar_continuum = stellar continuum spectrum resampled on dustem wavelengths [MJy/sr]
; OPTIONAL OUTPUT PARAMETERS:
; scope = if set, return the scope of the plugin
; paramdefault = default values of parameters
; paramtag = if set, return the return plugin parameter names as strings
; ACCEPTED KEY-WORDS:
; help = if set, print this help
; method = SSP used. can be EMILES or CB19. default='EMILES'
; COMMON BLOCKS:
; None
; SIDE EFFECTS:
; None
; RESTRICTIONS:
; The DustEMWrap IDL code must be installed
; PROCEDURE:
; This is a DustEMWrap plugin for phangs
; It differs from dustem_plugin_emiles_stellar_continuum.pro only by the way the extinction is applied to the output spectrum
; EXAMPLES
; dustem_init
; vec=dustem_plugin_phangs_stellar_continuum(scope=scope)
; MODIFICATION HISTORY:
; Written by JPB June 2023
; Evolution details on the DustEMWrap gitlab.
; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
;-
IF keyword_set(help) THEN BEGIN
doc_library,'dustem_plugin_phangs_stellar_continuum'
output=0.
goto,the_end
ENDIF
IF keyword_set(scope) THEN BEGIN
scope='ADD_SED'
out=0
goto, the_scope
ENDIF
;use_method='EMILES'
use_method='CB19'
age_values=[0.03, 0.05, 0.08, 0.15, 0.25, 0.40, 0.60, 1.0, 1.75, 3.0, 5.0, 8.5, 13.5]
metalicity_values=[-1.48630, -0.961400, -0.351200, +0.0600000, +0.255900, +0.397100]
;==== This does not work because the info files has many more entries than there are available ssps.
;age_values=age_values_emiles
;metalicity_values=met_values_emiles
;This does not work either
;st_templates=read_muse_templates_info(age_values=age_values,metalicity_values=metalicity_values,Nbins=Nbins,Nage=Nage,NZ=Nmetalicity,bins=bins)
;print,age_values
; 0.030000000 0.050000000 0.080000000 0.15000000 0.25000000 0.40000000 0.60000000 1.0000000 1.7500000 3.0000000 5.0000000 8.5000000 13.500000
;print,metalicity_values
; -1.49000 -0.960000 -0.350000 0.0600000 0.260000 0.400000
;stop
Nage=n_elements(age_values)
Nmetalicity=n_elements(metalicity_values)
Nparam=Nage*Nmetalicity+2 ;+2 is for Amplitude and E(B-V)
;==== define parameter tags
IF keyword_set(paramtag) THEN BEGIN
paramtag=strarr(Nparam)
paramtag[0]='Amplitude'
paramtag[1]='E(B-V)'
FOR i=2L,Nparam-1 DO BEGIN
ij=index2ij([i-2],[Nage,Nmetalicity])
paramtag[i]='Age'+strtrim(ij[0,0],2)+'Met'+strtrim(ij[0,1],2)
ENDFOR
GOTO,the_paramtag
ENDIF
paramdefault=dblarr(Nparam)
paramdefault[*]=0.D0 ;default parameter values is 0.
;==== decode key keyword.
;stop
paramvalues=paramdefault
IF keyword_set(key) THEN BEGIN
ind=where(key EQ 1,count)
IF count NE 0 THEN paramvalues[0]=val[ind[0]]
FOR i=1L,Nparam-1 DO BEGIN
ind=where(key EQ i+1,count)
IF count NE 0 THEN paramvalues[i]=val[ind[0]]
ENDFOR
ENDIF
;==== initialize templates if not already present
defsysv,'!dustem_plugin_phangs_stellar_continuum',exist=exist
IF exist EQ 0 THEN BEGIN
st={ssps:ptr_new(),wavs:ptr_new(),Mstars:ptr_new(),mstars_o_Lv:ptr_new(),Mstarprems:ptr_new()}
defsysv,'!dustem_plugin_phangs_stellar_continuum',st
ENDIF
IF not ptr_valid(!dustem_plugin_phangs_stellar_continuum.ssps) THEN BEGIN
message,'Initializing SSP templates',/info
CASE use_method of
'EMILES': BEGIN
template_flux=dustem_read_emiles_stellar_templates(age_values,metalicity_values,template_wav=template_wav,Mstars=Mstars,Ms_o_Lv=Ms_o_Lv,Mstarprems=Mstarprems)
END
'CB19': BEGIN
;== this is just to get the Mstars
bidon=dustem_read_emiles_stellar_templates(age_values,metalicity_values,template_wav=template_wav,Mstars=Mstars,Ms_o_Lv=Ms_o_Lv,Mstarprems=Mstarprems)
;== This is to read in the CB19 templates
;== Note, templates are already interpolated on dustem wavelenghts (!)
template_flux=dustem_read_cb19_stellar_templates(age_values=age_values,metalicity_values=metalicity_values,template_wav=template_wav);,Mstars=Mstars)
END
ENDCASE
!dustem_plugin_phangs_stellar_continuum.ssps=ptr_new(template_flux) ; in Lsun Msun^-1 A^-1
!dustem_plugin_phangs_stellar_continuum.wavs=ptr_new(template_wav)
!dustem_plugin_phangs_stellar_continuum.mstars=ptr_new(Mstars)
!dustem_plugin_phangs_stellar_continuum.mstars_o_Lv=ptr_new(Ms_o_Lv)
!dustem_plugin_phangs_stellar_continuum.mstarprems=ptr_new(Mstarprems)
ENDIF
;stop
use_Voronoi_Npix=1. ;deg^2
IF keyword_set(Voronoi_Npix) THEN use_Voronoi_Npix=Voronoi_Npix
;use_object_distance=1. ;MPc
;IF keyword_set(object_distance) THEN use_object_distance=object_distance
;bin_surface_pc2=use_bin_surface*!pi^2/180.^2*(use_object_distance*1.e6)^2 ;surface in pc^2
;Npixfact=(1./(0.2/60./60.))^2 ;how much Muse pixels of 0.2 arcsec in 1 deg^2
;bin_surface_Npixmuse=use_bin_surface*Npixfact
lambir=dustem_get_wavelengths()
Nwavs=n_elements(lambir)
output=dblarr(Nwavs,3)
wavs=*!dustem_plugin_phangs_stellar_continuum.wavs
Mstars=*!dustem_plugin_phangs_stellar_continuum.mstars
Ms_o_Lv=*!dustem_plugin_phangs_stellar_continuum.mstars_o_Lv
Mstarprems=*!dustem_plugin_phangs_stellar_continuum.Mstarprems
;==== sum up the stellar spectrum
;stop
FOR i=2L,Nparam-1 DO BEGIN
IF paramvalues[i] NE 0. THEN BEGIN
;stop
;only interpolate within the template wavelengths
;ind=where(lambir LE max(wavs) AND lambir GE min(wavs),count)
Mstar=Mstars[i-2]
Mstarprem=Mstarprems[i-2]
M_o_Lv=Ms_o_Lv[i-2]
ssp=*(*!dustem_plugin_phangs_stellar_continuum.ssps)[i-2]
;output[*,0]=output[*,0]+paramvalues[i]*Mstar*ssp
;==== below is trying to remove Mstar
;output[*,0]=output[*,0]+paramvalues[i]*ssp
;==== below is trying to multiply weights by L/M
;output[*,0]=output[*,0]+paramvalues[i]/M_o_Lv*ssp*use_bin_surface
;output[*,0]=output[*,0]+paramvalues[i]/M_o_Lv*ssp/use_bin_surface
;output[*,0]=output[*,0]+paramvalues[i]*ssp*Mstar/(bin_surface_pc2)
;output[*,0]=output[*,0]+paramvalues[i]*ssp*Mstar/(bin_surface_Npixmuse)
; output[*,0]=output[*,0]+paramvalues[i]*ssp*Mstar/(use_Voronoi_Npix)
output[*,0]=output[*,0]+paramvalues[i]*ssp*Mstarprem/(use_Voronoi_Npix)
message,'Mstar='+strtrim(Mstar,2)+' paramvalue='+strtrim(paramvalues[i],2),/info
ENDIF
ENDFOR
output[*,0]=output[*,0]*paramvalues[0]
;cgplot,wavs*1.e4,output[*,0],xrange=[4750,7000]
;stop
;flux=output[*,0]
;no polarization for now
;output[*,1]=0.
;output[*,2]=0.
;==== Caution:
;output is Flambda in Lsun/pc2/AA
;at this point, output is unredenned
;==== must be transformed into MJy/sr
Lsun2ergs=3.828e33
pc2_to_cm2=(3.0857e18)^2
fact1=Lsun2ergs/pc2_to_cm2
NHref=*!DUSTEM_HCD ;NH dustemwrap reference [H/cm^2]
fact2=lambir*1.e4 ;[AA] to go from F_lambda in ergs/s/cm2/AA to Fnu ergs/s/cm2
fact3=1./NHref ;[cm2/H] to go from ergs/s/cm2 to ergs/s/H
fact4=1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/lambir)*1.e20/1.e7 ;This goes from ergs/s/H to MJy/sr (see dustem_plot_dataset)
;fact=fact2*fact3*fact4
fact=fact1*fact2*fact3*fact4
;print,'flux(1mic)=',interpol(flux,lambir,1.),' ergs/s/cm2/AA'
;print,'fact2=',interpol(fact2,lambir,1.)
;print,'fact3=',fact3
;print,'fact4=',interpol(fact4,lambir,1.)
;print,'fact=',interpol(fact,lambir,1.)
;print,'flux(1mic)=',interpol(flux,lambir,1.)*interpol(fact,lambir,1.),' MJy/sr'
;goto,no_unred
;==== reden spectrum using Calzetti's extinction law (which is what Phangs does)
;cgplot,lambir,flux,/xlog,/ylog
;factor -1 is to reden, instead of unreden. factor 1e4 is to go from microns to Angstroem
;This uses the prescription of Calzetti et al 2000. It is applied from 912AA to 2.2 mic, no correction applied before and beyound that.
;Note that the supplied color excess should be that derived for the stellar continuum, EBV(stars), which (in Calzetti 2000) is related to the
;reddening derived from the gas, EBV(gas), via the Balmer decrement by EBV(stars) = 0.44*EBV(gas)
calz_unred, lambir*1.e4, output[*,0], -1.*paramvalues[1], flux_red, R_V = R_V
;Note: we could actually use our own dustemwrap extinction curve to do that.
output[*,0]=flux_red
no_unred:
output[*,0]=output[*,0]*fact
;output[*,1]=output[*,1]*fact
;output[*,2]=output[*,2]*fact
;print,'output=',interpol(output[*,0],lambir,1.),' MJy/sr'
;stop
the_scope:
the_paramtag:
the_end:
RETURN,output
END