cena_gfactor_ch3

import matplotlib
matplotlib.use('Agg')

from pylab import *
from irfpy.cena import energy
from irfpy.cena import gfactor

if __name__ == '__main__':
    energy_in_eV = energy.getEnergyE16()

    # Calculate the g-factor based on the E-dep g-factor from
    # Table 3.5 the eq. 3.19
    gf0 = gfactor.GFactorH_Table()
    g0 = gf0.geometricFactor()
    g0_ch3 = g0[:,3]

    # Calculate the g-factor based on the E-dep g-factor by
    # aggregation of parameters and the eq. 3.19
    gf1 = gfactor.GFactorH_Component()
    g1 = gf1.geometricFactor()
    g1_ch3 = g1[:,3]

    # G-factor of the CH-3 listed in Table 3.5
    g2 = np.array([-10, -10, -10, -10, 2.7, 3.1, 3.7, 4.5,
                5.8, 7.7, 11, 15, -10, -10, -10, -10]) * 1e-7
    g2_ch3 = ma.masked_less(g2, 0)
    g2_ch3

    # Calculate the difference between g0 and g2.

    diff = 100 * (g0_ch3 - g2_ch3) / g2_ch3

    figure()

    plot(energy_in_eV, g0_ch3, 'o', label='E-dep g-fact (Table 3.5) + Eq 3.19')
    plot(energy_in_eV, g1_ch3, label='E-dep g-fact (aggregation) + Eq 3.19')
    plot(energy_in_eV, g2_ch3, label='Table 3.5, last column')

    for i in range(4, 12):
        print(i, g0_ch3[i], diff[i])
        text(energy_in_eV[i], g0_ch3[i], '%.1f%% ' %
                diff[i], horizontalalignment='right', rotation=-10)

    yscale('log')
    xscale('log')
    xlabel('Energy [eV]')
    ylabel('G-factor of CH-3 [cm2 sr eV/eV]')
    legend(loc='upper left')

    savefig('cena_gfactor_ch3.png')