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()