<p>There are two issues here: one of definitions and one of code.</p>
<p>First, the value you're expecting to get here -- 1361 Watts per square meter, the solar constant -- is the total flux received from the Sun at 1 AU, i.e. on Earth. Here, you've set the distance at 1 parsec; since that's larger by a factor of <span class="math-container">$\sim10^5$</span>, this swap has decreased your answer by a factor of <span class="math-container">$\sim10^{10}$</span>!</p>
<p>The code-related problem likely has to do with how the wavelength integral<span class="math-container">$^{\dagger}$</span> is being computed by <code>integrate.quad</code>. <span class="math-container">$B_{\lambda}(\lambda, T)$</span> peaks at visible wavelengths (which you can show via <a href="https://en.wikipedia.org/wiki/Wien%27s_displacement_law" rel="nofollow noreferrer">Wien's law</a>), at a few hundred nanometers. The main contribution to the integral is within an order of magnitude or so from this range. However, you've set the upper and lower limits of the integral to <span class="math-container">$0$</span> and <span class="math-container">$2\times10^{10}$</span> meters, respectively. <code>integrate.quad</code> is evaluating <span class="math-container">$B_{\lambda}(\lambda, T)$</span> on some grid of wavelength values, and with the range you've picked, the vast majority of points on that grid are outside the range of relevant wavelengths, at which the Planck function is effectively 0.</p>
<p>(One way you can check this is to reduce the upper limit to, say, <span class="math-container">$10^5$</span> meters so that more of your grid points are at wavelengths closer to the range of interest. You should see that, counterintuitively, <em>decreasing</em> the upper limit causes the value your code spits out to <em>increase</em>.)</p>
<p>There are several ways you could go about addressing this (changing the limits, for example, or doing a logarithmic change-of-variables as suggested by David Z), but the simplest solution would be to not have your code do the integration at all. The integral of <span class="math-container">$B_{\lambda}(\lambda, T)$</span> over all wavelengths has an analytical solution; after we integrate over all wavelengths and solid angles, we get
<span class="math-container">$$L=\int d\Omega\;\cos(\theta)\int B_{\lambda}(\lambda, T)\;d\lambda=\sigma T^4$$</span>
with <span class="math-container">$\sigma$</span> the Stefan-Boltzmann constant. This avoids the numerical problems you've run into.</p>
<hr />
<p><span class="math-container">$^{\dagger}$</span>We also need to integrate <span class="math-container">$B_{\lambda}(\lambda, T)$</span> over all solid angles, but I see your code implicitly does this by multiplying by <span class="math-container">$4\pi^2$</span>, since the Planck function is isotropic.</p>

CC BY-SA 4.0
2023-11-27T20:01:28.083Z
Thank you very much. I changed the limits to (0, 1000) m and got 6.852257448717239e-15 W/m**2, which is still too far.
Anna-Kat
CC BY-SA 4.0
2023-11-27T20:17:21.543Z
@Anna-Kat I think that's still too high a limit -- reducing it by several more orders of magnitude than that did it for me. (I'm sure there's still a more elegant way of doing it than this, though!)
HDE 226868
CC BY-SA 4.0
2023-11-29T09:02:14.240Z
@HDE226868 [Adaptive quadrature](https://en.wikipedia.org/wiki/Adaptive_quadrature), probably, although of course even that won't help if the range of integration is large enough and the grid sparse enough to completely miss any interesting features in the function. Personally, for something like this which spans several orders of magnitude, I might try a logarithmic change of variables ($u = \ln(\lambda/L)$ for some arbitrarily chosen $L$) to hopefully increase the chance that some grid points land in the right spots.
David Z
CC BY-SA 4.0
2023-12-06T16:24:07.383Z
@DavidZ Yeah, logarithms might be the best way to go about it.
HDE 226868
CC BY-SA 4.0
2024-03-29T04:25:17.887Z
@HDE226868 Could you please find the mistake here https://astronomy.stackexchange.com/questions/57307/calculating-suns-j-band-mag-via-integrating-the-planck-law ?
questionhang
CC BY-SA 4.0
2023-11-27T18:53:13.437Z
2023-12-06T16:24:46.763Z
2023-12-06T16:24:46.763Z
HDE 226868
<p>The integral of the Planck function over all wavelengths gives the surface flux at the photosphere of the Sun divided by pi.</p>
<p>In order to get the solar constant, you should multiply this by <span class="math-container">$\pi \times 4\pi R_\odot^2$</span> and then divide by <span class="math-container">$4\pi d^2$</span>, where <span class="math-container">$d$</span> is an astronomical unit in metres. This looks correct.</p>
<p>The integration limits should be set to something sensible. The Sun isn't a blackbody, though the effective temperature should reasonably represent the total flux. I would use a lower limit of 100 nm and an upper limit of 5000 nm.</p>

CC BY-SA 4.0
2023-11-27T20:18:36.100Z
2023-11-27T20:18:36.100Z
ProfRob
<p>There is two issues:</p>
<p>a) you have to use Earth distance to get the Solar constant, but you use 1pc (<span class="math-container">$10^{16}$</span>m)instead of one astronomical unit (150 million kilometers, <span class="math-container">$1.5\cdot 10^{11}$</span>m). That's a factor of <span class="math-container">$10^{10}$</span> in your results.</p>
<p>b) Use more reasonable bounds like 10nm to 1cm. With the short wavelength bound reduced my python runs into an overflow in the exponential function. With the long-wavelength bound larger the integrator complains about not converging anymore.</p>
<p>Changing both, gives a result of 1356W/m^2</p>

