/src/php/Bairwell/Geocoder/LatLon.php

<?php
/**
 * Latitude/longitude spherical geodesy formulae & scripts
 *
 * @package Bairwell
 * @subpackage Geocoder
 * @author Richard Chiswell <richard@bairwell.com>
 * @copyright 2012 Bairwell Ltd
 * @license MIT
 */
namespace Bairwell\Geocoder;

/**
 * Latitude/longitude spherical geodesy formulae & scripts.
 *
 * Based on Javascript code written by and (c) Chris Veness 2002-2012
 * http://www.movable-type.co.uk/scripts/latlon.js
 * http://www.movable-type.co.uk/scripts/latlong.html
 */
class LatLon
{

    /**
     * Latitude in numeric degrees
     * @var float
     */
    private $lat;

    /**
     * Longitude in numeric degrees
     * @var float
     */
    private $lon;

    /**
     * Radius of earth
     * @var int
     */
    private $radius;

    /**
     * Creates a point on the earth's surface at the supplied latitude / longitude.
     *
     * @param float $lat latitude in numeric degrees
     * @param float $lon longitude in numeric degrees
     * @param float $rad radius of earth if different value is required from UGC standard 6,371km
     */
    public function __construct($lat = 0.0, $lon = 0.0, $rad = 6371.0)
    {
        if (TRUE === is_numeric($lat)) {
            $this->setLat(floatval($lat));
        }
        if (TRUE === is_numeric($lon)) {
            $this->setLon(floatval($lon));
        }
        if (TRUE === is_numeric($rad)) {
            $this->setRadius(floatval($rad));
        }
    }

    /**
     * Returns the distance from this point to the supplied point, in km.
     * (using Haversine formula)
     *
     * from: Haversine formula - R. W. Sinnott, "Virtues of the Haversine",
     *       Sky and Telescope, vol 68, no 2, 1984
     *
     * @param   LatLon $point Latitude/longitude of destination point
     * @param   int $precision no of significant digits to use for returned value.  4 sig figs reflects typical 0.3% accuracy of spherical model
     * @return string Distance in km between this point and destination point
     */
    public function distanceTo(LatLon $point, $precision = 4)
    {
        $R = $this->getRadius();
        $lat1 = deg2rad($this->getLat());
        $lon1 = deg2rad($this->getLon());
        $lat2 = deg2rad($point->getLat());
        $lon2 = deg2rad($point->getLon());
        $dLat = $lat2 - $lat1;
        $dLon = $lon2 - $lon1;

        $a = sin($dLat / 2) * sin($dLat / 2) +
            cos($lat1) * cos($lat2) *
                sin($dLon / 2) * sin($dLon / 2);
        $c = 2 * atan2(sqrt($a), sqrt(1 - $a));
        $d = $R * $c;
        return FixedPrecision::toPrecisionFixed($d, $precision);
    }

    /**
     * Returns the (initial) bearing from this point to the supplied point, in degrees.
     *
     *   see http://williams.best.vwh.net/avform.htm#Crs
     *
     * @param   LatLon $point Latitude/longitude of destination point
     * @return float Initial bearing in degrees from North
     */
    public function bearingTo(LatLon $point)
    {
        $lat1 = deg2rad($this->getLat());
        $lat2 = deg2rad($point->getLat());
        $dLon = deg2rad($point->getLon() - $this->getLon());

        $y = sin($dLon) * cos($lat2);
        $x = cos($lat1) * sin($lat2) -
            sin($lat1) * cos($lat2) * cos($dLon);
        $brng = atan2($y, $x);

        return fmod(rad2deg($brng) + 360, 360);
    }

    /**
     * Returns final bearing arriving at supplied destination point from this point.
     *
     * The final bearing will differ from the initial bearing by varying degrees
     * according to distance and latitude
     *
     * @param   LatLon $point Latitude/longitude of destination point
     * @return float Final bearing in degrees from North
     */
    public function finalBearingTo(LatLon $point)
    {
        // get initial bearing from supplied point back to this point...
        $lat1 = deg2rad($point->getLat());
        $lat2 = deg2rad($this->getLat());
        $dLon = deg2rad($this->getLon() - $point->getLon());

        $y = sin($dLon) * cos($lat2);
        $x = cos($lat1) * sin($lat2) -
            sin($lat1) * cos($lat2) * cos($dLon);
        $brng = atan2($y, $x);

        // ... & reverse it by adding 180°
        return fmod(rad2deg($brng) + 180, 360);
    }

