dustem_plugin_emiles_stellar_continuum.pro
5.99 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
FUNCTION dustem_plugin_emiles_stellar_continuum,key=key $
,val=val $
,scope=scope $
,paramtag=paramtag $
,paramdefault=paramdefault $
,help=help
;+
; NAME:
; dustem_plugin_emiles_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_emiles_stellar_continuum([,key=][,val=][,scope=][,paramtag=][paramdefault=][,/help])
; INPUTS:
; None
; OPTIONAL INPUT PARAMETERS:
; key = input parameter number
; First parameter: factor to apply to tau for redenning (1 means value derived from emission SED)
; 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
; EXAMPLES
; dustem_init
; vec=dustem_plugin_emiles_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_emiles_stellar_continuum'
output=0.
goto,the_end
ENDIF
;stop
;default values of input parameters
scope='ADD_SED'
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)
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]
Nage=n_elements(age_values)
Nmetalicity=n_elements(metalicity_values)
Nparam=Nage*Nmetalicity+1 ;+1 is for opacity
paramtag=strarr(Nparam)
paramdefault=dblarr(Nparam)
paramdefault[*]=0.D0
i=0L
paramtag[0]='tau_f'
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
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 present
defsysv,'!dustem_plugin_emiles_stellar_continuum',exist=exist
IF exist EQ 0 THEN BEGIN
st={ssps:ptr_new(),wavs:ptr_new(),Mstars:ptr_new(),Mstarprems:ptr_new()}
defsysv,'!dustem_plugin_emiles_stellar_continuum',st
ENDIF
IF not ptr_valid(!dustem_plugin_emiles_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,Mstarprems=Mstarprems)
!dustem_plugin_emiles_stellar_continuum.ssps=ptr_new(template_flux)
!dustem_plugin_emiles_stellar_continuum.wavs=ptr_new(template_wav)
!dustem_plugin_emiles_stellar_continuum.mstars=ptr_new(Mstars)
!dustem_plugin_emiles_stellar_continuum.mstarprems=ptr_new(Mstarprems)
ENDIF
lambir=dustem_get_wavelengths()
Nwavs=n_elements(lambir)
output=fltarr(Nwavs,3)
wavs=*!dustem_plugin_emiles_stellar_continuum.wavs
Mstars=*!dustem_plugin_emiles_stellar_continuum.mstars
Mstarprems=*!dustem_plugin_emiles_stellar_continuum.Mstarprems
;==== sum up the stellar spectrum
;stop
FOR i=1L,Nparam-1 DO BEGIN
;IF paramvalues[i] NE 0. AND ptr_valid((*!phangs_ssps)[i-1]) THEN BEGIN
IF paramvalues[i] NE 0. THEN BEGIN
;stop
;only interpolate within the template wavelengths
;ind=where(lambir LE max(*!phangs_wavs) AND lambir GE min(*!phangs_wavs),count)
ind=where(lambir LE max(wavs) AND lambir GE min(wavs),count)
;Mstar=(*!phangs_Mstars)[i]
Mstar=Mstars[i-1]
Mstarprem=Mstarprems[i-1]
IF count NE 0 THEN BEGIN
;output[ind,0]=output[ind,0]+paramvalues[i]*interpol(*(*!phangs_ssps)[i-1],*!phangs_wavs,lambir[ind])
;spectrum=*(*!phangs_ssps)[i-1]
;ssp=*(*!phangs_ssps)[i-1]
ssp=*(*!dustem_plugin_emiles_stellar_continuum.ssps)[i-1]
;stop
output[ind,0]=output[ind,0]+paramvalues[i]*Mstar*interpol(ssp,wavs,lambir[ind])
ENDIF
ENDIF
ENDFOR
;==== Caution: at this point, output is in Lsun/pc^2/Angstroem
;Caution: I am not sure about the nomalisation of the extinction curve that is in
IF ptr_valid(!dustem_current) THEN BEGIN
;stop
tau_fact=paramvalues[0]
extinction_curve=((*!dustem_current).ext).ext_tot*(*!dustem_HCD)/1.0e21 ;This is tau
;*(*!dustem_HCD)/1.0e21
;wav_ref=551.e-9*1.e6 ;visible 551 nm in microns
;tau_ref=interpol(extinction_curve,lambir,wav_ref)
;normalize extinction curve by given tau fact
extinction_curve=extinction_curve*tau_fact
;stop
output[*,0]=output[*,0]*exp(-1.*extinction_curve)
ENDIF
;no polarization for now
output[*,1]=0.
output[*,2]=0.
;at that stage, output is F_lambda in Lsun/A/pc^2
;goint to MJy/sr (1e-20 W/m2/sr/hz)
pc=3.085677581e16 ;m
;kpc=3.085677581e19 ;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
output=output*fact
the_end:
RETURN,output
END