appl121105_iotorus.neutralmodel

''' Neutral model.

An wrapper class for neutral model.
'''
import numpy as np

from pyana.pep import iotorus

from . import moon_io

class NeutralModelTwoPeak2:
    ''' A wrapping of neutral cloud model

    High performance version of :class:`NeutralModelTwoPeak`.

    >>> io = moon_io.Io()   # t=0, phase=0
    >>> nm = NeutralModelTwoPeak2(io)

    >>> print('%.2f' % nm.get_density(420000, 0, 0))
    10623.05

    If you want to know the different time, you
    can use ``t`` keyword. For example, after/before 3 hours,

    >>> n = nm.get_density([420000, 420000], [0, 0], [0, 0], tarr=[10800, -10800])
    >>> print('%.2f %.2f' % (n[0], n[1]))
    17.43 1849.77

    >>> n = nm.get_density([380000, 380000], [-175000, 175000], [0, 0])
    >>> print('%.2f %.2f' % (n[0], n[1]))
    16.53 1421.22

    '''
    def __init__(self, io):
        tpm = iotorus.TwoPeakModel()
        self.model = iotorus.RotatingNeutralCloud2(tpm)
        self.io = io

    def get_density(self, xarr, yarr, zarr, tarr=None):
        ''' Get density in cm^-3

        :param x: x coordinates in km/s.
        :returns: Density at (x, y, z) in cm^-3.
        '''
        rj = 71492.
        xrj = np.array(xarr) / rj
        yrj = np.array(yarr) / rj
        zrj = np.array(zarr) / rj

        if tarr == None:
            tarr = np.zeros_like(xrj)

        phase = np.deg2rad(self.io.get_phase(tarr))
        return self.model.get_density(xrj, yrj, zrj, iophase=phase)


class NeutralModelTwoPeak:
    ''' A wrapping (adapter) of neutral cloud model.

    Adapting :class:`pyana.pep.iotorus.RotatingNeutralCloud`
    with :class:`pyana.pep.iotorus.TwoPeakModel`.

    First, :class:`moon_io.Io` is needed to instance
    NeutralModel.  It is because neutral model depends on
    the Io position.

    >>> io = moon_io.Io()   # t=0, phase=0
    >>> nm = NeutralModelTwoPeak(io)

    Now you can get the density.
    The distance is approximately the Io position.

    >>> print('%.2f' % nm.get_density(420000, 0, 0))
    10623.05

    If you want to know the different time, you
    can use ``t`` keyword. For example, after/before 3 hours,

    >>> print('%.2f' % nm.get_density(420000, 0, 0, t=10800))
    17.43
    >>> print('%.2f' % nm.get_density(420000, 0, 0, t=-10800))
    1849.77

    After 3 hours, the position is behind Io.
    Before 3 hours, the position is in front of Io.
    The banana extends preferably toward the Io,
    the density profile looks ok.

    Also, you can check this values as follows.
    3 hours corresponds to 0.43 radiand (roughly).
    Thus, the corresponding position for t=10800 is
    (x, y, z) = (420000 * cos(0.43), -420000 * sin(0.43), 0)
    and t=-10800 is
    (x, y, z) = (420000 * cos(0.43), 420000 * sin(0.43), 0)

    >>> print('%.2f' % nm.get_density(380000, -175000, 0))
    16.53
    >>> print('%.2f' % nm.get_density(380000, 175000, 0))
    1421.22

    '''
    def __init__(self, io, rjlimit=10):
        tpm = iotorus.TwoPeakModel()
        self.model = iotorus.RotatingNeutralCloud(tpm)
        self.io = io
        self.rjlimit2 = rjlimit * rjlimit

    def get_density(self, x, y, z, t=0):
        ''' Get density in cm^-3

        :param x: x coordinates in km/s.
        :returns: Density at (x, y, z) in cm^-3.
        '''
        rj = 71492.
        xrj = np.array(x) / rj
        yrj = np.array(y) / rj
        zrj = np.array(z) / rj

        if xrj * xrj + yrj * yrj + zrj * zrj > self.rjlimit2:
            return 0

        phase = np.deg2rad(self.io.get_phase(t))
        self.model.set_iophase(phase)
        return self.model.get_density(xrj, yrj, zrj)

import unittest
import doctest


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