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?

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

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.

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?

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

Some typos have been found in our notebook tutorials, will take care of finding and correcting them.

We should try to crosscheck with different approaches the External Compton radiation:

• crosscheck against the literature;
• crosscheck against another code;
• crosscheck against back-of-the-envelope 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 EC computation is reliable.

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

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

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.

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

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'

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)

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

• setting the particle spectrum
  https://github.com/cosimoNigro/agnpy/blob/master/agnpy/emission_regions/blob.py#L200
• computing the particle distribution energy density
  https://github.com/cosimoNigro/agnpy/blob/master/agnpy/emission_regions/blob.py#L277
  and total energy content
  https://github.com/cosimoNigro/agnpy/blob/master/agnpy/emission_regions/blob.py#L286

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.

@jsitarek @dimaniad6

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/")

import agnpy classes

from agnpy.emission_regions import Blob
from agnpy.synchrotron import Synchrotron

z = 0.94 # redshift

delta_D=15
Gamma=15
theta=1./Gamma

far part

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)

let us compute synchrotron SED

synch_sed2 = synch2.sed_flux(nu)

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.

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).

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()

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.

@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?

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

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

1. create a new class (e.g. BLRShells) initialized only with disk luminosity
2. the new class will have store a list of SphericalShellBLR objects and we should add methods to add/remove a single line and add the whole list of lines (from lines_dictionary)
3. in the absorption.py we add opacity_shells_blr which simply make a loop over corresponding opacity_shell_blr (this function takes self as the first parameters, but I hope you can still run it over a specific target in an easy way)
4. similar things we can do also for IC

@cosimoNigro what do you think?
@pawel21 could do the actual implementation

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.

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.

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

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.

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

@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:

@cosimoNigro how were the fits files with the models generated?

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

We should try to crosscheck with different approaches the gamma-gamma optical depth calculation:

• crosscheck against the literature;
• crosscheck against another code;
• crosscheck against back-of-the-envelope 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.

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?

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 path
• from 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/[...]
• some issues with Quantities in Synchrotron.evaluate_sed_flux([....] in notebooks/experiments/basic/synchro_fit_trial.ipynb

Cheers,

Some SED tables we have present energy and some frequency columns.
I will make them uniform with frequency against energy flux.

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.

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.

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.

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:

1. 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

2. 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

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.

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

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

I have run ssc_gammapy_fit.ipynb in my jupyterlab, but I get the following error, I don't know how to solve this problem

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?

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:

• a citation.cff giving instructions on how to properly cite the software;
• a .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...

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.

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.

We should try to crosscheck with different approaches the Synchrotron Self Compton radiation:

• crosscheck against the literature;
• crosscheck against another code;
• crosscheck against back-of-the-envelope 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 SSC computation is reliable.

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 @classmethods, 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).

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?

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.