Fine tuning track distances #139

lib/track/calculator.dart

@@ -1,6 +1,7 @@
 import 'dart:math';
 
 import 'package:flutter/material.dart';
+import 'package:vector_math/vector_math.dart';
 
 import 'constants.dart';
 import 'track_descriptor.dart';
@@ -146,16 +147,99 @@ class TrackCalculator {
   }
 
   // https://stackoverflow.com/a/56499934/292502
-  static double distanceInKmBetweenEarthCoordinates(lat1, lon1, lat2, lon2) {
-    final dLat = (lat2 - lat1) * degreesToRadians;
-    final dLon = (lon2 - lon1) * degreesToRadians;
+  static double haversineDistance(double lat1, double lon1, double lat2, double lon2) {
+    final dLat = radians(lat2 - lat1);
+    final dLon = radians(lon2 - lon1);
 
     final latSin = sin(dLat / 2.0);
     final lonSin = sin(dLon / 2.0);
 
-    var a = latSin * latSin +
-        lonSin * lonSin * cos(lat1 * degreesToRadians) * cos(lat2 * degreesToRadians);
+    var a = latSin * latSin + lonSin * lonSin * cos(radians(lat1)) * cos(radians(lat2));
     var c = 2 * atan2(sqrt(a), sqrt(1 - a));
-    return earthRadiusKm * c;
+    const earthRadiusM = 6371008; // m
+    return earthRadiusM * c;
+  }
+
+  static double vincentyDistance(double lat1, double lon1, double lat2, double lon2) {
+    // Define WSG84 ellipsoid parameters
+    const a = 6378137.0;
+    const b = 6356752.314245;
+    const f = 1 / 298.257223563;
+
+    double L = radians(lon2 - lon1);
+    double u1 = atan((1 - f) * tan(radians(lat1)));
+    double u2 = atan((1 - f) * tan(radians(lat2)));
+    double sinU1 = sin(u1), cosU1 = cos(u1);
+    double sinU2 = sin(u2), cosU2 = cos(u2);
+
+    double lambda = L;
+    double lambdaP = 2 * pi;
+    int iterationLimit = 100;
+
+    double sinSigma = 0.0;
+    double cosSigma = 0.0;
+    double sigma = 0.0;
+    double cosSqAlpha = 0.0;
+    double cos2SigmaM = 0.0;
+
+    while ((lambda - lambdaP).abs() > 1e-12 && --iterationLimit > 0) {
+      double sinLambda = sin(lambda), cosLambda = cos(lambda);
+      sinSigma =
+          sqrt(pow(cosU2 * sinLambda, 2) + pow(cosU1 * sinU2 - sinU1 * cosU2 * cosLambda, 2));
+      if (sinSigma == 0) return 0; // co-incident points
+      cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
+      sigma = atan2(sinSigma, cosSigma);
+      double sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
+      cosSqAlpha = 1.0 - pow(sinAlpha, 2);
+      cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
+      if (cos2SigmaM.isNaN) cos2SigmaM = 0; // equatorial line: cosSqAlpha=0 (§6)
+      double C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
+      lambdaP = lambda;
+      lambda = L +
+          (1 - C) *
+              f *
+              sinAlpha *
+              (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * pow(cos2SigmaM, 2))));
+    }
+
+    // formula failed to converge
+    if (iterationLimit == 0) return double.nan;
+
+    double uSq = cosSqAlpha * (pow(a, 2) - pow(b, 2)) / pow(b, 2);
+    double A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
+    double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
+    double deltaSigma = B *
+        sinSigma *
+        (cos2SigmaM +
+            B /
+                4 *
+                (cosSigma * (-1 + 2 * pow(cos2SigmaM, 2)) -
+                    B /
+                        6 *
+                        cos2SigmaM *
+                        (-3 + 4 * pow(sinSigma, 2)) *
+                        (-3 + 4 * pow(cos2SigmaM, 2))));
+
+    return b * A * (sigma - deltaSigma);
+  }
+
+  static Offset degreesPerMeter(double latitude) {
+    // WGS-84 ellipsoid parameters
+    const a = 6378137.0;
+    const b = 6356752.314245;
+
+    double e2 = 1.0 - pow(b / a, 2);
+    double phi = radians(latitude);
+    double sinPhi = sin(phi);
+    double cosPhi = cos(phi);
+    double e2SinPhiSq = e2 * sinPhi * sinPhi;
+    double N = a / sqrt(1 - e2SinPhiSq);
+    double R = N * (1 - e2) / (1 - e2SinPhiSq);
+
+    double tmp = a * (1 - e2) / pow(1 - e2SinPhiSq, 1.5);
+    double D = sqrt(1 - e2SinPhiSq);
+    double S = tmp * D;
+
+    return Offset(degrees(1 / (R * cosPhi)), degrees(1 / S)); // lon, lat
   }
 }

lib/track/constants.dart

@@ -1,10 +1,6 @@
-import 'dart:math';
-
 const thick = 15.0;
 const trackLength = 400.0; // meters
 const fourHundredMTrackLengthFactor = 1.0;
 const fiveHundredMTrackLengthFactor = 1.25; // 500m, water
 const trackQuarter = trackLength / 4.0;
-const earthRadiusKm = 6371.008;
-const degreesToRadians = pi / 180.0;
 const trackPaintingRadiusBoost = 1.2;