vima_demo.plot_af_gfactor

A code to plot Andrei’s G-factor.

Code author: Moa Persson

""" A code to plot Andrei's G-factor.

.. codeauthor:: Moa Persson

.. image:: ../../../src/scripts/vima_demo/gf_Hp_pac6.png

.. image:: ../../../src/scripts/vima_demo/gf_Op_pac6.png
"""
import numpy as np
import matplotlib.pyplot as plt
import irfpy.vima.energy
import irfpy.vima.massring as mring
import irfpy.vima.gfactor as gfactor

def main():

    # Import tables
    etbl = irfpy.vima.energy.get_default_table_v3()
    species = [1, 16]  # H+ and O+
    name = ['Hp', 'Op']
    paccidx = 6

    for midx, mass in enumerate(species):

        fig, ax = plt.subplots(1, 1, figsize=[12, 10])

        # Plot Andrei Fedorov G-factor for the different cases of azimuth sector and elevation angle
        for elev in range(-45, 45, 5):
            for sect in range(16):
                Rm = [mring.rm(e_per_q=energy, m_per_q=mass, pi=paccidx) for energy in etbl]
                gfac = [gfactor.gf_af(energy, paccidx, rm, elev, sect) for (energy, rm) in zip(etbl, Rm)]
                # lab = 'AF, P: {}, S: {}'.format(name[midx], elev, sect)
                # plt.plot(etbl, gfac, label=lab)
                plt.plot(etbl, gfac)

        # Plot Marcus Franz G-factor
        g = np.array([gfactor.gf_mf(ee, midx) for ee in etbl])
        g *= 0.292
        plt.plot(etbl, g, 'ko-', lw=3, label='MF'.format(name[midx]))

        plt.title('Species: {}, AF g-fac for Elev -45 to 45 degrees, Azim sect 1-16, Pacc: {}'.format(name[midx], paccidx), fontsize=20)
        ax.tick_params(axis='both', direction='out', labelsize=18)

        # plt.xlim([1, 1e6])
        plt.xscale('log')
        plt.yscale('log')
        plt.xlabel('Energy [eV]', fontsize=18)
        plt.ylabel('G-factor', fontsize=18)
        plt.legend()
        plt.savefig('gf_{}_pac{}.png'.format(name[midx], paccidx))

#    plt.show()

if __name__ == '__main__':
    main()

# # Plot Marcus Franz G-factor
#     g1 = [gfactor.gf_mf(ee, 0) for ee in etbl]
#     g16 = [gfactor.gf_mf(ee, 1) for ee in etbl]
#     plt.plot(etbl, g1, 'o-', label='MF: Hp')
#     plt.plot(etbl, g16, 'x-', label='MF: Op')