Geopoint

/// Geopoint class
///
/// Used for working with latitude/longitude and angles/distances. One instance
/// represents single geographic location.
class GeoPoint
{
public:
    double latitude;        /// Latitude on the spherical surface
    double longitude;       /// Longitude on the spherical surface

    /// Empty constructor
    GeoPoint()
    {
    }

    /// Constructor from latitude and longitude
    GeoPoint(double _lat, double _long)
    {
        latitude = (_lat * 3.14159265f / 180.0f);
        longitude = (_long * 3.14159265f / 180.0f);
    }

    /// Calculates distance between two geographic location, uses spherical geometry
    friend double distance(const GeoPoint& a, const GeoPoint& b)
    {
        double d = 2.0 * 6372795.477598 * asin(sqrt(pow(sin((b.latitude - a.latitude) / 2.0), 2.0) +
            cos(a.latitude) * cos(b.latitude) * pow(sin((b.longitude - a.longitude) / 2.0), 2.0)));
        return d;
    }

    /// Calculates angle from a to b related to equator
    friend double angle(const GeoPoint& a, const GeoPoint& b)
    {
        double z = log(tan(b.latitude / 2.0 + 3.14159265 / 4.0) / tan(a.latitude / 2.0 + 3.14159265 / 4.0));
        double x = a.longitude - b.longitude;
        if (x < 0.0)
        {
            x = -x;
        }
        return atan2(z, x);
    }

    /// Calculates direction from a to b (z being north-south, x being east-west)
    friend float4 direction(const GeoPoint& a, const GeoPoint& b)
    {
        float4 result = float4(0.0f);
        result.z = (float)log(tan(b.latitude / 2.0 + (double)M_PI / 4.0) / tan(a.latitude / 2.0 + (double)M_PI / 4.0));
        result.x = (float)(a.longitude - b.longitude);
        return result;
    }

    /// Calculates angle from a to b related to north
    friend double bearing(const GeoPoint& a, const GeoPoint& b)
    {
        double longitudeDiffRad = b.longitude - a.longitude;

        double y = sin(longitudeDiffRad) * cos(b.latitude);
        double x = cos(a.latitude)*sin(b.latitude) -
            sin(a.latitude)*cos(b.latitude)*cos(longitudeDiffRad);
        return atan2(y, x);
    }
};