Fix the black body SEDs of the Disk and the Dust Torus #92

cosimoNigro · 2021-05-08T10:32:37Z

this PR implements a correct multi-temperature black body SED for the Disk and a (single-temperature) black body SED for the DT.
I renormalise these SEDs to the respective luminosities of the Disk and DT, i.e. I impose that the integral luminosity of the SED has to be equal to the luminosity of the disk (or the luminosity of the DT).
I added these checks of the luminosity to the tests.

This should solve #23.

Changes in files other than targets.py and test_targets.py are unrelated, and need no checking, they are mostly corrections of the docstrings.

Here a snippet to test the changes, developed on the notebook you posted in issue #23 (checkout fix_disk_and_dt_bb)

import numpy as np
import astropy.units as u
from astropy.constants import M_sun
from astropy.coordinates import Distance
import matplotlib.pyplot as plt
from agnpy.targets import SSDisk, RingDustTorus
from agnpy.utils.plot import load_mpl_rc, plot_sed

load_mpl_rc()

# create disks with same luminosity but different radii
L_disk = 6 * 1e44 * u.Unit("erg s-1")
M_BH = 1e8 * M_sun
z = 0.06
disk1 = SSDisk(M_BH, L_disk, 1 / 12, 6, 100, R_g_units=True)
disk2 = SSDisk(M_BH, L_disk, 1 / 12, 6, 1000, R_g_units=True)
disk3 = SSDisk(M_BH, L_disk, 1 / 12, 6, 10000, R_g_units=True)

nu = np.logspace(10, 20, 100) * u.Hz
sed_disk1 = disk1.sed_flux(nu, z)
sed_disk2 = disk2.sed_flux(nu, z)
sed_disk3 = disk3.sed_flux(nu, z)

plot_sed(nu, sed_disk1, ls="-", label=r"$R_{\rm out} = 10^2 R_{\rm g}$")
plot_sed(nu, sed_disk2, ls="-", label=r"$R_{\rm out} = 10^3 R_{\rm g}$")
plot_sed(nu, sed_disk3, ls="-", label=r"$R_{\rm out} = 10^4 R_{\rm g}$")
plt.ylim([1e-15, 1e-9])
plt.legend()
plt.show()

# compute back the luminosity
print(f"disk reference Luminosity: {L_disk:.2e}")
d_L = Distance(z=z).to("cm")
for sed in [sed_disk1, sed_disk2, sed_disk3]:
    F_nu = sed / nu
    L = 4 * np.pi * np.power(d_L, 2) * np.trapz(F_nu, nu, axis=0)
    print(f"SED luminosity: {L:.2e}")

# create dust toruses with same luminosity but different radii
xi_dt = 0.5
dt1 = RingDustTorus(L_disk, xi_dt, 1000 * u.K)
dt2 = RingDustTorus(L_disk, xi_dt, 2000 * u.K)
dt3 = RingDustTorus(L_disk, xi_dt, 5000 * u.K)

nu = np.logspace(10, 20, 100) * u.Hz
sed_dt1 = dt1.sed_flux(nu, z)
sed_dt2 = dt2.sed_flux(nu, z)
sed_dt3 = dt3.sed_flux(nu, z)
plot_sed(nu, sed_dt1, ls="-", label=r"$T_{\rm dt} = 10^3 \, {\rm K}$")
plot_sed(nu, sed_dt2, ls="-", label=r"$T_{\rm dt} = 2 \times 10^3 \, {\rm K}$")
plot_sed(nu, sed_dt3, ls="-", label=r"$T_{\rm dt} = 5 \times 10^3 \, {\rm K}$")
plt.ylim([1e-15, 1e-9])
plt.legend()
plt.show()

# compute back the luminosity
L_dt = xi_dt * L_disk
print(f"DT reference Luminosity: {L_dt:.2e}")
d_L = Distance(z=z).to("cm")
for sed in [sed_dt1, sed_dt2, sed_dt3]:
    F_nu = sed / nu
    L = 4 * np.pi * np.power(d_L, 2) * np.trapz(F_nu, nu, axis=0)
    print(f"SED luminosity: {L:.2e}")

