-
Notifications
You must be signed in to change notification settings - Fork 17
/
Copy pathspherical.js
51 lines (45 loc) · 1.5 KB
/
spherical.js
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
var pi = Math.PI,
tau = 2 * pi,
quarterPi = pi / 4,
radians = pi / 180,
abs = Math.abs,
atan2 = Math.atan2,
cos = Math.cos,
sin = Math.sin;
function halfArea(ring, closed) {
var i = 0,
n = ring.length,
sum = 0,
point = ring[closed ? i++ : n - 1],
lambda0, lambda1 = point[0] * radians,
phi1 = (point[1] * radians) / 2 + quarterPi,
cosPhi0, cosPhi1 = cos(phi1),
sinPhi0, sinPhi1 = sin(phi1);
for (; i < n; ++i) {
point = ring[i];
lambda0 = lambda1, lambda1 = point[0] * radians;
phi1 = (point[1] * radians) / 2 + quarterPi;
cosPhi0 = cosPhi1, cosPhi1 = cos(phi1);
sinPhi0 = sinPhi1, sinPhi1 = sin(phi1);
// Spherical excess E for a spherical triangle with vertices: south pole,
// previous point, current point. Uses a formula derived from Cagnoli’s
// theorem. See Todhunter, Spherical Trig. (1871), Sec. 103, Eq. (2).
// See https://github.com/d3/d3-geo/blob/master/README.md#geoArea
var dLambda = lambda1 - lambda0,
sdLambda = dLambda >= 0 ? 1 : -1,
adLambda = sdLambda * dLambda,
k = sinPhi0 * sinPhi1,
u = cosPhi0 * cosPhi1 + k * cos(adLambda),
v = k * sdLambda * sin(adLambda);
sum += atan2(v, u);
}
return sum;
}
export function sphericalRingArea(ring, interior) {
var sum = halfArea(ring, true);
if (interior) sum *= -1;
return (sum < 0 ? tau + sum : sum) * 2;
}
export function sphericalTriangleArea(t) {
return abs(halfArea(t, false)) * 2;
}