Skip to content

Instantly share code, notes, and snippets.

@Tythos
Last active August 12, 2021 18:23
Show Gist options
  • Save Tythos/885c2db3de71c6fb12aab159a61edf58 to your computer and use it in GitHub Desktop.
Save Tythos/885c2db3de71c6fb12aab159a61edf58 to your computer and use it in GitHub Desktop.
spheregeo
/**
*/
define(function(require, exports, module) {
let D2R = Math.PI / 180;
let sum = function(series) { return series.reduce(function(a=0, b) { return a + b; } ) };
/**
* Returns a unit vector (3-element Array) for the given phi and theta
* values.
*
* @param {Number} phi_rad - Angle from north pole, where equator will be at pi/2 [radians]
* @param {Number} theta_rad - Angle about north pole, from meridian [radians]
* @returns {Array} - Unit vector corresponding to given phi/theta
*/
exports.unit = function(phi_rad, theta_rad) {
return [
Math.cos(theta_rad) * Math.sin(phi_rad),
Math.sin(theta_rad) * Math.sin(phi_rad),
Math.cos(phi_rad)
];
};
/**
* Computes the "safe" arc-cosine of a ratio x, in which it is rounded to
* such degrees of precision as will prevent numerical noise from causing a
* domain error.
*
* @param {Number} x - Ratio of two trigonometric edges
* @returns {Number} - Arc-cosine of ratio given as input, guaranteed to be within [-1,1].
*/
exports.safeacos = function(x) {
if (x > 1) {
return 0;
} else if (x < -1) {
return Math.PI;
}
return Math.acos(x);
};
/**
* Returns the dot product of two identicall-dimensioned Array arguments.
*
* @param {Array} v1 - Vector as an n-element Array
* @param {Array} v2 - Vector as an n-element Array
* @returns {Number} - Dot product of two given input arrays
*/
exports.dot = function(v1, v2) {
/* Returns dot product of two identically-dimensioned Array arguments
*/
if (v1.length != v2.length) {
console.error("dot product must operate on identically-dimensioned Array arguments");
}
let sum = 0;
v1.forEach(function(v, i) {
sum += v * v2[i];
});
return sum;
};
/**
* Returns cross-product of two three-element vectors (Arrays). Unlike the
* dot product operator, of course, three-element Arrays are required here.
*
* @param {Array} xyz1 - Three-element vector; LHS operand of cross-product operator
* @param {Array} xyz2 - Three-element vector; RHS operand of cross-product operator
* @returns {Array} - Cross-product of two given vectors
*/
exports.cross = function(xyz1, xyz2) {
let x = xyz1[1] * xyz2[2] - xyz2[1] * xyz1[2];
let y = xyz1[2] * xyz2[0] - xyz2[2] * xyz1[0];
let z = xyz1[0] * xyz2[1] - xyz2[0] * xyz1[1];
return [x, y, z];
};
/**
* Returns the difference between two identically-dimensioned Array
* arguments. Result will have the same dimensions, so this will work on 2
* and 3 element vectors.
*
* @param {Array} v1 - LHS vector of difference operator
* @param {Array} v2 - RHS vector of difference operator
* @returns {Array} - Element-wise difference between v1 and v2, such that result + v2 = v1
*/
exports.sub = function(v1, v2) {
if (v1.length != v2.length) {
console.error("dot product must operate on identically-dimensioned Array arguments");
}
return v1.map(function(v, i) {
return v - v2[i];
});
};
/**
* Returns *true* if the lat/lon triangle defined by points ll1, ll2, and
* ll3 is "outward"--that is, a right-hand rotation through that sequence
* of points is outward from the surface of a sphere positioned at the
* origin.
*
* @param {Array} ll1 - Lat/lon coordinates at first point of triangle [degrees]
* @param {Array} ll2 - Lat/lon coordinates at second point of triangle [degrees]
* @param {Array} ll3 - Lat/lon coordinates at third point of triangle [degrees]
* @returns {Boolean} - Will have value *true* if triangle is right-handed "outward" from origin
*/
exports.isTriOut = function(ll1, ll2, ll3) {
let ll1r = [ll1[0] * D2R, ll1[1] * D2R];
let ll2r = [ll2[0] * D2R, ll2[1] * D2R];
let ll3r = [ll3[0] * D2R, ll3[1] * D2R];
let xyz1 = [
Math.cos(ll1r[1]) * Math.cos(ll1r[0]),
Math.sin(ll1r[1]) * Math.cos(ll1r[0]),
Math.sin(ll1r[0])
];
let xyz2 = [
Math.cos(ll2r[1]) * Math.cos(ll2r[0]),
Math.sin(ll2r[1]) * Math.cos(ll2r[0]),
Math.sin(ll2r[0])
];
let xyz3 = [
Math.cos(ll3r[1]) * Math.cos(ll3r[0]),
Math.sin(ll3r[1]) * Math.cos(ll3r[0]),
Math.sin(ll3r[0])
];
let r23 = exports.sub(xyz3, xyz2);
let r21 = exports.sub(xyz1, xyz2);
return exports.dot(exports.cross(r23, r21), xyz2) > 0;
};
/**
* Computes and returns the *signed* area under the given spherical segment
* defined by two lat/lon coordinates. Positive values indicate the
* triangle defined by the ll1-ll2-southpole sequence is *inward* by the
* right-hand rule. Summing areas for each polygone edge, then, results in
* the total area of a triangle because "southern" edges will result in
* negative areas.
*
* @param {Array} ll1 - Lat/lon coordinates of first point [degress]
* @param {Array} ll2 - Lat/lon coordinates of second point [degress]
* @returns {Number} - Area under southern triangle (third point defined by south pole of sphere) [kilometers]
*/
exports.southArea = function(ll1, ll2) {
let lat1 = ll1[0] * D2R;
let lon1 = ll1[1] * D2R;
let lat2 = ll2[0] * D2R;
let lon2 = ll2[1] * D2R;
let s1 = Math.pow(Math.sin(0.5 * (lat2 - lat1)), 2.0);
let s2 = Math.pow(Math.sin(0.5 * (lon2 - lon1)), 2.0);
let s2d2 = s1 + Math.cos(lat1) * Math.cos(lat2) * s2;
let d = 2 * Math.asin(Math.pow(s2d2, 0.5));
let s = 0.5 * (d + (0.5 * Math.PI + lat1) + (0.5 * Math.PI + lat2));
let t1 = Math.tan(0.5 * s);
let t2 = Math.tan(0.5 * (s - d));
let t3 = Math.tan(0.5 * (s - (0.5 * Math.PI + lat1)));
let t4 = Math.tan(0.5 * (s - (0.5 * Math.PI + lat2)));
let E = 4 * Math.atan(Math.pow(t1 * t2 * t3 * t4, 0.5));
let R = 6371;
let A = E * Math.pow(R, 2.0);
if (exports.isTriOut(ll1, ll2, [-90,0])) {
A = -1 * A;
}
return A;
};
/**
* Uses trapezoidal approximation of integration along latitude bands to
* estimate area of the given series of points.
*
* @param {Array} ll_deg - Array of n lat/lon points--e.g., Array of two-element Arrays [degrees]
* @returns {Number} - Estimation of area within polygon defined by an array of lat/lon points [kilometers]
*/
exports.getAreaAppx = function(ll_deg) {
let N = ll_deg.length;
let terms = (new Array(N)).fill(0);
for (let i = 0; i < N; i++) {
let iNext = i < N - 1 ? i + 1 : 0;
let iPrev = i > 0 ? i - 1 : N - 1;
let lat = ll_deg[i][0] * D2R;
let dlon = ((ll_deg[iNext][1] - ll_deg[iPrev][1]) % 360) * D2R;
terms[i] = dlon * Math.sin(lat);
}
let R = 6371;
let s = sum(terms);
return 0.5 * Math.pow(R, 2.0) * s;
};
/**
* Sums signed area calculation from *southArea()* to integrate area of
* spherical polygon defined by given series of lat/lon points.
*
* @param {Array} ll_deg - Array of n lat/lon points--e.g., Array of two-element Arrays [degrees]
* @returns {Number} - Area of given polygon, estimated from trapezoidal approximation [kilometers]
*/
exports.getArea = function(ll_deg) {
let N = ll_deg.length;
let part = (new Array(N)).fill(0);
for (let i = 0; i < N; i++) {
let i1 = i < N - 1 ? i + 1 : 0;
part[i] = exports.southArea(ll_deg[i], ll_deg[i1]);
}
return sum(part);
};
/**
* Given three lat/lon coordinates, computes and returns the surface angle
* defined by the spherical triangle ABC at the vertex B. Angle is positive
* if BA x BC is aligned with OB (i.e., right-hand of angle is "outward").
*
* @param {Array} llA - Two-element Array with lat/lon values of point A [degrees]
* @param {Array} llB - Two-element Array with lat/lon values of point B [degrees]
* @param {Array} llC - Two-element Array with lat/lon values of point C [degrees]
* @returns {Number} - Angle at surface of spherical triangle ABC [degrees]
*/
exports.surfAng = function(llA, llB, llC) {
let xyzA = [
Math.cos(llA[1] * D2R) * Math.cos(llA[0] * D2R),
Math.sin(llA[1] * D2R) * Math.cos(llA[0] * D2R),
Math.sin(llA[0] * D2R)
];
let xyzB = [
Math.cos(llB[1] * D2R) * Math.cos(llB[0] * D2R),
Math.sin(llB[1] * D2R) * Math.cos(llB[0] * D2R),
Math.sin(llB[0] * D2R)
];
let xyzC = [
Math.cos(llC[1] * D2R) * Math.cos(llC[0] * D2R),
Math.sin(llC[1] * D2R) * Math.cos(llC[0] * D2R),
Math.sin(llC[0] * D2R)
];
let a = exports.safeacos(exports.dot(xyzB, xyzC));
let b = exports.safeacos(exports.dot(xyzC, xyzA));
let c = exports.safeacos(exports.dot(xyzA, xyzB));
let B = exports.safeacos((Math.cos(b) - Math.cos(a) * Math.cos(c)) / (Math.sin(a) * Math.sin(c)));
if (exports.dot(xyzB, exports.cross(xyzA, xyzC)) < 0) {
B = -1 * B;
}
return B / D2R;
};
/**
* Returns *true* if the lat/lon point Q_deg is within the spherical
* polygon defined by the given array of lat/lon points. This is computed
* from integrating signed-turning angles between each edge endpoint and Q.
* (In case you're curious, *pip* is short for "point in polygon", which is
* apparently a sufficiently-common term to deserve it's own unexplained
* abbreviation in A LOT of literature. Not that I'm bitter.)
*
* @param {Array} Q_deg - Lat/lon coordinates of point in question [degrees]
* @param {Array} vertices_deg - Spherical polygon defined by array of n lat/lon points--e.g., Array of two-element Arrays [degrees]
* @returns {Boolean} - Is *true* if given point lies within given spherical polygon
*/
exports.pip = function(Q_deg, vertices_deg) {
let turn_rad = 0;
for (let i1 = 0; i1 < vertices_deg.length; i1++) {
let ll1 = vertices_deg[i1];
let i2 = i1 < vertices_deg.length - 1 ? i1 + 1 : 0;
let ll2 = vertices_deg[i2];
turn_rad += exports.surfAng(ll1, Q_deg, ll2);
}
return Math.abs(turn_rad) > Math.PI;
};
/**
* Computes altitude of zenith-staring (e.g., straight down) eyeball given
* a specific FOV and earth-centered angle/span at ground. Specific to
* earth equatorial radius. In other words, how high do you have to be if,
* for a given FOV, to see a given arc of the earth's surface while staring
* straight down?
*
* @param {Number} fov_rad - Angle of field of view of eyeball [radians]
* @param {Number} ground_rad - Angle (along spherical arc about central body) of FOV projection [radians]
* @returns {Number} - Altitude of eyeball to achieve given angle of given FOV projected onto earth's surface [meters]
*/
exports.groundStareHeight = function(fov_rad, ground_rad) {
let re_m = 6.371e6;
let beta_rad = Math.PI - 0.5 * fov_rad - 0.5 * ground_rad;
let alt_m = re_m * (Math.sin(beta_rad) / Math.sin(0.5 * fov_rad) - 1.0);
return alt_m;
};
/**
* Given two lat/lon points, returns the angle (e.g., length of surface
* arc) between them. Unlike great circle, this is directly computed from
* corresponding unit vectors, so it's really the length of the triangle
* (not the spherical arc length), though these will of course be close for
* small angles.
*
* @param {Array} ll1 - Lat/lon coordinates of first point [degrees]
* @param {Array} ll2 - Lat/lon coordinates of second point [degrees]
* @returns {Number} - Angle (length) of surface arc between two lat/lon points [degrees]
*/
exports.latLonAngle = function(ll1, ll2) {
let xyz1 = exports.unit(0.5 * Math.PI - ll1[0] * D2R, ll1[1] * D2R);
let xyz2 = exports.unit(0.5 * Math.PI - ll2[0] * D2R, ll2[1] * D2R);
let ang = exports.safeacos(exports.dot(xyz1, xyz2)); // will be in radians, initially
return ang / D2R;
};
/**
* Returns length (angle, in degrees) of the great-circle solution between
* two given lat/lon points. Unlike *latLonAngle()*, this is in fact the
* length of the spherical arc between these two points.
*
* @param {Array} ll0_deg - Lat/lon coordinates of first point [degrees]
* @param {Array} llF_deg - Lat/lon coordinates of second point [degrees]
* @returns {Number} - Length of great-circle arc [degrees]
*/
exports.getGreatCircleLength = function(ll0_deg, llF_deg) {
/* Returns the length (angle, in degrees) of the great-circle solution
between the two given lat/lon points.
*/
let ll0 = ll0_deg.map(function(deg) { return deg * Math.PI / 180; });
let llF = llF_deg.map(function(deg) { return deg * Math.PI / 180; });
let dlon = llF[1] - ll0[1];
let ss = Math.sin(ll0[0]) * Math.sin(llF[0]);
let cc = Math.cos(ll0[0]) * Math.cos(llF[0]);
let dang = Math.acos(ss + cc * Math.cos(dlon));
return dang * 180 / Math.PI;
};
/**
* Returns vector *v* rotated about unit vector *k* by angle *tht*, using
* Rodrigues' rotation formula.
*
* @param {Array} v - A three-element vector defining the vector to be rotated
* @param {Array} k - A three-element unit (length=1) vector defining the axis of rotation
* @param {Number} tht - The angle by which vector *v* will be rotated about vector *k* [radians]
* @returns {Array} - The result of rotating vector *v* about axis *k* by angle *tht*
*/
exports.rodrigues = function(v, k, tht) {
k = exports.normalize(k); // just in case caller didn't; weird things can happen
let t1 = v.map(r => r * Math.cos(tht));
let kxv = exports.cross(k, v);
let t2 = kxv.map(r => r * Math.sin(tht));
let kdv = exports.dot(k, v);
let t3 = k.map(r => r * kdv * (1 - Math.cos(tht)));
let vr = [
t1[0] + t2[0] + t3[0],
t1[1] + t2[1] + t3[1],
t1[2] + t2[2] + t3[2]
];
return vr;
};
/**
* Returns a normalized copy of the given three-element vector. No, it does
* not use the fast inverse square algorithm.
*
* @param {Array} xyz - A thre-element vector
* @returns {Array} - A normalized version of the given three-element vector
*/
exports.normalize = function(xyz) {
let n = Math.sqrt(xyz[0]**2 + xyz[1]**2 + xyz[2]**2);
return xyz.map(r => r / n);
};
/**
* Returns a great circle solution to the given initial and final lat/lon
* points. Solution is defined by an Array of nPoints elements, each of
* which is a two-element Array of lat/lon point values. The path is
* computed using rotations about the cross-product between unit vectors
* for initial and final coordinates.
*
* @param {Array} ll0_deg - Two-element Array defining lat/lon coordinates of initial point [degrees]
* @param {Array} llF_deg - Two-element Array defining lat/lon coordinates of final point [degrees]
* @param {Number} nPoints - Number of points sampled for great circle path; defaults to 100
* @returns {Array} - Array of lat/lon points (each of which is a two-element Array) defining a great circle path between given points [degrees]
*/
exports.getGreatCirclePath = function(ll0_deg, llF_deg, nPoints=100) {
let xyz0 = exports.unit(0.5 * Math.PI - ll0_deg[0] * D2R, ll0_deg[1] * D2R);
let xyzF = exports.unit(0.5 * Math.PI - llF_deg[0] * D2R, llF_deg[1] * D2R);
// axis of rotation will be cross-product from xyz0 to xyzF
let X0F = exports.normalize(exports.cross(xyz0, xyzF));
let lli = [];
let beta_rad = exports.latLonAngle(ll0_deg, llF_deg) * D2R;
for (let i = 0; i < nPoints; i++) {
let ang_rad = beta_rad * i / (nPoints - 1);
let rot = exports.rodrigues(xyz0, X0F, ang_rad);
let xy = Math.sqrt(rot[0]**2 + rot[1]**2);
lli.push([
Math.atan2(rot[2], xy) / D2R,
Math.atan2(rot[1], rot[0]) / D2R
]);
}
return lli;
};
Object.assign(exports, {
"__url__": "https://gist.github.com/885c2db3de71c6fb12aab159a61edf58.git",
"__semver__": "0.1.0",
"__author_": "code@tythos.net",
"__license__": "MIT", // SPDX Identifier
});
});
/**
* Module for spherical geometry calculations. Includes some basic linear
* algebra operations as well to avoid external calculations, mostly treating
* 3d vectors as 3-element Arrays of Number values. Uses the SFJM standard (see
* https://dev.tythos.net/?name=sfjm for details) for single-file JavaScript
* modules. Unless otherwise noted, all operators return copies of result and
* (unlike, say, THREE.js vector operations) do not modify operands.
*
* @author "Brian Kirkpatrick" <code@tythos.net>
*/
//if (typeof(define) != "function") { function define(callback) { return callback(require, exports, module); } } // apparently a "let" upgrade to RequireJS foils this polyfill; will need some new tricks
define(function(require, exports, module) {
let D2R = Math.PI / 180;
let sum = function(series) { return series.reduce(function(a=0, b) { return a + b; } ) };
/**
* Returns a unit vector (3-element Array) for the given phi and theta
* values.
*
* @param {Number} phi_rad - Angle from north pole, where equator will be at pi/2 [radians]
* @param {Number} theta_rad - Angle about north pole, from meridian [radians]
* @returns {Array} - Unit vector corresponding to given phi/theta
*/
exports.unit = function(phi_rad, theta_rad) {
return [
Math.cos(theta_rad) * Math.sin(phi_rad),
Math.sin(theta_rad) * Math.sin(phi_rad),
Math.cos(phi_rad)
];
};
/**
* Computes the "safe" arc-cosine of a ratio x, in which it is rounded to
* such degrees of precision as will prevent numerical noise from causing a
* domain error.
*
* @param {Number} x - Ratio of two trigonometric edges
* @returns {Number} - Arc-cosine of ratio given as input, guaranteed to be within [-1,1].
*/
exports.safeacos = function(x) {
if (x > 1) {
return 0;
} else if (x < -1) {
return Math.PI;
}
return Math.acos(x);
};
/**
* Returns the dot product of two identicall-dimensioned Array arguments.
*
* @param {Array} v1 - Vector as an n-element Array
* @param {Array} v2 - Vector as an n-element Array
* @returns {Number} - Dot product of two given input arrays
*/
exports.dot = function(v1, v2) {
/* Returns dot product of two identically-dimensioned Array arguments
*/
if (v1.length != v2.length) {
console.error("dot product must operate on identically-dimensioned Array arguments");
}
let sum = 0;
v1.forEach(function(v, i) {
sum += v * v2[i];
});
return sum;
};
/**
* Returns cross-product of two three-element vectors (Arrays). Unlike the
* dot product operator, of course, three-element Arrays are required here.
*
* @param {Array} xyz1 - Three-element vector; LHS operand of cross-product operator
* @param {Array} xyz2 - Three-element vector; RHS operand of cross-product operator
* @returns {Array} - Cross-product of two given vectors
*/
exports.cross = function(xyz1, xyz2) {
let x = xyz1[1] * xyz2[2] - xyz2[1] * xyz1[2];
let y = xyz1[2] * xyz2[0] - xyz2[2] * xyz1[0];
let z = xyz1[0] * xyz2[1] - xyz2[0] * xyz1[1];
return [x, y, z];
};
/**
* Returns the difference between two identically-dimensioned Array
* arguments. Result will have the same dimensions, so this will work on 2
* and 3 element vectors.
*
* @param {Array} v1 - LHS vector of difference operator
* @param {Array} v2 - RHS vector of difference operator
* @returns {Array} - Element-wise difference between v1 and v2, such that result + v2 = v1
*/
exports.sub = function(v1, v2) {
if (v1.length != v2.length) {
console.error("dot product must operate on identically-dimensioned Array arguments");
}
return v1.map(function(v, i) {
return v - v2[i];
});
};
/**
* Returns *true* if the lat/lon triangle defined by points ll1, ll2, and
* ll3 is "outward"--that is, a right-hand rotation through that sequence
* of points is outward from the surface of a sphere positioned at the
* origin.
*
* @param {Array} ll1 - Lat/lon coordinates at first point of triangle [degrees]
* @param {Array} ll2 - Lat/lon coordinates at second point of triangle [degrees]
* @param {Array} ll3 - Lat/lon coordinates at third point of triangle [degrees]
* @returns {Boolean} - Will have value *true* if triangle is right-handed "outward" from origin
*/
exports.isTriOut = function(ll1, ll2, ll3) {
let ll1r = [ll1[0] * D2R, ll1[1] * D2R];
let ll2r = [ll2[0] * D2R, ll2[1] * D2R];
let ll3r = [ll3[0] * D2R, ll3[1] * D2R];
let xyz1 = [
Math.cos(ll1r[1]) * Math.cos(ll1r[0]),
Math.sin(ll1r[1]) * Math.cos(ll1r[0]),
Math.sin(ll1r[0])
];
let xyz2 = [
Math.cos(ll2r[1]) * Math.cos(ll2r[0]),
Math.sin(ll2r[1]) * Math.cos(ll2r[0]),
Math.sin(ll2r[0])
];
let xyz3 = [
Math.cos(ll3r[1]) * Math.cos(ll3r[0]),
Math.sin(ll3r[1]) * Math.cos(ll3r[0]),
Math.sin(ll3r[0])
];
let r23 = exports.sub(xyz3, xyz2);
let r21 = exports.sub(xyz1, xyz2);
return exports.dot(exports.cross(r23, r21), xyz2) > 0;
};
/**
* Computes and returns the *signed* area under the given spherical segment
* defined by two lat/lon coordinates. Positive values indicate the
* triangle defined by the ll1-ll2-southpole sequence is *inward* by the
* right-hand rule. Summing areas for each polygone edge, then, results in
* the total area of a triangle because "southern" edges will result in
* negative areas.
*
* @param {Array} ll1 - Lat/lon coordinates of first point [degress]
* @param {Array} ll2 - Lat/lon coordinates of second point [degress]
* @returns {Number} - Area under southern triangle (third point defined by south pole of sphere) [kilometers]
*/
exports.southArea = function(ll1, ll2) {
let lat1 = ll1[0] * D2R;
let lon1 = ll1[1] * D2R;
let lat2 = ll2[0] * D2R;
let lon2 = ll2[1] * D2R;
let s1 = Math.pow(Math.sin(0.5 * (lat2 - lat1)), 2.0);
let s2 = Math.pow(Math.sin(0.5 * (lon2 - lon1)), 2.0);
let s2d2 = s1 + Math.cos(lat1) * Math.cos(lat2) * s2;
let d = 2 * Math.asin(Math.pow(s2d2, 0.5));
let s = 0.5 * (d + (0.5 * Math.PI + lat1) + (0.5 * Math.PI + lat2));
let t1 = Math.tan(0.5 * s);
let t2 = Math.tan(0.5 * (s - d));
let t3 = Math.tan(0.5 * (s - (0.5 * Math.PI + lat1)));
let t4 = Math.tan(0.5 * (s - (0.5 * Math.PI + lat2)));
let E = 4 * Math.atan(Math.pow(t1 * t2 * t3 * t4, 0.5));
let R = 6371;
let A = E * Math.pow(R, 2.0);
if (exports.isTriOut(ll1, ll2, [-90,0])) {
A = -1 * A;
}
return A;
};
/**
* Uses trapezoidal approximation of integration along latitude bands to
* estimate area of the given series of points.
*
* @param {Array} ll_deg - Array of n lat/lon points--e.g., Array of two-element Arrays [degrees]
* @returns {Number} - Estimation of area within polygon defined by an array of lat/lon points [kilometers]
*/
exports.getAreaAppx = function(ll_deg) {
let N = ll_deg.length;
let terms = (new Array(N)).fill(0);
for (let i = 0; i < N; i++) {
let iNext = i < N - 1 ? i + 1 : 0;
let iPrev = i > 0 ? i - 1 : N - 1;
let lat = ll_deg[i][0] * D2R;
let dlon = ((ll_deg[iNext][1] - ll_deg[iPrev][1]) % 360) * D2R;
terms[i] = dlon * Math.sin(lat);
}
let R = 6371;
let s = sum(terms);
return 0.5 * Math.pow(R, 2.0) * s;
};
/**
* Sums signed area calculation from *southArea()* to integrate area of
* spherical polygon defined by given series of lat/lon points.
*
* @param {Array} ll_deg - Array of n lat/lon points--e.g., Array of two-element Arrays [degrees]
* @returns {Number} - Area of given polygon, estimated from trapezoidal approximation [kilometers]
*/
exports.getArea = function(ll_deg) {
let N = ll_deg.length;
let part = (new Array(N)).fill(0);
for (let i = 0; i < N; i++) {
let i1 = i < N - 1 ? i + 1 : 0;
part[i] = exports.southArea(ll_deg[i], ll_deg[i1]);
}
return sum(part);
};
/**
* Given three lat/lon coordinates, computes and returns the surface angle
* defined by the spherical triangle ABC at the vertex B. Angle is positive
* if BA x BC is aligned with OB (i.e., right-hand of angle is "outward").
*
* @param {Array} llA - Two-element Array with lat/lon values of point A [degrees]
* @param {Array} llB - Two-element Array with lat/lon values of point B [degrees]
* @param {Array} llC - Two-element Array with lat/lon values of point C [degrees]
* @returns {Number} - Angle at surface of spherical triangle ABC [degrees]
*/
exports.surfAng = function(llA, llB, llC) {
let xyzA = [
Math.cos(llA[1] * D2R) * Math.cos(llA[0] * D2R),
Math.sin(llA[1] * D2R) * Math.cos(llA[0] * D2R),
Math.sin(llA[0] * D2R)
];
let xyzB = [
Math.cos(llB[1] * D2R) * Math.cos(llB[0] * D2R),
Math.sin(llB[1] * D2R) * Math.cos(llB[0] * D2R),
Math.sin(llB[0] * D2R)
];
let xyzC = [
Math.cos(llC[1] * D2R) * Math.cos(llC[0] * D2R),
Math.sin(llC[1] * D2R) * Math.cos(llC[0] * D2R),
Math.sin(llC[0] * D2R)
];
let a = exports.safeacos(exports.dot(xyzB, xyzC));
let b = exports.safeacos(exports.dot(xyzC, xyzA));
let c = exports.safeacos(exports.dot(xyzA, xyzB));
let B = exports.safeacos((Math.cos(b) - Math.cos(a) * Math.cos(c)) / (Math.sin(a) * Math.sin(c)));
if (exports.dot(xyzB, exports.cross(xyzA, xyzC)) < 0) {
B = -1 * B;
}
return B / D2R;
};
/**
* Returns *true* if the lat/lon point Q_deg is within the spherical
* polygon defined by the given array of lat/lon points. This is computed
* from integrating signed-turning angles between each edge endpoint and Q.
* (In case you're curious, *pip* is short for "point in polygon", which is
* apparently a sufficiently-common term to deserve it's own unexplained
* abbreviation in A LOT of literature. Not that I'm bitter.)
*
* @param {Array} Q_deg - Lat/lon coordinates of point in question [degrees]
* @param {Array} vertices_deg - Spherical polygon defined by array of n lat/lon points--e.g., Array of two-element Arrays [degrees]
* @returns {Boolean} - Is *true* if given point lies within given spherical polygon
*/
exports.pip = function(Q_deg, vertices_deg) {
let turn_rad = 0;
for (let i1 = 0; i1 < vertices_deg.length; i1++) {
let ll1 = vertices_deg[i1];
let i2 = i1 < vertices_deg.length - 1 ? i1 + 1 : 0;
let ll2 = vertices_deg[i2];
turn_rad += exports.surfAng(ll1, Q_deg, ll2);
}
return Math.abs(turn_rad) > Math.PI;
};
/**
* Computes altitude of zenith-staring (e.g., straight down) eyeball given
* a specific FOV and earth-centered angle/span at ground. Specific to
* earth equatorial radius. In other words, how high do you have to be if,
* for a given FOV, to see a given arc of the earth's surface while staring
* straight down?
*
* @param {Number} fov_rad - Angle of field of view of eyeball [radians]
* @param {Number} ground_rad - Angle (along spherical arc about central body) of FOV projection [radians]
* @returns {Number} - Altitude of eyeball to achieve given angle of given FOV projected onto earth's surface [meters]
*/
exports.groundStareHeight = function(fov_rad, ground_rad) {
let re_m = 6.371e6;
let beta_rad = Math.PI - 0.5 * fov_rad - 0.5 * ground_rad;
let alt_m = re_m * (Math.sin(beta_rad) / Math.sin(0.5 * fov_rad) - 1.0);
return alt_m;
};
/**
* Given two lat/lon points, returns the angle (e.g., length of surface
* arc) between them. Unlike great circle, this is directly computed from
* corresponding unit vectors, so it's really the length of the triangle
* (not the spherical arc length), though these will of course be close for
* small angles.
*
* @param {Array} ll1 - Lat/lon coordinates of first point [degrees]
* @param {Array} ll2 - Lat/lon coordinates of second point [degrees]
* @returns {Number} - Angle (length) of surface arc between two lat/lon points [degrees]
*/
exports.latLonAngle = function(ll1, ll2) {
let xyz1 = exports.unit(0.5 * Math.PI - ll1[0] * D2R, ll1[1] * D2R);
let xyz2 = exports.unit(0.5 * Math.PI - ll2[0] * D2R, ll2[1] * D2R);
let ang = exports.safeacos(exports.dot(xyz1, xyz2)); // will be in radians, initially
return ang / D2R;
};
/**
* Returns length (angle, in degrees) of the great-circle solution between
* two given lat/lon points. Unlike *latLonAngle()*, this is in fact the
* length of the spherical arc between these two points.
*
* @param {Array} ll0_deg - Lat/lon coordinates of first point [degrees]
* @param {Array} llF_deg - Lat/lon coordinates of second point [degrees]
* @returns {Number} - Length of great-circle arc [degrees]
*/
exports.getGreatCircleLength = function(ll0_deg, llF_deg) {
/* Returns the length (angle, in degrees) of the great-circle solution
between the two given lat/lon points.
*/
let ll0 = ll0_deg.map(function(deg) { return deg * Math.PI / 180; });
let llF = llF_deg.map(function(deg) { return deg * Math.PI / 180; });
let dlon = llF[1] - ll0[1];
let ss = Math.sin(ll0[0]) * Math.sin(llF[0]);
let cc = Math.cos(ll0[0]) * Math.cos(llF[0]);
let dang = Math.acos(ss + cc * Math.cos(dlon));
return dang * 180 / Math.PI;
};
/**
* Returns vector *v* rotated about unit vector *k* by angle *tht*, using
* Rodrigues' rotation formula.
*
* @param {Array} v - A three-element vector defining the vector to be rotated
* @param {Array} k - A three-element unit (length=1) vector defining the axis of rotation
* @param {Number} tht - The angle by which vector *v* will be rotated about vector *k* [radians]
* @returns {Array} - The result of rotating vector *v* about axis *k* by angle *tht*
*/
exports.rodrigues = function(v, k, tht) {
k = exports.normalize(k); // just in case caller didn't; weird things can happen
let t1 = v.map(r => r * Math.cos(tht));
let kxv = exports.cross(k, v);
let t2 = kxv.map(r => r * Math.sin(tht));
let kdv = exports.dot(k, v);
let t3 = k.map(r => r * kdv * (1 - Math.cos(tht)));
let vr = [
t1[0] + t2[0] + t3[0],
t1[1] + t2[1] + t3[1],
t1[2] + t2[2] + t3[2]
];
return vr;
};
/**
* Returns a normalized copy of the given three-element vector. No, it does
* not use the fast inverse square algorithm.
*
* @param {Array} xyz - A thre-element vector
* @returns {Array} - A normalized version of the given three-element vector
*/
exports.normalize = function(xyz) {
let n = Math.sqrt(xyz[0]**2 + xyz[1]**2 + xyz[2]**2);
return xyz.map(r => r / n);
};
/**
* Returns a great circle solution to the given initial and final lat/lon
* points. Solution is defined by an Array of nPoints elements, each of
* which is a two-element Array of lat/lon point values. The path is
* computed using rotations about the cross-product between unit vectors
* for initial and final coordinates.
*
* @param {Array} ll0_deg - Two-element Array defining lat/lon coordinates of initial point [degrees]
* @param {Array} llF_deg - Two-element Array defining lat/lon coordinates of final point [degrees]
* @param {Number} nPoints - Number of points sampled for great circle path; defaults to 100
* @returns {Array} - Array of lat/lon points (each of which is a two-element Array) defining a great circle path between given points [degrees]
*/
exports.getGreatCirclePath = function(ll0_deg, llF_deg, nPoints=100) {
let xyz0 = exports.unit(0.5 * Math.PI - ll0_deg[0] * D2R, ll0_deg[1] * D2R);
let xyzF = exports.unit(0.5 * Math.PI - llF_deg[0] * D2R, llF_deg[1] * D2R);
// axis of rotation will be cross-product from xyz0 to xyzF
let X0F = exports.normalize(exports.cross(xyz0, xyzF));
let lli = [];
let beta_rad = exports.latLonAngle(ll0_deg, llF_deg) * D2R;
for (let i = 0; i < nPoints; i++) {
let ang_rad = beta_rad * i / (nPoints - 1);
let rot = exports.rodrigues(xyz0, X0F, ang_rad);
let xy = Math.sqrt(rot[0]**2 + rot[1]**2);
lli.push([
Math.atan2(rot[2], xy) / D2R,
Math.atan2(rot[1], rot[0]) / D2R
]);
}
return lli;
};
/**
* Returns azimuth of heading that will define a great-circle route from
* ll0 to llF lat/lon points.
*
* @param {Array} ll0_deg - Two-element Array defining lat/lon coordinates of initial point [degrees]
* @param {Array} llF_deg - Two-element Array defining lat/lon coordinates of final point [degrees]
* @returns {Number} - Azimuthal heading (zero from north towards east) of route from initial point that each final point [degrees]
*/
exports.getGreatCircleHeading = function(ll0_deg, llF_deg) {
// determine ENU frame at initial point
let xyz0 = exports.unit(0.5 * Math.PI - ll0_deg[0] * D2R, ll0_deg[1] * D2R);
let xyzF = exports.unit(0.5 * Math.PI - llF_deg[0] * D2R, llF_deg[1] * D2R);
let node = exports.normalize(exports.cross(xyz0, xyzF));
let heading = exports.normalize(exports.cross(node, xyz0));
let up = xyz0;
let east = exports.normalize(exports.cross([0, 0, 1], up));
let north = exports.normalize(exports.cross(up, east));
// azimuth will be angle between north and heading, but cross-product is required to determine sign
let az_rad = Math.acos(exports.dot(north, heading));
let hxn = exports.normalize(exports.cross(heading, north));
if (exports.dot(hxn, up) < 0) {
az_rad = 2 * Math.PI - az_rad;
}
return az_rad * 180 / Math.PI;
};
/**
* Returns *true* if the two scalar values are within a percent value of
* one another. Primarily used by numerical assertions when numerical
* precision is not infinite/reliable.
*
* @param {Number} s1 - First scalar value, against which difference will be compared
* @param {Number} s2 - Second scalar value, used as comparison for difference
* @param {Number} pct - Percent (e.g,. 1.0 is 100%) against which relative difference will be compared
* @returns {Boolean} - Will be *true* if second value is within percent threshold of first value
*/
exports.scalarsWithinPct = function(s1, s2, pct) {
return Math.abs(s2 - s1) / s1 < pct;
};
/**
* Returns *true* if the error between two vectors (difference over first
* vector length) is within a given percent. Used primarily for test
* assertions. Deliberately avoids any internal dependencies, so it may
* look a little overcomplicated. *DOES* require (for the time being)
* three-element vectors for comparison, since all the operations are
* hard-coded, but this wouldn't be too hard to change in the near-future.
*
* @param {Array} v1 - Three-element vector to compare against
* @param {Array} v2 - Three-element vector subject to comparision
* @param {Number} pct - Percent value (e.g., 1.0 is 100%) within which deviation of vectors is acceptible
* @returns {Boolean} - Will be *true* if second vector is within percent threshold of first vector
*/
exports.vectorsWithinPct = function(v1, v2, pct) {
let dv = [ v2[0] - v1[0], v2[1] - v1[1], v2[2] - v1[2] ];
let l1 = Math.sqrt(Math.pow(v1[0], 2.0) + Math.pow(v1[1], 2.0) + Math.pow(v1[2], 2.0));
let lD = Math.sqrt(Math.pow(dv[0], 2.0) + Math.pow(dv[1], 2.0) + Math.pow(dv[2], 2.0));
return lD / l1 < pct;
};
Object.assign(exports, {
"__url__": "https://gist.github.com/885c2db3de71c6fb12aab159a61edf58.git",
"__semver__": "1.0.0",
"__author_": "code@tythos.net",
"__license__": "MIT", // SPDX Identifier
"__tests__": [
function(assert) { assert(exports.vectorsWithinPct(
exports.unit(0.1, 0.2),
[0.09784, 0.01983, 0.9950],
1e-3
)) },
function(assert) { assert(exports.safeacos(-1.00001) == Math.PI) },
function(assert) { assert(exports.safeacos(0.0) == 0.5 * Math.PI) },
function(assert) { assert(exports.safeacos(1.00001) == 0.0) },
function(assert) { assert(exports.dot([1, 2, 3], [4, 5, 6]) == 32) },
function(assert) { assert(exports.vectorsWithinPct(
exports.cross([1, 2, 3], [4, 5, 6]),
[-3, 6, -3],
1e-3
)) },
function(assert) { assert(exports.vectorsWithinPct(
exports.sub([1, 2, 3], [4, 5, 6]),
[-3, -3, -3],
1e-3
)) },
function(assert) { assert(exports.isTriOut([0, 0], [0, 1], [1, 0])) },
function(assert) { assert(!exports.isTriOut([0, 0], [1, 0], [0, 1])) },
function(assert) { assert(exports.scalarsWithinPct(
exports.southArea([-10,0], [-10,10]),
5.85102e6,
1e-3
)) },
function(assert) { assert(exports.scalarsWithinPct(
exports.getAreaAppx([ [-1,-1], [1,-1], [1,1], [-1,1] ]),
4.94547e4,
1e-3
)) },
function(assert) { assert(exports.scalarsWithinPct(
exports.getArea([ [-1,-1], [1,-1], [1,1], [-1,1] ]),
4.94547e4,
1e-3
) )},
function(assert) { assert(exports.scalarsWithinPct(
exports.surfAng([1,0], [0,0], [0,1]),
0.5 * Math.PI,
1e-3
)) },
function(assert) { assert(exports.pip([0,0], [ [-1,-1], [-1,1], [1,1], [1,-1] ])) },
function(assert) { assert(!exports.pip([2,2], [ [-1,-1], [-1,1], [1,1], [1,-1] ])) },
function(assert) { assert(exports.scalarsWithinPct(
exports.groundStareHeight(0.12, 0.34),
1.7851e7,
1e-3
)) },
function(assert) { assert(exports.scalarsWithinPct(
exports.groundStareHeight(0.43, 0.21),
3.0226e6,
1e-3
)) },
function(assert) { assert(exports.scalarsWithinPct(
exports.latLonAngle([0,0], [0,1]),
1.0,
1e-3
)) },
function(assert) { assert(exports.scalarsWithinPct(
exports.latLonAngle([0,0], [1,0]),
1.0,
1e-3
)) },
function(assert) { assert(exports.scalarsWithinPct(
exports.latLonAngle([0,0], [1,1]),
1.4142,
1e-3
)) },
function(assert) { assert(exports.scalarsWithinPct(
exports.getGreatCircleLength([0,0], [0,45]),
45.0,
1e-3
)) },
function(assert) { assert(exports.scalarsWithinPct(
exports.getGreatCircleLength([0,0], [45,45]),
60.0,
1e-3
)) },
function(assert) { assert(exports.vectorsWithinPct(
exports.rodrigues([1,0,0], [0,0,1], 0.5 * Math.PI),
[0, 1, 0],
1e-3
)) },
function(assert) { assert(exports.vectorsWithinPct(
exports.rodrigues([1,2,3], [0,0,1], 1.0),
[-1.1426, 1.9221, 3],
1e-3
)) },
function(assert) { assert(exports.vectorsWithinPct(
exports.rodrigues([0,0,1], [1,2,3], 1.0),
[0.54829, -0.027879, 0.83582],
1e-3
)) },
function(assert) { assert(exports.vectorsWithinPct(
exports.normalize([1, 2, 3]),
[0.26726, 0.53452, 0.80178],
1e-3
)) },
function(assert) { assert(exports.vectorsWithinPct(
exports.getGreatCirclePath([0,0], [0,30], 3)[1].concat([0]),
[0, 15, 0],
1e-3
)) },
function(assert) { assert(exports.vectorsWithinPct(
exports.getGreatCirclePath( // now here's a fun verification...
[33.9415, -118.4056], // LAX
[35.7719, +140.3928], // NRT
101)[50].concat([0]), // midpoint
[47.6540, -168.2303, 0],
1e-3
)) },
function(assert) { assert(exports.scalarsWithinPct(
exports.getGreatCircleHeading(
[33.9415, -118.4056], // LAX
[35.7719, +140.3928], // NRT
),
305.75,
1e-3
)) },
function(assert) { assert(exports.scalarsWithinPct(
exports.getGreatCircleHeading([0,0], [1,1]),
45.0,
1e-3
)) },
// yes, we're going to test our asserts, they're part of the module, too
function(assert) { assert(exports.scalarsWithinPct(1.0, 1.001, 1e-3)) },
function(assert) { assert(!exports.scalarsWithinPct(1.0, 1.01, 1e-3)) },
function(assert) { assert(exports.vectorsWithinPct([1,1,1], [1,1,1.001], 1e-3)) },
function(assert) { assert(!exports.vectorsWithinPct([1,1,1], [1,1,1.01], 1e-3)) }
]
});
});
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment