add geometry sub-package with convenience function (spherical coords)

also move true_angular_distance and great_circle_bearing from
pathprof.helper to geometry module

docs/geometry/index.rst

@@ -0,0 +1,27 @@
+.. pycraf-geometry:
+
+**************************************
+Geometry helpers (`pycraf.geometry`)
+**************************************
+
+.. currentmodule:: pycraf.geometry
+
+Introduction
+============
+
+The `~pycraf.geometry` sub-package just offers some convenience functions
+for various geometrical calculations, e.g., in spherical coordinates.
+
+
+See Also
+========
+
+- `Astropy Units and Quantities package <http://docs.astropy.org/en/stable/
+  units/index.html>`_, which is used extensively in pycraf.
+- `Spherical coordinates <https://en.wikipedia.org/wiki/Spherical_coordinate_system>`_
+
+Reference/API
+=============
+
+.. automodapi:: pycraf.geometry
+    :no-inheritance-diagram:

docs/index.rst

@@ -42,6 +42,7 @@ Available modules
    protection/index
    geospatial/index
    satellite/index
+   geometry/index
 
 .. automodapi:: pycraf
     :no-heading:

docs/satellite/index.rst

@@ -1,4 +1,4 @@
-.. pycraf-atm:
+.. pycraf-satellite:
 
 **************************************
 Satellites (`pycraf.satellite`)

notebooks/A01_case_study_fixed_link.ipynb

@@ -57,7 +57,7 @@
     "import matplotlib.pyplot as plt\n",
     "from astropy import units as u\n",
     "from pycraf import conversions as cnv\n",
