fritz2006.py 10.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
# -*- coding: utf-8 -*-
# Copyright (C) 2013, 2014 Department of Physics, University of Crete
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Laure Ciesla

"""
Fritz et al. (2006) AGN dust torus emission module
==================================================

10
This module implements the Fritz et al. (2006) models.
11
12

"""
13
from collections import OrderedDict
14

15
import numpy as np
16

17
from pcigale.data import Database
18
from . import SedModule
19

20
import scipy.constants as cst
Médéric Boquien's avatar
Médéric Boquien committed
21

22
class Fritz2006(SedModule):
23
24
    """Fritz et al. (2006) AGN dust torus emission

25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
    The AGN emission is computed from the library of Fritz et al. (2006) from
    which all of the models are available. They take into account two emission
    components linked to the AGN. The first one is the isotropic emission of
    the central source, which is assumed to be point-like. This emission is a
    composition of power laws with variable indices, in the wavelength range of
    0.001-20 microns. The second one is the thermal and scattering dust torus
    emission. The conservation of the energy is always verified within 1% for
    typical solutions, and up to 10% in the case of very high optical depth and
    non-constant dust density. We refer the reader to Fritz et al. (2006) for
    more information on the library.

    The relative normalization of these components is handled through a
    parameter which is the fraction of the total IR luminosity due to the AGN
    so that: L_AGN = fracAGN * L_IRTOT, where L_AGN is the AGN luminosity,
    fracAGN is the contribution of the AGN to the total IR luminosity
    (L_IRTOT), i.e. L_Starburst+L_AGN.
41
42
43

    """

44
    parameter_list = OrderedDict([
45
        ('r_ratio', (
46
            "cigale_list(options=10. & 30. & 60. & 100. & 150.)",
47
48
            "Ratio of the maximum to minimum radii of the dust torus. "
            "Possible values are: 10, 30, 60, 100, 150.",
49
            60.
50
51
        )),
        ('tau', (
52
53
            "cigale_list(options=0.1 & 0.3 & 0.6 & 1.0 & 2.0 & 3.0 & 6.0 & "
            "10.0)",
54
55
            "Optical depth at 9.7 microns. "
            "Possible values are: 0.1, 0.3, 0.6, 1.0, 2.0, 3.0, 6.0, 10.0.",
56
            1.0
57
58
        )),
        ('beta', (
59
60
            "cigale_list(options=-1.00 & -0.75 & -0.50 & -0.25 & 0.00)",
            "Beta. Possible values are: -1.00, -0.75, -0.50, -0.25, 0.00.",
61
            -0.50
62
63
        )),
        ('gamma', (
64
            'cigale_list(options=0.0 & 2.0 & 4.0 & 6.0)',
65
            "Gamma. Possible values are: 0.0, 2.0, 4.0, 6.0.",
66
            4.0
67
68
        )),
        ('opening_angle', (
69
            'cigale_list(options=60. & 100. & 140.)',
70
            "Full opening angle of the dust torus (Fig 1 of Fritz 2006). "
71
            "Possible values are: 60., 100., 140.",
72
            100.
73
74
        )),
        ('psy', (
75
76
            'cigale_list(options=0.001 & 10.100 & 20.100 & 30.100 & 40.100 & '
            '50.100 & 60.100 & 70.100 & 80.100 & 89.990)',
77
            "Angle between equatorial axis and line of sight. "
Médéric Boquien's avatar
Médéric Boquien committed
78
79
80
            "Psy = 90◦ for type 1 and Psy = 0° for type 2. Possible values "
            "are: 0.001, 10.100, 20.100, 30.100, 40.100, 50.100, 60.100, "
            "70.100, 80.100, 89.990.",
81
            50.100
82
83
        )),
        ('fracAGN', (
84
            'cigale_list(minvalue=0., maxvalue=1.)',
85
            "AGN fraction.",
86
            0.1
87
88
89
90
91
92
93
94
        )),
        ('EBV', (
            'cigale_list(minvalue=0.)',
            "E(B-V) for extinction in polar direction",
            0.1
        )),
        ('temperature', (
            'cigale_list(minvalue=0.)',
95
            "Temperature of the polar dust in K",
96
97
98
99
            100.
        )),
        ("emissivity", (
            "cigale_list(minvalue=0.)",
100
            "Emissivity index of the polar dust.",
101
            1.6
102
        ))
103
   ])
104
105
106

    def _init_code(self):
        """Get the template set out of the database"""
