Source code for irfpy.venus.icb

""" Collection of ICB 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
.. codeauthor:: Yoshifumi Futaana

"""
import logging as _logging

import numpy as _np

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

    - Dayside: Circle with the radius 1.109 Rv.
    - Nightside: :math:`y=kx+d`, where :math:`k=-0.097` and :math:`d=1.109`.

    :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}`.

    >>> from irfpy.venus import icb
    >>> x, r = icb.xr_martinecz08()
    >>> print(x.shape, r.shape)
    (590,) (590,)

    You can plot the boundary as

    .. code-block:: python

        import matplotlib.pyplot as plt
        plt.plot(x, r, 'r')
    """
    logger = _logging.getLogger(__name__)

    x_lims = _np.arange(0, -50., -0.1)
    k = -0.097
    d = 1.109 * scale

    if d <= 1.0:
        logger.warning('Given scale ({}) is so small that the boundary is inside the planet.'.format(scale))
        logger.warning('Nevertheless, the program continues.')

    icb_line = k * x_lims + d

    icb_circle_x = _np.cos(_np.deg2rad(_np.arange(0, 90))) * d
    icb_circle_y = _np.sin(_np.deg2rad(_np.arange(0, 90))) * d


    x = _np.concatenate([icb_circle_x, x_lims])
    r = _np.concatenate([icb_circle_y, icb_line])

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

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

    .. codeauthor:: Moa Persson

    >>> from irfpy.venus import icb
    >>> print(icb.inside_martinecz08(1, 0, 0))
    True
    >>> print(icb.inside_martinecz08(2, 0, 0))
    False
    >>> print(icb.inside_martinecz08(0, 1.1, 0))
    True
    >>> print(icb.inside_martinecz08(-1, 1.2, 0))
    True
    >>> print(icb.inside_martinecz08(-1, 1.208, 0))
    False
    """
    k = -0.097
    d = 1.109 * scale

    if d <= 1.0:
        logger.warning('Given scale ({}) is so small that the boundary is inside the planet.'.format(scale))
        logger.warning('Nevertheless, the program continues.')

    if x >= 0.:
        rr = _np.sqrt(x**2 + y**2 + z**2)
        return rr < d
    elif x < 0.:
        icb_line = k * x + d
        rr = _np.sqrt(y**2 + z**2)
        return rr < icb_line