1/**
2 * Mesh topology for a spherical (nlat, nphi) grid following shtns conventions:
3 * - latitudinal grid given as cos(theta) (e.g. Gauss nodes, poles not included)
4 * - phi equally spaced starting at 0, endpoint excluded
5 *
6 * The phi seam is stitched when the phi grid spans the full circle, and pole
7 * cap vertices are added when the grid does not reach the poles, so that the
8 * rendered surface is closed.
9 *
10 * Adapted from figpack's SphereEmbedding view (figpack_experimental).
11 */
13import type { ColormapFunc } from './colormaps.ts';
15export type SphereMeshTopology = {
16 nlat: number;
17 nphi: number;
18 wrapPhi: boolean;
19 // Cap adjacent to row 0 / row nlat-1 (extra vertex appended after the grid)
20 startCapIndex: number; // -1 if absent
21 endCapIndex: number; // -1 if absent
22 numVertices: number;
23 indices: Uint32Array;
24 // Unit-sphere positions, length numVertices * 3
25 sphereRef: Float32Array;
26};
28export const buildTopology = (
29 cosTheta: Float64Array | Float32Array,
30 phi: Float64Array | Float32Array,
31): SphereMeshTopology => {
32 const nlat = cosTheta.length;
33 const nphi = phi.length;
35 // Does the phi grid span the full circle (so the seam should be stitched)?
36 let wrapPhi = false;
37 if (nphi >= 3) {
38 const dphi = phi[1] - phi[0];
39 const gap = phi[0] + 2 * Math.PI - phi[nphi - 1];
40 wrapPhi = Math.abs(gap - dphi) < 0.25 * Math.abs(dphi);
41 }
43 // Add pole caps where the grid does not reach the pole (|cos_theta| < 1),
44 // only when the surface wraps in phi (otherwise there is no hole to close)
45 const poleEps = 1e-9;
46 const hasStartCap = wrapPhi && Math.abs(Math.abs(cosTheta[0]) - 1) > poleEps;
47 const hasEndCap =
48 wrapPhi && Math.abs(Math.abs(cosTheta[nlat - 1]) - 1) > poleEps;
50 const numGridVertices = nlat * nphi;
51 let numVertices = numGridVertices;
52 const startCapIndex = hasStartCap ? numVertices++ : -1;
53 const endCapIndex = hasEndCap ? numVertices++ : -1;
55 const numCols = wrapPhi ? nphi : nphi - 1;
56 let numTriangles = (nlat - 1) * numCols * 2;
57 if (hasStartCap) numTriangles += nphi;
58 if (hasEndCap) numTriangles += nphi;
60 const indices = new Uint32Array(numTriangles * 3);
61 let k = 0;
62 for (let i = 0; i < nlat - 1; i++) {
63 for (let j = 0; j < numCols; j++) {
64 const j2 = (j + 1) % nphi;
65 const a = i * nphi + j;
66 const b = i * nphi + j2;
67 const c = (i + 1) * nphi + j;
68 const d = (i + 1) * nphi + j2;
69 indices[k++] = a;
70 indices[k++] = c;
71 indices[k++] = b;
72 indices[k++] = b;
73 indices[k++] = c;
74 indices[k++] = d;
75 }
76 }
77 if (hasStartCap) {
78 for (let j = 0; j < nphi; j++) {
79 const j2 = (j + 1) % nphi;
80 indices[k++] = startCapIndex;
81 indices[k++] = j;
82 indices[k++] = j2;
83 }
84 }
85 if (hasEndCap) {
86 const rowOffset = (nlat - 1) * nphi;
87 for (let j = 0; j < nphi; j++) {
88 const j2 = (j + 1) % nphi;
89 indices[k++] = rowOffset + j;
90 indices[k++] = endCapIndex;
91 indices[k++] = rowOffset + j2;
92 }
93 }
95 // Unit-sphere positions (z along the polar axis)
96 const sphereRef = new Float32Array(numVertices * 3);
97 for (let i = 0; i < nlat; i++) {
98 const ct = cosTheta[i];
99 const st = Math.sqrt(Math.max(0, 1 - ct * ct));
100 for (let j = 0; j < nphi; j++) {
101 const p = (i * nphi + j) * 3;
102 sphereRef[p] = st * Math.cos(phi[j]);
103 sphereRef[p + 1] = st * Math.sin(phi[j]);
104 sphereRef[p + 2] = ct;
105 }
106 }
107 if (hasStartCap) {
108 const p = startCapIndex * 3;
109 sphereRef[p + 2] = cosTheta[0] >= 0 ? 1 : -1;
110 }
111 if (hasEndCap) {
112 const p = endCapIndex * 3;
113 sphereRef[p + 2] = cosTheta[nlat - 1] >= 0 ? 1 : -1;
114 }
116 return {
117 nlat,
118 nphi,
119 wrapPhi,
120 startCapIndex,
121 endCapIndex,
122 numVertices,
123 indices,
124 sphereRef,
125 };
126};
128/**
129 * Expand a field frame (nlat * nphi) to per-vertex values (numVertices),
130 * with cap values averaged from the adjacent ring.
131 */
132export const fillFieldValues = (
133 out: Float32Array,
134 fieldFrame: Float32Array | Float64Array,
135 topo: SphereMeshTopology,
136): void => {
137 const { nlat, nphi } = topo;
138 const n = nlat * nphi;
139 for (let p = 0; p < n; p++) {
140 out[p] = fieldFrame[p];
141 }
142 const ringMean = (rowIndex: number) => {
143 let sum = 0;
144 let count = 0;
145 for (let j = 0; j < nphi; j++) {
146 const v = fieldFrame[rowIndex * nphi + j];
147 if (!Number.isNaN(v)) {
148 sum += v;
149 count++;
150 }
151 }
152 return count > 0 ? sum / count : NaN;
153 };
154 if (topo.startCapIndex >= 0) out[topo.startCapIndex] = ringMean(0);
155 if (topo.endCapIndex >= 0) out[topo.endCapIndex] = ringMean(nlat - 1);
156};
158/**
159 * Fill the color buffer (numVertices * 3, floats in [0, 1]) from per-vertex
160 * field values using the given colormap and range. NaN values render gray.
161 */
162export const fillColors = (
163 out: Float32Array,
164 values: Float32Array,
165 valueMin: number,
166 valueMax: number,
167 cmap: ColormapFunc,
168): void => {
169 const span = valueMax - valueMin;
170 const invSpan = span !== 0 ? 1 / span : 0;
171 for (let i = 0; i < values.length; i++) {
172 const v = values[i];
173 const p = i * 3;
174 if (Number.isNaN(v)) {
175 out[p] = 0.35;
176 out[p + 1] = 0.35;
177 out[p + 2] = 0.35;
178 } else {
179 const t = span !== 0 ? (v - valueMin) * invSpan : 0.5;
180 const [r, g, b] = cmap(t);
181 out[p] = r / 255;
182 out[p + 1] = g / 255;
183 out[p + 2] = b / 255;
184 }
185 }
186};