snippet.solar_photoionization

Calculate the total flux of the photon flow

The photoionization occurs by the photons with energy more than the potential eneryg of the target.

For hydrogen atom, it is 13.598 eV and for oxygen, it is 13.618 eV as well. For sodium atom, it is 10.360 eV. See http://www5f.biglobe.ne.jp/~rokky/kaisetu/0/syuukihyou_pdf03.pdf for summary.

Now calculate the photon flux that is more than the energy above.

'''Calculate the total flux of the photon flow

The photoionization occurs by the photons with energy more than
the potential eneryg of the target.

For hydrogen atom, it is 13.598 eV and for oxygen, it is 13.618 eV as well.
For sodium atom, it is 10.360 eV.
See http://www5f.biglobe.ne.jp/~rokky/kaisetu/0/syuukihyou_pdf03.pdf for summary.

Now calculate the photon flux that is more than the energy above.
'''
import numpy as np
import matplotlib.pyplot as plt

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

def calculate_high_energy_photon_flux(solar_model, potential_energy_ev):

    wavelength_threshold = pl.energy2wavlen(potential_energy_ev * k.qe)
    irr_threshold = solar_model.spectral_irradiance(wavelength_threshold)
    fph_threshold = solar_model.photon_flux(wavelength_threshold)

    print('Wlen = %e m, Irr = %e W /m2 /m /s, Fph = %e photons /m2 /m /s' % (wavelength_threshold, irr_threshold, fph_threshold))

    ndiv = 1000

    wavelens = np.linspace(1e-30, wavelength_threshold, ndiv)
    dwavelen = wavelens[1] - wavelens[0]

    fph_spectra = solar_model.photon_flux(wavelens)
    fph_integral = (fph_spectra * dwavelen).sum()

    return fph_integral, wavelens, fph_spectra


def main():
    '''Main script'''
    ### Instance the solar spectra model.

    solar_models = [ssp.BlackBody(),   # Black body
                    ssp.BlackBodyModFuv(),
#                    ssp.AM0(),     # AM0 model
                    ]
    color_models = ['b', 'g', ]

    potential = {'H': 13.598, 'O': 13.618, 'S': 10.360, }

    ypos = 1e14
    for solar, color in zip(solar_models, color_models):
        ### For plotting
        f, w, fs = calculate_high_energy_photon_flux(solar, 5.0)
        plt.plot(w * 1e9, fs / 1e9, color, label=solar.name)

        for species in list(potential.keys()):
            p = potential[species]
            print(species, p)
            f, l, fph = calculate_high_energy_photon_flux(solar, p)
            print('%.3e' % f)

            wavlen = pl.energy2wavlen(p * k.qe) * 1e9  # in nm
            plt.plot([wavlen], [ypos], color + 'o')
            plt.text(wavlen, ypos, species + ': %.2e /m2 s' % f)
            ypos = ypos * 2

        

    plt.xlabel('Wavelength [nm]')
    plt.ylabel('Photon flux [/m2 nm s]')
    plt.yscale('log')
    plt.ylim(1e11, 1e17)
    plt.legend(loc='best')
    plt.savefig('solar_photoionization.png')



if __name__ == "__main__":
    main()