# snippet.sun_spectrum_earth_scalingΒΆ

Exercise for solar spectrum.

The aim is to understand the solar spectrum from black body assupmption.

''' Exercise for solar spectrum.

The aim is to understand the solar spectrum from black body assupmption.
'''
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

import irfpy.util.constant as k
import irfpy.util.planck as pl

def main():
'''Main script'''
temperature = 5800.
bb = pl.BlackBodyMksa(temperature)   # A black body.

wavelens = np.linspace(200, 4000, 3801) * 1e-9  # in m
dwl = wavelens[1] - wavelens[0]  # in m

spectral_radiance = bb.Blambda(wavelens) # in W/m2 m sr
plt.plot(wavelens / 1e-9, spectral_radiance * 1e-9, label='SUN')
plt.title('Black Body T=%g' % temperature)
plt.xlabel('Wavelength, nm')

print('###### At the Sun Surface')
### Numerical integral in energy direction
### In solid angle direction
print('Integral Raidance (numerical) = %e W/m2' % integral_radiance)

### The integral_radiance is the total power emitted from unit area.
sun_area = 6.0877e18   # in m2
print('total power = %e W' % total_power)

print('###### At the Earth Orbit')
### All the power that Earth will get is the fraction
### (pi * Re ** 2) / (4 * pi * AU ** 2)
re = 6.374e6   # in m
au = 1.496e11
fraction = re * re / (4 * au * au)
print('Fraction earth get = ', fraction)
power_at_earth = total_power * fraction
print('Radiation Power at earth = %e W' % power_at_earth)

### IF albedo is 0, what is the black-body earth temperature?
power_at_earth_m2 = power_at_earth / (4 * np.pi * re * re)
print('Power from earth per unit = %e W/m2' % power_at_earth_m2)
t = np.power(power_at_earth_m2 / k.stefan_boltzmann_constant, 0.25)
print('Black-body earth temperature = %f K' % t)

### Power at Earth's orbit
power_per_unit_are_at_1AU = total_power / (4 * np.pi * au * au)
print('Power at 1 AU per unit area = %e W/m2' % power_per_unit_are_at_1AU)
print('  This is the solar constant, 1.37 kW/m2 or 1.37 W/cm2')
print('  Fraction to get is 1 / 4pi AU^2 = %e' % (1 / (4 * np.pi * au * au)))

print('###### At the Earth Orbit for 500 nm spectrum irradiance')
### Spectrum calculation
spectrum500nm = bb.Blambda(500e-9) * np.pi   # np.pi is integral over solid angle
print('u(500nm) = %e W / m2 / m' % spectrum500nm)
spectrum500nm = spectrum500nm * 1e-9
print('         = %e W / m2 / nm' % spectrum500nm)

# Integral over all the Sun area
total_spectrum500nm = spectrum500nm * sun_area
print('F(500nm, emission from Sun) = %e W / nm' % total_spectrum500nm)
spectrum500nm_per_unitarea_1AU = total_spectrum500nm / (4 * np.pi * au * au)
print("u'(500nm, at 1AU) = %e W / m2 / nm" % spectrum500nm_per_unitarea_1AU)

### Number of photons.
energy_500nm = k.h * k.c / 500e-9
print('500 nm photon energy = %e J' % energy_500nm)
print('Photon Flux = %e # / m2 / nm / s' % (
spectrum500nm_per_unitarea_1AU / energy_500nm))

if __name__ == "__main__":
main()