cosimonigro / agnpy Goto Github PK
View Code? Open in Web Editor NEWModelling jetted Active Galactic Nuclei radiative processes with python
Home Page: https://agnpy.readthedocs.io/en/latest/
License: BSD 3-Clause "New" or "Revised" License
Modelling jetted Active Galactic Nuclei radiative processes with python
Home Page: https://agnpy.readthedocs.io/en/latest/
License: BSD 3-Clause "New" or "Revised" License
In the case of spectral index exactly equal to -2 the analytical normalization function results in division by zero. Either a dedicated analytical integral should be used in this case, or a simple fix it to substitute automatically -2 by a very close number e.g. -2.001
The current code applies the SSA only when the sed_flux is computed of synchrotron radiation, but non-absorbed synchrotron emission is used for the calculation of the IC photons from SSC, while it should be done consistently also with SSA taken into account
In https://github.com/jsitarek/agnpy/blob/master/tests/basic/simple_ec.ipynb I put a few tests of the EC code comparing it with back of the envelope calculations.
in [16] (test2) compares the flux of the EC spectrum for Thomson case without beaming
in [23] (test3) there is the same check done in K-N regime (still without beaming)
in [27] (test4) the same is done for Thomson regime but with strong beaming,
in all the 3 cases the obtained flux is about a factor 2 lower than what I would expect from the back on the envelope calculations
The emission_region
class is getting too long, we should consider moving the constraint on the maximum and break Lorentz factors in a new class.
This class should have the objective of returning a constraint on the physical parameters, how should we call it @jsitarek?
Some SED tables we have present energy and some frequency columns.
I will make them uniform with frequency against energy flux.
Hi,
I have tried use absorption example doc
and I got the following error:
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
<ipython-input-3-b07779d8d11b> in <module>
2 r = 1.1e16 * u.cm
3
----> 4 absorption_disk = Absorption(blob, disk, r=r)
5 absorption_blr = Absorption(blob, blr, r=r)
6 absorption_dt = Absorption(blob, dt, r=r)
~/anaconda3/envs/agn/lib/python3.7/site-packages/agnpy-0.0.7-py3.7.egg/agnpy/absorption.py in __init__(self, blob, target, r)
43 self.target = target
44 self.r = r
---> 45 self.set_mu()
46 self.set_phi()
47 self.set_l()
~/anaconda3/envs/agn/lib/python3.7/site-packages/agnpy-0.0.7-py3.7.egg/agnpy/absorption.py in set_mu(self, mu_size)
49 def set_mu(self, mu_size=100):
50 self.mu_size = mu_size
---> 51 if self.target.type == "SSDisk":
52 # in case of hte disk the mu interval does not go from -1 to 1
53 r_tilde = (self.r / self.target.R_g).to_value("")
AttributeError: 'SSDisk' object has no attribute 'type'
I guess this is due to disk not have attribute 'type', but 'name'.
Could you reploduce my error?
I think it would be good to improve how the integrals are being done in agnpy, for the moment integration is done with simple traperoid integration in each dimention one by one over a predefined grid. There is no check if the integration has converged, and the points are not selected more densely in the regions where the integrant is varying fast. It would be good to use a dedicated multidimentional integration function with a speficied accuracy.
Connected with it, for the case of BLR sometimes the integrant function goes singular. Since the default binning is 50 bins over 5 decades in distance when you select e.g. r = 0.1 * R_line one of the points falls on top of the BLR shell. Proper multidimentional integral should hopefully solve this automatically
I tried to plot EBL absorption from agnpy, however got really strange results.
See the attached notebook.
ebl_absorption_check.ipynb.gz
I think that the underlying fits files with the EBL models are buggy.
In addition there are two problems:
if one tries to interpolate beyond the z range (or energy range) of the model it tries either way.
I think that the extrapolation in energy should work such that it gives tau = 0 below the min energy and tau = infinity above the highest energy
the interpolation is done in the attenuation space, which can change very fast, I think it is much better choice to interpolate tau, and compute the absorption from it
Automatised tests have to be added.
Many of the quantities computed by the code can be analytically checked, they can be embedded in pytest suites.
As suggested by @jsitarek in issue #50, I am opening a different issue for each of the absorption crosschecks.
This one regards the absorption on the photon field of the BLR.
I obtained the opacity vs energy computed at several distances by Finke - before I was using values I had fetched with webplotdigitizer from the paper's figures.
I have made a comparison for all the absorption on the single Lyman alpha line for the 4 distances represented in Figure 14 of Finke (2016): r=10^(-1) R(Ly alpha)
, r=R(Ly alpha)
, r=10^(0.5) R(Ly alpha)
, r=10 R(Ly alpha)
.
The agreement is good only at small distances (within the BLR) and deteriorates at larger ones (outside the BLR).
SBPL normalization from energy density does not work (it gives an error because of wrong units).
BTW, according to the definition in the code, SmoothlyBrokenPowerLaw is not "smoothly" broken, there is a sharp break at gamma_b
In https://github.com/jsitarek/agnpy/blob/master/tests/basic/simple_ec.ipynb I put a few tests of the EC code comparing it with back of the envelope calculations.
in [15] (test2) compares the peak energy of the EC spectrum for Thomson case without beaming, there is a factor of 3.3 difference here, out of which I can understand a factor of 2
in [22] (test3) there is the same check done in K-N regime (still without beaming), here the values match (but they are also limitted by a somewhat different process)
in [26] (test4) the same is done for Thomson regime but with strong beaming, the difference is smaller in this case ~ 1.3, maybe good enough?
We should try to crosscheck with different approaches the Synchrotron Self Compton radiation:
Each of the tick box will be marked once the crosscheck is obtained.
After all of them are marked we will close the issue, and be sure that the SSC computation is reliable.
Hi,
I found and error when installing agnpy
from pip
, that brings the latests version of astropy
(v0.5.1).
When importing agnpy
it raises the following error (see below). I downgrade astropy
to 4.2 and agnpy
can be imported correctly.
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
/opt/conda/lib/python3.8/site-packages/astropy/units/_typing.py in <module>
8 try: # py 3.9+
----> 9 from typing import Annotated
10 except (ImportError, ModuleNotFoundError): # optional dependency
ImportError: cannot import name 'Annotated' from 'typing' (/opt/conda/lib/python3.8/typing.py)
During handling of the above exception, another exception occurred:
AttributeError Traceback (most recent call last)
<ipython-input-1-4024036b86cf> in <module>
----> 1 import agnpy
/opt/conda/lib/python3.8/site-packages/agnpy/__init__.py in <module>
----> 1 from .spectra import *
2 from .emission_regions import *
3 from .targets import *
4 from .synchrotron import *
5 from .compton import *
/opt/conda/lib/python3.8/site-packages/agnpy/spectra/__init__.py in <module>
----> 1 from .spectra import *
/opt/conda/lib/python3.8/site-packages/agnpy/spectra/spectra.py in <module>
1 # module containing the electron spectra
2 import numpy as np
----> 3 import astropy.units as u
4 from ..utils.math import trapz_loglog
5 from ..utils.conversion import mec2
/opt/conda/lib/python3.8/site-packages/astropy/units/__init__.py in <module>
38
39 from .structured import *
---> 40 from .decorators import *
41
42 del bases
/opt/conda/lib/python3.8/site-packages/astropy/units/decorators.py in <module>
11 import numpy as np
12
---> 13 from . import _typing as T
14 from .core import (Unit, UnitBase, UnitsError,
15 add_enabled_equivalencies, dimensionless_unscaled)
/opt/conda/lib/python3.8/site-packages/astropy/units/_typing.py in <module>
16
17 else:
---> 18 from typing_extensions import * # override typing
19
20 HAS_ANNOTATED = Annotated is not NotImplemented
AttributeError: module 'typing_extensions' has no attribute 'OrderedDictTypedDict'
@pawel21 saw that the EBL absorption is a suspiciously low energy (~3 orders of magnitude too low):
https://agnpy.readthedocs.io/en/latest/absorption.html#extragalactic-background-light
The code that is loading the EBL models says that the energies are in eV:
agnpy/agnpy/absorption/absorption.py
Line 785 in 075933e
@cosimoNigro how were the fits files with the models generated?
As we discussed in last week's call @pawel21 will take care of trying to speed-up the SED computations using numba
.
I suggested him to copy, in a script or in a notebook, one of the evaluate_sed_flux
(e.g. the synchrotron one which is the simplest) and try to decorate it with numba.
After my latest modifications (PR #52) functions are @classmethod
s, so you do not need any other agnpy
class (like the blob) to run them. They take directly as input all the physical parameters necessary to compute the SED (before some of them where provided by the blob - like spectral parameters and magnetic field).
the following code fails, while it was running fine in the past (Jun 2020 version).
the error message is:
UnitConversionError: 'cm(5/2) Fr3 g(1/2) / (cm2 J s2)' and 'erg / (cm2 s)' are not convertible
this is probably again a similar problem like we had with transform of B into cgs units, but this time with the electric charge
import numpy as np
import math
import astropy.units as u
import astropy.constants as const
from astropy.coordinates import Distance
import sys
sys.path.append("/home/jsitarek/zdalne/agnpy/agnpy_js/")
#sys.path.append("/home/jsitarek/zdalne/agnpy/agnpy_b2-1420/")
from agnpy.emission_regions import Blob
from agnpy.synchrotron import Synchrotron
z = 0.94 # redshift
delta_D=15
Gamma=15
theta=1./Gamma
spectrum_norm2=2.e-7u.Unit("erg cm-3")
B2=np.sqrt(spectrum_norm2.value8np.pi) u.G
gmin2=2. # minimum Lorentz factor
xi2=3.e-10 # acceleration parameter
dist2=100* 3.e18*u.cm # distance of the blob
p1_2=2.0
p2_2=p1_2+1
radius2=dist2*theta
gbreak=43.e3 # tentative values
gmax=43.e3 # tentative values
parameters2 = {
"p1": p1_2,
"p2": p2_2,
"gamma_b": gbreak,
"gamma_min": gmin2,
"gamma_max": gmax,
}
spectrum_dict2 = {"type": "BrokenPowerLaw", "parameters": parameters2}
blob2 = Blob(radius2, z, delta_D, Gamma, B2, spectrum_norm2, spectrum_dict2, xi=xi2)
nu = np.logspace(8, 28, 100) * u.Hz
synch2 = Synchrotron(blob2, ssa=True)
synch_sed2 = synch2.sed_flux(nu)
I had a look at the n_ph function in the Dust Torus, and I think there is a small bug here:
in this line
Line 391 in 02af6a5
I also need a similar function, but in the blob frame for computing break and max gamma factors for EC. It is basically this function multiplied by Gamma^2 (1-mu * beta)^2, do you prefer it in the target.py or in emission_region.py? It needs variables from both blob and DT.
Hi
First point is not an issue.
I just wanted to mention that I tried to run the published version of agnpy (on zenodo) in mybinder and that it worked fine. I thought you might be interested :)
https://mybinder.org/v2/zenodo/10.5281/zenodo.4311271
A couple of things that did not work:
sys.path.append("/home/jsitarek/zdalne/agnpy/agnpy/")
in notebooks/experiments/basic/SSDisk.ipynb
- you probably don't want to keep this pathfrom agnpy.synchrotron import R, epsilon_equivalency
in notebooks/experiments/basic/trapz_loglog_test.ipynb
does not work../../data/sampled_taus/tau_dt_figure_14_finke_2016.txt
in notebooks/experiments/basic/check_tau_dt.ipynb
should be replaced by ../../agnpy/data/[...]
Quantities
in Synchrotron.evaluate_sed_flux([....]
in notebooks/experiments/basic/synchro_fit_trial.ipynb
Cheers,
hi,
I am getting an error with the following output.
ValueError Traceback (most recent call last)
<ipython-input-94-beb33ca4751d> in <module>
1 r = np.logspace(15, 21) * u.cm
----> 2 u_disk = disk.u(r)
3
C:\Anaconda3\lib\site-packages\agnpy\targets\targets.py in u(self, r, blob)
273 """
274 r_tilde = (r / self.R_g).to_value("")
--> 275 mu = self.evaluate_mu_from_r_tilde(self.R_in_tilde, self.R_out_tilde, r_tilde)
276 prefactor = (
277 3
C:\Anaconda3\lib\site-packages\agnpy\targets\targets.py in evaluate_mu_from_r_tilde(R_in_tilde, R_out_tilde, r_tilde, size)
206 mu_min = 1 / np.sqrt(1 + np.power((R_out_tilde / r_tilde), 2))
207 mu_max = 1 / np.sqrt(1 + np.power((R_in_tilde / r_tilde), 2))
--> 208 return np.linspace(mu_min, mu_max, size)
209
210 @staticmethod
C:\Anaconda3\lib\site-packages\numpy\core\function_base.py in linspace(start, stop, num, endpoint, retstep, dtype)
122 if num > 1:
123 step = delta / div
--> 124 if step == 0:
125 # Special handling for denormal numbers, gh-5437
126 y /= div
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
One thing that we discussed yesterday during the call was to implement a BLR composed of multiple shells. I thought how to do it in the code in the most clean way, and I would suggest such a solution
@cosimoNigro what do you think?
@pawel21 could do the actual implementation
When I have tried to simulate observations with a SpectrumDatasetOnOff.fake() with a Agnpy Spectral Model defined as a wrapper of gammapy SpectralModel, an error has occurred :
ValueError (note: full exception trace is shown but execution is paused at: _run_module_as_main) operands could not be broadcast together with shapes (30,10) (300,)
Thus, the problem appears when program tries to compute
sed = (prefactor * epsilon * emissivity).to("erg cm-2 s-1")
Indeed epsilon is a (30,10) array when emissivity is a 300 one. This sizes come from the fact that self.geom.axes["energy_true"].edges is a (31,)
array.
I'm maybe wrong but I did not solve this error by modifiying the wrapper (cell 2 of https://agnpy.readthedocs.io/en/latest/tutorials/ssc_gammapy_fit.html).
I solved it by replacing the ligne 182 of the synchrotron module of agnpy:
emissivity = integrator(integrand, gamma, axis=0)
by :
emissivity = integrator(integrand, gamma, axis=0).reshape(epsilon.shape)
I am getting an error with the following output.
UnitConversionError Traceback (most recent call last)
<ipython-input-2-ee72776bd796> in <module>
13 # compute the SED over an array of frequencies
14 nu = np.logspace(8, 23) * u.Hz
---> 15 sed = synch.sed_flux(nu)
16
17 # plot it
C:\Anaconda3\lib\site-packages\agnpy\synchrotron\synchrotron.py in sed_flux(self, nu)
236 ssa=self.ssa,
237 integrator=self.integrator,
--> 238 gamma=self.blob.gamma,
239 )
240
C:\Anaconda3\lib\site-packages\agnpy\synchrotron\synchrotron.py in evaluate_sed_flux(nu, z, d_L, delta_D, B, R_b, n_e, ssa, integrator, gamma, *args)
186 emissivity = integrator(integrand, gamma, axis=0)
187 prefactor = np.power(delta_D, 4) / (4 * np.pi * np.power(d_L, 2))
--> 188 sed = (prefactor * epsilon * emissivity).to("erg cm-2 s-1")
189
190 if ssa:
C:\Anaconda3\lib\site-packages\astropy\units\quantity.py in to(self, unit, equivalencies)
667 # and don't want to slow down this method (esp. the scalar case).
668 unit = Unit(unit)
--> 669 return self._new_view(self._to_value(unit, equivalencies), unit)
670
671 def to_value(self, unit=None, equivalencies=[]):
C:\Anaconda3\lib\site-packages\astropy\units\quantity.py in _to_value(self, unit, equivalencies)
639 equivalencies = self._equivalencies
640 return self.unit.to(unit, self.view(np.ndarray),
--> 641 equivalencies=equivalencies)
642
643 def to(self, unit, equivalencies=[]):
C:\Anaconda3\lib\site-packages\astropy\units\core.py in to(self, other, value, equivalencies)
987 return UNITY
988 else:
--> 989 return self._get_converter(other, equivalencies=equivalencies)(value)
990
991 def in_units(self, other, value=1.0, equivalencies=[]):
C:\Anaconda3\lib\site-packages\astropy\units\core.py in _get_converter(self, other, equivalencies)
918 pass
919
--> 920 raise exc
921
922 def _to(self, other):
C:\Anaconda3\lib\site-packages\astropy\units\core.py in _get_converter(self, other, equivalencies)
904 try:
905 return self._apply_equivalencies(
--> 906 self, other, self._normalize_equivalencies(equivalencies))
907 except UnitsError as exc:
908 # Last hope: maybe other knows how to do it?
C:\Anaconda3\lib\site-packages\astropy\units\core.py in _apply_equivalencies(self, unit, other, equivalencies)
888 raise UnitConversionError(
889 "{0} and {1} are not convertible".format(
--> 890 unit_str, other_str))
891
892 def _get_converter(self, other, equivalencies=[]):
UnitConversionError: 'cm(1/2) Fr3 g(1/2) / (J s2)' and 'erg / (cm2 s)' are not convertible
If we were to follow the major.minor.patch
schema, it's been a while we are just "patching" agnpy
😁.
I think we have gathered a set of stable and cross-checked features and we can step-up to our first "minor" release.
What do you think @jsitarek?
Two things I would like to add to our first minor release:
.zenodo.json
file to directly produce the information displayed in the zenodo page (summary, authors, keywords, etc... - now I am adding them manually each time). I have found plenty of examples in github
but no clear instructions on how to actually generate this file, I will give it a try...We should try to crosscheck with different approaches the gamma-gamma optical depth calculation:
Each of the tick box will be marked once the crosscheck is obtained.
After all of them are marked we will close the issue, and be sure that the absorption computation is reliable.
We decided not to include any fit wrapper in the code but to simply show in a couple of tutorial notebooks how to implement one using sherpa
. Sherpa has several minimisation routines and data handling utilities.
We should have a working example with Synchrotron + SSC for a BL Lac, I started it in this notebook using the famous 2011 MWL observations of Mrk421.
The fit is not converging and it has to be fixed.
We discussed the best strategy to introduce hadronic processes and came to the conclusion that we will start with introducing a Particle
class describing either an electron or a proton energy distribution.
This will considerably simplify also the Blob
class. The latter now takes care e.g. of
these are clearly properties of the particle distribution and not of the emission region.
Will open a branch implementing this change and then we keep working on it.
Hi,
I deal with absorption from multi shell in BLR.
I found strange behaviour: when increasted r - distance of the blob from the Black Hole, absorption sometimes it increases sometimes it increases or sometimes go to very very low value.
Here plot absorption vs r for E = 1e5 MeV.
I am attaching notebook with my code.
absorption_blr.ipynb.zip
Now that we have both the release paper and the zenodo we should add clear citation instructions.
We might use this citation file format.
It seems a good option, but it appears also that this file is automatically parsed by zenodo, so we should check that there are no conflicts with the current .zenodo.json
file.
There seems to be some bug in the drawing of the SSDisk luminosity.
The actual luminosity depends on the assumed parameters (radii, mass of BH), while it should only depend on L_disk.
Attached a notebook showing the issue
SSDisk_test.ipynb.gz
We decided not to include any fit wrapper in the code but to simply show in a couple of tutorial notebooks how to implement one using sherpa
. Sherpa has several minimisation routines and data handling utilities.
We should have a working example with Synchrotron + EC + absorption for a FSRQ.
I asked @jsitarek to come up with a good MWL spectrum of a FSRQ to be fitted.
The latest class introduced by @jsitarek lacks documentation and pytest integration.
The tests are actually already in the scripts in the experiments directory, we just have to reformat them in a
agnpy/tests/test_spectral_constraints.py
script that can be run by pytest
.
Also the sphinx documentation has to be updated with them.
So far we had only numpy
, scipy
, astropy
and matplotlib
among the dependencies.
These are specified in setup.py
and automatically downloaded and installed when agnpy
is installed via pip
or conda
.
After #119, I have introduced gammapy
and sherpa
wrappers for the agnpy
radiative models.
This means that there are now classes in the new agnpy.fit
submodule that are built on sherpa
and gammapy
objects.
Trying to simulate an installation from pip, I created a new conda environment.
mamba create -n agnpy_0.2.0_test python=3.9
conda activate agnpy_0.2.0_test
I created a distribution (using the same command used by the github action automatically publishing our releases on PyPI) and installed it
python setup.py sdist bdist_wheel
cd dist
pip install agnpy-0.2.0.tar.gz
When I try to import any object I get the following error:
>>> from agnpy.emission_regions import Blob
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/Users/cosimo/software/miniconda3/envs/agnpy_0.2.0_test/lib/python3.9/site-packages/agnpy/__init__.py", line 9, in <module>
from .fit import *
File "/Users/cosimo/software/miniconda3/envs/agnpy_0.2.0_test/lib/python3.9/site-packages/agnpy/fit/__init__.py", line 1, in <module>
from .models import *
File "/Users/cosimo/software/miniconda3/envs/agnpy_0.2.0_test/lib/python3.9/site-packages/agnpy/fit/models.py", line 1, in <module>
from .gammapy_wrapper import (
File "/Users/cosimo/software/miniconda3/envs/agnpy_0.2.0_test/lib/python3.9/site-packages/agnpy/fit/gammapy_wrapper.py", line 6, in <module>
from gammapy.modeling import Parameter, Parameters
ModuleNotFoundError: No module named 'gammapy.modeling'
Curiously, just repeating the import a second time in the same python terminal does not raise any error
>>> from agnpy.emission_regions import Blob
>>> blob = Blob()
>>> print(blob)
* spherical emission region
- R_b (radius of the blob): 1.00e+16 cm
- t_var (variability time scale): 4.13e-01 d
- V_b (volume of the blob): 4.19e+48 cm3
- z (source redshift): 0.07
- d_L (source luminosity distance):9.92e+26 cm
- delta_D (blob Doppler factor): 1.00e+01
- Gamma (blob Lorentz factor): 1.00e+01
- Beta (blob relativistic velocity): 9.95e-01
- theta_s (jet viewing angle): 5.74e+00 deg
- B (magnetic field tangled to the jet): 1.00e+00 G
- xi (coefficient for 1st order Fermi acceleration) : 1.00e+00
* electron spectrum
- power law
- k_e: 9.27e+06 1 / cm3
- p: 2.80
- gamma_min: 1.00e+02
- gamma_max: 1.00e+07
I think the only solution to avoid the error is to add sherpa
and gammapy
to the install_requires
in setup.py
.
I do not think it is a big deal since these packages live basically all in the same environment.
On the other hand, if some users want to use agnpy
for modelling, then they might not want to install sherpa
and gammapy
as well.
@jsitarek, any opinion?
@adonath, is there a way to skip the import of certain submodules (in this case of agnpy.fit
)? If it is possible, would it be much more work than simply adding sherpa
and gammapy
to the dependencies?
Some typos have been found in our notebook tutorials, will take care of finding and correcting them.
As suggested by @jsitarek in issue #50, I am opening a different issue for each of the absorption crosschecks.
This one regards the absorption on the photon field of the DT.
I obtained the opacity vs energy computed at several distances by Finke - before I was using values I had fetched with webplotdigitizer from the paper's figures.
Here again the same problem of the BLR in issue #65 repeats: for small distances, within the disk, there is agreement, at larger distances r=10^2 R(Ly alpha)
there is mismatch.
On the other hand, at very large distances we have the consistency check from the approximation with the point-source behind the jet.
So I am a bit confused.
We should try to crosscheck with different approaches the External Compton radiation:
Each of the tick box will be marked once the crosscheck is obtained.
After all of them are marked we will close the issue, and be sure that the EC computation is reliable.
The changes made in the class since May somehow broke the way how synchrotron spectrum is calculated. Please find attached a notebook where I compare the old and new version and plot also a theoretical line
validate_synchr.ipynb.gz
The parameters of the emission region in Figure 7.4 of Dermer and Menon 2009 are repeated over and over in code, notebooks and tests.
I propose to hardcode its values in the default values of the Blob.
In this way
blob = Blob()
will automatically return a blob with the parameters in Figure 7.4 of Dermer Menon 2009.
Lines 31 to 45 in 489f1cc
this part of code is Eq 78 in Finke'16 paper (in the description there is also a small mistake mentioning Eq 76 instead),
however the formula in the paper is for mu_star^2, and the code is simply taking a sqrt of it, while the value of mu_star can be also negative for parts of the BLR shell. This most likely explains the difference that I was having with simple calculations of the energy threshold, and too good agreement between DT and BLR for a blob well inside the BLR sphere
We added the absorption on synchrotron photons to the code and the tests, but not to the docs.
A paragraph should be added in the Absorption section.
I have typed bullet points in several pages.
Despite they appear correctly when I build locally, they are not displayed on ReadTheDocs.
E.g. in https://agnpy.readthedocs.io/en/latest/spectra.html
The different particle distributions should have a bullet to the left of their names.
I tried the solution pointed out here
https://stackoverflow.com/questions/67542699/readthedocs-sphinx-not-rendering-bullet-list-from-rst-file
but some of them (e.g. fixing the sphinx and sphinx-rtd-theme versions) did not work.
Will try some of the other suggestions.
HI @cosimoNigro
I had a look into implementing a proper calculations of the absorption for a case of mu_s != 1 (moving at an angle to the jet axis)
I want to test is carefully because there is a lot of fun with geometry, and different "phi" angles, but you can have a look at the current code in this branch if you have any comments about the code structure:
https://github.com/jsitarek/agnpy/tree/absorption_mu_s_dependent
the idea for the code would be to leave the current absorption functions (we might remove the mu_s parameter from them, and also avoid "phi" integrals later on), and add new functions with ..._mu_s in the name that would calculate the more complicated integrals with those proper direction.
I think this can be also a very nice addition to the paper, I could write the corresponding formulas in a dedicated appendix.
attached below some results from comparing absorption at mu_s=1, with the current code (that we know that is buggy) and this new code. qualitatively it makes sense, but I want to have a look
if I can make also some quantitative checks.
The range of mu used for integration changes with the position of the blob above the BH (r) according to
mu_min = 1 / sqrt(1 + (R_out / r)^2)
mu_max = 1 / sqrt(1 + (R_in / r)^2)
In the absorption though an array of mu calculated at a single r (the one occupied by the blob) is used
https://github.com/cosimoNigro/agnpy/blob/master/agnpy/absorption.py#L97
when integrating from r to infinity for the optical depth, the array of mu has to be recomputed at each value of r and updated.
For the same reason the same problem might occur in the photon density:
https://github.com/cosimoNigro/agnpy/blob/master/agnpy/targets.py#L202
The integration has to be separated somehow, I'll take care of that.
As suggested by @jsitarek in issue #50, I am opening a different issue for each of the absorption crosschecks.
This one regards the absorption on the photon field of the SS Disk.
I obtained the opacity vs energy computed at several distances by Finke - before I was using values I had fetched with webplotdigitizer from the paper's figures.
I just compared two distances (r=10^(-1) * R(Ly alpha)
and r=R(Ly alpha)
, the same represented in Figure 14 of Finke 2016) and in both cases there is mismatch both in the energy threshold and in the maximum value of the absorption.
Let us continue the discussion started in issue #50 here.
Dear @jsitarek, @pawel21, @mwcraig,
I will be presenting agnpy
at the ADASS XXX with a poster.
You can find it at this link.
Let me know if you have any comment.
Small edits can be directly implemented in the overleaf, for major changes please leave a comment.
I listed all the contributors (i.e. the four of us) as authors.
Thanks
@jsitarek, @pawel21, let use this thread to keep discussing what I commented in PR #46.
what do you think about the mismatch of agnpy with Finke (2016) optical depth?
why my optical depths are consistent if I consider instead of the target a monochromatic point-like source approximating it.
Is the integration procedure wrong, so the two cases - full target and point-like source - match because I perform the same wrong integration in both cases?
do we have another code to crosscheck optical depths, I remember @jsitarek had available the code from
https://arxiv.org/abs/0807.4228?
Investigating the DT SED computation I realized that this part is quite problematic. If you compute DT SED for nu that are not sorted (e.g. if you generate lambdas and convert them to nu) the value will be wrong (I got negative). Also if the calculations are done only at course frequencies (or not covering the whole emission) this renormalization can screw up the integral emission. BB emission can be integrated analitically, so it would be better to use such analytical integration instead.
agnpy/agnpy/targets/targets.py
Line 596 in c7ede0c
Not all VHE instruments provide de-absorbed spectral points.
We need to add EBL absorption in the Absorption
class.
I will take care of this.
The calculation of the absorption (and I presume also of the External Compton) on the SphericalShellBLR
fails when r=R_line
.
This is due to the definition of x
(the distance from the reprocessing material to the blob), from Finke 2016:
x^2 = r^2 + R_line^2 - 2 * r * R_re * mu_re
where
r
= height of the blob above the BH;R_re
= radial coordinate of the reprocessing material;mu_re
= zenith angle of the reprocessing material.Clearly when r=R_re
, x=0
and the absorption (or EC) calculations gives nan
, as the optical depth (or the SED) is dependent on x^(-2)
.
Here a small snippet to trigger the issue:
import numpy as np
import astropy.units as u
import astropy.constants as const
from agnpy.targets import SphericalShellBLR
from agnpy.absorption import Absorption
# BLR definition
L_disk = 2 * 1e46 * u.Unit("erg s-1")
xi_line = 0.024
R_line = 1.1 * 1e17 * u.cm
blr = SphericalShellBLR(L_disk, xi_line, "Lyalpha", R_line)
# absorption
abs_blr = Absorption(blr, r=R_line)
# tau
nu = np.logspace(22, 28) * u.Hz
tau_blr = abs_blr.tau(nu)
print(tau_blr)
returning:
[nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan
nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan
nan nan nan nan nan nan nan nan nan nan nan nan nan nan]
@jsitarek, shall we put a check in the x
calculation?
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.