Calculate endpoint given distance, bearing, starting point

Question

I am trying to find the destination point, given a starting point lat/long, bearing & distance. The calculator from this website below gives me the desired results.

http://www.movable-type.co.uk/scripts/latlong.html

When I try to implement the same through code, I don't get the right results.

Below is my code -

    private  GLatLng pointRadialDistance(double lat1, double lon1,
               double radianBearing, double radialDistance)
    {
        double rEarth = 6371.01;
        lat1 = DegreeToRadian(lat1);
        lon1 = DegreeToRadian(lon1);
        radianBearing = DegreeToRadian(radianBearing);
        radialDistance = radialDistance / rEarth;
        double lat = Math.Asin(Math.Sin(lat1) * Math.Cos(radialDistance) + Math.Cos(lat1)
                        * Math.Sin(radialDistance) * Math.Cos(radianBearing));
        double lon;
        if (Math.Cos(lat) == 0)
        {  // Endpoint a pole 
            lon = lon1;
        }
        else
        {
            lon = ((lon1 - Math.Asin(Math.Sin(radianBearing) * Math.Sin(radialDistance) / Math.Cos(lat))
                            + Math.PI) % (2 * Math.PI)) - Math.PI;
        }
        lat = RadianToDegree(lat);
        lon = RadianToDegree(lon);
        GLatLng newLatLng = new GLatLng(lat, lon);
        return newLatLng;
    }

    public  double Bearing(double lat1, double long1, double lat2, double long2)
    {
        //Convert input values to radians   
        lat1 = DegreeToRadian(lat1);
        long1 = DegreeToRadian(long1);
        lat2 = DegreeToRadian(lat2);
        long2 = DegreeToRadian(long2);

        double deltaLong = long2 - long1;

        double y = Math.Sin(deltaLong) * Math.Cos(lat2);
        double x = Math.Cos(lat1) * Math.Sin(lat2) -
                Math.Sin(lat1) * Math.Cos(lat2) * Math.Cos(deltaLong);
        double bearing = Math.Atan2(y, x);
        return bearing;
    }   

   public double DegreeToRadian(double angle)
    {
    return Math.PI * angle / 180.0;
    }

    public double RadianToDegree(double angle)
    {
        return 180.0 * angle / Math.PI;
    }

From the main program, I call the sub procedures as follows -

double bearing = Bearing(-41.294444, 174.814444, -40.90521, 175.6604);
GLatLng endLatLng = pointRadialDistance(-41.294444, 174.814444, bearing, 80);

I get below results -

Bearing=1.02749621782165
endLatLng=-40.5751022737927,174.797458881699

The answer I expect is -40.939722,175.646389 (from website link above).

Can anyone suggest what mistake I am making in the code here?

Answer 1

Below is the implementation of the JavaScript code at http://www.movable-type.co.uk/scripts/latlong.html which I wrote for myself and use in my own projects. You can implement it to your project if you will.

Note: Coordinate is a class with X (longitude), Y (latitude), Z (altitude) properties. ToDegree() and ToRadian() are extensions for Double type. Finally, GetTarget() is an extension for a Coordinate instance.

/// <summary>Calculates the destination coordinate by given angle and distance.</summary>
/// <param name="origin">Origin.</param>
/// <param name="bearing">Azimuth.</param>
/// <param name="distance">Distance (km).</param>
/// <returns>Coordinate.</returns>
public static Coordinate GetTarget(
this Coordinate origin, double bearing, double distance, double altitude = 0)
{
    var d = distance / 6371;
    var rlat = origin.Y.ToRadian();
    var rlon = origin.X.ToRadian();
    var rbearing = bearing.ToRadian();
    var lat2 = rlat + (d * Math.Cos(rbearing));
    var dlat = lat2 - rlat;
    var dphi = Math.Log((Math.Tan((lat2 / 2) + (Math.PI / 4))) / (Math.Tan((rlat / 2) + (Math.PI / 4))));
    var q =
        Math.Abs(dlat) > 0.0000000001
        ? dlat / dphi
        : Math.Cos(rlat);
    var dlon = (d * Math.Sin(rbearing)) / q;

    if (Math.Abs(lat2) > Math.PI / 2)
    {
        lat2 = lat2 > 0 ? Math.PI : Math.PI - lat2;
    }

    var lon2 = (rlon + dlon + Math.PI) % (2 * Math.PI) - Math.PI;

    return new Coordinate
    {
        X = lon2.ToDegree(),
        Y = lat2.ToDegree(),
        Z = origin.Z
    };
}