diff --git a/MapControl/Shared/TransverseMercatorProjection.cs b/MapControl/Shared/TransverseMercatorProjection.cs index 6b03f70e..a4e1eeca 100644 --- a/MapControl/Shared/TransverseMercatorProjection.cs +++ b/MapControl/Shared/TransverseMercatorProjection.cs @@ -14,9 +14,15 @@ namespace MapControl /// public class TransverseMercatorProjection : MapProjection { - private readonly double[] alpha = new double[3]; // α_j - private readonly double[] beta = new double[3]; // β_j - private readonly double[] delta = new double[3]; // δ_j + private double a1; // α1 + private double a2; // α2 + private double a3; // α3 + private double b1; // β1 + private double b2; // β2 + private double b3; // β3 + private double d1; // δ1 + private double d2; // δ2 + private double d3; // δ3 private double f1; // A / a private double f2; // 2 * sqrt(n) / (1+n) @@ -26,26 +32,20 @@ namespace MapControl set { field = value; - // n, n^2, n^3 var n = field / (2d - field); - var n2 = n * n; - var n3 = n * n2; - // A / a - f1 = (1d + n2 / 4d + n2 * n2 / 64d) / (1d + n); - // 2 * sqrt(n) / (1+n) + var nn = n * n; + var nnn = n * nn; + a1 = n / 2d - nn * 2d / 3d + nnn * 5d / 16d; + a2 = nn * 13d / 48d - nnn * 3d / 5d; + a3 = nnn * 61d / 240d; + b1 = n / 2d - nn * 2d / 3d + nnn * 37d / 96d; + b2 = nn / 48d + nnn / 15d; + b3 = nnn * 17d / 480d; + d1 = n * 2d - nn * 2d / 3d - nnn * 2d; + d2 = nn * 7d / 3d - nnn * 8d / 5d; + d3 = nnn * 56d / 15d; + f1 = (1d + nn / 4d + nn * nn / 64d) / (1d + n); f2 = 2d * Math.Sqrt(n) / (1d + n); - // α_j - alpha[0] = n / 2d - n2 * 2d / 3d + n3 * 5d / 16d; - alpha[1] = n2 * 13d / 48d - n3 * 3d / 5d; - alpha[2] = n3 * 61d / 240d; - // β_j - beta[0] = n / 2d - n2 * 2d / 3d + n3 * 37d / 96d; - beta[1] = n2 / 48d + n3 / 15d; - beta[2] = n3 * 17d / 480d; - // δ_j - delta[0] = n * 2d - n2 * 2d / 3d - n3 * 2d; - delta[1] = n2 * 7d / 3d - n3 * 8d / 5d; - delta[2] = n3 * 56d / 15d; } } @@ -63,7 +63,7 @@ namespace MapControl public override Point RelativeScale(double latitude, double longitude) { - return new Point(ScaleFactor, ScaleFactor); + return new Point(ScaleFactor, ScaleFactor); // sufficiently precise } public override Point? LocationToMap(double latitude, double longitude) @@ -73,7 +73,7 @@ namespace MapControl #else static double Atanh(double x) => Math.Atanh(x); #endif - // k_0 * A + // k0 * A var k0A = ScaleFactor * EquatorialRadius * f1; // φ var phi = latitude * Math.PI / 180d; @@ -88,14 +88,14 @@ namespace MapControl var eta_ = Atanh(Math.Sin(lambda) / Math.Sqrt(1d + t * t)); // ξ var xi = xi_ - + alpha[0] * Math.Sin(2d * xi_) * Math.Cosh(2d * eta_) - + alpha[1] * Math.Sin(4d * xi_) * Math.Cosh(4d * eta_) - + alpha[2] * Math.Sin(6d * xi_) * Math.Cosh(6d * eta_); + + 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_); // η var eta = eta_ - + alpha[0] * Math.Cos(2d * xi_) * Math.Sinh(2d * eta_) - + alpha[1] * Math.Cos(4d * xi_) * Math.Sinh(4d * eta_) - + alpha[2] * Math.Cos(6d * xi_) * Math.Sinh(6d * 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_); return new Point( k0A * eta + FalseEasting, @@ -104,7 +104,7 @@ namespace MapControl public override Location MapToLocation(double x, double y) { - // k_0 * A + // k0 * A var k0A = ScaleFactor * EquatorialRadius * f1; // ξ var xi = (y - FalseNorthing) / k0A; @@ -112,21 +112,21 @@ namespace MapControl var eta = (x - FalseEasting) / k0A; // ξ' var xi_ = xi - - beta[0] * Math.Sin(2d * xi) * Math.Cosh(2d * eta) - - beta[1] * Math.Sin(4d * xi) * Math.Cosh(4d * eta) - - beta[2] * Math.Sin(6d * xi) * Math.Cosh(6d * eta); + - b1 * Math.Sin(2d * xi) * Math.Cosh(2d * eta) + - b2 * Math.Sin(4d * xi) * Math.Cosh(4d * eta) + - b3 * Math.Sin(6d * xi) * Math.Cosh(6d * eta); // η' var eta_ = eta - - beta[0] * Math.Cos(2d * xi) * Math.Sinh(2d * eta) - - beta[1] * Math.Cos(4d * xi) * Math.Sinh(4d * eta) - - beta[2] * Math.Cos(6d * xi) * Math.Sinh(6d * eta); + - b1 * Math.Cos(2d * xi) * Math.Sinh(2d * eta) + - b2 * Math.Cos(4d * xi) * Math.Sinh(4d * eta) + - b3 * Math.Cos(6d * xi) * Math.Sinh(6d * eta); // χ var chi = Math.Asin(Math.Sin(xi_) / Math.Cosh(eta_)); // φ var phi = chi - + delta[0] * Math.Sin(2d * chi) - + delta[1] * Math.Sin(4d * chi) - + delta[2] * Math.Sin(6d * chi); + + d1 * Math.Sin(2d * chi) + + d2 * Math.Sin(4d * chi) + + d3 * Math.Sin(6d * chi); // λ - λ0 var lambda = Math.Atan(Math.Sinh(eta_) / Math.Cos(xi_));