Last active
August 12, 2021 18:23
-
-
Save Tythos/885c2db3de71c6fb12aab159a61edf58 to your computer and use it in GitHub Desktop.
spheregeo
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/** | |
*/ | |
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_": "[email protected]", | |
"__license__": "MIT", // SPDX Identifier | |
}); | |
}); |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/** | |
* 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" <[email protected]> | |
*/ | |
//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_": "[email protected]", | |
"__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