cena_f2c_sample

A sample code to convert flux 2 count rate for CENA

''' A sample code to convert flux 2 count rate for CENA
'''

import os
import sys
import logging
logging.basicConfig()
import datetime
import math

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp

import irfpy.cena.cena_flux
import irfpy.cena.energy

import irfpy.cena.empirical

def main():
    ''' There are several ways of getting "counts" from flux.

    Below shows the sample.
    '''
    # To plot the original differential flux.
    fig = plt.figure()
    # ax1 is for differential flxu domain
    ax1 = fig.add_subplot(121)
    # ax2 is for count rate domain
    ax2 = fig.add_subplot(122)


    # For dense E axis
    e0 = np.logspace(0, 4, num=1000)  # 1eV to 100 keV

    # First of all, make the flux to be observed.

    # Assuming Maxwell, T=100 eV. n=10 /cm3.
    # The complete expression is
    #       2 n / sqrt(m) (2 pi k T) ^ (-3/2) E exp (-E/kT).
    # Roughly 1e2 * E * (-E / 100)  where E is in eV and
    # J is in cm-2 sr-1 eV-1 s-1.
    temp = 100.
#    temp = 120.
#    temp = 60.
    j = lambda E: 1e2 * E * np.exp(-E/temp)

    # Simplest j for debugging is to set all 1.
#    j = lambda E: E ** 0

    p0 = ax1.plot(e0, j(e0))


    # Ok, now just calculate the count rate in the simplest way.
    f2c_simple = irfpy.cena.cena_flux.Flux2Count_Simple()
    # To get the count rate, the flux at the energy step should be specified.
    sect = 3
    estep = irfpy.cena.energy.getEnergyE16()
    flux = j(estep)
    count_simple = f2c_simple.get_counts(flux, 3)
    p1 = ax2.plot(estep, count_simple, '+-', label='Simple')

    # Then, Matrix way.
    f2c_matrix = irfpy.cena.cena_flux.Flux2Count_Matrix()
    count_matrix = f2c_matrix.get_counts(flux, 3)
    p2 = ax2.plot(estep, count_matrix, 'o-', label='Matrix')

    # Then, integral way.
    type='simulated'
#    type='simple_triangle'
#    type='delta'
    f2c_integral = irfpy.cena.cena_flux.Flux2Count_Integral(energy_response_type=type)
    count_integral = [f2c_integral.getCount(j, ie, 3) for ie in range(16)]
    p3 = ax2.plot(estep, count_integral, 'v-', label='Integral')

    # The counts are considered as observation.
    # Then, they are used simple c->j model and fit them.
    # Now you may get the parameter estimation.


    m_simple = irfpy.cena.empirical.data_maxwellfit(count_simple, 3)
    m_matrix = irfpy.cena.empirical.data_maxwellfit(count_matrix, 3)
    m_integral = irfpy.cena.empirical.data_maxwellfit(count_integral, 3)
    print('Simple16', m_simple[0], irfpy.cena.empirical.flux_half_maxwell(m_simple[0][0], m_simple[0][1]))
    print('Matrix16', m_matrix[0], irfpy.cena.empirical.flux_half_maxwell(m_matrix[0][0], m_matrix[0][1]))
    print('Integral16', m_integral[0], irfpy.cena.empirical.flux_half_maxwell(m_integral[0][0], m_integral[0][1]))

#    mask = [1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1,]
#    m2_simple = irfpy.cena.empirical.data_maxwellfit(np.ma.masked_array(count_simple, mask), 3)
#    m2_matrix = irfpy.cena.empirical.data_maxwellfit(np.ma.masked_array(count_matrix, mask), 3)
#    m2_integral = irfpy.cena.empirical.data_maxwellfit(np.ma.masked_array(count_integral, mask), 3)
#    print 'Simple8', m2_simple[0], irfpy.cena.empirical.flux_half_maxwell(m2_simple[0][0], m2_simple[0][1])
#    print 'Matrix8', m2_matrix[0], irfpy.cena.empirical.flux_half_maxwell(m2_matrix[0][0], m2_matrix[0][1])
#    print 'Integral8', m2_integral[0], irfpy.cena.empirical.flux_half_maxwell(m2_integral[0][0], m2_integral[0][1])

    ax1.set_yscale('log')
    ax1.set_xscale('log')
    ax1.set_ylim(1e-2, 1e4)

    ax2.set_yscale('log')
    ax2.set_xscale('log')
    ax2.set_ylim(1e-5, 1e1)

    return j(e0)


if __name__ == '__main__':
    main()