pyrocko/src/orthodrome.py

# http://pyrocko.org - GPLv3
#
# The Pyrocko Developers, 21st Century
# ---|P------/S----------~Lg----------
from __future__ import division, absolute_import

import math
import numpy as num

from .moment_tensor import euler_to_matrix
from .config import config
from .plot.beachball import spoly_cut

from matplotlib.path import Path

d2r = math.pi/180.
r2d = 1./d2r
earth_oblateness = 1./298.257223563
earthradius_equator = 6378.14 * 1000.
earthradius = config().earthradius
d2m = earthradius_equator*math.pi/180.
m2d = 1./d2m

_testpath = Path([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)], closed=True)

raise_if_slow_path_contains_points = False


class Slow(Exception):
    pass


if hasattr(_testpath, 'contains_points') and num.all(
        _testpath.contains_points([(0.5, 0.5), (1.5, 0.5)]) == [True, False]):

    def path_contains_points(verts, points):
        p = Path(verts, closed=True)
        return p.contains_points(points).astype(num.bool)

else:
    # work around missing contains_points and bug in matplotlib ~ v1.2.0

    def path_contains_points(verts, points):
        if raise_if_slow_path_contains_points:
            # used by unit test to skip slow gshhg_example.py
            raise Slow()

        p = Path(verts, closed=True)
        result = num.zeros(points.shape[0], dtype=num.bool)
        for i in range(result.size):
            result[i] = p.contains_point(points[i, :])

        return result


try:
    cbrt = num.cbrt
except AttributeError:
    def cbrt(x):
        return x**(1./3.)


def float_array_broadcast(*args):
    return num.broadcast_arrays(*[
        num.asarray(x, dtype=float) for x in args])


class Loc(object):
    '''
    Simple location representation.

    :attrib lat: Latitude in [deg].
    :attrib lon: Longitude in [deg].
    '''
    def __init__(self, lat, lon):
        self.lat = lat
        self.lon = lon


def clip(x, mi, ma):
    '''
    Clip values of an array.

    :param x: Continunous data to be clipped.
    :param mi: Clip minimum.
    :param ma: Clip maximum.
    :type x: :py:class:`numpy.ndarray`
    :type mi: float
    :type ma: float

    :return: Clipped data.
    :rtype: :py:class:`numpy.ndarray`
    '''
    return num.minimum(num.maximum(mi, x), ma)


def wrap(x, mi, ma):
    '''
    Wrapping continuous data to fundamental phase values.

    .. math::
        x_{\\mathrm{wrapped}} = x_{\\mathrm{cont},i} -
            \\frac{ x_{\\mathrm{cont},i} - r_{\\mathrm{min}} }
                  { r_{\\mathrm{max}} -  r_{\\mathrm{min}}}
            \\cdot ( r_{\\mathrm{max}} -  r_{\\mathrm{min}}),\\quad
        x_{\\mathrm{wrapped}}\\; \\in
            \\;[ r_{\\mathrm{min}},\\, r_{\\mathrm{max}}].

    :param x: Continunous data to be wrapped.
    :param mi: Minimum value of wrapped data.
    :param ma: Maximum value of wrapped data.
    :type x: :py:class:`numpy.ndarray`
    :type mi: float
    :type ma: float

    :return: Wrapped data.
    :rtype: :py:class:`numpy.ndarray`
    '''
    return x - num.floor((x-mi)/(ma-mi)) * (ma-mi)


def _latlon_pair(args):
    if len(args) == 2:
        a, b = args
        return a.lat, a.lon, b.lat, b.lon

    elif len(args) == 4:
        return args


def cosdelta(*args):
    '''
    Cosine of the angular distance between two points ``a`` and ``b`` on a
    sphere.

    This function (find implementation below) returns the cosine of the
    distance angle 'delta' between two points ``a`` and ``b``, coordinates of
    which are expected to be given in geographical coordinates and in degrees.
    For numerical stability a maximum of 1.0 is enforced.

    .. math::

        A_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot A_{lat}, \\quad
        A_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot A_{lon}, \\quad
        B_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot B_{lat}, \\quad
        B_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot B_{lon}\\\\[0.5cm]

        \\cos(\\Delta) = \\min( 1.0, \\quad  \\sin( A_{\\mathrm{lat'}})
                     \\sin( B_{\\mathrm{lat'}} ) +
                     \\cos(A_{\\mathrm{lat'}})  \\cos( B_{\\mathrm{lat'}} )
                     \\cos( B_{\\mathrm{lon'}} - A_{\\mathrm{lon'}} )

    :param a: Location point A.
    :type a: :py:class:`pyrocko.orthodrome.Loc`
    :param b: Location point B.
    :type b: :py:class:`pyrocko.orthodrome.Loc`

    :return: Cosdelta.
    :rtype: float
    '''

    alat, alon, blat, blon = _latlon_pair(args)

    return min(
        1.0,
        math.sin(alat*d2r) * math.sin(blat*d2r) +
        math.cos(alat*d2r) * math.cos(blat*d2r) *
        math.cos(d2r*(blon-alon)))


def cosdelta_numpy(a_lats, a_lons, b_lats, b_lons):
    '''
    Cosine of the angular distance between two points ``a`` and ``b`` on a
    sphere.

    This function returns the cosines of the distance
    angles *delta* between two points ``a`` and ``b`` given as
    :py:class:`numpy.ndarray`.
    The coordinates are expected to be given in geographical coordinates
    and in degrees. For numerical stability a maximum of ``1.0`` is enforced.

    Please find the details of the implementation in the documentation of
    the function :py:func:`pyrocko.orthodrome.cosdelta` above.

    :param a_lats: Latitudes in [deg] point A.
    :param a_lons: Longitudes in [deg] point A.
    :param b_lats: Latitudes in [deg] point B.
    :param b_lons: Longitudes in [deg] point B.
    :type a_lats: :py:class:`numpy.ndarray`
    :type a_lons: :py:class:`numpy.ndarray`
    :type b_lats: :py:class:`numpy.ndarray`
    :type b_lons: :py:class:`numpy.ndarray`

    :return: Cosdelta.
    :type b_lons: :py:class:`numpy.ndarray`, ``(N)``
    '''
    return num.minimum(
        1.0,
        num.sin(a_lats*d2r) * num.sin(b_lats*d2r) +
        num.cos(a_lats*d2r) * num.cos(b_lats*d2r) *
        num.cos(d2r*(b_lons-a_lons)))


def azimuth(*args):
    '''
    Azimuth calculation.

    This function (find implementation below) returns azimuth ...
    between points ``a`` and ``b``, coordinates of
    which are expected to be given in geographical coordinates and in degrees.

    .. math::

        A_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot A_{lat}, \\quad
        A_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot A_{lon}, \\quad
        B_{\\mathrm{lat'}} = \\frac{ \\pi}{180} \\cdot B_{lat}, \\quad
        B_{\\mathrm{lon'}} = \\frac{ \\pi}{180} \\cdot B_{lon}\\\\

        \\varphi_{\\mathrm{azi},AB} = \\frac{180}{\\pi} \\arctan \\left[
            \\frac{
                \\cos( A_{\\mathrm{lat'}}) \\cos( B_{\\mathrm{lat'}} )
                \\sin(B_{\\mathrm{lon'}} - A_{\\mathrm{lon'}} )}
                {\\sin ( B_{\\mathrm{lat'}} ) - \\sin( A_{\\mathrm{lat'}}
                  cosdelta) } \\right]

    :param a: Location point A.
    :type a: :py:class:`pyrocko.orthodrome.Loc`
    :param b: Location point B.
    :type b: :py:class:`pyrocko.orthodrome.Loc`

    :return: Azimuth in degree
    '''

    alat, alon, blat, blon = _latlon_pair(args)

    return r2d*math.atan2(
        math.cos(alat*d2r) * math.cos(blat*d2r) *
        math.sin(d2r*(blon-alon)),
        math.sin(d2r*blat) - math.sin(d2r*alat) * cosdelta(
            alat, alon, blat, blon))


def azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta=None):
    '''
    Calculation of the azimuth (*track angle*) from a location A towards B.

    This function returns azimuths (*track angles*) from locations A towards B
    given in :py:class:`numpy.ndarray`. Coordinates are expected to be given in
    geographical coordinates and in degrees.

    Please find the details of the implementation in the documentation of the
    function :py:func:`pyrocko.orthodrome.azimuth`.


    :param a_lats: Latitudes in [deg] point A.
    :param a_lons: Longitudes in [deg] point A.
    :param b_lats: Latitudes in [deg] point B.
    :param b_lons: Longitudes in [deg] point B.
    :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
    :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
    :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
    :type b_lons: :py:class:`numpy.ndarray`, ``(N)``

    :return: Azimuths in [deg].
    :rtype: :py:class:`numpy.ndarray`, ``(N)``
    '''
    if _cosdelta is None:
        _cosdelta = cosdelta_numpy(a_lats, a_lons, b_lats, b_lons)

    return r2d*num.arctan2(
        num.cos(a_lats*d2r) * num.cos(b_lats*d2r) *
        num.sin(d2r*(b_lons-a_lons)),
        num.sin(d2r*b_lats) - num.sin(d2r*a_lats) * _cosdelta)


def azibazi(*args, **kwargs):
    '''
    Azimuth and backazimuth from location A towards B and back.

    :returns: Azimuth in [deg] from A to B, back azimuth in [deg] from B to A.
    :rtype: tuple[float, float]
    '''

    alat, alon, blat, blon = _latlon_pair(args)
    if alat == blat and alon == blon:
        return 0., 180.

    implementation = kwargs.get('implementation', 'c')
    assert implementation in ('c', 'python')
    if implementation == 'c':
        from pyrocko import orthodrome_ext
        return orthodrome_ext.azibazi(alat, alon, blat, blon)

    cd = cosdelta(alat, alon, blat, blon)
    azi = r2d*math.atan2(
        math.cos(alat*d2r) * math.cos(blat*d2r) *
        math.sin(d2r*(blon-alon)),
        math.sin(d2r*blat) - math.sin(d2r*alat) * cd)
    bazi = r2d*math.atan2(
        math.cos(blat*d2r) * math.cos(alat*d2r) *
        math.sin(d2r*(alon-blon)),
        math.sin(d2r*alat) - math.sin(d2r*blat) * cd)

    return azi, bazi


def azibazi_numpy(a_lats, a_lons, b_lats, b_lons, implementation='c'):
    '''
    Azimuth and backazimuth from location A towards B and back.

    Arguments are given as :py:class:`numpy.ndarray`.

    :param a_lats: Latitude(s) in [deg] of point A.
    :type a_lats: :py:class:`numpy.ndarray`
    :param a_lons: Longitude(s) in [deg] of point A.
    :type a_lons: :py:class:`numpy.ndarray`
    :param b_lats: Latitude(s) in [deg] of point B.
    :type b_lats: :py:class:`numpy.ndarray`
    :param b_lons: Longitude(s) in [deg] of point B.
    :type b_lons: :py:class:`numpy.ndarray`

    :returns: Azimuth(s) in [deg] from A to B,
        back azimuth(s) in [deg] from B to A.
    :rtype: :py:class:`numpy.ndarray`, :py:class:`numpy.ndarray`
    '''

    a_lats, a_lons, b_lats, b_lons = float_array_broadcast(
        a_lats, a_lons, b_lats, b_lons)

    assert implementation in ('c', 'python')
    if implementation == 'c':
        from pyrocko import orthodrome_ext
        return orthodrome_ext.azibazi_numpy(a_lats, a_lons, b_lats, b_lons)

    _cosdelta = cosdelta_numpy(a_lats, a_lons, b_lats, b_lons)
    azis = azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta)
    bazis = azimuth_numpy(b_lats, b_lons, a_lats, a_lons, _cosdelta)

    eq = num.logical_and(a_lats == b_lats, a_lons == b_lons)
    ii_eq = num.where(eq)[0]
    azis[ii_eq] = 0.0
    bazis[ii_eq] = 180.0
    return azis, bazis


def azidist_numpy(*args):
    '''
    Calculation of the azimuth (*track angle*) and the distance from locations
    A towards B on a sphere.

    The assisting functions used are :py:func:`pyrocko.orthodrome.cosdelta` and
    :py:func:`pyrocko.orthodrome.azimuth`

    :param a_lats: Latitudes in [deg] point A.
    :param a_lons: Longitudes in [deg] point A.
    :param b_lats: Latitudes in [deg] point B.
    :param b_lons: Longitudes in [deg] point B.
    :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
    :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
    :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
    :type b_lons: :py:class:`numpy.ndarray`, ``(N)``

    :return: Azimuths in [deg], distances in [deg].
    :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
    '''
    _cosdelta = cosdelta_numpy(*args)
    _azimuths = azimuth_numpy(_cosdelta=_cosdelta, *args)
    return _azimuths, r2d*num.arccos(_cosdelta)