107
108
109
110
111
112
113
        self.r_ratio = float(self.parameters["r_ratio"])
        self.tau = float(self.parameters["tau"])
        self.beta = float(self.parameters["beta"])
        self.gamma = float(self.parameters["gamma"])
        self.opening_angle = (180. - self.parameters["opening_angle"]) / 2.
        self.psy = float(self.parameters["psy"])
        self.fracAGN = float(self.parameters["fracAGN"])
114
115
116
        self.EBV = float(self.parameters["EBV"])
        self.temperature = float(self.parameters["temperature"])
        self.emissivity = float(self.parameters["emissivity"])
117
118

        with Database() as base:
119
120
121
            self.fritz2006 = base.get_fritz2006(self.r_ratio, self.tau,
                                                self.beta, self.gamma,
                                                self.opening_angle, self.psy)
122
123
124
125
        self.l_agn_scatt = np.trapz(self.fritz2006.lumin_scatt,
                                    x=self.fritz2006.wave)
        self.l_agn_agn = np.trapz(self.fritz2006.lumin_agn,
                                  x=self.fritz2006.wave)
126

Guang's avatar
Guang committed
127
        # Apply wavelength cut to avoid X-ray wavelength
128
        lam_cut = 10**0.9
Guang's avatar
Guang committed
129
130
131
132
133
134
135
136
137
138
        lam_idxs = self.fritz2006.wave>=lam_cut
        # Calculate the re-normalization factor to keep energy conservation
        norm_fac = np.trapz(self.fritz2006.lumin_intrin_agn, x=self.fritz2006.wave) /\
            np.trapz(self.fritz2006.lumin_intrin_agn[lam_idxs], x=self.fritz2006.wave[lam_idxs])
        # Perform the cut
        self.fritz2006.wave = self.fritz2006.wave[lam_idxs]
        self.fritz2006.lumin_agn = self.fritz2006.lumin_agn[lam_idxs]*norm_fac
        self.fritz2006.lumin_scatt = self.fritz2006.lumin_scatt[lam_idxs]*norm_fac
        self.fritz2006.lumin_therm  = self.fritz2006.lumin_therm[lam_idxs]
        self.fritz2006.lumin_intrin_agn = self.fritz2006.lumin_intrin_agn[lam_idxs]*norm_fac
Guang's avatar
Guang committed
139

Guang's avatar
Guang committed
140
        # Apply polar-dust obscuration
141
        # We define various constants necessary to compute the model
Guang's avatar
Guang committed
142
        self.c = cst.c * 1e9
143
144
145
146
147
        lambda_0 = 200e3
        # Calculate the extinction (SMC)
        k_lam = 1.39*(self.fritz2006.wave/1e3)**-1.2-0.38
        A_lam = k_lam * self.EBV
        # Set negative value to zero
Guang's avatar
Guang committed
148
        A_lam[A_lam<0] = 0
149
150
151
        # The extinction factor, flux_1/flux_0
        ext_fac = 10**(A_lam/-2.5)
        # Calculate the new AGN SED shape after extinction
Guang's avatar
Guang committed
152
153
154
155
156
157
158
159
        if self.psy > (90-self.opening_angle):
            # The direct and scattered components (line-of-sight) are extincted for type-1 AGN
            lumin_agn_new = self.fritz2006.lumin_agn * ext_fac
            lumin_scatt_new = self.fritz2006.lumin_scatt * ext_fac
        else:
            # Keep the direct and scatter components for type-2
            lumin_agn_new = self.fritz2006.lumin_agn
            lumin_scatt_new = self.fritz2006.lumin_scatt
160
        # Calculate the total extincted luminosity averaged over all directions
Guang's avatar
Guang committed
161
        # note that self.opening_angle is different from Fritz's definiting!
162
163
        l_ext = np.trapz(self.fritz2006.lumin_intrin_agn * (1-ext_fac),
                         x=self.fritz2006.wave) * \
Guang's avatar
Guang committed
164
                        (1 - np.cos( np.deg2rad(self.opening_angle) ))
165
        # Casey (2012) modified black body model
Guang's avatar
Guang committed
166
        conv = self.c / (self.fritz2006.wave * self.fritz2006.wave)
167
        lumin_blackbody = (conv * (1. - np.exp(-(lambda_0 / self.fritz2006.wave)
Guang's avatar
Guang committed
168
169
                                ** self.emissivity)) * (self.c / self.fritz2006.wave) ** 3. / (np.exp(
                                cst.h * self.c / (self.fritz2006.wave * cst.k * self.temperature)) - 1.))
