calc_azimuth.js

/**
 * 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));