def distance_accurate50m(*args, **kwargs):
    '''
    Accurate distance calculation based on a spheroid of rotation.

    Function returns distance in meter between points A and B, coordinates of
    which must be given in geographical coordinates and in degrees.
    The returned distance should be accurate to 50 m using WGS84.
    Values for the Earth's equator radius and the Earth's oblateness
    (``f_oblate``) are defined in the pyrocko configuration file
    :py:class:`pyrocko.config`.

    From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:

    ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell,
    Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1``

    .. math::

        F = \\frac{\\pi}{180}
                            \\frac{(A_{lat} + B_{lat})}{2}, \\quad
        G = \\frac{\\pi}{180}
                            \\frac{(A_{lat} - B_{lat})}{2}, \\quad
        l = \\frac{\\pi}{180}
                            \\frac{(A_{lon} - B_{lon})}{2} \\quad
        \\\\[0.5cm]
        S = \\sin^2(G) \\cdot \\cos^2(l) +
              \\cos^2(F) \\cdot \\sin^2(l), \\quad \\quad
        C = \\cos^2(G) \\cdot \\cos^2(l) +
                  \\sin^2(F) \\cdot \\sin^2(l)

    .. math::

        w = \\arctan \\left( \\sqrt{ \\frac{S}{C}} \\right) , \\quad
        r =  \\sqrt{\\frac{S}{C} }

    The spherical-earth distance D between A and B, can be given with:

    .. math::

        D_{sphere} = 2w \\cdot R_{equator}

    The oblateness of the Earth requires some correction with
    correction factors h1 and h2:

     .. math::

        h_1 = \\frac{3r - 1}{2C}, \\quad
        h_2 = \\frac{3r +1 }{2S}\\\\[0.5cm]

        D = D_{\\mathrm{sphere}} \\cdot [ 1 + h_1 \\,f_{\\mathrm{oblate}}
                 \\cdot \\sin^2(F)
                 \\cos^2(G) - h_2\\, f_{\\mathrm{oblate}}
                 \\cdot \\cos^2(F) \\sin^2(G)]


    :param a: Location point A.
    :type a: :py:class:`pyrocko.orthodrome.Loc`
    :param b: Location point B.
    :type b: :py:class:`pyrocko.orthodrome.Loc`

    :return: Distance in [m].
    :rtype: float
    '''

    alat, alon, blat, blon = _latlon_pair(args)

    implementation = kwargs.get('implementation', 'c')
    assert implementation in ('c', 'python')
    if implementation == 'c':
        from pyrocko import orthodrome_ext
        return orthodrome_ext.distance_accurate50m(alat, alon, blat, blon)

    f = (alat + blat)*d2r / 2.
    g = (alat - blat)*d2r / 2.
    h = (alon - blon)*d2r / 2.

    s = math.sin(g)**2 * math.cos(h)**2 + math.cos(f)**2 * math.sin(h)**2
    c = math.cos(g)**2 * math.cos(h)**2 + math.sin(f)**2 * math.sin(h)**2

    w = math.atan(math.sqrt(s/c))

    if w == 0.0:
        return 0.0

    r = math.sqrt(s*c)/w
    d = 2.*w*earthradius_equator
    h1 = (3.*r-1.)/(2.*c)
    h2 = (3.*r+1.)/(2.*s)

    return d * (1. +
                earth_oblateness * h1 * math.sin(f)**2 * math.cos(g)**2 -
                earth_oblateness * h2 * math.cos(f)**2 * math.sin(g)**2)


