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 * Fill the position buffer (numVertices * 3) from the surface's coordinates
130 * (nlat * nphi * 3), interpolating toward the reference unit sphere.
131 * morph = 1 gives the surface itself; morph = 0 pulls it back to the sphere,
132 * which is the mesh the solver's parametrization actually lives on. Sweeping
133 * between them shows which points went where.
134 *
135 * The pole caps are not on the grid, so they take the mean of the adjacent
136 * ring — for a surface that is smooth at the pole, where the ring is a small
137 * circle around it, that is the pole to the accuracy the ring resolves.
138 */
139export const fillPositions = (
140 out: Float32Array,
141 coords: Float32Array | Float64Array,
142 topo: SphereMeshTopology,
143 morph: number,
144): void => {
145 const { nlat, nphi, sphereRef } = topo;
146 const n = nlat * nphi * 3;
147 for (let p = 0; p < n; p++) {
148 out[p] = (1 - morph) * sphereRef[p] + morph * coords[p];
149 }
150 const fillCap = (capIndex: number, rowIndex: number): void => {
151 let x = 0;
152 let y = 0;
153 let z = 0;
154 const rowOffset = rowIndex * nphi * 3;
155 for (let j = 0; j < nphi; j++) {
156 x += coords[rowOffset + j * 3];
157 y += coords[rowOffset + j * 3 + 1];
158 z += coords[rowOffset + j * 3 + 2];
159 }
160 const p = capIndex * 3;
161 out[p] = (1 - morph) * sphereRef[p] + (morph * x) / nphi;
162 out[p + 1] = (1 - morph) * sphereRef[p + 1] + (morph * y) / nphi;
163 out[p + 2] = (1 - morph) * sphereRef[p + 2] + (morph * z) / nphi;
164 };
165 if (topo.startCapIndex >= 0) fillCap(topo.startCapIndex, 0);
166 if (topo.endCapIndex >= 0) fillCap(topo.endCapIndex, nlat - 1);
167};
169/**
170 * Expand a field frame (nlat * nphi) to per-vertex values (numVertices),
171 * with cap values averaged from the adjacent ring.
172 */
173export const fillFieldValues = (
174 out: Float32Array,
175 fieldFrame: Float32Array | Float64Array,
176 topo: SphereMeshTopology,
177): void => {
178 const { nlat, nphi } = topo;
179 const n = nlat * nphi;
180 for (let p = 0; p < n; p++) {
181 out[p] = fieldFrame[p];
182 }
183 const ringMean = (rowIndex: number) => {
184 let sum = 0;
185 let count = 0;
186 for (let j = 0; j < nphi; j++) {
187 const v = fieldFrame[rowIndex * nphi + j];
188 if (!Number.isNaN(v)) {
189 sum += v;
190 count++;
191 }
192 }
193 return count > 0 ? sum / count : NaN;
194 };
195 if (topo.startCapIndex >= 0) out[topo.startCapIndex] = ringMean(0);
196 if (topo.endCapIndex >= 0) out[topo.endCapIndex] = ringMean(nlat - 1);
197};
199/**
200 * Fill the color buffer (numVertices * 3, floats in [0, 1]) from per-vertex
201 * field values using the given colormap and range. NaN values render gray.
202 */
203export const fillColors = (
204 out: Float32Array,
205 values: Float32Array,
206 valueMin: number,
207 valueMax: number,
208 cmap: ColormapFunc,
209): void => {
210 const span = valueMax - valueMin;
211 const invSpan = span !== 0 ? 1 / span : 0;
212 for (let i = 0; i < values.length; i++) {
213 const v = values[i];
214 const p = i * 3;
215 if (Number.isNaN(v)) {
216 out[p] = 0.35;
217 out[p + 1] = 0.35;
218 out[p + 2] = 0.35;
219 } else {
220 const t = span !== 0 ? (v - valueMin) * invSpan : 0.5;
221 const [r, g, b] = cmap(t);
222 out[p] = r / 255;
223 out[p + 1] = g / 255;
224 out[p + 2] = b / 255;
225 }
226 }
227};