Source code for irfpy.swim.flux

''' A collection of generic functions to convert from count rate to flux.
'''

from numpy import array
import irfpy.util.maxwell as maxwell
import math
import logging

MP_KG=1.67262158e-27
EV2JOULE=1.60217646e-19

def simple_c2j(count_, dutytime_s, gfactor_cm2sreV_eV, efficiency_, energy_eV):
    ''' A simple converter from counts observed to differential flux.  Proton is assumed.

    :param count_:   Observed count [particles/sample].  It can be a count rate [particles/sec], in this case dutytime_s should be set to 1.
    :param dutytime_s:   Duty time of the observation [sec/sample]
    :param gfactor_cm2sreV_eV:   G-factor [cm^2 sr eV/eV].  Efficiency can be multipled in this variable, in this case *efficiency_* should be set to 1.
    :param efficiency_:    Efficiency [unitless].  G-factor can include efficiency (multiply). In that case, specify one for *efficiency_*.
    :param energy_eV:   Energy of the particles [eV].
    :returns: Differential flux [particles / cm^2 sr eV s]
    '''

    J = count_ / dutytime_s / gfactor_cm2sreV_eV / efficiency_ / energy_eV

    return J

def simple_j2f(J__cm2sreVs, energy_eV):
    ''' A simple converter from the differential flux [#/cm^2 sr eV s] to the velocity distribution functions [s^3 / m^6].  Proton is assumed.

    :param J__cm2sreVs:    Differential flux [#/cm^2 sr eV s]
    :param energy_eV:    Energy [eV]
    :return: Velocity distribution function [s^3 / m^6]
    '''
    m_kg = MP_KG
    energy_J = energy_eV * EV2JOULE
    J__m2srJs = J__cm2sreVs / 1e-4 / EV2JOULE
    f = m_kg**2 / 2 / energy_J * J__m2srJs
    return f

def simple_moment(velocitytbl_km_s, fslice_s3_m6):
    ''' From a velocity distribution function, the density [/m^3]is simply calculated.

    Distribution function of Maxwell distribution is expressed by
    f(v) = n0 ( 1 / 2 pi vth^2 ) ^ 3/2  exp [ - (v-v0)^2 / 2 vth^2 ].
    If you can know max(f(v)) and vth, the density is calculated as
    '''

    vtable_m_s = array(velocitytbl_km_s).copy()*1e3
    logging.debug('V=%s'%str(vtable_m_s))

    ## Now fitting.
    (a,b,c), s = maxwell.fit1d(vtable_m_s, fslice_s3_m6)

    logging.debug('a,b,c=%g,%g,%g'%(a,b,c))

    # The result is a * exp ( - (v-b)^2 / 2c^2)
    nsw = a * (2 * math.pi * c**2)**(1.5)
    vsw = b/1e3
    vth = c/1e3

    return nsw, vsw, vth