CC BY-SA 4.0
2023-11-27T20:31:51.140Z
2023-11-27T20:38:23.560Z
2023-11-27T20:38:23.560Z
planetmaker
<p>Why I did not get the value of 1361 W m^{-2}, please? I am getting 2.0125538614437745e-47.</p>
<pre><code>from scipy import integrate
import numpy as np
k = 1.38e-23  # J K^-1
h = 6.63e-34  # J s
c = 299792458.  # m s^-1
pc = 3.08567758e16  # m
R_S = 6.957e8  # m
pc = 3.08567758e16  # m
R_S = 6.957e8  # m


R = 1 * R_S

d = 1 * pc

T = 5780.  # K

def B_lambda(lam, T):
    return (2. * h * c**2) / (lam**5) * 1. / (np.exp(h * c / (lam * k * T)) - 1.)

lambda_min = 0  # Lower limit for integration
lambda_max = 2e+10  # Upper limit for integration

result, error = integrate.quad(B_lambda, lambda_min, lambda_max, args=(T,))
print(result*np.pi*R**2/d**2, error)
</code></pre>
<p><a href="https://i.sstatic.net/ON9dn.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/ON9dn.png" alt="enter image description here" /></a></p>

CC BY-SA 4.0
2023-11-27T20:14:20.027Z
Please update your question if you have made changes to the code/results.
ProfRob
<the-sun><flux>
Integration of Plank curve of the Sun
Anna-Kat
# Integration of Plank curve of the Sun

(Asked by: username_0 on 2023-11-27)

Why I did not get the value of 1361 W m^{-2}, please? I am getting 2.0125538614437745e-47.
from scipy import integrate
import numpy as np
k = 1.38e-23  # J K^-1
h = 6.63e-34  # J s
c = 299792458.  # m s^-1
pc = 3.08567758e16  # m
R_S = 6.957e8  # m
pc = 3.08567758e16  # m
R_S = 6.957e8  # m


R = 1 * R_S

d = 1 * pc

T = 5780.  # K

def B_lambda(lam, T):
    return (2. * h * c**2) / (lam**5) * 1. / (np.exp(h * c / (lam * k * T)) - 1.)

lambda_min = 0  # Lower limit for integration
lambda_max = 2e+10  # Upper limit for integration

result, error = integrate.quad(B_lambda, lambda_min, lambda_max, args=(T,))
print(result*np.pi*R**2/d**2, error)

**Tags:** the-sun, flux
> **Comments:**
> **username_1**: Please update your question if you have made changes to the code/results.  

---

## 3 answers

### Answer by username_5 on 2023-11-27. Score: 8 | Accepted answer

There is two issues:

