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')
    plt.ylabel('Spectral radiance, W/m2 nm sr')

    print('###### At the Sun Surface')
    ### Numerical integral in energy direction
    integral_radiance = (spectral_radiance * dwl).sum()
    ### In solid angle direction
    integral_radiance = np.pi * integral_radiance
    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
    total_power = integral_radiance * sun_area
    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()