fuzzyworld

shadow-zone / files / src / rayTracer.ts

rayTracer.ts

Blame: the commit that last changed each line · file · history

0ed88387 2026-06-12 Pavel Safronov1import type { EarthModel } from './earth';
2import { R_EARTH } from './earth';
3
4export interface RayPoint {
5 r: number;
6 theta: number;
7}
8
9export interface RayPath {
10 points: RayPoint[];
11 rayParameter: number;
12 takeoffAngle: number;
13 epicentralDistance: number | null;
14 travelTimeSec: number | null; // total travel time in seconds; null if ray never resurfaces
15 wave: 'P' | 'S';
16 blocked: boolean;
17}
18
19const STEP = 5; // km
20const MAX_STEPS = 60000;
21
22export function traceRay(
23 model: EarthModel,
24 rayParameter: number,
25 wave: 'P' | 'S',
26 takeoffAngle: number,
27 sourceDepthKm = 0,
28): RayPath {
29 const p = rayParameter;
30 const points: RayPoint[] = [];
31 let blocked = false;
32 let travelTimeSec = 0;
33
34 // Start slightly inside the surface (or at source depth)
35 let r = R_EARTH - sourceDepthKm - 0.1;
36 let theta = 0;
37 let sign = -1; // -1 = going down
38
e86a2c69 2026-06-13 Pavel Safronov39 // Impossible-core modifier state
40 // thetaSign: +1 normal angular direction, -1 for negative-refraction inside core
41 // forceRadial: true inside zero-index core — ray travels radially, no θ accumulation
42 let thetaSign = 1;
43 let forceRadial = false;
44
0ed88387 2026-06-12 Pavel Safronov45 const v0 = model.vAt(r, wave);
46 if (v0 === 0) {
47 return { points, rayParameter: p, takeoffAngle, epicentralDistance: null, travelTimeSec: null, wave, blocked: false };
48 }
49 const sinI0 = p * v0 / r;
50 if (sinI0 > 1) {
51 return { points, rayParameter: p, takeoffAngle, epicentralDistance: null, travelTimeSec: null, wave, blocked: false };
52 }
53
54 // For downgoing ray from a subsurface source, we also need to trace the upgoing
55 // leg to the surface on the source side (theta = 0 side). We start downward only —
56 // the upgoing leg is symmetric and the full path from source to receiver is traced
57 // as one continuous arc.
58
59 let epicentralDistance: number | null = null;
60
61 for (let step = 0; step < MAX_STEPS; step++) {
62 const v = model.vAt(r, wave);
63
64 if (wave === 'S' && v === 0) {
65 blocked = true;
66 break;
67 }
68
69 const sinI = p * v / r;
70
e86a2c69 2026-06-13 Pavel Safronov71 // For zero-index core the ray is forced radial — skip the Bouguer turning check.
72 if (!forceRadial && sinI >= 1) {
0ed88387 2026-06-12 Pavel Safronov73 sign = 1;
74 }
75
76 const rawCosI = Math.sqrt(Math.max(0, 1 - sinI * sinI));
e86a2c69 2026-06-13 Pavel Safronov77 // Zero-index core: ray is perfectly radial regardless of ray parameter.
78 const cosI = forceRadial ? 1 : (sinI >= 1 ? Math.max(rawCosI, 0.001) : rawCosI);
0ed88387 2026-06-12 Pavel Safronov79 const dr = sign * cosI * STEP;
e86a2c69 2026-06-13 Pavel Safronov80 // Zero-index: no angular accumulation; negative-refraction: reverse angular direction.
81 const dtheta = forceRadial ? 0 : thetaSign * (p * v / (r * r)) * STEP;
0ed88387 2026-06-12 Pavel Safronov82
83 const rNext = r + dr;
84 const layerNow = model.layerAt(r);
85 const layerNext = model.layerAt(rNext);
86
87 let newR = rNext;
88 let newSign = sign;
89
90 if (layerNow !== layerNext) {
91 const bR = sign < 0 ? layerNow.rMin : layerNow.rMax;
e86a2c69 2026-06-13 Pavel Safronov92 const isCMB = model.cmbRadius > 0 && Math.abs(bR - model.cmbRadius) < 2;
0ed88387 2026-06-12 Pavel Safronov93
94 const vOld = model.vAt(bR - 0.001 * sign, wave);
95 const vNew = model.vAt(bR + 0.001 * sign, wave);
96
e86a2c69 2026-06-13 Pavel Safronov97 // S-wave absorbed by liquid outer core
0ed88387 2026-06-12 Pavel Safronov98 if (wave === 'S' && vNew === 0) {
99 const fracStep = Math.abs(bR - r) / (Math.abs(dr) + 1e-10);
100 travelTimeSec += fracStep * STEP / v;
101 points.push({ r: bR, theta: theta + dtheta * fracStep });
102 blocked = true;
103 break;
104 }
105
e86a2c69 2026-06-13 Pavel Safronov106 // Perfect absorber: all waves stopped at CMB
107 if (isCMB && sign < 0 && model.coreModifier === 'absorber') {
108 const fracStep = Math.abs(bR - r) / (Math.abs(dr) + 1e-10);
109 travelTimeSec += fracStep * STEP / v;
110 points.push({ r: bR, theta: theta + dtheta * fracStep });
111 blocked = true;
112 break;
113 }
114
115 // Normal Snell's law TIR check — skipped inside zero-index core since the
116 // ray is perfectly radial (angle = 0) and passes through any inner boundary.
117 let didTIR = false;
118 if (!forceRadial) {
119 const sinIBound = p * vOld / bR;
120 const sinINew = sinIBound * (vNew / vOld);
121 if (Math.abs(sinINew) > 1) {
122 newSign = -sign;
123 newR = r + newSign * cosI * STEP;
124 didTIR = true;
125 }
126 }
127
128 // CMB impossible-core modifiers (only when ray actually crosses, not TIR'd)
129 if (isCMB && !didTIR) {
130 if (model.coreModifier === 'negative-refraction') {
131 // Negate tangential component: flip angular direction on each CMB crossing.
132 thetaSign = sign < 0 ? -1 : 1;
133 }
134 if (model.coreModifier === 'zero-index') {
135 // Inside core: radial-only travel. Restore normal Bouguer on exit.
136 forceRadial = sign < 0;
137 }
0ed88387 2026-06-12 Pavel Safronov138 }
139 }
140
141 points.push({ r, theta });
142 travelTimeSec += STEP / v;
143 r = newR;
144 theta += dtheta;
145 sign = newSign;
146
147 if (r >= R_EARTH) {
148 r = R_EARTH;
149 points.push({ r, theta });
e86a2c69 2026-06-13 Pavel Safronov150 // Use absolute value: negative theta means the ray curved backward
151 // (negative refraction), but distance is always a positive angle.
152 epicentralDistance = Math.abs((theta * 180) / Math.PI);
0ed88387 2026-06-12 Pavel Safronov153 break;
154 }
155
e86a2c69 2026-06-13 Pavel Safronov156 if (r <= 1) {
157 if (forceRadial) {
158 // Zero-index core: bounce off centre and head back out radially.
159 sign = 1;
160 r = 2;
161 continue;
162 }
163 // Spiral-trap and other rays that reach the centre are absorbed.
164 blocked = true;
165 break;
166 }
0ed88387 2026-06-12 Pavel Safronov167 }
168
169 return {
170 points,
171 rayParameter: p,
172 takeoffAngle,
173 epicentralDistance,
174 travelTimeSec: epicentralDistance !== null ? travelTimeSec : null,
175 wave,
176 blocked,
177 };
178}
179
180export function traceAllRays(
181 model: EarthModel,
182 numRays: number,
183 waveType: 'P' | 'S' | 'both',
184 sourceDepthKm = 0,
185): RayPath[] {
186 const results: RayPath[] = [];
187 const waves: Array<'P' | 'S'> = waveType === 'both' ? ['P', 'S'] : [waveType];
188 const sourceR = R_EARTH - sourceDepthKm;
189
190 for (const wave of waves) {
191 const vSource = model.vAt(sourceR - 0.1, wave);
192 if (vSource === 0) continue;
193
194 for (let i = 0; i < numRays; i++) {
195 // Non-linear distribution: t^0.6 biases toward larger takeoff angles
196 // (deeper, more diverse paths), reducing the dense cluster of near-identical
197 // shallow arcs near the epicenter.
198 const t = numRays > 1 ? i / (numRays - 1) : 0;
199 const angleDeg = 5 + Math.pow(t, 0.6) * 84; // 5° – 89°
200 const angleRad = (angleDeg * Math.PI) / 180;
201 const p = (sourceR * Math.sin(angleRad)) / vSource;
202 results.push(traceRay(model, p, wave, angleDeg, sourceDepthKm));
203 }
204 }
205
206 return results;
207}
208
209/**
210 * Given a model, trace enough rays to build a travel-time lookup:
211 * epicentral distance (degrees) → travel time (seconds).
212 * Used for validation against real earthquake observations.
213 */
214export function buildTravelTimeCurve(
215 model: EarthModel,
216 wave: 'P' | 'S',
217 sourceDepthKm: number,
218 numRays = 720,
219): Array<{ distanceDeg: number; travelTimeSec: number }> {
220 const curve: Array<{ distanceDeg: number; travelTimeSec: number }> = [];
221 const sourceR = R_EARTH - sourceDepthKm;
222 const vSource = model.vAt(sourceR - 0.1, wave);
223 if (vSource === 0) return curve;
224
225 for (let i = 0; i < numRays; i++) {
226 const angleDeg = 0.5 + (i / (numRays - 1)) * 88.5;
227 const angleRad = (angleDeg * Math.PI) / 180;
228 const p = (sourceR * Math.sin(angleRad)) / vSource;
229 const ray = traceRay(model, p, wave, angleDeg, sourceDepthKm);
230 if (ray.epicentralDistance !== null && ray.travelTimeSec !== null) {
231 curve.push({ distanceDeg: ray.epicentralDistance, travelTimeSec: ray.travelTimeSec });
232 }
233 }
234
235 curve.sort((a, b) => a.distanceDeg - b.distanceDeg);
236 return curve;
237}