sphere-hello.xtl

1⍝!/usr/bin/env xetal
2⍝│ Sphere hello: a handful of points on the unit sphere. Two ways to
3⍝│ measure how far apart two of them are, the straight line through
4⍝│ the sphere and the great circle along it, order every point's
5⍝│ neighbors the same way; and the plain average of the points is not
6⍝│ on the sphere, where the mean along it is.
7"mf:" u̲se< "Manifold"
8"ch:" u̲se< "Check"
9
10⍝━━ The points
11
12⍝│ Six directions in three dimensions, chosen by hand (nothing is
13⍝│ random), each put onto the unit sphere: one point per row.
14P ← ᵐᶠs̲phereProject 6 3 r̲eshape 1.0 0.2 0.1 0.9 0.5 0.0 0.3 1.0 0.2 -0.4 0.8 0.5 0.1 0.3 1.0 -0.6 -0.5 0.6
15P
16"every point is on the sphere" ᶜʰt̲est (6 r̲eshape 0.0) ᶜʰn̲ear ᵐᶠs̲phereOff P
17
18⍝━━ Two distances
19
20⍝│ The straight line through the sphere (the chord) between every
21⍝│ pair: the length of the difference, summed one coordinate at a time.
22ᵘc̲hords ← { Q →
23  n ← f̲irst s̲hape Q
24  (D, _) ← (2 s̲elect s̲hape Q) '{ (d, k) → c ← k s̲elect₂ Q◆ (d + (c '− t̲able c) ^ 2.0, k + 1) } p̲ower ((n c̲at n) r̲eshape 0.0, 1)
25  s̲qrt D
26}
27chord ← ᵘc̲hords P
28chord
29
30⍝│ The great circle along the sphere between every pair: the angle,
31⍝│ in radians (Manifold).
32arc ← ᵐᶠs̲phereDistances P
33arc
34
35⍝│ The two agree on a formula: a chord is 2 sin of half its arc, so
36⍝│ the arc can be read off the chord and back. Checked on every pair.
37"chord = 2 sin (arc / 2)" ᶜʰt̲est chord ᶜʰn̲ear 2.0 × s̲in arc ÷ 2.0
38⍝│ And the arc is never the shorter: along the sphere is the long way
39⍝│ round, by a little for near points and by more for far ones.
40"arc >= chord" ᶜʰt̲est 1 ᶜʰn̲ear '∧ r̲/ r̲avel arc ≥ chord
41⍝│ The widest gap between the two, over all the pairs, and the widest arc.
42'm̲ax r̲/ r̲avel arc − chord
43'm̲ax r̲/ r̲avel arc
44
45⍝━━ The same neighbors
46
47⍝│ Each point's neighbors, nearest first (itself first of all, at
48⍝│ distance 0): one row per point, the points' numbers in order of
49⍝│ that row of the distance matrix.
50ᵘo̲rders ← { D →
51  n ← f̲irst s̲hape D
52  (R, _) ← n '{ (acc, i) → (acc c̲at (1 c̲at n) r̲eshape g̲rade i s̲elect D, i + 1) } p̲ower ((0 c̲at n) r̲eshape 0, 1)
53  R
54}
55byChord ← ᵘo̲rders chord
56byChord
57byArc ← ᵘo̲rders arc
58"the chord and the arc order every point's neighbors the same way" ᶜʰt̲est byChord ᶜʰi̲s byArc
59
60⍝━━ Two means
61
62⍝│ The plain mean of the points, coordinate by coordinate: it lies
63⍝│ inside the sphere, and its length says how far in.
64plain ← ('+ r̲/ P) ÷ 6.0
65plain
66s̲qrt '+ r̲/ plain × plain
67"the plain mean is off the sphere" ᶜʰt̲est 1 ᶜʰi̲s 0.0 < ᵐᶠs̲phereOff plain
68
69⍝│ The mean on the sphere (Manifold): the point along the sphere
70⍝│ whose summed squared arcs to the points is least. It is on the
71⍝│ sphere, and it is not the plain mean put onto it.
72onSphere ← ᵐᶠs̲phereMean P
73onSphere
74"the mean on the sphere is on the sphere" ᶜʰt̲est 0.0 ᶜʰn̲ear ᵐᶠs̲phereOff onSphere
75(ᵐᶠs̲phereProject plain) ᵐᶠs̲phereDistance onSphere