dustem_plugin_phangs_stellar_continuum.pro
6.1 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
FUNCTION dustem_plugin_phangs_stellar_continuum,key=key $
,val=val $
,scope=scope $
,paramtag=paramtag $
,paramdefault=paramdefault $
,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_brightness=dustem_plugin_phangs_stellar_continuum([,key=][,val=][,scope=][,paramtag=][paramdefault=][,/help])
; INPUTS:
; None
; OPTIONAL INPUT PARAMETERS:
; key = input parameter number
; First parameter: E(B-V) for extinction 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_brightness = stellar spectrum spectrum resampled on dustem wavelengths [MJy/sr]
; OPTIONAL OUTPUT PARAMETERS:
; scope = scope of the plugin
; paramdefault = default values of parameters
; paramtag = plugin parameter names as strings
; ACCEPTED KEY-WORDS:
; help = if set, print this help
; 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
out=0
goto, the_scope
ENDIF
;IF keyword_set(paramtag) THEN BEGIN
; out=0
; goto, the_paramtag
;ENDIF
;stop
;default values of input parameters
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]
;read template info
info_st=dustem_read_emiles_stellar_templates(age_values,metalicity_values,/info_only)
;stop
age_values_emiles=info_st.age
order=sort(age_values_emiles)
age_values_emiles=age_values_emiles[order]
un=uniq(age_values_emiles)
age_values_emiles=age_values_emiles[un]
met_values_emiles=info_st.z
order=sort(met_values_emiles)
met_values_emiles=met_values_emiles[order]
un=uniq(met_values_emiles)
met_values_emiles=met_values_emiles[un]
;==== 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
Nage=n_elements(age_values)
Nmetalicity=n_elements(metalicity_values)
Nparam=Nage*Nmetalicity+1 ;+1 is for opacity
paramdefault=dblarr(Nparam)
paramdefault[*]=0.D0 ;default parameter values is 0.
;==== define parameter tags
;If the test below is activated (ie paramtags computed only when requested) the pluggin does not work.
;IF keyword_set(paramtag) THEN BEGIN
paramtag=strarr(Nparam)
paramtag[0]='E(B-V)'
ii=1L
FOR i=1L,Nparam-1 DO BEGIN
ij=index2ij([i-1],[Nage,Nmetalicity])
paramtag[ii]='Age'+strtrim(ij[0,0],2)+'Met'+strtrim(ij[0,1],2)
ii=ii+1
ENDFOR
;out=0.
;GOTO,the_paramtag
;ENDIF
;==== 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()}
defsysv,'!dustem_plugin_phangs_stellar_continuum',st
ENDIF
IF not ptr_valid(!dustem_plugin_phangs_stellar_continuum.ssps) THEN BEGIN
message,'Initializing SSP templates',/continue
template_flux=dustem_read_emiles_stellar_templates(age_values,metalicity_values,template_wav=template_wav,Mstars=Mstars)
!dustem_plugin_phangs_stellar_continuum.ssps=ptr_new(template_flux)
!dustem_plugin_phangs_stellar_continuum.wavs=ptr_new(template_wav)
!dustem_plugin_phangs_stellar_continuum.mstars=ptr_new(Mstars)
ENDIF
;stop
lambir=dustem_get_wavelengths()
Nwavs=n_elements(lambir)
output=fltarr(Nwavs,3)
wavs=*!dustem_plugin_phangs_stellar_continuum.wavs
Mstars=*!dustem_plugin_phangs_stellar_continuum.mstars
;==== sum up the stellar spectrum
;stop
FOR i=1L,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-1]
IF count NE 0 THEN BEGIN
ssp=*(*!dustem_plugin_phangs_stellar_continuum.ssps)[i-1]
output[ind,0]=output[ind,0]+paramvalues[i]*Mstar*interpol(ssp,wavs,lambir[ind])
message,'Mstar='+strtrim(Mstar,2)+' paramvalue='+strtrim(paramvalues[i],2),/continue
ENDIF
ENDIF
ENDFOR
;==== Caution: at this point, output is unredenned Flambda in Lsun/pc^2/Angstroem
stop
;==== reden spectrum using Calzetti's extinction law (which is what Phangs does)
flux=output[*,0]
;factor -1 is to reden, instead of unreden. factor 1e4 is to go from microns to Angstroem
calz_unred, lambir*1.e4, flux, -1.*paramvalues[0], flux_red, R_V = R_V
output[*,0]=flux_red
;no polarization for now
output[*,1]=0.
output[*,2]=0.
;at that stage, output is F_lambda in Lsun/A/pc^2
;going to W/m2/sr/hz
pc=3.085677581e16 ;m
Lsun=3.828e26 ;Watts
c=3.e5*1e3 ;m/s
;Fact=Lsun ;Watts/A
;fact=fact*c ;Watts/hz
;fact=fact/4./!pi ;Watts/Hz/sr
;fact=fact/pc^2 ;Watts/m2/Hz/sr
fact=Lsun/pc*c/4./!pi/pc ;to go from F_lambda in Lsun/A/pc^2 to Watts/m2/Hz/sr
output=output*fact
;print,fact
; 9.59803
the_scope:
scope='ADD_SED'
the_paramtag:
the_end:
RETURN,output
END