170
171
172
173
174
175
176
177
        lumin_blackbody *= l_ext / np.trapz(lumin_blackbody, x=self.fritz2006.wave)
        # Add the black body to dust thermal emission
        lumin_therm_new = self.fritz2006.lumin_therm + lumin_blackbody
        # Normalize direct, scatter, and thermal components
        norm = np.trapz(lumin_therm_new, x=self.fritz2006.wave)
        lumin_therm_new /= norm
        lumin_scatt_new /= norm
        lumin_agn_new /= norm
Guang's avatar
Guang committed
178
179
180
181
        # Update fritz2006 lumin
        self.fritz2006.lumin_therm = lumin_therm_new
        self.fritz2006.lumin_scatt = lumin_scatt_new
        self.fritz2006.lumin_agn = lumin_agn_new
182
        self.fritz2006.lumin_intrin_agn /= norm
Guang's avatar
Guang committed
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

        # Integrate AGN luminosity for different components
        self.l_agn_scatt = np.trapz(self.fritz2006.lumin_scatt, x=self.fritz2006.wave)
        self.l_agn_agn = np.trapz(self.fritz2006.lumin_agn, x=self.fritz2006.wave)
        # Intrinsic (de-reddened) AGN luminosity from the central source
        self.l_agn_intrin_agn = np.trapz(self.fritz2006.lumin_intrin_agn, x=self.fritz2006.wave)
        # Calculate L_lam(2500A)
        self.l_agn_2500A = np.interp(250, self.fritz2006.wave, self.fritz2006.lumin_intrin_agn)
        # Convert L_lam to L_nu
        self.l_agn_2500A *= 250**2/self.c

    def process(self, sed):
        """Add the IR re-emission contributions

        Parameters
        ----------
        sed: pcigale.sed.SED object
        parameters: dictionary containing the parameters

        """

        if 'dust.luminosity' not in sed.info:
            sed.add_info('dust.luminosity', 1., True)
        luminosity = sed.info['dust.luminosity']

        sed.add_module(self.name, self.parameters)
        sed.add_info('agn.r_ratio', self.r_ratio)
        sed.add_info('agn.tau', self.tau)
        sed.add_info('agn.beta', self.beta)
        sed.add_info('agn.gamma', self.gamma)
        sed.add_info('agn.opening_angle', self.parameters["opening_angle"])
        sed.add_info('agn.psy', self.psy)
        sed.add_info('agn.fracAGN', self.fracAGN)
        sed.add_info('agn.EBV', self.EBV)
        sed.add_info('agn.temperature', self.temperature)
        sed.add_info('agn.emissivity', self.emissivity)
219

220
        # Compute the AGN luminosity
221
222
        if self.fracAGN < 1.:
            agn_power = luminosity * (1./(1.-self.fracAGN) - 1.)
223
            l_agn_therm = agn_power
224
225
            l_agn_scatt = agn_power * self.l_agn_scatt
            l_agn_agn = agn_power * self.l_agn_agn
226
            l_agn_total = l_agn_therm + l_agn_scatt + l_agn_agn
Guang's avatar
Guang committed
227
228
            l_agn_intrin_agn = agn_power * self.l_agn_intrin_agn
            l_agn_2500A = agn_power * self.l_agn_2500A
229
230
231
        else:
            raise Exception("AGN fraction is exactly 1. Behaviour "
                            "undefined.")
232

233
234
        sed.add_info('agn.therm_luminosity', l_agn_therm, True)
        sed.add_info('agn.scatt_luminosity', l_agn_scatt, True)
Guang's avatar
Guang committed
235
        sed.add_info('agn.disk_luminosity', l_agn_agn, True)
236
        sed.add_info('agn.luminosity', l_agn_total, True)
237
        sed.add_info('agn.accretion_power', l_agn_intrin_agn, True)
Guang's avatar
Guang committed
238
        sed.add_info('agn.intrin_Lnu_2500A', l_agn_2500A, True)
239
240

        sed.add_contribution('agn.fritz2006_therm', self.fritz2006.wave,
Guang's avatar
Guang committed
241
                             agn_power * self.fritz2006.lumin_therm)
242
        sed.add_contribution('agn.fritz2006_scatt', self.fritz2006.wave,
Guang's avatar
Guang committed
243
                             agn_power * self.fritz2006.lumin_scatt)
244
        sed.add_contribution('agn.fritz2006_agn', self.fritz2006.wave,
Guang's avatar
Guang committed
245
                             agn_power * self.fritz2006.lumin_agn)
246

247

248
# SedModule to be returned by get_module
249
Module = Fritz2006