forked from cosimoNigro/agnpy_paper
-
Notifications
You must be signed in to change notification settings - Fork 0
/
appendix_B_sed_resolution.py
138 lines (133 loc) · 4.71 KB
/
appendix_B_sed_resolution.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
130
131
132
133
134
135
136
137
138
# script to study the effect of the integration grid on the final SED resolution
import sys
import numpy as np
import astropy.units as u
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
from agnpy.emission_regions import Blob
from agnpy.targets import RingDustTorus
from agnpy.compton import ExternalCompton
from agnpy.utils.plot import load_mpl_rc, sed_x_label, sed_y_label
from pathlib import Path
from utils import logging
# blob
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)
# dust torus definition
L_disk = 2 * 1e46 * u.Unit("erg s-1")
T_dt = 1e3 * u.K
xi_dt = 0.1
dt = RingDustTorus(L_disk, xi_dt, T_dt)
# different resolution
n_gamma_list = [100, 200, 300, 400, 500]
n_phi_list = [50, 100, 200, 300, 400]
colors = ["crimson", "dodgerblue", "goldenrod", "lightseagreen", "k"]
# same frequency grid
nu = np.logspace(15, 29, 50) * u.Hz
# show the effect of changing the Lorentz factor and frequency grid
seds_variable_gamma = []
labels_variable_gamma = []
for n_gamma in n_gamma_list:
logging.info(f"computing SED with {n_gamma} Lorentz factor points")
blob.set_gamma_size(n_gamma)
ec = ExternalCompton(blob, dt, r=1e21 * u.cm)
sed = ec.sed_flux(nu)
label = r"$n_{\gamma}=$" + f"{n_gamma}, " + r"$n_{\phi}=50,\,n_{\nu}=100$"
seds_variable_gamma.append(sed)
labels_variable_gamma.append(label)
# show the effect of changing the azimuth angle grid
seds_variable_phi = []
labels_variable_phi = []
for n_phi in n_phi_list:
logging.info(f"computing SED with {n_phi} azimuth points")
blob.set_gamma_size(500)
ec = ExternalCompton(blob, dt, r=1e21 * u.cm)
ec.set_phi(n_phi)
sed = ec.sed_flux(nu)
label = r"$n_{\gamma}=500,\,n_{\phi}=$" + f"{n_phi}" + r"$,\,n_{\nu}=100$"
seds_variable_phi.append(sed)
labels_variable_phi.append(label)
# figure
load_mpl_rc()
plt.rcParams["text.usetex"] = True
# gridspec plot setting
fig = plt.figure(figsize=(12, 6), tight_layout=True)
spec = gridspec.GridSpec(ncols=2, nrows=2, height_ratios=[2, 1], figure=fig)
ax1 = fig.add_subplot(spec[0, 0])
ax2 = fig.add_subplot(spec[0, 1])
ax3 = fig.add_subplot(spec[1, 0], sharex=ax1)
ax4 = fig.add_subplot(spec[1, 1], sharex=ax2, sharey=ax3)
# changing Lorentz factor and frequency grid
for sed, label, color in zip(
seds_variable_gamma[:-1], labels_variable_gamma[:-1], colors[:-1]
):
ax1.loglog(nu, sed, ls="-", color=color, label=label)
# plot the last one as reference
ax1.loglog(
nu,
seds_variable_gamma[-1],
ls="--",
color=colors[-1],
label=labels_variable_gamma[-1],
)
ax1.set_ylabel(sed_y_label)
ax1.legend(loc="best", fontsize=10)
ax1.set_title("EC on ring DT, " + r"$r=10^{21}\,{\rm cm}$")
# changing azimuth angle grid
for sed, label, color in zip(
seds_variable_phi[:-1], labels_variable_phi[:-1], colors[:-1]
):
ax2.loglog(nu, sed, ls="-", color=color, label=label)
# plot the last one as reference
ax2.loglog(
nu, seds_variable_phi[-1], ls="--", color=colors[-1], label=labels_variable_phi[-1]
)
ax2.legend(loc="best", fontsize=10)
ax2.set_title("EC on ring DT, " + r"$r=10^{21}\,{\rm cm}$")
# plot the deviation from the denser SED in the bottom panel
for sed, label, color in zip(
seds_variable_gamma[:-1], labels_variable_gamma[:-1], colors[:-1]
):
deviation = (sed / seds_variable_gamma[-1]) - 1
ax3.semilogx(nu, deviation, ls="-", color=color, label=label)
ax3.grid(False)
ax3.axhline(0, ls="-", color="darkgray")
ax3.axhline(0.2, ls="--", color="darkgray")
ax3.axhline(-0.2, ls="--", color="darkgray")
ax3.axhline(0.3, ls=":", color="darkgray")
ax3.axhline(-0.3, ls=":", color="darkgray")
ax3.set_ylim([-0.5, 0.5])
ax3.set_yticks([-0.4, -0.2, 0.0, 0.2, 0.4])
ax3.set_xlabel(sed_x_label)
ax3.set_ylabel(r"$\frac{\nu F_{\nu, \rm agnpy}}{\nu F_{\nu, \rm ref}} - 1$")
# plot the deviation from the denser SED in the bottom panel
for sed, label, color in zip(
seds_variable_phi[:-1], labels_variable_phi[:-1], colors[:-1]
):
deviation = (sed / seds_variable_phi[-1]) - 1
ax4.semilogx(nu, deviation, ls="-", color=color, label=label)
ax4.grid(False)
ax4.axhline(0, ls="-", color="darkgray")
ax4.axhline(0.2, ls="--", color="darkgray")
ax4.axhline(-0.2, ls="--", color="darkgray")
ax4.axhline(0.3, ls=":", color="darkgray")
ax4.axhline(-0.3, ls=":", color="darkgray")
ax4.set_ylim([-0.5, 0.5])
ax4.set_yticks([-0.4, -0.2, 0.0, 0.2, 0.4])
ax4.set_xlabel(sed_x_label)
Path("figures").mkdir(exist_ok=True)
fig.savefig(f"figures/figure_appendix_B.png")
fig.savefig(f"figures/figure_appendix_B.pdf")