class Location
{
    protected $longitude;
    protected $latitude;

    /**
     * Constructs a new location.
     *
     * @param float $lat
     * @param float $lon
     * @throws InvalidArgumentException
     */
    public function __construct($lat, $lon)
    {
        $this->setLatitude($lat);
        $this->setLongitude($lon);
    }

    public function bearingTo(Location $location)
    {
        $result = $this->distanceAndBearing($location);
        return $result['bearing'];
    }

    /**
     * Returns the approximate distance in meters between this location and the given location.
     *
     * @returns float
     */
    public function distanceTo(Location $location)
    {
        $result = $this->distanceAndBearing($location);
        return $result['distance'];
    }

    /**
     * Returns the latitude for this location.
     *
     * @return float
     */
    public function getLatitude()
    {
        return $this->latitude;
    }

    /**
     * Returns the longitude for this location.
     *
     * @return float
     */
    public function getLongitude()
    {
        return $this->longitude;
    }

    /**
     * Sets the latitude for this location.
     *
     * @param float $latitude
     * @throws InvalidArgumentException
     */
    public function setLatitude($latitude)
    {
        if (!is_numeric($latitude))
        {
            throw new InvalidArgumentException("Argument must be numeric.");
        }
        $this->latitude = (float) $latitude;
    }

    /**
     * Sets the longitude for this location.
     *
     * @param float $longitude
     * @throws InvalidArgumentException
     */
    public function setLongitude($longitude)
    {
        if (!is_numeric($longitude))
        {
            throw new InvalidArgumentException("Argument must be numeric.");
        }
        $this->longitude = (float) $longitude;
    }

    private function distanceAndBearing(Location $location)
    {
        // Based on http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
        // using the "Inverse Formula" (section 4)
        $lat1 = deg2rad($this->getLatitude());
        $lat2 = deg2rad($location->getLatitude());
        $lon1 = deg2rad($this->getLongitude());
        $lon2 = deg2rad($location->getLongitude());

        $major = 6378137.0;
        $minor = 6356752.3142;
        $f = ($major-$minor)/$major;

        $aSqMinusBSqOverBSq = (pow($major, 2) - pow($minor, 2))/($minor*$minor);
        $L = $lon2 - $lon1;
        $U1 = atan((1.0 - $f) * tan($lat1));
        $U2 = atan((1.0 - $f) * tan($lat2));

        $cosU1 = cos($U1);
        $cosU2 = cos($U2);
        $sinU1 = sin($U1);
        $sinU2 = sin($U2);
        $sinU1sinU2 = $sinU1 * $sinU2;
        $cosU1cosU2 = $cosU1 * $cosU2;

        $A = 0.0;
        $cos2SM = 0.0;
        $sinSigma = 0.0;
        $cosSigma = 0.0;
        $deltaSigma = 0.0;
        $sinLambda = 0.0;
        $cosLambda = 0.0;
        $sigma = 0.0;
        $lambda = $L;

        for ($i = 0; $i < 20; $i++)
        {
            $lambdaOrig = $lambda;
            $cosLambda = cos($lambda);
            $sinLambda = sin($lambda);
            $t1 = $cosU1 * $sinLambda;
            $t2 = $cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda;
            $sinSqSigma = pow($t1, 2) + pow($t2, 2);
            $sinSigma = sqrt($sinSqSigma);
            $sinAlpha = ($sinSigma == 0) ? 0.0 : $cosU1cosU2 * $sinLambda / $sinSigma;
            $cosSqAlpha = 1.0 - pow($sinAlpha, 2);
            $uSq = $cosSqAlpha * $aSqMinusBSqOverBSq;

            $A = 1 + ($uSq / 16384.0) * (4096.0 + $uSq * (-768 + $uSq * (320.0 - 175.0 * $uSq)));
            $B = ($uSq / 1024.0) * (256.0 + $uSq * (-128.0 + $uSq * (74.0 - 47.0 * $uSq)));
            $C = ($f / 16.0) * $cosSqAlpha * (4.0 + $f * (4.0 - 3.0 * $cosSqAlpha));

            $cosSigma = $sinU1sinU2 + $cosU1cosU2 * $cosLambda;
            $cos2SM = ($cosSqAlpha == 0) ? 0.0 : $cosSigma - 2.0 * $sinU1sinU2 / $cosSqAlpha;

            $sigma = atan2($sinSigma, $cosSigma);
            $cos2SMSq = pow($cos2SM, 2);
            $deltaSigma = $B * $sinSigma * ($cos2SM + ($B / 4.0) * ($cosSigma * (-1.0 + 2.0 * $cos2SMSq) - ($B / 6.0) * $cos2SM * (-3.0 + 4.0 * pow($sinSigma, 2)) * (-3.0 + 4.0 * $cos2SMSq)));

            $lambda = $L + (1.0 - $C) * $f * $sinAlpha * ($sigma + $C * $sinSigma * ($cos2SM + $C * $cosSigma * (-1.0 + 2.0 * $cos2SM * $cos2SM)));
            $delta = ($lambda - $lambdaOrig) / $lambda;
            if (abs($delta) < 1.0e-12) break;
        }
        $distance = (float) ($minor * $A * ($sigma - $deltaSigma));
        $bearing = (float) rad2deg(atan2($cosU2 * $sinLambda, $cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda));

        return array(
            'distance' => $distance,
            'bearing' => $bearing
        );
    }
}