disk reference Luminosity: 6.00e+44 erg / s
SED luminosity: 6.00e+44 erg / s
SED luminosity: 6.00e+44 erg / s
SED luminosity: 6.00e+44 erg / s

DT reference Luminosity: 3.00e+44 erg / s
SED luminosity: 3.00e+44 erg / s
SED luminosity: 3.00e+44 erg / s
SED luminosity: 3.00e+44 erg / s

codecov · 2021-05-08T10:36:52Z

Codecov Report

Merging #92 (80cca68) into master (8c89b47) will increase coverage by 1.02%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master      #92      +/-   ##
==========================================
+ Coverage   94.46%   95.48%   +1.02%     
==========================================
  Files          30       30              
  Lines        2094     2127      +33     
==========================================
+ Hits         1978     2031      +53     
+ Misses        116       96      -20

Flag	Coverage Δ
unittests	`95.48% <100.00%> (+1.02%)`	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
agnpy/compton/synchrotron_self_compton.py	`87.80% <ø> (ø)`
agnpy/synchrotron/synchrotron.py	`94.11% <ø> (ø)`
agnpy/compton/external_compton.py	`87.50% <100.00%> (ø)`
agnpy/targets/targets.py	`90.62% <100.00%> (+12.09%)`	⬆️
agnpy/tests/test_targets.py	`100.00% <100.00%> (ø)`
agnpy/utils/math.py	`100.00% <100.00%> (ø)`

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 8c89b47...80cca68. Read the comment docs.

jsitarek

sorry for a delay on checking this code, I spotted an issue with lack of observation angle dependence on the disk luminosity.

jsitarek · 2021-05-12T13:28:42Z

agnpy/targets/targets.py

@@ -343,17 +336,50 @@ def phi_disk_mu(self, mu, r_tilde):
    def phi_disk(self, R_tilde):
        return 1 - np.sqrt(self.R_in_tilde / R_tilde)

