snippet.solar_photon

Plot the sun photon spectum at the Earth orbit.

'''Plot the sun photon spectum at the Earth orbit.

.. image:: ../../../src/scripts/snippet/solar_photon.png
'''
import numpy as np
import matplotlib.pyplot as plt
import irfpy.util.solarspectra
import irfpy.util.planck as pl
import irfpy.util.constant as k

def main():
    '''Main script'''
    wavelength = np.linspace(200, 4000, 3801) * 1e-9  # Wavelength to examine
    dwl = wavelength[1] - wavelength[0]

    bb = irfpy.util.solarspectra.BlackBody(temperature=6050)
    plt.figure()

    ### Plot number of photons
    plt.subplot(211)

    fph = bb.photon_flux(wavelength)   # in /m2 /m /s
    totfph = (fph * dwl).sum()
    print('Total flux = %e /m2 /s' % (fph * dwl).sum())
    
    plt.plot(wavelength / 1e-9, fph * 1e-9, 'b-', label='Black body 6050 K (TOT=%.2e)' % totfph)   # For converting nm and /m2/nm /s
    plt.xlabel('Wavelength [nm]')
    plt.ylabel('Photon flux [1/m^2/nm/s]')
    plt.legend()

    ### In the energy domain
    wavelength0 = wavelength - dwl / 2.
    wavelength1 = wavelength + dwl / 2.

    energy_eV = pl.wavlen2energy(wavelength) / k.qe  # Center energy
    energy_eV0 = pl.wavlen2energy(wavelength1) / k.qe  # Lower energy bound
    energy_eV1 = pl.wavlen2energy(wavelength0) / k.qe  # Upper energy bound
    denergy = energy_eV1 - energy_eV0  # Energy pass band.
    fph_eV = fph * dwl / denergy  # #/m^2/eV/s

    plt.subplot(212)
    totfph_eV = (fph_eV * denergy).sum()
    print('Total flux = %e /m2 /s' % totfph_eV)
    
    plt.plot(energy_eV, fph_eV, 'b-', label='Black body 6050 K (TOT=%.2e)' % totfph)   # For converting nm and /m2/nm /s
    plt.xlabel('Energy [eV]')
    plt.ylabel('Photon flux [1/m^2/eV/s]')
    plt.legend()


    plt.savefig('solar_photon.png')


if __name__ == "__main__":
    main()