def distance_accurate50m_numpy(
        a_lats, a_lons, b_lats, b_lons, implementation='c'):

    '''
    Accurate distance calculation based on a spheroid of rotation.

    Function returns distance in meter between points ``a`` and ``b``,
    coordinates of which must be given in geographical coordinates and in
    degrees.
    The returned distance should be accurate to 50 m using WGS84.
    Values for the Earth's equator radius and the Earth's oblateness
    (``f_oblate``) are defined in the pyrocko configuration file
    :py:class:`pyrocko.config`.

    From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:

    ``Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell,
    Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1``

    .. math::

        F_i = \\frac{\\pi}{180}
                           \\frac{(a_{lat,i} + a_{lat,i})}{2}, \\quad
        G_i = \\frac{\\pi}{180}
                           \\frac{(a_{lat,i} - b_{lat,i})}{2}, \\quad
        l_i= \\frac{\\pi}{180}
                            \\frac{(a_{lon,i} - b_{lon,i})}{2} \\\\[0.5cm]
        S_i = \\sin^2(G_i) \\cdot \\cos^2(l_i) +
              \\cos^2(F_i) \\cdot \\sin^2(l_i), \\quad \\quad
        C_i = \\cos^2(G_i) \\cdot \\cos^2(l_i) +
                  \\sin^2(F_i) \\cdot \\sin^2(l_i)

    .. math::

        w_i = \\arctan \\left( \\sqrt{\\frac{S_i}{C_i}} \\right), \\quad
        r_i =  \\sqrt{\\frac{S_i}{C_i} }

    The spherical-earth distance ``D`` between ``a`` and ``b``,
    can be given with:

    .. math::

        D_{\\mathrm{sphere},i} = 2w_i \\cdot R_{\\mathrm{equator}}

    The oblateness of the Earth requires some correction with
    correction factors ``h1`` and ``h2``:

     .. math::

        h_{1.i} = \\frac{3r - 1}{2C_i}, \\quad
        h_{2,i} = \\frac{3r +1 }{2S_i}\\\\[0.5cm]

        D_{AB,i} = D_{\\mathrm{sphere},i} \\cdot [1 + h_{1,i}
            \\,f_{\\mathrm{oblate}}
            \\cdot \\sin^2(F_i)
            \\cos^2(G_i) - h_{2,i}\\, f_{\\mathrm{oblate}}
            \\cdot \\cos^2(F_i) \\sin^2(G_i)]

    :param a_lats: Latitudes in [deg] point A.
    :param a_lons: Longitudes in [deg] point A.
    :param b_lats: Latitudes in [deg] point B.
    :param b_lons: Longitudes in [deg] point B.
    :type a_lats: :py:class:`numpy.ndarray`, ``(N)``
    :type a_lons: :py:class:`numpy.ndarray`, ``(N)``
    :type b_lats: :py:class:`numpy.ndarray`, ``(N)``
    :type b_lons: :py:class:`numpy.ndarray`, ``(N)``

    :return: Distances in [m].
    :rtype: :py:class:`numpy.ndarray`, ``(N)``
    '''

    a_lats, a_lons, b_lats, b_lons = float_array_broadcast(
        a_lats, a_lons, b_lats, b_lons)

    assert implementation in ('c', 'python')
    if implementation == 'c':
        from pyrocko import orthodrome_ext
        return orthodrome_ext.distance_accurate50m_numpy(
            a_lats, a_lons, b_lats, b_lons)

    eq = num.logical_and(a_lats == b_lats, a_lons == b_lons)
    ii_neq = num.where(num.logical_not(eq))[0]

    if num.all(eq):
        return num.zeros_like(eq, dtype=float)

    def extr(x):
        if isinstance(x, num.ndarray) and x.size > 1:
            return x[ii_neq]
        else:
            return x

    a_lats = extr(a_lats)
    a_lons = extr(a_lons)
    b_lats = extr(b_lats)
    b_lons = extr(b_lons)

    f = (a_lats + b_lats)*d2r / 2.
    g = (a_lats - b_lats)*d2r / 2.
    h = (a_lons - b_lons)*d2r / 2.

    s = num.sin(g)**2 * num.cos(h)**2 + num.cos(f)**2 * num.sin(h)**2
    c = num.cos(g)**2 * num.cos(h)**2 + num.sin(f)**2 * num.sin(h)**2

    w = num.arctan(num.sqrt(s/c))

    r = num.sqrt(s*c)/w

    d = 2.*w*earthradius_equator
    h1 = (3.*r-1.)/(2.*c)
    h2 = (3.*r+1.)/(2.*s)

    dists = num.zeros(eq.size, dtype=float)
    dists[ii_neq] = d * (
        1. +
        earth_oblateness * h1 * num.sin(f)**2 * num.cos(g)**2 -
        earth_oblateness * h2 * num.cos(f)**2 * num.sin(g)**2)

    return dists