    /**
     * Returns the midpoint between this point and the supplied point.
     *
     *   see http://mathforum.org/library/drmath/view/51822.html for derivation
     *
     * @param   LatLon $point Latitude/longitude of destination point
     * @return LatLon Midpoint between this point and the supplied point
     */
    public function midpointTo(LatLon $point)
    {
        $lat1 = deg2rad($this->getLat());
        $lon1 = deg2rad($this->getLon());
        $lat2 = deg2rad($point->getLat());
        $dLon = deg2rad($point->getLon() - $this->getLon());

        $Bx = cos($lat2) * cos($dLon);
        $By = cos($lat2) * sin($dLon);

        $lat3 = atan2(
            sin($lat1) + sin($lat2),
            sqrt((cos($lat1) + $Bx) * (cos($lat1) + $Bx) + $By * $By)
        );

        $lon3 = $lon1 + atan2($By, cos($lat1) + $Bx);
        $lon3 = fmod(($lon3 + 3 * pi()), (2 * pi()) - pi()); // normalise to -180..+180º
        return new \Bairwell\Geocoder\LatLon(rad2deg($lat3), rad2deg($lon3));
    }

    /**
     * Returns the destination point from this point having travelled the given distance.
     *
     * (in km) on the given initial bearing (bearing may vary before destination is reached)
     *
     *   see http://williams.best.vwh.net/avform.htm#LL
     *
     * @param   float $brng Initial bearing in degrees
     * @param   float $dist Distance in km
     * @return LatLon Destination point
     */
    public function destinationPoint($brng, $dist)
    {
        $dist = $dist / $this->getRadius(); // convert dist to angular distance in radians
        $brng = deg2rad($brng);
        $lat1 = deg2rad($this->getLat());
        $lon1 = deg2rad($this->getLon());

        $lat2 = asin(
            sin($lat1) * cos($dist) +
                cos($lat1) * sin($dist) * cos($brng)
        );
        $lon2 = $lon1 + atan2(
            sin($brng) * sin($dist) * cos($lat1),
            cos($dist) - sin($lat1) * sin($lat2)
        );
        // normalise to -180..+180º
        $lon2 = fmod(($lon2 + 3 * pi()), (2 * pi()) - pi());

        return new \Bairwell\Geocoder\LatLon(rad2deg($lat2), rad2deg($lon2));
    }

    /**
     * Returns the point of intersection of two paths defined by point and bearing.
     *
     *   see http://williams.best.vwh.net/avform.htm#Intersection
     *
     * @param   LatLon $p1 First point
     * @param   float $brng1 Initial bearing from first point
     * @param   LatLon $p2 Second point
     * @param   float $brng2 Initial bearing from second point
     * @return LatLon Destination point (null if no unique intersection defined)
     */
    public function intersection(LatLon $p1, $brng1, LatLon $p2, $brng2)
    {
        $lat1 = deg2rad($p1->getLat());
        $lon1 = deg2rad($p1->getLon());
        $lat2 = deg2rad($p2->getLat());
        $lon2 = deg2rad($p2->getLon());
        $brng13 = deg2rad($brng1);
        $brng23 = deg2rad($brng2);
        $dLat = $lat2 - $lat1;
        $dLon = $lon2 - $lon1;

        $dist12 = 2 * asin(
            sqrt(
                sin($dLat / 2) * sin($dLat / 2) +
                    cos($lat1) * cos($lat2) * sin($dLon / 2) * sin($dLon / 2)
            )
        );
        if ($dist12 === 0) {
            return NULL;
        }

        // initial/final bearings between points
        $brngA = acos(
            (sin($lat2) - sin($lat1) * cos($dist12)) /
                (sin($dist12) * cos($lat1))
        );
        if (FALSE === is_numeric($brngA)) {
            // protect against rounding
            $brngA = 0;
        }
        $brngB = acos(
            (sin($lat1) - sin($lat2) * cos($dist12)) /
                (sin($dist12) * cos($lat2))
        );

        if (sin($lon2 - $lon1) > 0) {
            $brng12 = $brngA;
            $brng21 = 2 * pi() - $brngB;
        } else {
            $brng12 = 2 * pi() - $brngA;
            $brng21 = $brngB;
        }
        $alpha1 = fmod(($brng13 - $brng12 + pi()), (2 * pi())) - pi(); // angle 2-1-3
        $alpha2 = fmod(($brng21 - $brng23 + pi()), (2 * pi())) - pi(); // angle 1-2-3

        if (sin($alpha1) === 0 && sin($alpha2) === 0) {
            // infinite intersections
            return NULL;
        }
        if (sin($alpha1) * sin($alpha2) < 0) {
            // ambiguous intersection
            return NULL;
        }

        //alpha1 = Math.abs(alpha1);
        //alpha2 = Math.abs(alpha2);
        // ... Ed Williams takes abs of alpha1/alpha2, but seems to break calculation?

        $alpha3 = acos(
            -cos($alpha1) * cos($alpha2) +
                sin($alpha1) * sin($alpha2) * cos($dist12)
        );
        $dist13 = atan2(
            sin($dist12) * sin($alpha1) * sin($alpha2),
            cos($alpha2) + cos($alpha1) * cos($alpha3)
        );
        $lat3 = asin(
            sin($lat1) * cos($dist13) +
                cos($lat1) * sin($dist13) * cos($brng13)
        );
        $dLon13 = atan2(
            sin($brng13) * sin($dist13) * cos($lat1),
            cos($dist13) - sin($lat1) * sin($lat3)
        );
        $lon3 = $lon1 + $dLon13;
        $lon3 = fmod(($lon3 + 3 * pi()), (2 * pi())) - pi(); // normalise to -180..+180º

        return new \Bairwell\Geocoder\LatLon(rad2deg($lat3), rad2deg($lon3));
    }

