diff --git a/src/client/types/geoms.ts b/src/client/types/geoms.ts new file mode 100644 index 0000000..628c61b --- /dev/null +++ b/src/client/types/geoms.ts @@ -0,0 +1,177 @@ +// import L from "leaflet" + +const ellipsoid = { + a: 6378137, + b: 6356752.3142, + f: 1 / 298.257223563 +}; + +function mod(n: number, p: number): number { + const r = n % p; + return r < 0 ? r + p : r; +} + +function wrap(degrees: number, max = 360) { + if (-max <= degrees && degrees <= max) { + return degrees; + } else { + return mod(degrees + max, 2 * max) - max; + } +} + +function dist(src, dst, itr = 100, mit = true): number { + const p1 = src, + p2 = dst; + const φ1 = toRadians(p1.lat), + λ1 = toRadians(p1.lng); + const φ2 = toRadians(p2.lat), + λ2 = toRadians(p2.lng); + const π = Math.PI; + const ε = Number.EPSILON; + + // allow alternative ellipsoid to be specified + const { a, b, f } = ellipsoid; + + const dL = λ2 - λ1; // L = difference in longitude, U = reduced latitude, defined by tan U = (1-f)·tanφ. + const tanU1 = (1 - f) * Math.tan(φ1), + cosU1 = 1 / Math.sqrt(1 + tanU1 * tanU1), + sinU1 = tanU1 * cosU1; + const tanU2 = (1 - f) * Math.tan(φ2), + cosU2 = 1 / Math.sqrt(1 + tanU2 * tanU2), + sinU2 = tanU2 * cosU2; + + const antipodal = Math.abs(dL) > π / 2 || Math.abs(φ2 - φ1) > π / 2; + + let λ = dL, + sinλ: number | null = null, + cosλ: number | null = null; // λ = difference in longitude on an auxiliary sphere + let σ = antipodal ? π : 0, + sinσ = 0, + cosσ = antipodal ? -1 : 1, + sinSqσ: number | null = null; // σ = angular distance P₁ P₂ on the sphere + let cos2σₘ = 1; // σₘ = angular distance on the sphere from the equator to the midpoint of the line + let sinα: number | null = null, + cosSqα = 1; // α = azimuth of the geodesic at the equator + let C: number | null = null; + + let λʹ: number | null = null, + iterations = 0; + do { + sinλ = Math.sin(λ); + cosλ = Math.cos(λ); + sinSqσ = + cosU2 * sinλ * (cosU2 * sinλ) + + (cosU1 * sinU2 - sinU1 * cosU2 * cosλ) * (cosU1 * sinU2 - sinU1 * cosU2 * cosλ); + if (Math.abs(sinSqσ) < ε) { + break; // co-incident/antipodal points (falls back on λ/σ = L) + } + sinσ = Math.sqrt(sinSqσ); + cosσ = sinU1 * sinU2 + cosU1 * cosU2 * cosλ; + σ = Math.atan2(sinσ, cosσ); + sinα = (cosU1 * cosU2 * sinλ) / sinσ; + cosSqα = 1 - sinα * sinα; + cos2σₘ = cosSqα !== 0 ? cosσ - (2 * sinU1 * sinU2) / cosSqα : 0; // on equatorial line cos²α = 0 (§6) + C = (f / 16) * cosSqα * (4 + f * (4 - 3 * cosSqα)); + λʹ = λ; + λ = dL + (1 - C) * f * sinα * (σ + C * sinσ * (cos2σₘ + C * cosσ * (-1 + 2 * cos2σₘ * cos2σₘ))); + const iterationCheck = antipodal ? Math.abs(λ) - π : Math.abs(λ); + if (iterationCheck > π) { + throw new EvalError("λ > π"); + } + } while (Math.abs(λ - λʹ) > 1e-12 && ++iterations < itr); + + if (iterations >= itr) { + if (mit) { + return dist( + src, + new L.LatLng(dst.lat, dst.lng - 0.01), + itr, + mit + ); + } else { + throw new EvalError(`Inverse vincenty formula failed to converge after ${itr} iterations + (start=${src.lat}/${src.lng}; dest=${dst.lat}/${dst.lng})`); + } + } + const uSq = (cosSqα * (a * a - b * b)) / (b * b); + const A = 1 + (uSq / 16384) * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq))); + const B = (uSq / 1024) * (256 + uSq * (-128 + uSq * (74 - 47 * uSq))); + const Δσ = + B * + sinσ * + (cos2σₘ + + (B / 4) * + (cosσ * (-1 + 2 * cos2σₘ * cos2σₘ) - + (B / 6) * cos2σₘ * (-3 + 4 * sinσ * sinσ) * (-3 + 4 * cos2σₘ * cos2σₘ))); + + const s = b * A * (σ - Δσ); // s = length of the geodesic + return s + +} + +function pointDistance(src: L.LatLng, dst: L.LatLng): number { + return dist( + new L.LatLng(src.lat, wrap(src.lng, 180)), + new L.LatLng(dst.lat, wrap(dst.lng, 180)) + ); +} + +function toRadians(degree: number): number { + return (degree * Math.PI) / 180; +} + +function toDegrees(radians: number): number { + return (radians * 180) / Math.PI; +} + +function midpoint(src: L.LatLng, dst: L.LatLng): L.LatLng { + // φm = atan2( sinφ1 + sinφ2, √( (cosφ1 + cosφ2⋅cosΔλ)² + cos²φ2⋅sin²Δλ ) ) + // λm = λ1 + atan2(cosφ2⋅sinΔλ, cosφ1 + cosφ2⋅cosΔλ) + // midpoint is sum of vectors to two points: mathforum.org/library/drmath/view/51822.html + + const φ1 = toRadians(src.lat); + const λ1 = toRadians(src.lng); + const φ2 = toRadians(dst.lat); + const Δλ = toRadians(dst.lng - src.lng); + + // get cartesian coordinates for the two points + const A = { x: Math.cos(φ1), y: 0, z: Math.sin(φ1) }; // place point A on prime meridian y=0 + const B = { x: Math.cos(φ2) * Math.cos(Δλ), y: Math.cos(φ2) * Math.sin(Δλ), z: Math.sin(φ2) }; + + // vector to midpoint is sum of vectors to two points (no need to normalise) + const C = { x: A.x + B.x, y: A.y + B.y, z: A.z + B.z }; + + const φm = Math.atan2(C.z, Math.sqrt(C.x * C.x + C.y * C.y)); + const λm = λ1 + Math.atan2(C.y, C.x); + + return new L.LatLng(toDegrees(φm), toDegrees(λm)); +} + +function recursiveMidPoint(src: L.LatLng, dst: L.LatLng, opt: { step?: number, dist?: number }): L.LatLng[] { + const geom: L.LatLng[] = [src, dst]; + const mp = midpoint(src, dst); + + if (opt.step != undefined) { + if (opt.step > 0) { + geom.splice(0, 1, ...recursiveMidPoint(src, mp, { step: opt.step - 1 })); + geom.splice(geom.length - 2, 2, ...recursiveMidPoint(mp, dst, { step: opt.step - 1 })); + } else { + geom.splice(1, 0, mp); + } + } else if (opt.dist != undefined) { + // console.log(src, dst, pointDistance(src, dst)) + if (pointDistance(src, dst) > opt.dist) { + geom.splice(0, 1, ...recursiveMidPoint(src, mp, { dist: opt.dist })); + geom.splice(geom.length - 2, 2, ...recursiveMidPoint(mp, dst, { dist: opt.dist })); + } else { + geom.splice(1, 0, mp); + } + } + console.log(geom) + return geom; +} + +export function getGeoLine(src: L.LatLng, dst: L.LatLng, opt: { step: number, dist: number }) { + + return recursiveMidPoint(src, dst, opt) +} \ No newline at end of file