Source code for irfpy.util.lorentztrans

''' Functions for Lorenz transformation

.. codeauthor:: Yoshifumi Futaana

.. todo::

    Consider importing to physics module
'''
import numpy as np

from irfpy.util import constant as k
from irfpy.util import unitq as u

c = k.c
c_u = u.c.simplified

def beta(v):
    r''' Return :math:`\beta`.

    Return the parameter :math:`\beta`.

    .. math::

        \beta = \frac{v}{c}

    where :math:`v` is the velocity of the object
    and :math:`c` is the velocity of light.

    :param v: Velocity (in meter) of the object.
    :type v: float or np.array.

    >>> print('%.5f' % beta(3e7))  # About 0.1
    0.10007
    '''
    return v / c


beta_u = u.make_unitful_function(beta, [u.m], 1)
''' Unitful version of :func:`beta`.
'''


def gamma(v):
    r''' Return :math:`\gamma`.

    Return the parameter :math:`\gamma`.

    .. math::

        \gamma = \frac{1}{\sqrt{1 - \bigl(\frac{v}{c}\bigl)^2}}

    where :math:`v` is the velocity of the object
    and :math:`c` is the velocity of light.

    :param v: Velocity of the object.
    :type v: float or np.array.

    >>> print('%.6f' % gamma(3e7))
    1.005045
    '''
    p = 1 + v / c
    m = 1 - v / c
    return 1 / np.sqrt(p * m)


gamma_u = u.make_unitful_function(gamma, [u.m], 1)
''' Unitful version of :func:`gamma`.
'''

def gamma2vel(gamma):
    r''' Inverse function of :func:`gamma`.

    :param gamma: The parameter :math:`\gamma`.
    :returns: :math:`v` in ``m/s``.

    The relation is
    
    .. math::

        \gamma = \frac{1}{\sqrt{1 - \bigl(\frac{v}{c}\bigl)^2}}

    so that

    .. math::

        v = c \sqrt{1 - \frac{1}{\gamma^2}}

    >>> print('%.3e' % gamma2vel(1.01))
    4.208e+07

    You may use unitful version:

    >>> print('%.3e' % gamma2vel_u(1.01).rescale(u.km / u.s))
    4.208e+04
    '''
    return np.sqrt(1 - 1 / (gamma ** 2)) * c


gamma2vel_u = u.make_unitful_function(gamma2vel, [], u.m / u.s, offset=1)


def transformx(x, t, v):
    r""" Lorenz transform along x.

    An inertia frame, *K(x, y, z)*, is Lorentz-transformed to
    *K'(x', y', z')*, another inertia frame.
    *K'* is moving relative to *K* with velocity *v*
    along *x* axis.

    The formulation is

    .. math::

        x' &= \gamma(x-vt)    \\
        t' &= \gamma(t-\frac{v}{c^2}x)

    :param x: The *x* coordinate in K frame.
    :param t: The *t* in K frame.
    :param v: The velocity of K' relative to K.
    :return: A tuple, (x', t').

    A simple example for x=3e8, t=10, v=2e8.
    If you consider Galileo transform, the time
    does **not** change (t=10) and the position
    is :math:`x-vt=3e8 - 10\times 2e8=-17e9`.
    However, if you use the Lorentz transform,
    these values differ.

    >>> x1, t1 = transformx(3e8, 10, 2e8)
    >>> print('%.3e' % x1)
    -2.282e+09
    >>> print('%.3f' % t1)
    12.528

    Unitful version is also supported.
    Numpy array argument is also supported.

    >>> x1u, t1u = transformx_u(np.array([3e8, 3e8]) * u.m, np.array([5, 10]) * u.s, np.array([2e8, 2e8]) * u.m/u.s)
    >>> print(x1u)   # doctest: +NORMALIZE_WHITESPACE
    [-9.39669291e+08  -2.28205399e+09] m
    >>> print(t1u)   # doctest: +NORMALIZE_WHITESPACE
    [  5.81576086  12.52768436] s

    If velocity is very small, (almost) identical solution
    to Galilei transform can be obtained.
    Example is x=5, t=10, v=3. In this case
    time does **not** change (t=10).
    The position is :math:`5 - 3 \times 10 = -25`.
    
    >>> x2, t2 = transformx(5, 10, 3)
    >>> print('{:.2f}'.format(x2))
    -25.00
    >>> print('{:.2f}'.format(t2))
    10.00
    """
    gm = gamma(v)
    x1 = gm * (x - v * t)
    t1 = gm * (t - (v / c) * (x / c))
    return (x1, t1)


transformx_u = u.make_unitful_function2(transformx, [u.m, u.s, u.m/u.s], [u.m, u.s])
''' Unitful version of :func:`transformx`
'''

import unittest
import doctest


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