OV-GSO-DC / dustem-wrapper_idl

Blame view

src/idl/dustem_plugin_phangs_stellar_continuum.pro 5.32 KB
Edit Raw Normal View History
b58e9569   Jean-Philippe Bernard   First commit 
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
 
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

;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]='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
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_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])
		ENDIF
	ENDIF
ENDFOR
;==== Caution: at this point, output is unredenned Flambda in Lsun/pc^2/Angstroem

;==== 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
;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