Integration of Plank curve of the Sun

Question

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)

Answer 1

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

Answer 2

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