I have made a sample project using dot42 and have taken this C# code and placed it inside a library:
var gcd = GetGreatCircleDistanceKm(52.0, -1.90, 21.0, 39.0); // returns NaN should be ~4915
public double GetGreatCircleDistanceKm(double startLat, double startLong, double endLat, double endLong)
{
var earthRadius = Constants.EarthRadiusKm;
var φ1 = DegreesToRadians(startLat);
var φ2 = DegreesToRadians(endLat);
var Δφ = DegreesToRadians(endLat - startLat);
var Δλ = DegreesToRadians(endLong - startLong);
var a = Math.Sin(Δφ / 2) * Math.Sin(Δφ / 2) + Math.Cos(φ1) * Math.Cos(φ2) * Math.Sin(Δλ / 2) * Math.Sin(Δλ / 2);
var c = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));
return earthRadius * c;
}
public static double DegreesToRadians(double angle)
{
return Math.PI * angle / 180.0;
}
problem is however that c
always returns NaN
and thus the function returns NaN
. This same method works fine in a normal .net project.
I have tried many versions of the same haversine
formula thinking I may have mistyped something, but no, all return NaN
.
Source of formula : http://www.movable-type.co.uk/scripts/latlong.html
The solution was (strangely) to do :
var a1 = Δφ/2;
var a = Math.Sin(a1) * Math.Sin(a1) + Math.Cos(φ1) * Math.Cos(φ2) * Math.Sin(Δλ / 2) * Math.Sin(Δλ / 2);
No idea how it works, but it works. Please if someone could shed light.
Full fixed version:
/// <summary>
/// Gets the shortest possible distance between two points on earth
/// </summary>
/// <param name="startLat"></param>
/// <param name="startLong"></param>
/// <param name="endLat"></param>
/// <param name="endLong"></param>
/// <returns></returns>
public double GetGreatCircleDistanceKm(double startLat, double startLong, double endLat, double endLong)
{
var startLatRad = DegreesToRadians(startLat);
var endLatRad = DegreesToRadians(endLat);
var latDiffRad = DegreesToRadians(endLat - startLat);
var longDiffRad = DegreesToRadians(endLong - startLong);
var halfLatDiff = latDiffRad/2;
var a = Math.Sin(halfLatDiff) * Math.Sin(halfLatDiff) + Math.Cos(startLatRad) * Math.Cos(endLatRad) * Math.Sin(longDiffRad / 2) * Math.Sin(longDiffRad / 2);
var c = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));
return Constants.EarthRadiusKm * c;
}