diff --git a/MapControl/Shared/PolarStereographicProjection.cs b/MapControl/Shared/PolarStereographicProjection.cs index c3a9ee34..0f461601 100644 --- a/MapControl/Shared/PolarStereographicProjection.cs +++ b/MapControl/Shared/PolarStereographicProjection.cs @@ -32,6 +32,7 @@ namespace MapControl { var sign = Math.Sign(LatitudeOfOrigin); var phi = sign * latitude * Math.PI / 180d; + var e = Math.Sqrt((2d - Flattening) * Flattening); var eSinPhi = e * Math.Sin(phi); var t = Math.Tan(Math.PI / 4d - phi / 2d) @@ -55,7 +56,6 @@ namespace MapControl var e = Math.Sqrt((2d - Flattening) * Flattening); var eSinPhi = e * Math.Sin(phi); - var t = Math.Tan(Math.PI / 4d - phi / 2d) / Math.Pow((1d - eSinPhi) / (1d + eSinPhi), e / 2d); // p.161 (15-9) // ρ diff --git a/MapControl/Shared/TransverseMercatorProjection.cs b/MapControl/Shared/TransverseMercatorProjection.cs index 98a4dea4..c1db6dad 100644 --- a/MapControl/Shared/TransverseMercatorProjection.cs +++ b/MapControl/Shared/TransverseMercatorProjection.cs @@ -17,6 +17,7 @@ namespace MapControl public class TransverseMercatorProjection : MapProjection { private readonly double n; + private readonly double m; // 2*sqrt(n)/(1+n) private readonly double a1; // α1 private readonly double a2; // α2 private readonly double a3; // α3 @@ -26,8 +27,7 @@ namespace MapControl private readonly double d1; // δ1 private readonly double d2; // δ2 private readonly double d3; // δ3 - private readonly double f1; // A/a - private readonly double f2; // 2*sqrt(n)/(1+n) + private readonly double A; public TransverseMercatorProjection( double equatorialRadius, double flattening, @@ -42,8 +42,10 @@ namespace MapControl FalseNorthing = falseNorthing; n = flattening / (2d - flattening); + m = 2d * Math.Sqrt(n) / (1d + n); var n2 = n * n; var n3 = n * n2; + a1 = n / 2d - n2 * 2d / 3d + n3 * 5d / 16d; a2 = n2 * 13d / 48d - n3 * 3d / 5d; a3 = n3 * 61d / 240d; @@ -53,8 +55,7 @@ namespace MapControl d1 = n * 2d - n2 * 2d / 3d - n3 * 2d; d2 = n2 * 7d / 3d - n3 * 8d / 5d; d3 = n3 * 56d / 15d; - f1 = (1d + n2 / 4d + n2 * n2 / 64d) / (1d + n); - f2 = 2d * Math.Sqrt(n) / (1d + n); + A = equatorialRadius / (1d + n) * (1d + n2 / 4d + n2 * n2 / 64d); } public TransverseMercatorProjection( @@ -87,7 +88,7 @@ namespace MapControl var cosLambda = Math.Cos(lambda); var tanLambda = Math.Tan(lambda); // t - var t = Math.Sinh(Atanh(sinPhi) - f2 * Atanh(f2 * sinPhi)); + var t = Math.Sinh(Atanh(sinPhi) - m * Atanh(m * sinPhi)); var u = Math.Sqrt(1d + t * t); // ξ' var xi_ = Math.Atan2(t, cosLambda); @@ -104,8 +105,9 @@ namespace MapControl 4d * a2 * Math.Sin(4d * xi_) * Math.Sinh(4d * eta_) + 6d * a3 * Math.Sin(6d * xi_) * Math.Sinh(6d * eta_); - var m = (1d - n) / (1d + n) * Math.Tan(phi); - var k = ScaleFactor * f1 * Math.Sqrt((1d + m * m) * (sigma * sigma + tau * tau) / (t * t + cosLambda * cosLambda)); + var q = (1d - n) / (1d + n) * Math.Tan(phi); + var k = ScaleFactor * A / EquatorialRadius + * Math.Sqrt((1d + q * q) * (sigma * sigma + tau * tau) / (t * t + cosLambda * cosLambda)); // γ, grid convergence var gamma = Math.Atan2(tau * u + sigma * t * tanLambda, sigma * u - tau * t * tanLambda); @@ -122,22 +124,20 @@ namespace MapControl var phi = latitude * Math.PI / 180d; var sinPhi = Math.Sin(phi); // t - var t = Math.Sinh(Atanh(sinPhi) - f2 * Atanh(f2 * sinPhi)); + var t = Math.Sinh(Atanh(sinPhi) - m * Atanh(m * sinPhi)); // λ - λ0 var lambda = (longitude - CentralMeridian) * Math.PI / 180d; // ξ' var xi_ = Math.Atan2(t, Math.Cos(lambda)); // η' var eta_ = Atanh(Math.Sin(lambda) / Math.Sqrt(1d + t * t)); - // k0 * A - var k0A = ScaleFactor * EquatorialRadius * f1; - var x = FalseEasting + k0A * (eta_ + + var x = FalseEasting + ScaleFactor * A * (eta_ + a1 * Math.Cos(2d * xi_) * Math.Sinh(2d * eta_) + a2 * Math.Cos(4d * xi_) * Math.Sinh(4d * eta_) + a3 * Math.Cos(6d * xi_) * Math.Sinh(6d * eta_)); - var y = FalseNorthing + k0A * (xi_ + + var y = FalseNorthing + ScaleFactor * A * (xi_ + a1 * Math.Sin(2d * xi_) * Math.Cosh(2d * eta_) + a2 * Math.Sin(4d * xi_) * Math.Cosh(4d * eta_) + a3 * Math.Sin(6d * xi_) * Math.Cosh(6d * eta_)); @@ -147,12 +147,10 @@ namespace MapControl public override Location MapToLocation(double x, double y) { - // k0 * A - var k0A = ScaleFactor * EquatorialRadius * f1; // ξ - var xi = (y - FalseNorthing) / k0A; + var xi = (y - FalseNorthing) / (ScaleFactor * A); // η - var eta = (x - FalseEasting) / k0A; + var eta = (x - FalseEasting) / (ScaleFactor * A); // ξ' var xi_ = xi - b1 * Math.Sin(2d * xi) * Math.Cosh(2d * eta) -