Source code for irfpy.util.kappa

''' Kappa distribution.

.. codeauthor:: Yoshifumi Futaana
'''
from irfpy.util import constant as k

import numpy as np
from scipy.special import gamma

class kappa_distribution:
    r''' Kappa distribution expressed by

    .. math::

        f_i^\kappa(r,v) = \frac{n_i}{2\pi(\kappa w_{\kappa i}^2)^{3/2}}
                \frac{\Gamma(\kappa+1)}{\Gamma(\kappa-1/2)\Gamma(3/2)}
                (1+\frac{v^2}{\kappa w_{\kappa i}^2})^{-(\kappa+1)}

    Here :math:`w_{\kappa i}^2=(2\kappa-3)kT_i/\kappa mi`

    In this class, use of MKSA unit system is implicitly assumed.

    >>> fk2 = kappa_distribution(2)    # Returned is the function.  f(v, k=2)
    >>> fk2_at0 = fk2(0)     # f(0, k=2) is returned.

    .. warning::

        Not yet fully tested.

    .. todo::

        Validate the formulation.
    '''

    def __init__(self, kappa, ni=1, Ti=11604., mi=k.m_proton):
        if kappa <= 1.5:
            raise ValueError('Given kappa={} not acceptable (>1.5)'.format(kappa))
        self.kappa = kappa
        self.wki2 = (2.0 - 3.0 / kappa) * k.kB * Ti / mi
#        print(self.wki2)
        t1 = ni / (2 * np.pi * np.power(kappa * self.wki2, 1.5))
        t2 = gamma(kappa + 1) / gamma(kappa - 0.5) / gamma(1.5)
        self.coeff = t1 * t2

    def __call__(self, xarr):
        kappa = self.kappa
        t3 = np.power((1 + xarr ** 2 / (kappa * self.wki2)), -kappa - 1)
        return self.coeff * t3
        
import unittest
import doctest


def doctests():
    return unittest.TestSuite((
        doctest.DocTestSuite(),
        ))
if __name__ == '__main__':
    unittest.main(defaultTest='doctests')