a) you have to use Earth distance to get the Solar constant, but you use 1pc ($10^{16}$m)instead of one astronomical unit (150 million kilometers, $1.5\cdot 10^{11}$m). That's a factor of $10^{10}$ in your results.

b) Use more reasonable bounds like 10nm to 1cm. With the short wavelength bound reduced my python runs into an overflow in the exponential function. With the long-wavelength bound larger the integrator complains about not converging anymore.

Changing both, gives a result of 1356W/m^2


---

### Answer by username_2 on 2023-11-27. Score: 9

There are two issues here: one of definitions and one of code.

First, the value you're expecting to get here -- 1361 Watts per square meter, the solar constant -- is the total flux received from the Sun at 1 AU, i.e. on Earth. Here, you've set the distance at 1 parsec; since that's larger by a factor of $\sim10^5$, this swap has decreased your answer by a factor of $\sim10^{10}$!

The code-related problem likely has to do with how the wavelength integral$^{\dagger}$ is being computed by `integrate.quad`. $B\_{\lambda}(\lambda, T)$ peaks at visible wavelengths (which you can show via [Wien's law](https://en.wikipedia.org/wiki/Wien%27s_displacement_law)), at a few hundred nanometers. The main contribution to the integral is within an order of magnitude or so from this range. However, you've set the upper and lower limits of the integral to $0$ and $2\times10^{10}$ meters, respectively. `integrate.quad` is evaluating $B\_{\lambda}(\lambda, T)$ on some grid of wavelength values, and with the range you've picked, the vast majority of points on that grid are outside the range of relevant wavelengths, at which the Planck function is effectively 0.

(One way you can check this is to reduce the upper limit to, say, $10^5$ meters so that more of your grid points are at wavelengths closer to the range of interest. You should see that, counterintuitively, *decreasing* the upper limit causes the value your code spits out to *increase*.)

There are several ways you could go about addressing this (changing the limits, for example, or doing a logarithmic change-of-variables as suggested by David Z), but the simplest solution would be to not have your code do the integration at all. The integral of $B\_{\lambda}(\lambda, T)$ over all wavelengths has an analytical solution; after we integrate over all wavelengths and solid angles, we get
$$L=\int d\Omega\;\cos(\theta)\int B\_{\lambda}(\lambda, T)\;d\lambda=\sigma T^4$$
with $\sigma$ the Stefan-Boltzmann constant. This avoids the numerical problems you've run into.

---

$^{\dagger}$We also need to integrate $B\_{\lambda}(\lambda, T)$ over all solid angles, but I see your code implicitly does this by multiplying by $4\pi^2$, since the Planck function is isotropic.

> **Comments:**
> **username_0**: Thank you very much. I changed the limits to (0, 1000) m and got 6.852257448717239e-15 W/m**2, which is still too far.  
> **username_2**: @Anna-Kat I think that's still too high a limit -- reducing it by several more orders of magnitude than that did it for me. (I'm sure there's still a more elegant way of doing it than this, though!)  
> **username_3**: @HDE226868 [Adaptive quadrature](https://en.wikipedia.org/wiki/Adaptive_quadrature), probably, although of course even that won't help if the range of integration is large enough and the grid sparse enough to completely miss any interesting features in the function. Personally, for something like this which spans several orders of magnitude, I might try a logarithmic change of variables ($u = \ln(\lambda/L)$ for some arbitrarily chosen $L$) to hopefully increase the chance that some grid points land in the right spots.  
> **username_2**: @DavidZ Yeah, logarithms might be the best way to go about it.  
> **username_4**: @HDE226868 Could you please find the mistake here https://astronomy.stackexchange.com/questions/57307/calculating-suns-j-band-mag-via-integrating-the-planck-law ?  

---

### Answer by username_1 on 2023-11-27. Score: 3

The integral of the Planck function over all wavelengths gives the surface flux at the photosphere of the Sun divided by pi.

In order to get the solar constant, you should multiply this by $\pi \times 4\pi R\_\odot^2$ and then divide by $4\pi d^2$, where $d$ is an astronomical unit in metres. This looks correct.

The integration limits should be set to something sensible. The Sun isn't a blackbody, though the effective temperature should reasonably represent the total flux. I would use a lower limit of 100 nm and an upper limit of 5000 nm.


---

