<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
Step1:  Day 3 Planck problem
Verify using np.sum that
Step2:  so good agreement with 10000 points -- make a plot as well
<ASSISTANT_TASK:>
Python Code:
import numpy as np
from matplotlib import pyplot as plt
#
# get Stull's c_1 and c_2 from fundamental constants
#
c=2.99792458e+08  #m/s -- speed of light in vacuum
h=6.62606876e-34  #J s  -- Planck's constant
kb=1.3806503e-23  # J/K  -- Boltzman's constant
c=3.e8  #speed of light in vacuum (m/s)
c1=2.*h*c**2.*np.pi
c2=h*c/kb
sigma=2.*np.pi**5.*kb**4./(15*h**3.*c**2.)
def planckwavelen(wavel,Temp):
    
    Calculate the blackbody radiant exitence (Stull 2.13)
    
    Parameters
    ----------
    
      wavel: float or array
           wavelength (meters)
           
      Temp: float
           temperature (K)
           
    Returns
    -------
    
    Elambda:  float or arr
           monochromatic radiant exitence (W/m^2/m)
    
    Elambda=c1/(wavel**5.*(np.exp(c2/(wavel*Temp)) -1))
    return Elambda
Temp = 300.
SBflux = sigma*Temp**4.
print('stefan-boltzman says flux is {:8.3f} W/m^2'.format(SBflux))
#
# get wavelength in meters between 0.1 and 500 microns
#
npoints = 10000
wavelengths=np.linspace(0.1,500.,npoints)*1.e-6
flux = planckwavelen(wavelengths,Temp)
avg_flux = (flux[1:] + flux[:-1])/2.
integral = np.sum(avg_flux*np.diff(wavelengths))
print('numerically integrated flux is: {:8.3f} W/m^2'.format(integral))
Explanation: Day 3 Planck problem
Verify using np.sum that:
$$\int_0^\infty  E^*\lambda(T)d\lambda   =\int_0^\infty \frac{\pi h c^2}{\lambda^5 \left [ \exp (h c/(\lambda k_B T )) -1 \right ] } d\lambda = \int_0^\infty \frac{c_1}{\lambda^5 \left [ \exp (c_2/(\lambda T )) -1 \right ] } d\lambda  = \sigma{SB} T^4$$
For a temperature of 300 K.
Answer -- check for Temp=300 K
End of explanation
%matplotlib inline
fig,ax = plt.subplots(1,1)
#
# change units from meters to microns to get prettier numbers
#
ax.semilogx(wavelengths*1.e6,flux*1.e-6)
_=ax.set(xlim=(1,100),xlabel="wavelength ($\mu m$)",ylabel=r"$E^*_\lambda\ (W/m^2/\mu m$)",
        title='Radiant exitence for Temp = 300 K')
Explanation: so good agreement with 10000 points -- make a plot as well
End of explanation
<END_TASK>