def ne_to_latlon(lat0, lon0, north_m, east_m):
    '''
    Transform local cartesian coordinates to latitude and longitude.

    From east and north coordinates (``x`` and ``y`` coordinate
    :py:class:`numpy.ndarray`)  relative to a reference differences in
    longitude and latitude are calculated, which are effectively changes in
    azimuth and distance, respectively:

    .. math::

        \\text{distance change:}\\; \\Delta {\\bf{a}} &= \\sqrt{{\\bf{y}}^2 +
                                           {\\bf{x}}^2 }/ \\mathrm{R_E},

        \\text{azimuth change:}\\; \\Delta \\bf{\\gamma} &= \\arctan( \\bf{x}
                                                         / \\bf{y}).

    The projection used preserves the azimuths of the input points.

    :param lat0: Latitude origin of the cartesian coordinate system in [deg].
    :param lon0: Longitude origin of the cartesian coordinate system in [deg].
    :param north_m: Northing distances from origin in [m].
    :param east_m: Easting distances from origin in [m].
    :type north_m: :py:class:`numpy.ndarray`, ``(N)``
    :type east_m: :py:class:`numpy.ndarray`, ``(N)``
    :type lat0: float
    :type lon0: float

    :return: Array with latitudes and longitudes in [deg].
    :rtype: :py:class:`numpy.ndarray`, ``(2xN)``

    '''

    a = num.sqrt(north_m**2+east_m**2)/earthradius
    gamma = num.arctan2(east_m, north_m)

    return azidist_to_latlon_rad(lat0, lon0, gamma, a)


def azidist_to_latlon(lat0, lon0, azimuth_deg, distance_deg):
    '''
    Absolute latitudes and longitudes are calculated from relative changes.

    Convenience wrapper to :py:func:`azidist_to_latlon_rad` with azimuth and
    distance given in degrees.

    :param lat0: Latitude origin of the cartesian coordinate system in [deg].
    :type lat0: float
    :param lon0: Longitude origin of the cartesian coordinate system in [deg].
    :type lon0: float
    :param azimuth_deg: Azimuth from origin in [deg].
    :type azimuth_deg: :py:class:`numpy.ndarray`, ``(N)``
    :param distance_deg: Distances from origin in [deg].
    :type distance_deg: :py:class:`numpy.ndarray`, ``(N)``

    :return: Array with latitudes and longitudes in [deg].
    :rtype: :py:class:`numpy.ndarray`, ``(2xN)``
    '''

    return azidist_to_latlon_rad(
        lat0, lon0, azimuth_deg/180.*num.pi, distance_deg/180.*num.pi)


def azidist_to_latlon_rad(lat0, lon0, azimuth_rad, distance_rad):
    '''
    Absolute latitudes and longitudes are calculated from relative changes.

    For numerical stability a range between of ``-1.0`` and ``1.0`` is
    enforced for ``c`` and ``alpha``.

    .. math::

        \\Delta {\\bf a}_i \\; \\text{and} \\; \\Delta \\gamma_i \\;
            \\text{are relative distances and azimuths from lat0 and lon0 for
                \\textit{i} source points of a finite source.}

    .. math::

        \\mathrm{b} &= \\frac{\\pi}{2} -\\frac{\\pi}{180}\\;\\mathrm{lat_0}\\\\
        {\\bf c}_i &=\\arccos[\\; \\cos(\\Delta {\\bf{a}}_i)
            \\cos(\\mathrm{b}) + |\\Delta \\gamma_i| \\,
            \\sin(\\Delta {\\bf a}_i)
                \\sin(\\mathrm{b})\\; ] \\\\
        \\mathrm{lat}_i &=  \\frac{180}{\\pi}
                          \\left(\\frac{\\pi}{2} - {\\bf c}_i \\right)

    .. math::

         \\alpha_i &= \\arcsin \\left[ \\; \\frac{ \\sin(\\Delta {\\bf a}_i )
                         \\sin(|\\Delta \\gamma_i|)}{\\sin({\\bf c}_i)}\\;
                         \\right] \\\\
        \\alpha_i &= \\begin{cases}
                 \\alpha_i, &\\text{if}  \\; \\cos(\\Delta {\\bf a}_i) -
                 \\cos(\\mathrm{b}) \\cos({\\bf{c}}_i) > 0, \\;
                 \\text{else} \\\\
                 \\pi - \\alpha_i, & \\text{if} \\; \\alpha_i > 0,\\;
                 \\text{else}\\\\
                -\\pi - \\alpha_i, & \\text{if} \\; \\alpha_i < 0.
                \\end{cases} \\\\
        \\mathrm{lon}_i &=   \\mathrm<