-    def T(self, R_tilde):
+    def T(self, R):
        r"""temperature of the disk at radius :math:`\tilde{R}`. 


\tilde(R) ==> R in docstring

jsitarek · 2021-05-12T13:56:43Z

agnpy/targets/targets.py

+        sed_disk = SSDisk.evaluate_multi_T_bb_sed(nu, z, M_BH, m_dot, R_in, R_out, d_L)
+        # renormalise
+        L = (np.trapz(sed_disk / nu, nu) * 4 * np.pi * d_L ** 2).to("erg s-1")
+        norm = (L_disk / L).to_value("")


there is a small issue here, for the renormalization you assume that the emission is isotropic (hence 4 pi), but since it is a black body flat surface (for a thin disk approximation) it should emit ~ cos (observation angle)
if you integrate it over the sphere I think what you are missing in norm is an additional factor: 2* cos (observation_angle).
you would also need to update the test accordingly.

Note that this is only for the disk, the dust torus emission is reprocessed so it should be basically isotropic at high distance

done, I have included the cosine in the SED calculation (in order to obtain the correct nuFnu depending on the viewing angle) and added the factor 2 in the luminosity calculations

jsitarek · 2021-05-12T13:59:34Z

agnpy/tests/test_targets.py

+        d_L = Distance(z=z).to("cm")
+        F_nu = sed / nu
+        L = 4 * np.pi * np.power(d_L, 2) * np.trapz(F_nu, nu, axis=0)
+        assert u.isclose(L, L_disk_test, atol=0 * u.Unit("erg s-1"), rtol=1e-2)


here you will also need to multiply L_disk_test times 2*cos (observation_angle)

the test is done for observation angle=0, do you think it is work to include in the parametrizaiton also a second case with observation at an angle (to check those cos factors), or would it be overdoing the testing?

Thanks for asking to add it, made me realise that the new variable mu_s was not propagated to all the functions computing the SEDs. Now it should be ok. Added three values of mu_s to the tests.

…its luminosities in renormalisation and tests

Fix the black body SEDs of the Disk and the Dust Torus (Sourcery refactored)

cosimoNigro · 2021-06-25T11:18:31Z

Hello @jsitarek,

so sorry for the huge delay in this, June has been full of meetings!

I took into account in your suggestions, now I multiply the nuFnu SED by the cosine of the viewing angle (newly introduced mu_s parameter with default 1). This way we get the correct SED value for an inclined observer.
This is how the docstring reads:

I am not totally convinced it is sufficient to multiply by cosine to get the correct integral, maybe it works only when you simplify R^2/d_L^2 to 1?
To be honest I am not sure how to go from the mus in blue to the green ones in the figure I sketched

When you incline the observer you change not only the viewing angles but also the distances to the disk surface...
Only reference I found giving this formula is this one, Eq. 10.

Beside the mathematical discussion, are you satisfied with the changes now?
Thanks!
Cosimo

jsitarek · 2021-06-25T11:34:01Z

thx @cosimoNigro
I'm not sure what do you mean by "when you simplify R^2/d_L^2 to 1?"

the function is used to determine the SED of the whole disk at the observer, i.e. you integrate the full disk, and then I think cos should work fine (it is also there in the formula in the reference that you provided. As far as I can understand it , the cos is there due to the black body property, it simply emits more perpendicularly to the surface than at an angle to it. And since the solid angle in which you see the disk is very small, also the change of distances along it does not matter.

The situation gets more complicated when you are close to the disk, and see the emission from its different points at different angles, than you also need to take into account that with a given small solid angle you will also catch a different part of the disk surface depending on this viewing angle. This needs to be taken into account when computing EC and absorption from disk, but for just showing the disk SED at the observer is not important

cosimoNigro · 2021-06-25T11:47:06Z

I'm not sure what do you mean by "when you simplify R^2/d_L^2 to 1?"

if you compare Eq. 10 of the paper and the formula in the docstrings I wrote it's clear they are simplifying the factor (1 + R^2/d_L^2) by assuming R << d_L, right? I thought that only in that case you could transform everything multiplying just by the cosine of the viewing angle. But I think I was wrong.

Let me know if you are satisfied with the changes (you can merge if you are).

Sorry again for taking so long.

jsitarek · 2021-06-25T12:02:53Z

Hi @cosimoNigro
Now I understand what you mean, the formula is only valid for distant observer, so R<< d_L.

cos angle to the surface normal is always there, but only if the observer is far away from the disk it is the same angle as the observation angle. If you need to calculate it close to the disk, not only the distances will vary, but also each part of the disk will be observed at a different angle. I think this function is only meant for plotting SED for distant observer, then I suggest that the doc string is revised to take this into account, and give the cosine explicitely. If you need a similar function for computing the density of the photons from the disk either for IC or for pair production, this is more complicated, but since those functions were implemented from Dermen's and Finke's paper they should already have this cosine hidden somewhere.

cosimoNigro · 2021-06-25T14:46:37Z

What about now?

About this:

this is more complicated, but since those functions were implemented from Dermen's and Finke's paper they should already have this cosine hidden somewhere.

looking at the formulas that we re-written and the schemes we draw I think it is always assumed that the disk axis matches the jet axis, which I think it's fine for most calculations. It's the observer that is in turn inclined to an angle \theta_s to the jet (and disk) axes.

jsitarek · 2021-06-25T15:01:28Z

Sorry for being picky, but I think that the part 1+ R^2 /dL^2 is really confusing here. The formula is correct only if R<< dL, so I would suggest to skip this (1+...).
Otherwise it gives an impression that it might work also close to the disk (even if the docstring says overwise)

sourcery-ai · 2021-06-25T15:26:22Z

Sourcery Code Quality Report

❌ Merging this PR will decrease code quality in the affected files by 0.24%.

Quality metrics	Before	After	Change
Complexity	0.91 ⭐	0.87 ⭐	-0.04 👍
Method Length	62.17 🙂	62.82 🙂	0.65 👎
Working memory	10.94 😞	11.04 😞	0.10 👎
Quality	69.88% 🙂	69.64% 🙂	-0.24% 👎

Other metrics	Before	After	Change
Lines	1999	2015	16

Changed files	Quality Before	Quality After	Quality Change
agnpy/compton/external_compton.py	59.59% 🙂	59.61% 🙂	0.02% 👍
agnpy/compton/synchrotron_self_compton.py	67.77% 🙂	67.78% 🙂	0.01% 👍
agnpy/synchrotron/synchrotron.py	72.06% 🙂	72.06% 🙂	0.00%
agnpy/targets/targets.py	68.79% 🙂	68.90% 🙂	0.11% 👍
agnpy/tests/test_targets.py	82.83% ⭐	80.05% ⭐	-2.78% 👎
agnpy/utils/math.py	63.55% 🙂	63.55% 🙂	0.00%
docs/conf.py	80.51% ⭐	80.38% ⭐	-0.13% 👎

Here are some functions in these files that still need a tune-up:

File	Function	Complexity	Length	Working Memory	Quality	Recommendation
agnpy/compton/external_compton.py	ExternalCompton.evaluate_sed_flux_ss_disk	0 ⭐	261 ⛔	25 ⛔	40.26% 😞	Try splitting into smaller methods. Extract out complex expressions
agnpy/targets/targets.py	SSDisk.__init__	9 🙂	208 ⛔	13 😞	43.83% 😞	Try splitting into smaller methods. Extract out complex expressions
agnpy/compton/external_compton.py	ExternalCompton.evaluate_sed_flux_blr	0 ⭐	197 😞	23 ⛔	44.85% 😞	Try splitting into smaller methods. Extract out complex expressions
agnpy/utils/math.py	trapz_loglog	3 ⭐	238 ⛔	15 😞	45.05% 😞	Try splitting into smaller methods. Extract out complex expressions
agnpy/compton/external_compton.py	ExternalCompton.evaluate_sed_flux_dt	0 ⭐	182 😞	22 ⛔	46.50% 😞	Try splitting into smaller methods. Extract out complex expressions

Legend and Explanation

The emojis denote the absolute quality of the code:

⭐ excellent
🙂 good
😞 poor
⛔ very poor

The 👍 and 👎 indicate whether the quality has improved or gotten worse with this pull request.

Please see our documentation here for details on how these metrics are calculated.

We are actively working on this report - lots more documentation and extra metrics to come!

Help us improve this quality report!

cosimoNigro · 2021-06-25T15:27:11Z

No problem, removed that factor from the docs, and - of course - from the actual formula in the code!

jsitarek · 2021-06-25T15:33:37Z

I think it is fine now, thx

cosimoNigro · 2021-06-25T15:45:54Z

I think you have to accept my changes or accept the review, it still says "changes requested".
After that you can merge yourself.

jsitarek · 2021-06-25T15:56:10Z

right, sorry, I just accepted the review and merged

cosimoNigro added 3 commits May 8, 2021 11:39

created a proper multi-temperature black body SED for the Disk

e5ee6de

minor fix in docs

cd9a13a

added also normalised function to compute the black body SED of the DT

d6215c3

cosimoNigro requested a review from jsitarek May 8, 2021 10:32

cosimoNigro added 2 commits May 8, 2021 15:26

added PKS1510-089 2012 MWL SED to agnpy/data

d17ff41

added low state PKS1510-089 SED in data

16001fe

jsitarek requested changes May 12, 2021

View reviewed changes

cosimoNigro and others added 3 commits June 25, 2021 12:33

fixed typo in disk T docstring

92a465a

added viewing angle correction to SS disk SED, multiplied by 2 sides …

95f1a85

…its luminosities in renormalisation and tests

'Refactored by Sourcery'

b3edb05

sourcery-ai bot mentioned this pull request Jun 25, 2021

Fix the black body SEDs of the Disk and the Dust Torus (Sourcery refactored) #95

Merged

Merge pull request #95 from cosimoNigro/sourcery/fix_disk_and_dt_bb

daf9a3b

Fix the black body SEDs of the Disk and the Dust Torus (Sourcery refactored)

added tests on SS Disk nuFnu SED including viewing angle

e7a8f69

expanded docstring of multi-T BB SED for the SS Disk

c641866

removed factor [1 + (R/d_L)^2] from docs and code in BB SS disk SED

80cca68

jsitarek approved these changes Jun 25, 2021

View reviewed changes

jsitarek merged commit bec5989 into master Jun 25, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix the black body SEDs of the Disk and the Dust Torus #92

Fix the black body SEDs of the Disk and the Dust Torus #92

cosimoNigro commented May 8, 2021 •

edited

Loading

codecov bot commented May 8, 2021 •

edited

Loading

jsitarek left a comment

jsitarek May 12, 2021

cosimoNigro Jun 25, 2021

jsitarek May 12, 2021

cosimoNigro Jun 25, 2021

jsitarek May 12, 2021

cosimoNigro Jun 25, 2021

jsitarek Jun 25, 2021

cosimoNigro Jun 25, 2021 •

edited

Loading

cosimoNigro commented Jun 25, 2021

jsitarek commented Jun 25, 2021

cosimoNigro commented Jun 25, 2021

jsitarek commented Jun 25, 2021

cosimoNigro commented Jun 25, 2021 •

edited

Loading

jsitarek commented Jun 25, 2021

sourcery-ai bot commented Jun 25, 2021

cosimoNigro commented Jun 25, 2021 •

edited

Loading

jsitarek commented Jun 25, 2021

cosimoNigro commented Jun 25, 2021 •

edited

Loading

jsitarek commented Jun 25, 2021

Fix the black body SEDs of the Disk and the Dust Torus #92

Fix the black body SEDs of the Disk and the Dust Torus #92

Conversation

cosimoNigro commented May 8, 2021 • edited Loading

codecov bot commented May 8, 2021 • edited Loading

Codecov Report

jsitarek left a comment

Choose a reason for hiding this comment

jsitarek May 12, 2021

Choose a reason for hiding this comment

cosimoNigro Jun 25, 2021

Choose a reason for hiding this comment

jsitarek May 12, 2021

Choose a reason for hiding this comment

cosimoNigro Jun 25, 2021

Choose a reason for hiding this comment

jsitarek May 12, 2021

Choose a reason for hiding this comment

cosimoNigro Jun 25, 2021

Choose a reason for hiding this comment

jsitarek Jun 25, 2021

Choose a reason for hiding this comment

cosimoNigro Jun 25, 2021 • edited Loading

Choose a reason for hiding this comment

cosimoNigro commented Jun 25, 2021

jsitarek commented Jun 25, 2021

cosimoNigro commented Jun 25, 2021

jsitarek commented Jun 25, 2021

cosimoNigro commented Jun 25, 2021 • edited Loading

jsitarek commented Jun 25, 2021

sourcery-ai bot commented Jun 25, 2021

Sourcery Code Quality Report

Legend and Explanation

cosimoNigro commented Jun 25, 2021 • edited Loading

jsitarek commented Jun 25, 2021

cosimoNigro commented Jun 25, 2021 • edited Loading

jsitarek commented Jun 25, 2021

cosimoNigro commented May 8, 2021 •

edited

Loading

codecov bot commented May 8, 2021 •

edited

Loading

cosimoNigro Jun 25, 2021 •

edited

Loading

cosimoNigro commented Jun 25, 2021 •

edited

Loading

cosimoNigro commented Jun 25, 2021 •

edited

Loading

cosimoNigro commented Jun 25, 2021 •

edited

Loading