-    "from pycraf import pathprof, protection, antenna"
+    "from pycraf import pathprof, protection, antenna, geometry"
    ]
   },
   {
@@ -480,7 +480,7 @@
     }
    ],
    "source": [
-    "ang_dist = pathprof.true_angular_distance(\n",
+    "ang_dist = geometry.true_angular_distance(\n",
     "    pprop_fl.alpha_tr, pprop_fl.eps_pt,\n",
     "    pprop_fl_ras.alpha_tr, pprop_fl_ras.eps_pt,\n",
     "    )\n",

pycraf/__init__.py

@@ -19,6 +19,7 @@
     from . import antenna
     from . import atm
     from . import conversions
+    from . import geometry
     from . import geospatial
     from . import utils
     from . import pathprof

pycraf/geometry/__init__.py

@@ -0,0 +1,7 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+'''
+Contains convenience functions for geometrical calculations on the sphere.
+'''
+
+from .geometry import *

pycraf/geometry/geometry.py

@@ -0,0 +1,173 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+
+from __future__ import (
+    absolute_import, unicode_literals, division, print_function
+    )
+
+from astropy import units as apu
+import numpy as np
+from .. import utils
+
+
+__all__ = [
+    'true_angular_distance', 'great_circle_bearing',
+    'cart_to_sphere', 'sphere_to_cart',
+    ]
+
+
+@utils.ranged_quantity_input(
+    l1=(None, None, apu.deg),
+    b1=(-90, 90, apu.deg),
+    l2=(None, None, apu.deg),
+    b2=(-90, 90, apu.deg),
+    strip_input_units=True, output_unit=apu.deg,
+    )
+def true_angular_distance(l1, b1, l2, b2):
+    '''
+    True angular distance between points (l1, b1) and (l2, b2).
+
+    Based on Vincenty formula
+    (http://en.wikipedia.org/wiki/Great-circle_distance).
+    This was spotted in astropy source code.
+
+    Parameters
+    ----------
+    l1, b1 : `~astropy.units.Quantity`
+        Longitude/Latitude of point 1 [deg]
+    l2, b2 : `~astropy.units.Quantity`
+        Longitude/Latitude of point 2 [deg]
+
+    Returns
+    -------
+    adist : `~astropy.units.Quantity`
+        True angular distance [deg]
+    '''
+
+    sin_diff_lon = np.sin(np.radians(l2 - l1))
+    cos_diff_lon = np.cos(np.radians(l2 - l1))
+    sin_lat1 = np.sin(np.radians(b1))
+    sin_lat2 = np.sin(np.radians(b2))
+    cos_lat1 = np.cos(np.radians(b1))
+    cos_lat2 = np.cos(np.radians(b2))
+
+    num1 = cos_lat2 * sin_diff_lon
+    num2 = cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_diff_lon
+    denominator = sin_lat1 * sin_lat2 + cos_lat1 * cos_lat2 * cos_diff_lon
+
+    return np.degrees(np.arctan2(
+        np.sqrt(num1 ** 2 + num2 ** 2), denominator
+        ))
+
+
+@utils.ranged_quantity_input(
+    l1=(None, None, apu.deg),
+    b1=(-90, 90, apu.deg),
+    l2=(None, None, apu.deg),
+    b2=(-90, 90, apu.deg),
+    strip_input_units=True, output_unit=apu.deg,
+    )
+def great_circle_bearing(l1, b1, l2, b2):
+    '''
+    Great circle bearing between points (l1, b1) and (l2, b2).
+
+    Parameters
+    ----------
+    l1, b1 : `~astropy.units.Quantity`
+        Longitude/Latitude of point 1 [deg]
+    l2, b2 : `~astropy.units.Quantity`
+        Longitude/Latitude of point 2 [deg]
+
+    Returns
+    -------
+    bearing : `~astropy.units.Quantity`
+        Great circle bearing [deg]
+    '''
+
+    sin_diff_lon = np.sin(np.radians(l2 - l1))
+    cos_diff_lon = np.cos(np.radians(l2 - l1))
+    sin_lat1 = np.sin(np.radians(b1))
+    sin_lat2 = np.sin(np.radians(b2))
+    cos_lat1 = np.cos(np.radians(b1))
+    cos_lat2 = np.cos(np.radians(b2))
+
+    a = cos_lat2 * sin_diff_lon
+    b = cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_diff_lon
+
+    return np.degrees(np.arctan2(a, b))
+
+
+@utils.ranged_quantity_input(
+    x=(None, None, apu.m),
+    y=(None, None, apu.m),
+    z=(None, None, apu.m),
+    strip_input_units=True, output_unit=(apu.m, apu.deg, apu.deg)
+    )
+def cart_to_sphere(x, y, z):
+    '''
+    Spherical coordinates from Cartesian representation.
+
+    Parameters
+    ----------
+    x, y, z : `~astropy.units.Quantity`
+        Cartesian position [m]
+
+    Returns
+    -------
+    r : `~astropy.units.Quantity`
+        Radial distance [m]
+    theta : `~astropy.units.Quantity`
+        Elevation [deg]
+    phi : `~astropy.units.Quantity`
+        Azimuth [deg]
+
+    Notes
+    -----
+    Unlike with the mathematical definition, `theta` is not the angle
+    to the (positive) `z` axis, but the elevation above the `x`-`y` plane.
+    '''
+
+    r = np.sqrt(x ** 2 + y ** 2 + z ** 2)
+    theta = 90 - np.degrees(np.arccos(z / r))
+    phi = np.degrees(np.arctan2(y, x))
+
+    return r, theta, phi
+
+
+@utils.ranged_quantity_input(
+    r=(None, None, apu.m),
+    theta=(-90, 90, apu.deg),
+    phi=(None, None, apu.deg),
+    strip_input_units=True, output_unit=(apu.m, apu.m, apu.m)
+    )
+def sphere_to_cart(r, theta, phi):
+    '''
+    Spherical coordinates from Cartesian representation.
+
+    Parameters
+    ----------
+    r : `~astropy.units.Quantity`
+        Radial distance [m]
+    theta : `~astropy.units.Quantity`
+        Elevation [deg]
+    phi : `~astropy.units.Quantity`
+        Azimuth [deg]
+
+    Returns
+    -------
+    x, y, z : `~astropy.units.Quantity`
+        Cartesian position [m]
+
+    Notes
+    -----
+    Unlike with the mathematical definition, `theta` is not the angle
+    to the (positive) `z` axis, but the elevation above the `x`-`y` plane.
+    '''
+
+    theta = 90. - theta
+
+    x = r * np.sin(np.radians(theta)) * np.cos(np.radians(phi))
+    y = r * np.sin(np.radians(theta)) * np.sin(np.radians(phi))
+    z = r * np.cos(np.radians(theta))
+
+    return x, y, z

pycraf/geometry/tests/test_geometry.py

@@ -0,0 +1,137 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+# from __future__ import absolute_import
+# from __future__ import division
+# from __future__ import print_function
+# from __future__ import unicode_literals
+
+import pytest
+from functools import partial
+import numpy as np
+from numpy.testing import assert_equal, assert_allclose
+from astropy.tests.helper import assert_quantity_allclose
+from astropy import units as apu
+from astropy.units import Quantity
+from ... import conversions as cnv
+from ... import geometry
+from ...utils import check_astro_quantities
+from astropy.utils.misc import NumpyRNGContext
+
+
+TOL_KWARGS = {'atol': 1.e-4, 'rtol': 1.e-4}
+
+
+def test_true_angular_distance():
+
+    pfunc = geometry.true_angular_distance
+    args_list = [
+        (None, None, apu.deg),
+        (-90, 90, apu.deg),
+        (None, None, apu.deg),
+        (-90, 90, apu.deg),
+        ]
+    check_astro_quantities(pfunc, args_list)
+
+    l1 = Quantity([10., 20., 30., 140., 350.], apu.deg)
+    b1 = Quantity([0., 20., 60., 40., 80.], apu.deg)
+    l2 = Quantity([-40., 200., 30.1, 130., 320.], apu.deg)
+    b2 = Quantity([0., -30., 60.05, 30., 10.], apu.deg)
+
+    adist = Quantity(
+        [
+            5.00000000e+01, 1.70000000e+02, 7.06839454e-02,
+            1.29082587e+01, 7.13909421e+01
+            ], apu.deg
+        )
+    assert_quantity_allclose(
+        pfunc(l1, b1, l2, b2),
+        adist,
+        )
+
+
+def test_great_circle_bearing():
+
+    pfunc = geometry.great_circle_bearing
+    args_list = [
+        (None, None, apu.deg),
+        (-90, 90, apu.deg),
+        (None, None, apu.deg),
+        (-90, 90, apu.deg),
+        ]
+    check_astro_quantities(pfunc, args_list)
+
+    l1 = Quantity([10., 20., 30., 140., 350.], apu.deg)
+    b1 = Quantity([0., 20., 60., 40., 80.], apu.deg)
+    l2 = Quantity([-40., 200., 30.1, 130., 320.], apu.deg)
+    b2 = Quantity([0., -30., 60.05, 30., 10.], apu.deg)
+
+    adist = Quantity(
+        [
+            -90., 180., 44.935034, -137.686456, -148.696726
+            ], apu.deg
+        )
+    assert_quantity_allclose(
+        pfunc(l1, b1, l2, b2),
+        adist,
+        )
+
+
+def test_cart_to_sphere():
+
+    pfunc = geometry.cart_to_sphere
+    args_list = [
+        (None, None, apu.m),
+        (None, None, apu.m),
+        (None, None, apu.m),
+        ]
+    check_astro_quantities(pfunc, args_list)
+
+    x = Quantity([0., 20., 60., 40., 80.], apu.m)
+    y = Quantity([10., 20., 30., 140., 350.], apu.m)
+    z = Quantity([-40., 200., 30.1, 130., 320.], apu.m)
+
+    r = Quantity([
+        41.23105626, 201.99009877, 73.52557378, 195.19221296,
+        480.93658626
+        ], apu.m)
+    theta = Quantity([
+        -75.96375653, 81.95053302, 24.16597925, 41.75987004,
+        41.71059314
+        ], apu.deg)
+    phi = Quantity([90., 45., 26.56505118, 74.0546041, 77.12499844], apu.deg)
+
+    _r, _theta, _phi = pfunc(x, y, z)
+    assert_quantity_allclose(_r, r)
+    assert_quantity_allclose(_theta, theta)
+    assert_quantity_allclose(_phi, phi)
+
+
+def test_sphere_to_cart():
+
+    pfunc = geometry.sphere_to_cart
+    args_list = [
+        (None, None, apu.m),
+        (-90, 90, apu.deg),
+        (None, None, apu.deg),
+        ]
+    check_astro_quantities(pfunc, args_list)
+
+    x = Quantity([0., 20., 60., 40., 80.], apu.m)
+    y = Quantity([10., 20., 30., 140., 350.], apu.m)
+    z = Quantity([-40., 200., 30.1, 130., 320.], apu.m)
+
+    r = Quantity([
+        41.23105626, 201.99009877, 73.52557378, 195.19221296,
+        480.93658626
+        ], apu.m)
+    theta = Quantity([
+        -75.96375653, 81.95053302, 24.16597925, 41.75987004,
+        41.71059314
+        ], apu.deg)
+    phi = Quantity([90., 45., 26.56505118, 74.0546041, 77.12499844], apu.deg)
+
+    _x, _y, _z = pfunc(r, theta, phi)
+    assert_quantity_allclose(_x, x, atol=1.e-8 * apu.m)
+    assert_quantity_allclose(_y, y)
+    assert_quantity_allclose(_z, z)