appl121105_iotorus.frames

Provide frame information.

Internally, most of the calculation is done in JEqS system. But in some cases, you may need System-III, for example for the Io plasma torus model.

''' Provide frame information.

Internally, most of the calculation is done in JEqS system.
But in some cases, you may need System-III, for example
for the Io plasma torus model.
'''
from . import parameters as p
import numpy as np

from irfpy.util import ringcache


__sys3tojeqs_cache = ringcache.SimpleRingCache(1000000)   # Will need a few 10s MB.

def system3_to_jeqs(sun_longitude_s3_deg):
    ''' Return 3x3 matrix converting system3 to jeqs.

    >>> matx = system3_to_jeqs(0)
    >>> print(matx[0, 0] == matx[1, 1] == matx[2, 2] == 1.)
    True
    >>> print(matx[0, 1] == matx[0, 2] == matx[1, 0] == 0.)
    True
    >>> print(matx[2, 0] == matx[1, 2] == matx[2, 1] == 0.)
    True

    If System-III longitude of sun (subsolar point) is 60 degrees,

    >>> matx = system3_to_jeqs(60.)

    The z axis is not changed.

    >>> print(matx.dot([0, 0, 1]))
    [ 0.  0.  1.]

    The x axis of JEqS is the sun direction, so that the
    System-III longitude is 60 degrees. In other words,
    (1, 0, 0) in JEqS is (cos 60, sin 60, 0) in System III.
    matx provides the conversion from System III to JEqS,
    the opposite will be provided its invert (identical to
    transpose in this case).

    >>> print(matx.T.dot([1, 0, 0]))
    [ 0.5        0.8660254  0.       ]
    '''

    if __sys3tojeqs_cache.hasKey(sun_longitude_s3_deg):
        return __sys3tojeqs_cache.get(sun_longitude_s3_deg)

    sunlon = np.deg2rad(sun_longitude_s3_deg)
    c = np.cos(sunlon)
    s = np.sin(sunlon)

    matx = np.array([[c, s, 0],
                [-s, c, 0],
                [0, 0, 1]])

    __sys3tojeqs_cache.add(sun_longitude_s3_deg, matx)

    return matx


__sys3matcache = {}

def system3_to_magnetic(maglon=p.maglon, maglat=p.maglat):
    ''' Return 3x3 matrix converting the System III to magnetic 

    Magnetic frame is defined as:
    z-axis along the magnetic dipole,
    x-axis as being in the System III x-z plane.

    :keyword maglon: Magnetic longitude (east). Default is :attr:`parameters.maglon`.
    :keyword maglat: Magnetic latitude (north). Default is :attr:`parameters.maglat`.

    z-axis in magnetic coordinate should be converted to
    the magnetic dipole vector in System3 (=[cos(mlat)*cos(mlon),
    cos(mlat)*sin(mlon), sin(mlat)], see :func:`magnetic_dipole_in_system3`).

    >>> matx = system3_to_magnetic()
    >>> print(matx.T.dot([0, 0, 1]))
    [-0.15495027  0.0616622   0.98599604]

    x-axis in magnetic coordinate should be converted to
    the vector in system 3 whose y component is 0.

    >>> print(matx.T.dot([1, 0, 0]))
    [ 0.98599604  0.          0.15495027]

    Check if the obtaind matrix is rotational matrix.

    >>> print(np.round(matx.dot(matx.T), 2))
    [[ 1.  0.  0.]
     [ 0.  1.  0.]
     [ 0.  0.  1.]]
    '''
    if not (maglon, maglat) in __sys3matcache:
        clon = np.cos(np.deg2rad(maglon))
        slon = np.sin(np.deg2rad(maglon))
        clat = np.cos(np.deg2rad(maglat))
        slat = np.sin(np.deg2rad(maglat))

        z = [clat * clon, clat * slon, slat]
        x = [slat, 0, -clat * clon]
        y = np.cross(z, x)

        matx = np.array([[x[0], x[1], x[2]], [y[0],  y[1], y[2]], [z[0], z[1], z[2]]])

        __sys3matcache[(maglon, maglat)] = matx

    matx = __sys3matcache[(maglon, maglat)]

    return matx


__jeqs_to_magnetic_cache = ringcache.SimpleRingCache(100000)

def jeqs_to_magnetic(sun_longitude_s3_deg, maglon=p.maglon, maglat=p.maglat):
    key = (sun_longitude_s3_deg, maglon, maglat)
    if __jeqs_to_magnetic_cache.hasKey(key):
        return __jeqs_to_magnetic_cache.get(key)

    jeqs2sys3 = system3_to_jeqs(sun_longitude_s3_deg).T
    sys32mag = system3_to_magnetic(maglon=maglon, maglat=maglat)
    jeqs2mag = sys32mag.dot(jeqs2sys3)
    __jeqs_to_magnetic_cache.add(key, jeqs2mag)

    return jeqs2mag

def magnetic_dipole_in_system3(maglon=p.maglon, maglat=p.maglat):
    mlat = np.deg2rad(maglat)
    mlon = np.deg2rad(maglon)
    return np.array([np.cos(mlat) * np.cos(mlon),
                np.cos(mlat) * np.sin(mlon),
                np.sin(mlat)])



import unittest
import doctest


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