    /**
     * Returns the distance from this point to the supplied point, in km, travelling along a rhumb line
     *
     *   see http://williams.best.vwh.net/avform.htm#Rhumb
     *
     * @param   LatLon $point Latitude/longitude of destination point
     * @return float Distance in km between this point and destination point
     */
    public function rhumbDistanceTo(LatLon $point)
    {
        $R = $this->getRadius();
        $lat1 = deg2rad($this->getLat());
        $lat2 = deg2rad($point->getLat());
        $dLat = deg2rad($point->getLat() - $this->getLat());
        $dLon = deg2rad(abs($point->getLon() - $this->getLon()));
        $dPhi = log(tan($lat2 / 2 + pi() / 4) / tan($lat1 / 2 + pi() / 4));
        // E-W line gives dPhi=0
        if (FALSE === is_nan($dLat / $dPhi)) {
            $q = $dLat / $dPhi;
        } else {
            $q = cos($lat1);
        }
        // if dLon over 180° take shorter rhumb across 180° meridian:
        if ($dLon > pi()) {
            $dLon = 2 * pi() - $dLon;
        }
        $dist = sqrt($dLat * $dLat + $q * $q * $dLon * $dLon) * $R;
        return FixedPrecision::toPrecisionFixed($dist, 4); // 4 sig figs reflects typical 0.3% accuracy of spherical model
    }

    /**
     * Returns the bearing from this point to the supplied point along a rhumb line, in degrees
     *
     * @param   LatLon $point Latitude/longitude of destination point
     * @return float Bearing in degrees from North
     */
    public function rhumbBearingTo(LatLon $point)
    {
        $lat1 = deg2rad($this->getLat());
        $lat2 = deg2rad($point->getLat());
        $dLon = deg2rad($point->getLon() - $this->getLon());

        $dPhi = log(tan($lat2 / 2 + pi() / 4) / tan($lat1 / 2 + pi() / 4));
        if (abs($dLon) > pi()) {
            if ($dLon>0) {
                $dLon=-(2 * pi() - $dLon);
            } else {
                $dLon=(2 * pi() + $dLon);
            }
        }
        $brng = atan2($dLon, $dPhi);

        return fmod(rad2deg($brng) + 360, 360);
    }

    /**
     * Returns the destination point from this point having travelled the given distance (in km) on the
     * given bearing along a rhumb line
     *
     * @param   float $brng Bearing in degrees from North
     * @param   float $dist Distance in km
     * @return LatLon Destination point
     */
    public function rhumbDestinationPoint($brng, $dist)
    {
        $R = $this->getRadius();
        $d = floatval($dist) / $R; // d = angular distance covered on earth's surface
        $lat1 = deg2rad($this->getLat());
        $lon1 = deg2rad($this->getLon());
        $brng = deg2rad($brng);

        $lat2 = $lat1 + $d * cos($brng);
        $dLat = $lat2 - $lat1;
        $dPhi = log(tan($lat2 / 2 + pi() / 4) / tan($lat1 / 2 + pi() / 4));
        // E-W line gives dPhi=0
        if (FALSE === is_nan($dLat / $dPhi)) {
            $q = $dLat / $dPhi;
        } else {
            $q = cos($lat1);
        }
        $dLon = $d * sin($brng) / $q;
        // check for some daft bugger going past the pole
        if (abs($lat2) > pi() / 2) {
            if ($lat2>0) {
                $lat2=pi() - $lat2;
            } else {
                $lat2=-pi() - $lat2;
            }
        }
        $lon2 = fmod(($lon1 + $dLon + 3 * pi()), (2 * pi())) - pi();

        return new \Bairwell\Geocoder\LatLon(rad2deg($lat2), rad2deg($lon2));
    }

    /**
     * Sets latitude in a fluent manner
     *
     * @param float $lat
     * @return LatLon
     */
    public function setLat($lat)
    {
        $this->lat = $lat;
        return $this;
    }

    /**
     * Gets the latitude
     * @return float
     */
    public function getLat()
    {
        return $this->lat;
    }

    /**
     * Sets longitude in a fluent manner
     *
     * @param float $lon
     * @return LatLon
     */
    public function setLon($lon)
    {
        $this->lon = $lon;
        return $this;
    }

    /**
     * Gets the longitude
     * @return float
     */
    public function getLon()
    {
        return $this->lon;
    }

    /**
     * Sets radius in a fluent manner
     *
     * @param int $radius
     * @return LatLon
     */
    public function setRadius($radius)
    {
        $this->radius = $radius;
        return $this;
    }

    /**
     * Gets the radius
     * @return int
     */
    public function getRadius()
    {
        return $this->radius;
    }


}