Skip to content

Instantly share code, notes, and snippets.

@adaburrows
Last active May 17, 2020 00:23
Show Gist options
  • Save adaburrows/72cc44f7b82522fb15e41e726b510f38 to your computer and use it in GitHub Desktop.
Save adaburrows/72cc44f7b82522fb15e41e726b510f38 to your computer and use it in GitHub Desktop.
/**
* Example code for computing azimuths and distance.
*
* This code is good for approximating distance, but some of the code is not
* gauranteed to converge. There are better algorithms located at:
* [Geodesics on an Ellipsoid](https://geographiclib.sourceforge.io/geod.html)
*
* For more background on why this is actually a more complicated problem, see:
* * [Most Changes in Earth's Shape Are Due to Changes in Climate](https://www.nasa.gov/centers/goddard/earthandsun/earthshape.html)
* * [THE EGM96 GEOID UNDULATION WITH RESPECT TO THE WGS84 ELLIPSOID](https://cddis.nasa.gov/926/egm96/doc/S11.HTML)
*/
// WGS84
const Earth = {
a: 6378137.0, // Major axis
f: 1 / 298.257222100882711243, // inverse flattening
b: 6356752.314140347, // b = (1 − ƒ) a
GM: 3.986004418, // x10^14 m^3 s^-2
};
/**
* Math utility functions
*/
const Util = {
deg2rad: (deg) => deg*Math.PI/180,
rad2deg: (rad) => (360 + rad * 180 / Math.PI) % 360,
sin2: (x) => (1 - Math.cos(2 * x)) / 2,
cos2: (x) => (1 + Math.cos(2 * x)) / 2,
hav: (x) => Util.sin2(x / 2),
};
/**
* GPSCoord implements spherical math
*/
class GPSCoord {
/**
* Constructs a coord in degrees
*/
constructor(lat, lng) {
this.lat = lat;
this.lng = lng;
}
/**
* Returns latitude in rads
*/
get radLat() {
return Util.deg2rad(this.lat);
}
/**
* Returns longitude in rads
*/
get radLng() {
return Util.deg2rad(this.lng);
}
/**
* Computes the difference in longitude
*/
L(that) {
// L = L2 - L1
return Util.deg2rad(that.lng - this.lng);
}
/**
* Calculates the azimuth at U_1 pointing towards U_2.
*/
alpha1(U_2, U_1, λ) {
return Math.atan2(
Math.sin(λ) * Math.cos(U_2),
Math.cos(U_1) * Math.sin(U_2) - Math.sin(U_1) * Math.cos(U_2) * Math.cos(λ)
);
}
/**
* Calculates the azimuth at U_2 after arriving along the arc from U_1
*/
alpha2(U_2, U_1, λ) {
return Math.atan2(
Math.sin(λ) * Math.cos(U_1),
-1 * Math.sin(U_1) * Math.cos(U_2) + Math.cos(U_1) * Math.sin(U_2) * Math.cos(λ)
);
}
/**
* Calculates the azimuth at this, looking at that.
*/
azimuthTo(that) {
let λ = this.L(that);
let U_2 = that.radLat;
let U_1 = this.radLat;
let a1 = this.alpha1(U_2, U_1, λ);
return Util.rad2deg(a1);
}
/**
* Calculates the azimuth at that after arriving along the arc from this
*/
azimuthAt(that) {
let λ = this.L(that);
let U_2 = that.radLat;
let U_1 = this.radLat;
let a2 = this.alpha2(U_2, U_1, λ);
return Util.rad2deg(a2);
}
/**
* Computes the distance between this and that using haversine
*/
distanceTo(that) {
let R = (Earth.a + Earth.b) / 2;
let λ = this.L(that);
let U_2 = that.radLat;
let U_1 = this.radLat;
let dU = U_2 - U_1;
let a = Util.hav(dU) + Math.cos(U_1) * Math.cos(U_2) * Util.hav(λ);
let d = 2 * R * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a))
return d;
}
}
/**
* GPSVCoord implements Vincenty's formulae:
* https://en.wikipedia.org/wiki/Vincenty%27s_formulae
*/
class GPSVCoord extends GPSCoord {
/**
* Returns latitude on the auxilary sphere adjusted for oblateness
*/
get U() {
return Math.atan((1 - Earth.f) * Math.tan(this.radLat));
}
sinSigma(U_2, U_1, λ) {
return Math.hypot(
Math.cos(U_2) * Math.sin(λ),
Math.cos(U_1) * Math.sin(U_2) -
Math.sin(U_1) * Math.cos(U_2) * Math.cos(λ)
);
}
cosSigma(U_2, U_1, λ) {
return (
Math.sin(U_1) * Math.sin(U_2) +
Math.cos(U_1) * Math.cos(U_2) * Math.cos(λ)
);
}
sigma(U_2, U_1, λ) {
return Math.atan2(
this.sinSigma(U_2, U_1, λ),
this.cosSigma(U_2, U_1, λ)
);
}
sinAlpha(U_2, U_1, λ) {
return (
Math.cos(U_1) * Math.cos(U_2) * Math.sin(λ) /
this.sinSigma(U_2, U_1, λ)
);
}
alpha(U_2, U_1, λ) {
return Math.asin(this.sinAlpha(U_2, U_1, λ));
}
cos2sig_m(U_2, U_1, λ) {
return this.cosSigma(U_2, U_1, λ) - (
2 * Math.sin(U_1) * Math.sin(U_2) /
Util.cos2(this.alpha(U_2, U_1, λ))
);
}
cos2alpha(U_2, U_1, λ) {
return Util.cos2(this.alpha(U_2, U_1, λ));
}
u2(U_2, U_1, λ) {
let a = Earth.a;
let b = Earth.b;
let cos2alpha = this.cos2alpha(U_2, U_1, λ);
return cos2alpha * (a*a - b*b) / (b*b);
}
A(U_2, U_1, λ) {
let u2 = this.u2(U_2, U_1, λ);
return 1 + u2 / 16384 * (4094 + u2 * (-768 + u2 * (320 - 175 * u2)))
}
B(U_2, U_1, λ) {
let u2 = this.u2(U_2, U_1, λ);
return u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2)));
}
C(U_2, U_1, λ) {
let f = Earth.f;
let cos2alpha = this.cos2alpha(U_2, U_1, λ);
return (f/16 * cos2alpha * (4 + f * (4 - 3 * cos2alpha)));
}
cos22sig_m(U_2, U_1, λ) {
let _2sigM = Math.acos(this.cos2sig_m(U_2, U_1, λ));
return Util.cos2(_2sigM);
}
/**
* Computes an update to lamba
*/
lambda(L, U_2, U_1, λ) {
let C = this.C(U_2, U_1, λ);
let sinA = this.sinAlpha(U_2, U_1, λ);
let sigma = this.sigma(U_2, U_1, λ);
let sinS = this.sinSigma(U_2, U_1, λ);
let c2sm = this.cos2sig_m(U_2, U_1, λ);
let cosS = this.cosSigma(U_2, U_1, λ);
let c22sm = this.cos22sig_m(U_2, U_1, λ);
return L + (
(1 - C) * Earth.f * sinA * (
sigma + C * sinS * (c2sm + C * cosS * (-1 + 2 * c22sm))
)
);
}
dSigma(U_2, U_1, λ) {
let B = this.B(U_2, U_1, λ);
let c2sm = this.cos2sig_m(U_2, U_1, λ);
let c22sm = this.cos22sig_m(U_2, U_1, λ);
let sinS = this.sinSigma(U_2, U_1, λ);
let cosS = this.cosSigma(U_2, U_1, λ);
let s2s = Util.sin2(this.sigma(U_2, U_1, λ));
return B * sinS * (
c2sm + B / 4 * (
(cosS * (-1 + 2 * c22sm)) -
(B / 6 * c2sm * (-3 + 4 * s2s) * (-3 + 4 * c22sm))
)
);
}
/**
* Computes the length of the line between the two points
*/
s(U_2, U_1, λ) {
let sigma = this.sigma(U_2, U_1, λ);
let dSigma = this.dSigma(U_2, U_1, λ);
let A = this.A(U_2, U_1, λ);
return Earth.b * A * (sigma - dSigma);
}
/**
* Solves for λ at the defined precision. May not coverge for certain values.
*/
solveLambda(L, U_2, U_1) {
let λ = L;
let λ_p = Number.POSITIVE_INFINITY; // an absurd value that will never match
let e = 0.000000000000001;
while (λ_p - λ > e) {
λ_p = λ;
λ = this.lambda(L, U_2, U_1, λ);
}
return λ;
}
/**
* Calculates the azimuth at this, looking at that.
*/
azimuthTo(that) {
let L = this.L(that);
let U_2 = that.U;
let U_1 = this.U;
let λ = this.solveLambda(L, U_2, U_1);
let a1 = this.alpha1(U_2, U_1, λ);
return Util.rad2deg(a1);
}
/**
* Calculates the azimuth at that after arriving along the arc from this
*/
azimuthAt(that) {
let L = this.L(that);
let U_2 = that.U;
let U_1 = this.U;
let λ = this.solveLambda(L, U_2, U_1);
let a2 = this.alpha2(U_2, U_1, λ);
return Util.rad2deg(a2);
}
/**
* Computes the distance between this and that using iterative solver
*/
distanceTo(that) {
let L = this.L(that);
let U_2 = that.U;
let U_1 = this.U;
let λ = this.solveLambda(L, U_2, U_1);
return this.s(U_2, U_1, λ);
}
}
let tower = new GPSCoord(45.034441, -68.682151);
let bsr = new GPSCoord(45.061331, -68.601991);
console.log("Spherical", bsr.distanceTo(tower), bsr.azimuthTo(tower), bsr.azimuthAt(tower))
let towerV = new GPSVCoord(45.034441, -68.682151);
let bsrV = new GPSVCoord(45.061331, -68.601991);
console.log("Vincenty's", bsrV.distanceTo(towerV), bsrV.azimuthTo(towerV), bsrV.azimuthAt(towerV));
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment