Source code for irfpy.venus.bowshock

""" Collection of bow shock model

Reference:

  Martinecz et al. (2008): "Location of the bow shock and ion composition boundaries
  at Venus--initial determinations from Venus Express ASPERA-4" PSS 56, 780-784.

.. codeauthor:: Moa Persson
"""
import numpy as _np

def xr_martinecz08(scale=1.0):
    r""" Return the (x,r) coordinates of Martinecz 08 model.

    :returns: A tuple, (x, r), where x is the array of x coordinates
        and r is the array of :math:`\sqrt{y^2+z^2}`. The unit is Rv.

    The parameter is :math:`L = 1.303, \epsilon=1.056, x_0=0.788`.
    Then, :math:`r=L/(1+\epsilon\cos\theta)`, where :math:`(r, \theta)` are
    measured at :math:`(x_0, 0, 0)`.

    >>> from irfpy.venus import bowshock
    >>> x, r = bowshock.xr_martinecz08()
    >>> print(x.shape, r.shape)
    (161,) (161,)

    You can plot the boundary as

    .. code-block:: python

        import matplotlib.pyplot as plt
        plt.plot(x, r, 'r')
    """
    l = 1.303 * scale
    e = 1.056
    x0 = 0.788
    angle = _np.arange(0, 180, 1)*_np.pi/180

    BS_radius = l/(1 + e*_np.cos(angle))
    BS_x = BS_radius*_np.cos(angle) + x0
    BS_y = BS_radius*_np.sin(angle)

    return BS_x[BS_x < BS_radius+x0], BS_y[BS_x < BS_radius+x0]


def inside_martinecz08(x, y, z, scale=1.0):
    """ Return True if the given position is inside Martinecz shock model.

    :param x, y, z: Position.  The unit should be Rv.
    :returns: True if the given position is inside the Martinecz model.

    >>> inside_martinecz08(1, 0, 0)
    True
    >>> inside_martinecz08(2, 0, 0)
    False
    >>> inside_martinecz08(1, 1, 0)
    True
    >>> inside_martinecz08(1, 1, 1)
    False
    """
    l = 1.303 * scale
    e = 1.056
    x0 = 0.788

    xx = x - x0
    yy = _np.sqrt(y ** 2 + z ** 2)

    theta = _np.arctan2(yy, xx)

    rr2 = xx ** 2 + yy ** 2

    rb = l / (1 + e * _np.cos(theta))
    rb2 = rb ** 2

    return rr2 < rb2