-
Notifications
You must be signed in to change notification settings - Fork 2
/
figure_4_absorption_example.py
129 lines (109 loc) · 3.47 KB
/
figure_4_absorption_example.py
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
import numpy as np
import astropy.units as u
from agnpy.emission_regions import Blob
from agnpy.targets import lines_dictionary, SphericalShellBLR, RingDustTorus
from agnpy.absorption import Absorption
import matplotlib.pyplot as plt
from pathlib import Path
from utils import time_tau
# distance of the source
z = 1
# dictionary of the shells we want to use
ly_alpha_line = lines_dictionary["Lyalpha"]
H_alpha_line = lines_dictionary["Halpha"]
# build a BLR composed of two spherical shells emitting two different lines: Ly alpha and H alpha
L_disk = 2 * 1e46 * u.Unit("erg s-1")
# the luminosity of the shell formed only by the Lyman Alpha line is a fraction 0.1 of the disk luminosity
xi_ly_alpha = 0.024
L_ly_alpha = xi_ly_alpha * L_disk
R_ly_alpha = 1.1e17 * u.cm
# all the shells radiuses and luminosities are given as a function of the Lybeta shell
# obtain its radius and luminosity from the Ly alpha dictionary
R_ly_beta = R_ly_alpha / ly_alpha_line["R_Hbeta_ratio"]
L_ly_beta = L_ly_alpha / ly_alpha_line["L_Hbeta_ratio"]
# obtain H alpha radius and luminosity from the
L_H_alpha = L_ly_beta * H_alpha_line["L_Hbeta_ratio"]
R_H_alpha = R_ly_beta * H_alpha_line["R_Hbeta_ratio"]
# spherical shells emitting the single lines
blr_ly_alpha = SphericalShellBLR(
L_disk, (L_ly_alpha / L_disk).to_value(""), "Lyalpha", R_ly_alpha
)
blr_H_alpha = SphericalShellBLR(
L_disk, (L_H_alpha / L_disk).to_value(""), "Halpha", R_H_alpha
)
# dust torus
dt = RingDustTorus(L_disk, 0.1, 1000 * u.K)
# print targets to check their parameters
print(blr_ly_alpha)
print(blr_H_alpha)
print(dt)
# distance from the central sources
r = 1e16 * u.cm
# add absorption on synchrotron photons
# blob definition
spectrum_norm = 6e42 * u.erg
parameters = {
"p1": 2.0,
"p2": 3.5,
"gamma_b": 1e4,
"gamma_min": 20,
"gamma_max": 5e7,
}
spectrum_dict = {"type": "BrokenPowerLaw", "parameters": parameters}
R_b = 1e16 * u.cm
B = 0.56 * u.G
z = 1
delta_D = 40
Gamma = 40
blob = Blob(R_b, z, delta_D, Gamma, B, spectrum_norm, spectrum_dict)
# absorptions
abs_blr_ly_alpha = Absorption(blr_ly_alpha, r=r, z=z)
abs_blr_H_alpha = Absorption(blr_H_alpha, r=r, z=z)
abs_dt = Absorption(dt, r=r, z=z)
abs_synch = Absorption(blob)
# array of energies to compute the absorption
nu = np.logspace(24, 30, 100) * u.Hz
# opacities
tau_blr_ly_alpha = time_tau(abs_blr_ly_alpha, nu)
tau_blr_H_alpha = time_tau(abs_blr_H_alpha, nu)
tau_dt = time_tau(abs_dt, nu)
tau_synch = time_tau(abs_synch, nu)
tau_ext = tau_blr_ly_alpha + tau_blr_H_alpha + tau_dt
# make figure 4
fig, ax = plt.subplots()
ax.loglog(
nu,
tau_blr_ly_alpha,
lw=2,
ls="--",
label="BLR, " + r"${\rm Ly\alpha}$" + " shell",
color="crimson",
)
ax.loglog(
nu,
tau_blr_H_alpha,
lw=2,
ls="--",
label="BLR, " + r"${\rm H\alpha}$" + " shell",
color="dodgerblue",
)
ax.loglog(nu, tau_dt, lw=2, ls="--", label="DT", color="goldenrod")
ax.loglog(nu, tau_ext, lw=2, ls="-", label="external", color="k")
ax.fill_between(nu, np.zeros(len(nu)), tau_ext, alpha=0.6, color="darkgray")
ax.loglog(
nu,
1e3 * tau_synch,
lw=2,
color="lightseagreen",
label=r"$10^3 \times$" + " synchrotron",
)
ax.fill_between(
nu, np.zeros(len(nu)), 1e3 * tau_synch, alpha=0.6, color="lightseagreen"
)
ax.legend(loc="best", fontsize=12)
ax.set_xlabel(r"$\nu\,/\,{\rm Hz}$")
ax.set_ylabel(r"$\tau_{\gamma \gamma}$")
ax.set_ylim([1e-1, 1e3])
Path("figures").mkdir(exist_ok=True)
fig.savefig("figures/figure_4.pdf")
fig.savefig("figures/figure_4.png")