# Copyright (c) 2011-2014, B.I.Stepanov Institute of Physics, National Academy
# of Sciences of Belarus.
#
# This program is free software; you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free Software
# Foundation; either version 2 of the License, or (at your option) any later
# version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# this program; if not, write to the Free Software Foundation, Inc., 51
# Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

import math

__all__ = ['GeodeticPoint']

# *****************************************************************************
class GeodeticPoint:
    """Geodetic coordinates of a point with respect to WGS84 ellipsoid."""

    # Major radius of WGS84 ellipsoid in meters.
    ellipsoidMajorRadius = 6378137.0
    # Minor radius of WGS84 ellipsoid in meters.
    ellipsoidMinorRadius = 6356752.314245

    # ---- Public methods -----------------------------------------------------
    def __init__(self, latitude, longitude, altitude = 0.0):
        """Parameters:
          - 'latitide': geodetic latitude, in degrees (angle between the
            equatorial plane and a normal to the reference ellipsoid passing
            through the point, positive direction northward).
          - 'longitude': geodetic longitude, in degrees (angle between the
            meridional plane of the reference meridian and a normal to the
            reference ellipsoid passing through the point, positive direction
            eastward).
          - 'altitude': geodetic altitude, in meters (normal distance to the
            point from the reference ellipsoid, positive direction upward)."""

        self.latitude = latitude
        self.longitude = longitude
        self.altitude = altitude

    def zeroAltitudePoint(self):
        """Return a copy of 'self' with the same latitude and longitude, but
        zero altitude."""

        return GeoPoint(self.latitude, self.longitude)

    def cartesianCoords(self):
        """Calculate geocentric cartesian coordinates of the point, in meters,
        and return them as a 3-tuple '(x, y, z)'.

        'z' axis points to the North pole, 'x' axis is aligned with the
        reference meridian, whereas 'y' axis complements the system to a
        right-handed one."""

        latRadians = math.radians(self.latitude)
        longRadians = math.radians(self.longitude)

        # First consider the vertical cross section, passing through 'self',
        # the center of the ellipsoid, and the north pole. Let 'xy' and 'z' be
        # cartesian coordinates of 'self' with respect to that cross section.

        denominator = math.sqrt(
            self.ellipsoidMajorRadius**2 * math.cos(latRadians)**2 +
            self.ellipsoidMinorRadius**2 * math.sin(latRadians)**2)

        # Get coordinates without altitude taken into account.
        xy = self.ellipsoidMajorRadius**2 * math.cos(latRadians) / denominator
        z  = self.ellipsoidMinorRadius**2 * math.sin(latRadians) / denominator

        # Shift the point by 'self.altitude' meters along the normal.
        xy += self.altitude * math.cos(latRadians)
        z += self.altitude * math.sin(latRadians)

        # Rotate the horizontal component of the point by 'self.long'.
        x  = xy * math.cos(longRadians)
        y  = xy * math.sin(longRadians)

        return (x, y, z)

    def distToPoint(self, other):
        """Return the distance between 'self' and 'other' points, measured
        along a straight line, in meters.

        Warning: Altitudes of both points are taken into account."""

        selfCoords = self.cartesianCoords()
        otherCoords = other.cartesianCoords()

        return math.sqrt(sum((selfCoords[i] - otherCoords[i])**2
            for i in range(3)))

    def distFromLine(self, pt1, pt2):
        """Return the distance between 'self' point and a plane passing through
        'pt1', 'pt2', and Earth's center of mass, measured along the normal to
        the plane, in meters.

        The return value is positive if 'self' belongs to right-hand hemisphere
        with respect to the vector 'pt1'->'pt2' and negative otherwise.
        Exception is raised if points 'pt1' and 'pt2' coincide or are opposite
        to each other with respect to Earth's center of mass."""

        (x1, y1, z1) = pt1.cartesianCoords()
        (x2, y2, z2) = pt2.cartesianCoords()

        # Calculate the normal vector for the plane passing through 'pt1',
        # 'pt2', and the datum origin (don't normalize the vector length).
        xNorm = y2 * z1 - y1 * z2
        yNorm = z2 * x1 - z1 * x2
        zNorm = x2 * y1 - x1 * y2

        (x0, y0, z0) = self.cartesianCoords()

        # Return the distance between 'self' and the plane (this will raise
        # 'ZeroDivisionError' if 'pt1' and 'pt2' coincide or are opposite to
        # each other).
        return ((x0 * xNorm + y0 * yNorm + z0 * zNorm) /
            math.sqrt(xNorm**2 + yNorm**2 + zNorm**2))

    # ---- Public utility methods ---------------------------------------------
    def __eq__(self, other):
        return (self.latitide == other.latitude and
            self.longitude == other.longitude and
            self.altitude == other.altitude)

    def __hash__(self):
        return hash(self.latitude) ^ hash(self.longitude) ^ hash(self.altitude)