diff --git a/MapControl/Shared/AzimuthalEquidistantProjection.cs b/MapControl/Shared/AzimuthalEquidistantProjection.cs
index 67c9d458..86d62c9b 100644
--- a/MapControl/Shared/AzimuthalEquidistantProjection.cs
+++ b/MapControl/Shared/AzimuthalEquidistantProjection.cs
@@ -31,7 +31,7 @@ namespace MapControl
return new Point();
}
- GetAzimuthDistance(Center, location, out double azimuth, out double distance);
+ Center.GetAzimuthDistance(location, out double azimuth, out double distance);
var mapDistance = distance * Wgs84EquatorialRadius;
@@ -50,7 +50,7 @@ namespace MapControl
var distance = mapDistance / Wgs84EquatorialRadius;
- return GetLocation(Center, azimuth, distance);
+ return Center.GetLocation(azimuth, distance);
}
}
}
diff --git a/MapControl/Shared/AzimuthalProjection.cs b/MapControl/Shared/AzimuthalProjection.cs
index 62a39f19..e63ee905 100644
--- a/MapControl/Shared/AzimuthalProjection.cs
+++ b/MapControl/Shared/AzimuthalProjection.cs
@@ -1,5 +1,4 @@
-using System;
-#if WPF
+#if WPF
using System.Windows;
#endif
@@ -46,55 +45,5 @@ namespace MapControl
return boundingBox;
}
-
- ///
- /// Calculates azimuth and spherical distance in radians from location1 to location2.
- /// The returned distance has to be multiplied with an appropriate earth radius.
- ///
- public static void GetAzimuthDistance(Location location1, Location location2, out double azimuth, out double distance)
- {
- var lat1 = location1.Latitude * Math.PI / 180d;
- var lon1 = location1.Longitude * Math.PI / 180d;
- var lat2 = location2.Latitude * Math.PI / 180d;
- var lon2 = location2.Longitude * Math.PI / 180d;
- var cosLat1 = Math.Cos(lat1);
- var sinLat1 = Math.Sin(lat1);
- var cosLat2 = Math.Cos(lat2);
- var sinLat2 = Math.Sin(lat2);
- var cosLon12 = Math.Cos(lon2 - lon1);
- var sinLon12 = Math.Sin(lon2 - lon1);
- var cosDistance = sinLat1 * sinLat2 + cosLat1 * cosLat2 * cosLon12;
-
- azimuth = Math.Atan2(sinLon12, cosLat1 * sinLat2 / cosLat2 - sinLat1 * cosLon12);
- distance = Math.Acos(Math.Min(Math.Max(cosDistance, -1d), 1d));
- }
-
- ///
- /// Calculates the Location of the point given by azimuth and spherical distance in radians from location.
- ///
- public static Location GetLocation(Location location, double azimuth, double distance)
- {
- var lat = location.Latitude;
- var lon = location.Longitude;
-
- if (distance > 0d)
- {
- var lat1 = lat * Math.PI / 180d;
- var sinDistance = Math.Sin(distance);
- var cosDistance = Math.Cos(distance);
- var cosAzimuth = Math.Cos(azimuth);
- var sinAzimuth = Math.Sin(azimuth);
- var cosLat1 = Math.Cos(lat1);
- var sinLat1 = Math.Sin(lat1);
- var sinLat2 = sinLat1 * cosDistance + cosLat1 * sinDistance * cosAzimuth;
- var lat2 = Math.Asin(Math.Min(Math.Max(sinLat2, -1d), 1d));
- var dLon = Math.Atan2(sinDistance * sinAzimuth, cosLat1 * cosDistance - sinLat1 * sinDistance * cosAzimuth);
-
- lat = lat2 * 180d / Math.PI;
- lon += dLon * 180d / Math.PI;
- }
-
- return new Location(lat, lon);
- }
}
}
diff --git a/MapControl/Shared/GnomonicProjection.cs b/MapControl/Shared/GnomonicProjection.cs
index 3e826200..26781707 100644
--- a/MapControl/Shared/GnomonicProjection.cs
+++ b/MapControl/Shared/GnomonicProjection.cs
@@ -31,7 +31,7 @@ namespace MapControl
return new Point();
}
- GetAzimuthDistance(Center, location, out double azimuth, out double distance);
+ Center.GetAzimuthDistance(location, out double azimuth, out double distance);
var mapDistance = distance < Math.PI / 2d
? Math.Tan(distance) * Wgs84EquatorialRadius
@@ -52,7 +52,7 @@ namespace MapControl
var distance = Math.Atan(mapDistance / Wgs84EquatorialRadius);
- return GetLocation(Center, azimuth, distance);
+ return Center.GetLocation(azimuth, distance);
}
}
}
diff --git a/MapControl/Shared/Location.cs b/MapControl/Shared/Location.cs
index 29ab2386..cc1475c0 100644
--- a/MapControl/Shared/Location.cs
+++ b/MapControl/Shared/Location.cs
@@ -5,6 +5,10 @@ namespace MapControl
{
///
/// A geographic location with latitude and longitude values in degrees.
+ /// For calculations with azimuth and distance on great circles, see
+ /// https://en.wikipedia.org/wiki/Great_circle
+ /// https://en.wikipedia.org/wiki/Great-circle_distance
+ /// https://en.wikipedia.org/wiki/Great-circle_navigation
///
#if UWP || WINUI
[Windows.Foundation.Metadata.CreateFromString(MethodName = "Parse")]
@@ -81,54 +85,63 @@ namespace MapControl
}
///
- /// Calculates the great circle distance between this and the specified Location.
- /// https://en.wikipedia.org/wiki/Great_circle
- /// https://en.wikipedia.org/wiki/Great-circle_distance
- /// https://en.wikipedia.org/wiki/Great-circle_navigation
+ /// Calculates great circle azimuth and distance in radians between this and the specified Location.
///
- public double GetDistance(
- Location location, double earthRadius = MapProjection.Wgs84EquatorialRadius)
+ public void GetAzimuthDistance(Location location, out double azimuth, out double distance)
{
var lat1 = Latitude * Math.PI / 180d;
var lon1 = Longitude * Math.PI / 180d;
var lat2 = location.Latitude * Math.PI / 180d;
var lon2 = location.Longitude * Math.PI / 180d;
- var sinLat1 = Math.Sin(lat1);
var cosLat1 = Math.Cos(lat1);
- var sinLat2 = Math.Sin(lat2);
+ var sinLat1 = Math.Sin(lat1);
var cosLat2 = Math.Cos(lat2);
- var sinLon12 = Math.Sin(lon2 - lon1);
+ var sinLat2 = Math.Sin(lat2);
var cosLon12 = Math.Cos(lon2 - lon1);
- var a = cosLat1 * sinLat2 - sinLat1 * cosLat2 * cosLon12;
- var b = cosLat2 * sinLon12;
- var c = sinLat1 * sinLat2 + cosLat1 * cosLat2 * cosLon12;
- var s12 = Math.Atan2(Math.Sqrt(a * a + b * b), c);
-
- return earthRadius * s12;
+ var sinLon12 = Math.Sin(lon2 - lon1);
+ var a = cosLat2 * sinLon12;
+ var b = cosLat1 * sinLat2 - sinLat1 * cosLat2 * cosLon12;
+ // α1
+ azimuth = Math.Atan2(a, b);
+ // σ12
+ distance = Math.Atan2(Math.Sqrt(a * a + b * b), sinLat1 * sinLat2 + cosLat1 * cosLat2 * cosLon12);
}
///
- /// Calculates the Location on a great circle at the specified azimuth angle and distance from this Location.
- /// https://en.wikipedia.org/wiki/Great_circle
- /// https://en.wikipedia.org/wiki/Great-circle_navigation
+ /// Calculates the great circle distance in meters between this and the specified Location.
///
- public Location GetLocation(
- double azimuth, double distance, double earthRadius = MapProjection.Wgs84EquatorialRadius)
+ public double GetDistance(Location location, double earthRadius = MapProjection.Wgs84EquatorialRadius)
+ {
+ GetAzimuthDistance(location, out double _, out double distance);
+
+ return earthRadius * distance;
+ }
+
+ ///
+ /// Calculates the Location on a great circle at the specified azimuth and distance in radians from this Location.
+ ///
+ public Location GetLocation(double azimuth, double distance)
{
- var s12 = distance / earthRadius;
- var az1 = azimuth * Math.PI / 180d;
var lat1 = Latitude * Math.PI / 180d;
var lon1 = Longitude * Math.PI / 180d;
- var sinS12 = Math.Sin(s12);
- var cosS12 = Math.Cos(s12);
- var sinAz1 = Math.Sin(az1);
- var cosAz1 = Math.Cos(az1);
- var sinLat1 = Math.Sin(lat1);
+ var cosD = Math.Cos(distance);
+ var sinD = Math.Sin(distance);
+ var cosA = Math.Cos(azimuth);
+ var sinA = Math.Sin(azimuth);
var cosLat1 = Math.Cos(lat1);
- var lat2 = Math.Asin(sinLat1 * cosS12 + cosLat1 * sinS12 * cosAz1);
- var lon2 = lon1 + Math.Atan2(sinS12 * sinAz1, cosLat1 * cosS12 - sinLat1 * sinS12 * cosAz1);
+ var sinLat1 = Math.Sin(lat1);
+ var lat2 = Math.Asin(sinLat1 * cosD + cosLat1 * sinD * cosA);
+ var lon2 = lon1 + Math.Atan2(sinD * sinA, cosLat1 * cosD - sinLat1 * sinD * cosA);
return new Location(lat2 * 180d / Math.PI, lon2 * 180d / Math.PI);
}
+
+ ///
+ /// Calculates the Location on a great circle at the specified azimuth in degrees and distance in meters from this Location.
+ ///
+ public Location GetLocation(double azimuth, double distance, double earthRadius = MapProjection.Wgs84EquatorialRadius)
+ {
+ return GetLocation(azimuth * Math.PI / 180d, distance / earthRadius);
+ }
}
}
diff --git a/MapControl/Shared/LocationCollection.cs b/MapControl/Shared/LocationCollection.cs
index 8ffc5f6d..90c89274 100644
--- a/MapControl/Shared/LocationCollection.cs
+++ b/MapControl/Shared/LocationCollection.cs
@@ -67,31 +67,29 @@ namespace MapControl
///
/// Calculates a series of Locations on a great circle, or orthodrome, that connects the two specified Locations,
/// with an optional angular resolution specified in degrees.
- ///
- /// See https://en.wikipedia.org/wiki/Great-circle_navigation
+ /// https://en.wikipedia.org/wiki/Great-circle_navigation
///
public static LocationCollection OrthodromeLocations(Location location1, Location location2, double resolution = 1d)
{
if (resolution <= 0d)
{
throw new ArgumentOutOfRangeException(
- nameof(resolution), "The resolution argument must be greater than zero.");
+ nameof(resolution), $"The {nameof(resolution)} argument must be greater than zero.");
}
var lat1 = location1.Latitude * Math.PI / 180d;
var lon1 = location1.Longitude * Math.PI / 180d;
var lat2 = location2.Latitude * Math.PI / 180d;
var lon2 = location2.Longitude * Math.PI / 180d;
-
var cosLat1 = Math.Cos(lat1);
var sinLat1 = Math.Sin(lat1);
var cosLat2 = Math.Cos(lat2);
var sinLat2 = Math.Sin(lat2);
var cosLon12 = Math.Cos(lon2 - lon1);
var sinLon12 = Math.Sin(lon2 - lon1);
-
- var a = cosLat1 * sinLat2 - sinLat1 * cosLat2 * cosLon12;
- var b = cosLat2 * sinLon12;
+ var a = cosLat2 * sinLon12;
+ var b = cosLat1 * sinLat2 - sinLat1 * cosLat2 * cosLon12;
+ // σ12
var s12 = Math.Atan2(Math.Sqrt(a * a + b * b), sinLat1 * sinLat2 + cosLat1 * cosLat2 * cosLon12);
var n = (int)Math.Ceiling(s12 / resolution * 180d / Math.PI); // s12 in radians
@@ -100,15 +98,18 @@ namespace MapControl
if (n > 1)
{
- var az1 = Math.Atan2(sinLon12, cosLat1 * sinLat2 / cosLat2 - sinLat1 * cosLon12);
+ // https://en.wikipedia.org/wiki/Great-circle_navigation#Finding_way-points
+ // α1
+ var az1 = Math.Atan2(a, b);
var cosAz1 = Math.Cos(az1);
var sinAz1 = Math.Sin(az1);
-
+ // α0
var az0 = Math.Atan2(sinAz1 * cosLat1, Math.Sqrt(cosAz1 * cosAz1 + sinAz1 * sinAz1 * sinLat1 * sinLat1));
- var sinAz0 = Math.Sin(az0);
var cosAz0 = Math.Cos(az0);
-
+ var sinAz0 = Math.Sin(az0);
+ // σ01
var s01 = Math.Atan2(sinLat1, cosLat1 * cosAz1);
+ // λ0
var lon0 = lon1 - Math.Atan2(sinAz0 * Math.Sin(s01), Math.Cos(s01));
for (var i = 1; i < n; i++)
@@ -139,7 +140,7 @@ namespace MapControl
if (resolution <= 0d)
{
throw new ArgumentOutOfRangeException(
- nameof(resolution), "The resolution argument must be greater than zero.");
+ nameof(resolution), $"The {nameof(resolution)} argument must be greater than zero.");
}
var lat1 = location1.Latitude;
@@ -153,13 +154,13 @@ namespace MapControl
if (double.IsInfinity(y1))
{
throw new ArgumentOutOfRangeException(
- nameof(location1), "The location1 argument must have an absolute latitude value of less than 90.");
+ nameof(location1), $"The {nameof(location1)} argument must have an absolute latitude value of less than 90.");
}
if (double.IsInfinity(y2))
{
throw new ArgumentOutOfRangeException(
- nameof(location2), "The location2 argument must have an absolute latitude value of less than 90.");
+ nameof(location2), $"The {nameof(location2)} argument must have an absolute latitude value of less than 90.");
}
var dlat = lat2 - lat1;
diff --git a/MapControl/Shared/StereographicProjection.cs b/MapControl/Shared/StereographicProjection.cs
index 59f40528..7b680417 100644
--- a/MapControl/Shared/StereographicProjection.cs
+++ b/MapControl/Shared/StereographicProjection.cs
@@ -31,7 +31,7 @@ namespace MapControl
return new Point();
}
- GetAzimuthDistance(Center, location, out double azimuth, out double distance);
+ Center.GetAzimuthDistance(location, out double azimuth, out double distance);
var mapDistance = Math.Tan(distance / 2d) * 2d * Wgs84EquatorialRadius;
@@ -50,7 +50,7 @@ namespace MapControl
var distance = 2d * Math.Atan(mapDistance / (2d * Wgs84EquatorialRadius));
- return GetLocation(Center, azimuth, distance);
+ return Center.GetLocation(azimuth, distance);
}
}
}