diff --git a/Makefile b/Makefile index feb5618..7fb16f3 100644 --- a/Makefile +++ b/Makefile @@ -106,6 +106,8 @@ endif # built for a different backend. TEST_OUT_DIR := $(OBJECT_DIR)/$(HDR_BUILD_TAG) TEST_TARGET := $(TEST_OUT_DIR)/test_geodesic +ASYMPTOTIC_TEST_TARGET := $(TEST_OUT_DIR)/test_asymptotic +ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET := $(TEST_OUT_DIR)/test_asymptotic_schwarzschild FRAME_TEST_TARGET := $(TEST_OUT_DIR)/test_frame SCHWARZSCHILD_TEST_TARGET := $(TEST_OUT_DIR)/test_schwarzschild ALCUBIERRE_TEST_TARGET := $(TEST_OUT_DIR)/test_alcubierre @@ -204,6 +206,12 @@ $(TEST_OUT_DIR): | $(BUILD_DIR) $(TEST_TARGET): tests/test_geodesic.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ +$(ASYMPTOTIC_TEST_TARGET): tests/test_asymptotic.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) + $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ + +$(ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET): tests/test_asymptotic_schwarzschild.c $(COMMON_SOURCES) src/spacetime_schwarzschild.c $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) + $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -DSPACETIME_SCHWARZSCHILD -Isrc $^ $(LDLIBS) -o $@ + $(FRAME_TEST_TARGET): tests/test_frame.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ @@ -259,10 +267,12 @@ FAST_PSF_FFTW_TEST_DEP := FAST_PSF_FFTW_TEST_RUN := endif -test: $(CAMERA_TEST_TARGETS) $(TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(ALCUBIERRE_TEST_TARGET) $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET) $(FAST_PSF_FFTW_TEST_DEP) $(TONE_MAP_TEST_TARGET) $(MOVIE_OUTPUT_TEST_TARGET) $(SENSOR_BLOOM_TEST_TARGET) +test: $(CAMERA_TEST_TARGETS) $(TEST_TARGET) $(ASYMPTOTIC_TEST_TARGET) $(ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(ALCUBIERRE_TEST_TARGET) $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET) $(FAST_PSF_FFTW_TEST_DEP) $(TONE_MAP_TEST_TARGET) $(MOVIE_OUTPUT_TEST_TARGET) $(SENSOR_BLOOM_TEST_TARGET) $(TEST_OUT_DIR)/test_observer_minkowski $(TEST_OUT_DIR)/test_observer_schwarzschild $(TEST_TARGET) + $(ASYMPTOTIC_TEST_TARGET) + $(ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(ALCUBIERRE_TEST_TARGET) diff --git a/nr_spacetime_movie_renderer_design.md b/nr_spacetime_movie_renderer_design.md index 8764aa0..2e5efa9 100644 --- a/nr_spacetime_movie_renderer_design.md +++ b/nr_spacetime_movie_renderer_design.md @@ -732,6 +732,288 @@ production renderer 不希望依赖每次 NR run 都开启昂贵的 AH finder。 --- +# 18A. 渐近外区、escape worldtube 与 endpoint 协议 + +本节冻结“到达 escape 区域就终止”这一旧行为被替换后的职责边界，是该协议的权威 +约定与唯一长期记录。 + +## 18A.1 问题 + +旧判定只看光线当前位置是否在某个 escape 半径之外，不看传播方向。因此当相机 +本身位于 escape 球外时，所有光线在初始化后立即被判为逃逸，连本应进入强场区的 +光线也不会积分。旧实现还把“外推到无穷远”简化为“在 escape 球处删除”，频移因此 +带有 $O(M/R)$ 的误差。 + +## 18A.2 ray 生命周期三段 + +1. **相机位于 escape worldtube 外**：由公共渐近外区模块判断光线是否会与 + worldtube 相交。 + - 相交：把光线外推到第一次由外向内穿越，并从该 entry event 开始交给 backend + 内区积分； + - 不相交：直接把光线外推到对应无穷远天球，写 `ESCAPED` endpoint。 +2. **光线在 backend 内区积分**：不再因为“当前位置处于 escape 区域”立即终止； + 只有沿 ray 的过去传播方向发生**有方向的 inside -> outside 穿越**时才进入外区 + 收尾。 +3. **穿越后**：公共渐近外区模块把有限半径处的 canonical photon state 推到无穷远， + 得到 `n_infinity` 和 `frequency_ratio`。 + +一个 backend 可以声明多个渐近远端；当前实现只暴露一个 `end_id`，但 endpoint 与 +ray 状态中不得把“整个时空只有一个无穷远”写死。 + +## 18A.3 职责边界 + +backend 负责声明： + +- `end_id` 及其稳定编号； +- 外区模型种类：`MINKOWSKI` 或 `SCHWARZSCHILD_MONOPOLE`； +- 从该远端看见的渐近质量 `mass`（允许为零）； +- 渐近参考系的 origin 与空间基（在 backend 坐标中表达）； +- 给定 coordinate time 的 escape worldtube 球心 `center`、速度 `velocity`、 + 半径 `radius`、半径变化率 `radius_rate`； +- worldtube 描述有效的时间区间与运动分段边界； +- backend 内一点属于哪个候选 end 的 outer region，或当前无法分类。 + +backend **不**负责：球外传播、entry/miss 判定、无穷远方向、pre-route 或 endpoint +写回。 + +公共渐近模块负责： + +- 将 backend photon state 与统一 canonical state 双向转换； +- 球外相机的 entry/miss 判定； +- 对 miss 光线直接生成 infinity endpoint； +- 把 entry 光线传播到 worldtube 的第一次由外向内穿越； +- 把内区积分产生的由内向外穿越传播到无穷远； +- $M=0$ 使用精确 Minkowski 几何；$M>0$ 固定同心球使用内建 monopole 近似； +- 返回明确状态码，而不是用 NaN 或任意 fallback 掩盖适用域错误。 + +geodesic/ray 生命周期层负责： + +- 初始化时调用 pre-route； +- 保存 entry event，并在 slab sweep 到达 entry time 时激活内区积分； +- 每个 accepted ODE step 后检测有方向的 crossing 并局部化第一次根； +- 调用公共外区模块完成 endpoint。 + +## 18A.4 canonical photon state 与时间约定 + +canonical state 至少包含：coordinate time $t$、渐近参考系中的位置、传播方向与 +能量/频移所需的 photon momentum 信息、以及 `end_id`。 + +renderer 沿过去方向积分。文档中使用的空间单位方向唯一约定为**过去传播方向** +$\mathbf w$：令 $s=t_\text{camera}-t\ge 0$，则局部轨迹满足 +$\mathbf x(s)=\mathbf x_0+s\,\mathbf w+\dots$。`n_infinity` 是光线在无穷远处的 +来向，即 $\mathbf w$ 在 $s\to\infty$ 的极限，因此与旧实现中 +`normalize(-gamma^{ij} Pi_j)` 的符号约定一致。 + +在静止时空（Minkowski 与 Kerr–Schild Schwarzschild）中，沿测地线守恒的 photon +能量为 +\[ +E_\infty=-p_t=\alpha p^0\left(\alpha-\beta^i\Pi_i\right), +\] +其中 $p_i$（即 `Pi` 的协变版本）满足 $p_i=\alpha p^0\,\Pi_i$。相机归一化取 +$E_\text{camera}=1$，故 +\[ +g=\frac{E_\text{camera}}{E_\infty}=\frac{1}{\alpha p^0(\alpha-\beta^i\Pi_i)}. +\] +旧实现对 $M>0$ 在 escape 球处直接返回 $\exp(-\log(\alpha p^0))$，是上式在 +$\beta^i\Pi_i\to0,\alpha\to1$ 下的近似。 + +## 18A.5 worldtube 与穿越方向 + +球面 worldtube： +\[ +F(t,\mathbf x)=|\mathbf x-\mathbf c(t)|^2-R(t)^2. +\] +$F>0$ 外、$F<0$ 内、$F=0$ 边界。边界点（$F=0$）必须结合过去传播方向的斜率 +$\mathrm dF/\mathrm ds$ 分类：$\mathrm dF/\mathrm ds<0$ 视为即将进入、$\ge0$ 视为 +向外或切触；因此相机恰在 $F=0$ 且 past-inward 才按 INSIDE 处理，past-outward +与 tangent 都按外层 route 处理。inside -> outside crossing 要求 $\mathrm dF/ +\mathrm ds>0$ 的严格符号变化（$F_{\rm before}\le0$ 且 $F_{\rm after}>0$）； +仅有 $F_{\rm after}=0$ 的单点切触不算 crossing，需等下一步是否真正到 $F>0$。 +必须区分方向： + +- camera pre-route 的 entry 是沿过去传播方向第一次 outside -> inside； +- 内区 escape 是沿过去传播方向第一次 inside -> outside； +- 某次采样发现 $F\ge0$ 不能独立构成 escape。 + +对步进端点接近零、切触和跨越 motion-segment 边界，使用显式容差和有界 root +localization；不得用固定位置 epsilon 把 tangent 误判成 crossing。 + +## 18A.6 $M=0$ 外区 + +渐近惯性系中为解析直线传播。固定球用 ray-sphere 二次方程取沿过去传播方向最早 +的合法根；匀速移动球在分段内把球心写成 $\mathbf c(t)=\mathbf c(t_0)+\mathbf v(t-t_0)$， +令 $\mathbf d=\mathbf x_0-\mathbf c(t_0)$、$\mathbf q=\mathbf w+\mathbf v$，entry 满足 +$|\mathbf d+s\mathbf q|^2=R^2$（半径线性变化时右端为 $(R_0-R_\text{rate}s)^2$）。 +任意加速球的未来接口使用分段 bracketed root driver；若 backend 历史在判定完成前 +结束，返回 `TIME_RANGE_EXHAUSTED`，不得武断判为 miss。 + +$M=0$ 的 finish 是平凡的：$\mathbf n_\infty=\mathbf w$、 +$g=\exp(-\log(\alpha p^0))$（flat 中守恒）。 + +## 18A.7 $M>0$ Schwarzschild-like 外区 + +首版严格限制：球心固定、escape 球与 monopole 同心、$R/M\ge64$、$M>0$、相机与 +worldtube 位于该外区。不满足则返回明确的 unsupported/domain 状态；不静默退回 +Minkowski，也不把一般移动 Schwarzschild 球解释成瞬时静态球。 + +无量纲量 $\rho=r/M$、$\beta=b/M$， +\[ +Q(\rho,\beta)=1-\beta^2\frac{1-2/\rho}{\rho^2}. +\] +escape 球处切触阈值 $\beta_R=\rho_R/\sqrt{1-2/\rho_R}$。entry/miss 的拓扑分类优先 +使用解析阈值和方向信息，不由低精度查表决定。 + +角度 primitive 为过去传播方向从半径 $\rho$ 到无穷远扫过的单调外向方位角 +\[ +\Phi(\rho,\beta)=\int_\rho^\infty +\frac{\beta}{\rho'^2\sqrt{Q(\rho',\beta)}}\,\mathrm d\rho' +=\int_0^{1/\rho}\frac{\beta\,\mathrm du}{\sqrt{1-\beta^2u^2+2\beta^2u^3}}. +\] +turning radius 满足 $\beta^2=\rho_\text{turn}^3/(\rho_\text{turn}-2)$。守恒的 +impact parameter 与角动量满足 +\[ +\beta=\frac{|x\times\Pi|}{\alpha-\beta^i\Pi_i},\qquad +\mathbf N=\widehat{x\times\Pi}, +\] +无穷远方向由 $\hat{\mathbf r}=x/|x|$ 绕 $\mathbf N$ 旋转 $\Phi(\rho,\beta)$ 得到。 +turning map 与 Schwarzschild coordinate-time transfer 使用离线验证过的有界 +residual、Chebyshev 表或解析主项；运行期不得建表。 + +若 time-transfer 无法在声明域内满足误差标准，则保留 $M=0$ 实现和接口，不把未经 +验证的时间公式写入生产代码，也不得降低验收标准。 + +## 18A.8 nmesh outer-shell 约定（仅约定，不实现） + +- 最外层必须是有明确六个面的 cubed-sphere shell； +- outer boundary 在渐近 frame 中是固定中心、固定半径球面； +- 提供 $R$、对应 end 的质量和 frame metadata； +- Schwarzschild monopole 模式要求 $R/M\ge64$；$64M$ 只是拒绝更靠内 junction + 的硬下限。nmesh 有 AMR，生产数据应把 outer shell 放到尽可能大的半径，建议以 + 至少接近解析 backend 当前的 $256M$ 为目标； +- DG element 边界本来允许场跳变，故 worldtube junction 不要求两侧 metric + pointwise 连续； +- 穿越时匹配 boundary local tetrad 中的 photon direction/energy，并用分辨率与 + outer-radius convergence test 验证，而不是强行匹配坐标分量。 + +## 18A.9 失败语义 + +遇到以下情况必须返回明确状态并上报，不得 fallback： + +- worldtube 历史不足，无法判断 first entry； +- Schwarzschild 外区不是固定同心球； +- $R/M<64$； +- canonical state 无法保持 null constraint 或 round-trip 精度； +- time coordinate 约定不明确； +- 高精度 reference evaluator 在域内不收敛； +- 表在 seam 或 grazing 区域超过误差限。 + +## 18A.10 当前实现状态 + +- 公共接口（end descriptor、escape worldtube sample、canonical photon state 与 + endpoint `end_id`）已实现； +- $M=0$ 固定球与匀速移动球的 camera pre-route、directed inside -> outside + crossing、`PENDING_ENTRY` 生命周期已实现，并接入 Minkowski 与 Alcubierre； +- 单帧与 movie 共用同一 pre-route/生命周期实现；`make test` 全绿； +- $M>0$ 固定同心 Schwarzschild monopole 外区以**解析 Carlson 椭圆积分**实现 + （见 18A.11），接入解析 Kerr–Schild backend；相机位于 escape 球外的 + pre-route、entry 传播、directed crossing 与 infinity endpoint 均可用； +- 任意加速 worldtube 的 bracketed root driver 已实现并由 synthetic accelerated + worldtube 测试覆盖（含跨 motion-segment 与 history-exhausted failure + semantics），但尚无生产 backend 使用该路径。该 driver 只检查离散端点的 + $F$ 符号，尚不能保证捕获一个步长内的窄进入；任意加速 backend 落地前需加入 + 段内速度/加速度界或自适应子区间搜索。 +- 多渐近远端路由尚未实现：`end_id` 已进入 endpoint、`LensVertex`、 + `terminal_mismatch` 与 `discrete_jacobian`，避免跨 end 插值；但 + `asymptotic_route_camera` 目前对多个 end 只处理第一个，lens-map 文件格式也 + 未序列化 `end_id`。虫洞/最大延拓接入前需要补这些。 +- 生命周期为三态：`NO_ENDS`（才允许 legacy）、`DIRECTED_READY`、 + `PROTOCOL_ERROR`（descriptor 读取失败或不支持的 exterior，显式 + `INVALID/UNSUPPORTED`，绝不退回 legacy）。pre-route 也先校验所有 descriptor + 与 exterior kind，再判定 worldtube 内外，因此“声明了不受支持 exterior 而相机 + 恰在球内”不会被静默接受。 +- 内区积分时若 worldtube sample 变为 `valid == 0`（历史耗尽），作为一等 + terminal reason 返回 `RAY_ENDPOINT_TIME_RANGE_EXHAUSTED`、保留 `end_id`，并由 + `RayPool` 记为 `RAY_POOL_TERMINATED`，与 pre-route 的耗尽语义一致。generic + accelerated driver 与内区 crossing localizer 的二分过程中任何 sample 失败都 + 直接传播该具体状态，不返回一个正常 entry 或普通 integration failure。 +- Minkowski entry quadratic 用稳定根公式（$q=-\tfrac12(b+\mathrm{copysign} + (\sqrt\Delta,b))$，取最小正根）；$c=0$（相机在边界）时按 $b=\mathrm dF/ + \mathrm ds$ 分类。worldtube sample 必须有限且 $R>0$，否则 `INVALID`。 +- **构造期验证优先**：`spacetime_create_*()` 成功即承诺该 source 已可安全光追。 + 每个 constructor 在安装 ops/context 后调用公共 `spacetime_source_finalize()`： + 检查 ops/context 完整、`end_id` 唯一且非 `NONE`、exterior kind 受支持、 + mass/frame 有限。**完整 worldtube 历史（所有 segment 边界、中心/半径连续性、 + 全程 $R>0$）由 backend constructor 负责**，普通解析 backend 在参数校验中完成； + 失败时 constructor 销毁 context 并返回错误，不存在半构造可用的 source。 +- **motion-segment 有效域**：候选二次根必须满足 $0\le s\le s_{\rm segment}$； + 本段无有效根时推进到下一段边界（用 `nextafter(boundary,-∞)` 进入下一段）重新 + 求根；只有最后一个向过去开放的恒速段无根时才判 `ESCAPED`；segment 预算耗尽 + 属于内部失败，返回 `INVALID`。因为 constructor 已保证全程 $R>0$，正常 + routing 不再做 segment-collapse 判定；但根处仍保留一次 + $R_0-\dot R\,\sigma>0$ 检查（一次乘减），用于防御绕过 constructor 的 backend， + 避免负半径伪 entry。 +- **运行期防御**：所有 backend（含 Schwarzschild）都经公共 + `worldtube_sample()` 读取，并保留一个便宜的 callback trust-boundary 检查 + `isfinite(radius) && radius > 0`：callback 失败 → `INVALID`、`valid=0` → + `TIME_RANGE_EXHAUSTED`、NaN/非正半径 → `INVALID`。这是防止第三方 backend 或 + 测试绕过 constructor 的安全网，不再承担正常配置验证；containment 循环原样 + 传播这些状态并保留 `end_id`。 +- turning radius 用 bracketed bisection + Newton polishing。相机紧贴大球时的 + entry 精度受 $|d|^2-R^2$ 输入条件数限制，由 entry-time 预算覆盖。 +- generic accelerated driver 的初始条件同样用 $F<0$，或 $F=0$ 且 + $\mathrm dF/\mathrm ds<0$ 才计为 entry；搜索只用严格 $F<0$ 的采样点确认 + entry，单点 $F=0$ 切触不算。 + +## 18A.11 $M>0$ 解析外区：闭式约化与验证 + +阶段 C 采用的不是 Chebyshev 表，而是把外区积分闭式约化到 Carlson 对称积分， +并在 double 下验证。记 $\rho=r/M$、$u=1/\rho$，$P(u)=1-\beta^2u^2+2\beta^2u^3$。 + +- **角度 primitive.** 对三次 $P$ 的实根分支用实 Legendre 形式 + \[ + \Phi=\frac{\sqrt2}{\sqrt{C-A}}\left[F(\varphi(u),\kappa) + -F(\varphi(0),\kappa)\right],\qquad + \sin^2\varphi=\frac{u-A}{B-A},\quad \kappa^2=\frac{B-A}{C-A}, + \] + $A +#include +#include + +static double dot3(const double a[3], const double b[3]) { + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; +} + +static double normalize3(double v[3]) { + const double length = sqrt(dot3(v, v)); + if (length > 0.0) + for (int i = 0; i < 3; ++i) + v[i] /= length; + return length; +} + +static int invert3(double a[3][3], double inv[3][3]) { + const double det = + a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1]) - + a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0]) + + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]); + if (!isfinite(det) || fabs(det) < 1e-300) + return -1; + inv[0][0] = (a[1][1] * a[2][2] - a[1][2] * a[2][1]) / det; + inv[0][1] = (a[0][2] * a[2][1] - a[0][1] * a[2][2]) / det; + inv[0][2] = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) / det; + inv[1][0] = (a[1][2] * a[2][0] - a[1][0] * a[2][2]) / det; + inv[1][1] = (a[0][0] * a[2][2] - a[0][2] * a[2][0]) / det; + inv[1][2] = (a[0][2] * a[1][0] - a[0][0] * a[1][2]) / det; + inv[2][0] = (a[1][0] * a[2][1] - a[1][1] * a[2][0]) / det; + inv[2][1] = (a[0][1] * a[2][0] - a[0][0] * a[2][1]) / det; + inv[2][2] = (a[0][0] * a[1][1] - a[0][1] * a[1][0]) / det; + return 0; +} + +static int find_end(const SpacetimeSource *source, SpacetimeEndId end_id, + SpacetimeAsymptoticEnd *out) { + const size_t count = spacetime_asymptotic_end_count(source); + for (size_t i = 0; i < count; ++i) { + SpacetimeAsymptoticEnd end; + if (spacetime_asymptotic_end(source, i, &end)) + continue; + if (end.end_id == end_id) { + if (out != NULL) + *out = end; + return 0; + } + } + return -1; +} + +/* Backend vector -> asymptotic-frame vector, where the frame axes are the + * columns of end->frame_axes expressed in backend coordinates. */ +static void backend_vector_to_frame(const SpacetimeAsymptoticEnd *end, + const double a[3], double out[3]) { + for (int j = 0; j < 3; ++j) + out[j] = end->frame_axes[0][j] * a[0] + end->frame_axes[1][j] * a[1] + + end->frame_axes[2][j] * a[2]; +} + +static void frame_vector_to_backend(const SpacetimeAsymptoticEnd *end, + const double a[3], double out[3]) { + for (int i = 0; i < 3; ++i) + out[i] = end->frame_axes[i][0] * a[0] + end->frame_axes[i][1] * a[1] + + end->frame_axes[i][2] * a[2]; +} + +static void backend_position_to_frame(const SpacetimeAsymptoticEnd *end, + const double x[3], double out[3]) { + const double shifted[3] = {x[0] - end->frame_origin[0], + x[1] - end->frame_origin[1], + x[2] - end->frame_origin[2]}; + backend_vector_to_frame(end, shifted, out); +} + +static void frame_position_to_backend(const SpacetimeAsymptoticEnd *end, + const double x[3], double out[3]) { + double rotated[3]; + frame_vector_to_backend(end, x, rotated); + for (int i = 0; i < 3; ++i) + out[i] = rotated[i] + end->frame_origin[i]; +} + +static AsymptoticStatus worldtube_sample( + const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end, double t, + SpacetimeEscapeWorldtubeSample *out) { + if (spacetime_escape_worldtube_sample(source, end->end_id, t, out)) + return ASYMPTOTIC_INVALID; + if (!out->valid) + return ASYMPTOTIC_TIME_RANGE_EXHAUSTED; + /* A declared worldtube must be a finite, positive-radius sphere. */ + if (!(out->radius > 0.0) || !isfinite(out->radius) || + !isfinite(out->radius_rate)) + return ASYMPTOTIC_INVALID; + for (int i = 0; i < 3; ++i) + if (!isfinite(out->center[i]) || !isfinite(out->velocity[i])) + return ASYMPTOTIC_INVALID; + return ASYMPTOTIC_OK; +} + +/* Worldtube value F and its derivative dF/ds along the past direction `w` + * (unit past spatial velocity, s = t0 - t). On the boundary F == 0 the sign + * of the slope decides inside vs outside. */ +static AsymptoticStatus worldtube_value_and_slope( + const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end, double t, + const double x[3], const double w[3], double *value, double *slope) { + SpacetimeEscapeWorldtubeSample sample; + const AsymptoticStatus status = worldtube_sample(source, end, t, &sample); + if (status != ASYMPTOTIC_OK) + return status; + double d[3], q[3]; + for (int i = 0; i < 3; ++i) { + d[i] = x[i] - sample.center[i]; + q[i] = w[i] + sample.velocity[i]; + } + *value = dot3(d, d) - sample.radius * sample.radius; + *slope = 2.0 * dot3(d, q) + 2.0 * sample.radius * sample.radius_rate; + return ASYMPTOTIC_OK; +} + +int asymptotic_worldtube_value(const SpacetimeSource *source, + SpacetimeEndId end_id, double t, + const double x[3], double *value) { + SpacetimeAsymptoticEnd end; + if (value == NULL || find_end(source, end_id, &end)) + return ASYMPTOTIC_INVALID; + double slope; + return worldtube_value_and_slope(source, &end, t, x, + (const double[3]){0.0, 0.0, 0.0}, value, + &slope); +} + +int asymptotic_canonical_from_backend(const SpacetimeSource *source, + SpacetimeEndId end_id, + const MetricData *metric, double t, + const double x[3], const double Pi[3], + double log_alpha_p0, + AsymptoticPhotonState *out) { + SpacetimeAsymptoticEnd end; + double inv[3][3], gamma[3][3]; + if (out == NULL || metric == NULL || find_end(source, end_id, &end)) + return -1; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + gamma[i][j] = metric->gamma[i][j]; + if (invert3(gamma, inv)) + return -1; + if (end.exterior_kind != ASYMPTOTIC_EXTERIOR_MINKOWSKI) + return -1; + double w_backend[3] = {0.0, 0.0, 0.0}; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + w_backend[i] -= inv[i][j] * Pi[j]; + if (normalize3(w_backend) <= 0.0) + return -1; + out->end_id = end_id; + out->t = t; + out->log_alpha_p0 = log_alpha_p0; + backend_position_to_frame(&end, x, out->x); + backend_vector_to_frame(&end, w_backend, out->w); + return 0; +} + +int asymptotic_backend_from_canonical(const SpacetimeSource *source, + const MetricData *metric, + const AsymptoticPhotonState *canonical, + double x[3], double Pi[3], + double *log_alpha_p0) { + SpacetimeAsymptoticEnd end; + (void)metric; + if (canonical == NULL || x == NULL || Pi == NULL || + find_end(source, canonical->end_id, &end)) + return -1; + if (end.exterior_kind != ASYMPTOTIC_EXTERIOR_MINKOWSKI) + return -1; + double w_backend[3]; + frame_position_to_backend(&end, canonical->x, x); + frame_vector_to_backend(&end, canonical->w, w_backend); + /* Flat exterior: the covariant momentum is the unit past direction negated. */ + for (int i = 0; i < 3; ++i) + Pi[i] = -w_backend[i]; + if (log_alpha_p0 != NULL) + *log_alpha_p0 = canonical->log_alpha_p0; + return 0; +} + +/* Solve |d + q s|^2 = (R0 - rr s)^2 for the smallest s >= 0 with outside -> + * inside crossing. Returns 1 on entry (sets s), 0 on miss, -1 on error. */ +static int solve_entry_quadratic(const double d[3], const double q[3], + double R0, double rr, double *s_out) { + const double qq = dot3(q, q); + const double a = qq - rr * rr; + const double b = 2.0 * (dot3(d, q) + R0 * rr); + const double c = dot3(d, d) - R0 * R0; + if (c < 0.0) { + /* Strictly inside; the lifecycle normally handles this as INSIDE. */ + *s_out = 0.0; + return 1; + } + if (c == 0.0) { + /* On the boundary: classify by dF/ds = b. Past-inward enters at once; + * outward/tangent rays may still re-enter later when the sphere shrinks + * (a < 0), so do not declare a permanent miss on the zero root. */ + if (b < 0.0) { + *s_out = 0.0; + return 1; + } + if (b == 0.0) { + if (a < 0.0) { + *s_out = 0.0; + return 1; + } + return 0; + } + if (a < 0.0) { + *s_out = -b / a; + return 1; + } + return 0; + } + /* Compare the quadratic coefficient against the velocity-squared scale it + * is built from; mixing in R0^2 would let a large radius misclassify a + * genuinely quadratic entry as linear. */ + const double scale = qq + rr * rr; + if (fabs(a) <= 32.0 * DBL_EPSILON * scale) { + if (!isfinite(b) || b >= 0.0) + return 0; + const double s = -c / b; + if (s <= 0.0) + return 0; + *s_out = s; + return 1; + } + const double disc = b * b - 4.0 * a * c; + if (!isfinite(disc) || disc <= 0.0) + return 0; + /* Numerically stable quadratic roots: q avoids cancellation in the root + * with the same sign as b, which is exactly the small entry root when the + * camera sits just outside a large sphere. */ + const double root = sqrt(disc); + const double qq2 = -0.5 * (b + copysign(root, b)); + const double r1 = qq2 / a; + const double r2 = c / qq2; + /* The first outside->inside crossing is the smallest positive root. */ + double s = INFINITY; + if (r1 > 0.0) + s = r1; + if (r2 > 0.0 && r2 < s) + s = r2; + if (!(s < INFINITY)) + return 0; + *s_out = s; + return 1; +} + +static double worldtube_F_frame(const double c_frame[3], double radius, + const double x_frame[3]) { + const double d[3] = {x_frame[0] - c_frame[0], x_frame[1] - c_frame[1], + x_frame[2] - c_frame[2]}; + return dot3(d, d) - radius * radius; +} + +/* Bracketed first-entry search for a worldtube whose motion is not constant. + * This is the future interface for accelerated worldtubes; the current + * backends always take the closed quadratic path above. */ +static AsymptoticStatus minkowski_generic_first_entry( + const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end, + double t0, const double x_frame[3], const double w_frame[3], + double *s_out) { + double s = 0.0; + SpacetimeEscapeWorldtubeSample sample; + AsymptoticStatus status = worldtube_sample(source, end, t0, &sample); + if (status != ASYMPTOTIC_OK) + return status; + double c_frame[3]; + backend_position_to_frame(end, sample.center, c_frame); + const double f_start = worldtube_F_frame(c_frame, sample.radius, x_frame); + if (f_start < 0.0) { + *s_out = 0.0; + return ASYMPTOTIC_OK; + } + if (f_start == 0.0) { + /* On the boundary only an inward slope is an entry; outward and tangent + * rays keep searching. */ + double v_frame[3], d0[3], q0[3]; + backend_vector_to_frame(end, sample.velocity, v_frame); + for (int i = 0; i < 3; ++i) { + d0[i] = x_frame[i] - c_frame[i]; + q0[i] = w_frame[i] + v_frame[i]; + } + const double slope0 = + 2.0 * dot3(d0, q0) + 2.0 * sample.radius * sample.radius_rate; + if (slope0 < 0.0) { + *s_out = 0.0; + return ASYMPTOTIC_OK; + } + } + const double speed = fabs(sample.velocity[0]) + fabs(sample.velocity[1]) + + fabs(sample.velocity[2]) + fabs(sample.radius_rate) + 1.0; + const double base_step = 0.5 * fmax(sample.radius, 1.0) / speed; + for (int iteration = 0; iteration < 1000000; ++iteration) { + double step = base_step; + const double boundary = spacetime_escape_worldtube_next_segment( + source, end->end_id, t0 - s); + if (isfinite(boundary)) { + /* The backward-integration distance to a past segment boundary. */ + const double to_boundary = (t0 - boundary) - s; + if (to_boundary > 0.0) + step = fmin(step, to_boundary); + } + const double s_next = s + step; + const double t_next = t0 - s_next; + SpacetimeEscapeWorldtubeSample next; + status = worldtube_sample(source, end, t_next, &next); + if (status != ASYMPTOTIC_OK) + return status; + double c_next[3], ray_next[3]; + backend_position_to_frame(end, next.center, c_next); + for (int i = 0; i < 3; ++i) + ray_next[i] = x_frame[i] + w_frame[i] * s_next; + const double f_next = worldtube_F_frame(c_next, next.radius, ray_next); + /* Only a strictly negative sample is an entry; a single touch at F == 0 + * (tangent) is not. */ + if (f_next < 0.0) { + double lo = s, hi = s_next; + AsymptoticStatus bisect_status = ASYMPTOTIC_OK; + for (int bisect = 0; bisect < 80; ++bisect) { + const double mid = 0.5 * (lo + hi); + SpacetimeEscapeWorldtubeSample mid_sample; + bisect_status = worldtube_sample(source, end, t0 - mid, &mid_sample); + if (bisect_status != ASYMPTOTIC_OK) + break; + double c_mid[3], ray_mid[3]; + backend_position_to_frame(end, mid_sample.center, c_mid); + for (int i = 0; i < 3; ++i) + ray_mid[i] = x_frame[i] + w_frame[i] * mid; + const double f_mid = + worldtube_F_frame(c_mid, mid_sample.radius, ray_mid); + if (f_mid <= 0.0) + hi = mid; + else + lo = mid; + } + /* A sample failure inside the bracket must not be disguised as a + * normal entry. */ + if (bisect_status != ASYMPTOTIC_OK) + return bisect_status; + *s_out = hi; + return ASYMPTOTIC_OK; + } + s = s_next; + } + return ASYMPTOTIC_INVALID; +} + +static void minkowski_route_escaped(const SpacetimeAsymptoticEnd *end, + const double w_frame[3], + SpacetimeEndId end_id, + AsymptoticRoute *route) { + route->kind = ASYMPTOTIC_ROUTE_ESCAPED; + route->end_id = end_id; + double w_backend[3]; + frame_vector_to_backend(end, w_frame, w_backend); + for (int i = 0; i < 3; ++i) + route->n_infinity[i] = w_backend[i]; +} + +static void minkowski_route_entry(const SpacetimeAsymptoticEnd *end, double t0, + const double x_frame[3], + const double w_frame[3], double s_entry, + SpacetimeEndId end_id, + AsymptoticRoute *route) { + double x_entry_frame[3]; + for (int i = 0; i < 3; ++i) + x_entry_frame[i] = x_frame[i] + w_frame[i] * s_entry; + route->kind = ASYMPTOTIC_ROUTE_ENTRY; + route->end_id = end_id; + route->activate_t = t0 - s_entry; + frame_position_to_backend(end, x_entry_frame, route->x); + double w_backend[3]; + frame_vector_to_backend(end, w_frame, w_backend); + for (int i = 0; i < 3; ++i) + route->Pi[i] = -w_backend[i]; +} + +static AsymptoticStatus minkowski_preroute( + const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end, double t0, + const double x_frame[3], const double w_frame[3], SpacetimeEndId end_id, + AsymptoticRoute *route) { + /* Walk constant-velocity motion segments. A quadratic root is only valid + * inside the current segment and while the radius stays positive; otherwise + * advance to the next segment boundary and re-sample. */ + double s = 0.0; + for (int segment = 0; segment < 1000000; ++segment) { + const double t = t0 - s; + SpacetimeEscapeWorldtubeSample sample; + AsymptoticStatus status = worldtube_sample(source, end, t, &sample); + if (status != ASYMPTOTIC_OK) + return status; + double x_cur[3]; + for (int i = 0; i < 3; ++i) + x_cur[i] = x_frame[i] + w_frame[i] * s; + if (!sample.velocity_constant) { + double s_rel; + status = minkowski_generic_first_entry(source, end, t, x_cur, w_frame, + &s_rel); + if (status != ASYMPTOTIC_OK) + return status; + minkowski_route_entry(end, t0, x_frame, w_frame, s + s_rel, end_id, + route); + return ASYMPTOTIC_OK; + } + double c_frame[3], v_frame[3], d[3], q[3]; + backend_position_to_frame(end, sample.center, c_frame); + backend_vector_to_frame(end, sample.velocity, v_frame); + for (int i = 0; i < 3; ++i) { + d[i] = x_cur[i] - c_frame[i]; + q[i] = w_frame[i] + v_frame[i]; + } + const double boundary = spacetime_escape_worldtube_next_segment( + source, end->end_id, t); + const double s_segment = isfinite(boundary) ? (t - boundary) : INFINITY; + if (!(s_segment >= 0.0)) + return ASYMPTOTIC_INVALID; + double sigma; + const int hit = solve_entry_quadratic(d, q, sample.radius, + sample.radius_rate, &sigma); + /* The backend constructor guarantees R > 0 throughout every segment, so + * a root inside the segment is a real entry. A root past the segment + * boundary is not adopted here; the next segment is sampled instead. + * The cheap R > 0 test at the root guards against a backend that + * bypasses its constructor. */ + if (hit && sigma >= 0.0 && sigma <= s_segment) { + if (sample.radius - sample.radius_rate * sigma <= 0.0) + return ASYMPTOTIC_INVALID; + minkowski_route_entry(end, t0, x_frame, w_frame, s + sigma, end_id, + route); + return ASYMPTOTIC_OK; + } + if (!isfinite(s_segment)) { + /* Open final segment with no entry: a genuine miss. */ + minkowski_route_escaped(end, w_frame, end_id, route); + return ASYMPTOTIC_OK; + } + const double next_s = t0 - nextafter(boundary, -INFINITY); + if (!(next_s > s)) + return ASYMPTOTIC_INVALID; + s = next_s; + } + /* Segment budget exhausted without a conclusion: never disguise this as an + * escape. */ + return ASYMPTOTIC_INVALID; +} + +static AsymptoticStatus schwarzschild_route( + const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end, + const MetricData *metric, const GeodesicRayState *state, + AsymptoticRoute *route) { + SpacetimeEscapeWorldtubeSample sample; + const AsymptoticStatus sample_status = worldtube_sample( + source, end, state->coordinate_time, &sample); + if (sample_status != ASYMPTOTIC_OK) + return sample_status; + /* The analytic monopole exterior only covers a fixed, concentric sphere. */ + const double center_tol = 1e-12 * fmax(1.0, sample.radius); + for (int i = 0; i < 3; ++i) { + if (fabs(sample.center[i] - end->frame_origin[i]) > center_tol || + fabs(sample.velocity[i]) > 1e-12 || + !(sample.radius > 2.0 * end->mass)) + return ASYMPTOTIC_UNSUPPORTED; + } + if (fabs(sample.radius_rate) > 1e-12) + return ASYMPTOTIC_UNSUPPORTED; + if (sample.radius / end->mass < 64.0) + return ASYMPTOTIC_UNSUPPORTED; + SchwarzschildCanonical camera; + if (asymptotic_schwarzschild_canonical_from_state( + end, metric, state->coordinate_time, state->x, state->Pi, + state->log_alpha_p0, &camera)) + return ASYMPTOTIC_INVALID; + SchwarzschildRouteKind kind = SCH_ROUTE_UNSUPPORTED; + double activate_t = 0.0, x[3] = {0.0, 0.0, 0.0}, Pi[3] = {0.0, 0.0, 0.0}; + double log_alpha_p0 = 0.0, n_inf[3] = {0.0, 0.0, 0.0}, frequency = 0.0; + if (asymptotic_schwarzschild_preroute( + end, sample.radius / end->mass, &camera, &kind, &activate_t, x, Pi, + &log_alpha_p0, n_inf, &frequency)) + return ASYMPTOTIC_INVALID; + if (kind == SCH_ROUTE_UNSUPPORTED) + return ASYMPTOTIC_UNSUPPORTED; + if (kind == SCH_ROUTE_TIME_RANGE_EXHAUSTED) { + route->kind = ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED; + route->end_id = end->end_id; + return ASYMPTOTIC_TIME_RANGE_EXHAUSTED; + } + route->end_id = end->end_id; + if (kind == SCH_ROUTE_ENTRY) { + route->kind = ASYMPTOTIC_ROUTE_ENTRY; + route->activate_t = activate_t; + for (int i = 0; i < 3; ++i) { + route->x[i] = x[i]; + route->Pi[i] = Pi[i]; + } + route->log_alpha_p0 = log_alpha_p0; + return ASYMPTOTIC_OK; + } + route->kind = ASYMPTOTIC_ROUTE_ESCAPED; + for (int i = 0; i < 3; ++i) + route->n_infinity[i] = n_inf[i]; + route->frequency_ratio = frequency; + return ASYMPTOTIC_OK; +} + +AsymptoticStatus asymptotic_route_camera(const SpacetimeSource *source, + const ObserverState *observer, + const double direction[3], + AsymptoticRoute *route) { + if (source == NULL || observer == NULL || direction == NULL || route == NULL) + return ASYMPTOTIC_INVALID; + *route = (AsymptoticRoute){.kind = ASYMPTOTIC_ROUTE_INVALID, + .end_id = SPACETIME_END_NONE}; + MetricData metric; + if (spacetime_eval(source, observer->coordinate_time, + observer->coordinate_position, &metric)) + return ASYMPTOTIC_INVALID; + GeodesicRayState state; + if (geodesic_initialize_past_ray_metric(&metric, observer, direction, &state)) + return ASYMPTOTIC_INVALID; + const size_t count = spacetime_asymptotic_end_count(source); + if (count == 0) { + route->kind = ASYMPTOTIC_ROUTE_INSIDE; + route->end_id = SPACETIME_END_NONE; + route->activate_t = state.coordinate_time; + for (int i = 0; i < 3; ++i) { + route->x[i] = state.x[i]; + route->Pi[i] = state.Pi[i]; + } + route->log_alpha_p0 = state.log_alpha_p0; + return ASYMPTOTIC_OK; + } + /* A backend that declares ends must describe them consistently and use a + * supported exterior; otherwise the protocol is broken and no route may be + * fabricated (not even an "inside" one). */ + for (size_t i = 0; i < count; ++i) { + SpacetimeAsymptoticEnd end; + if (spacetime_asymptotic_end(source, i, &end)) + return ASYMPTOTIC_INVALID; + if (end.exterior_kind != ASYMPTOTIC_EXTERIOR_MINKOWSKI && + end.exterior_kind != ASYMPTOTIC_EXTERIOR_SCHWARZSCHILD_MONOPOLE) + return ASYMPTOTIC_UNSUPPORTED; + } + double past_w[3]; + { + double inv[3][3]; + if (invert3(metric.gamma, inv)) + return ASYMPTOTIC_INVALID; + for (int i = 0; i < 3; ++i) { + double dxdt = -metric.beta[i]; + for (int j = 0; j < 3; ++j) + dxdt += metric.alpha * inv[i][j] * state.Pi[j]; + past_w[i] = -dxdt; + } + } + for (size_t i = 0; i < count; ++i) { + SpacetimeAsymptoticEnd end; + if (spacetime_asymptotic_end(source, i, &end)) + return ASYMPTOTIC_INVALID; + double value, slope; + const AsymptoticStatus ws = worldtube_value_and_slope( + source, &end, state.coordinate_time, state.x, past_w, &value, &slope); + if (ws == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) { + route->kind = ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED; + route->end_id = end.end_id; + return ASYMPTOTIC_TIME_RANGE_EXHAUSTED; + } + if (ws != ASYMPTOTIC_OK) + return ASYMPTOTIC_INVALID; + if (value < 0.0 || (value == 0.0 && slope < 0.0)) { + route->kind = ASYMPTOTIC_ROUTE_INSIDE; + route->end_id = end.end_id; + route->activate_t = state.coordinate_time; + for (int k = 0; k < 3; ++k) { + route->x[k] = state.x[k]; + route->Pi[k] = state.Pi[k]; + } + route->log_alpha_p0 = state.log_alpha_p0; + return ASYMPTOTIC_OK; + } + } + + AsymptoticPhotonState canonical; + int have_entry = 0, have_miss = 0; + double best_s = INFINITY; + AsymptoticRoute best = {.kind = ASYMPTOTIC_ROUTE_INVALID}; + SpacetimeEndId first_end = SPACETIME_END_NONE; + double miss_n_inf[3] = {0.0, 0.0, 0.0}; + for (size_t i = 0; i < count; ++i) { + SpacetimeAsymptoticEnd end; + if (spacetime_asymptotic_end(source, i, &end)) + return ASYMPTOTIC_INVALID; + if (first_end == SPACETIME_END_NONE) + first_end = end.end_id; + if (end.exterior_kind == ASYMPTOTIC_EXTERIOR_SCHWARZSCHILD_MONOPOLE) + return schwarzschild_route(source, &end, &metric, &state, route); + if (end.exterior_kind != ASYMPTOTIC_EXTERIOR_MINKOWSKI) + return ASYMPTOTIC_UNSUPPORTED; + if (asymptotic_canonical_from_backend(source, end.end_id, &metric, + state.coordinate_time, state.x, + state.Pi, state.log_alpha_p0, + &canonical)) + return ASYMPTOTIC_INVALID; + AsymptoticRoute candidate = {.kind = ASYMPTOTIC_ROUTE_INVALID}; + const AsymptoticStatus status = minkowski_preroute( + source, &end, state.coordinate_time, canonical.x, canonical.w, + end.end_id, &candidate); + if (status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) { + route->kind = ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED; + route->end_id = end.end_id; + return ASYMPTOTIC_TIME_RANGE_EXHAUSTED; + } + if (status != ASYMPTOTIC_OK) + return status; + if (candidate.kind == ASYMPTOTIC_ROUTE_ENTRY) { + const double s = state.coordinate_time - candidate.activate_t; + if (!have_entry || s < best_s) { + have_entry = 1; + best_s = s; + best = candidate; + } + } else if (!have_miss) { + have_miss = 1; + for (int k = 0; k < 3; ++k) + miss_n_inf[k] = candidate.n_infinity[k]; + } + } + if (have_entry) { + *route = best; + route->log_alpha_p0 = state.log_alpha_p0; + return ASYMPTOTIC_OK; + } + if (!have_miss) { + /* Declared ends exist but none produced a route: broken protocol. */ + return ASYMPTOTIC_INVALID; + } + route->kind = ASYMPTOTIC_ROUTE_ESCAPED; + route->end_id = first_end; + for (int k = 0; k < 3; ++k) + route->n_infinity[k] = miss_n_inf[k]; + route->frequency_ratio = exp(-state.log_alpha_p0); + return ASYMPTOTIC_OK; +} + +AsymptoticStatus asymptotic_finish_escape(const SpacetimeSource *source, + SpacetimeEndId end_id, double t, + const double x[3], const double Pi[3], + double log_alpha_p0, + RayEndpoint *endpoint) { + if (source == NULL || x == NULL || Pi == NULL || endpoint == NULL) + return ASYMPTOTIC_INVALID; + SpacetimeAsymptoticEnd end; + if (find_end(source, end_id, &end)) + return ASYMPTOTIC_INVALID; + if (end.exterior_kind == ASYMPTOTIC_EXTERIOR_SCHWARZSCHILD_MONOPOLE) { + MetricData sch_metric; + if (spacetime_eval(source, t, x, &sch_metric)) + return ASYMPTOTIC_INVALID; + SchwarzschildCanonical canonical; + if (asymptotic_schwarzschild_canonical_from_state( + &end, &sch_metric, t, x, Pi, log_alpha_p0, &canonical)) + return ASYMPTOTIC_INVALID; + double n_inf[3], frequency; + if (asymptotic_schwarzschild_finish(&end, &canonical, n_inf, &frequency)) + return ASYMPTOTIC_INVALID; + for (int i = 0; i < 3; ++i) + endpoint->n_infinity[i] = n_inf[i]; + endpoint->frequency_ratio = frequency; + endpoint->end_id = end_id; + endpoint->status = RAY_ENDPOINT_ESCAPED; + endpoint->magnification = 1.0; + return ASYMPTOTIC_OK; + } + if (end.exterior_kind != ASYMPTOTIC_EXTERIOR_MINKOWSKI) + return ASYMPTOTIC_UNSUPPORTED; + MetricData metric; + double inv[3][3]; + if (spacetime_eval(source, t, x, &metric) || invert3(metric.gamma, inv)) + return ASYMPTOTIC_INVALID; + double n[3] = {0.0, 0.0, 0.0}; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + n[i] -= inv[i][j] * Pi[j]; + if (normalize3(n) <= 0.0) + return ASYMPTOTIC_INVALID; + double beta_dot_pi = 0.0; + for (int i = 0; i < 3; ++i) + beta_dot_pi += metric.beta[i] * Pi[i]; + const double energy = exp(log_alpha_p0) * (metric.alpha - beta_dot_pi); + if (!isfinite(energy) || energy <= 0.0) + return ASYMPTOTIC_INVALID; + for (int i = 0; i < 3; ++i) + endpoint->n_infinity[i] = n[i]; + endpoint->frequency_ratio = 1.0 / energy; + endpoint->end_id = end_id; + endpoint->status = RAY_ENDPOINT_ESCAPED; + endpoint->magnification = 1.0; + return ASYMPTOTIC_OK; +} diff --git a/src/asymptotic.h b/src/asymptotic.h new file mode 100644 index 0000000..c63ca39 --- /dev/null +++ b/src/asymptotic.h @@ -0,0 +1,82 @@ +#ifndef ASYMPTOTIC_H +#define ASYMPTOTIC_H + +#include "geodesic.h" +#include "spacetime.h" + +/* Status codes for the common asymptotic-exterior module. Unsupported and + * exhausted are reported explicitly; callers must not turn them into a + * plausible-looking escape. */ +typedef enum { + ASYMPTOTIC_INVALID = -1, + ASYMPTOTIC_OK = 0, + ASYMPTOTIC_UNSUPPORTED = 1, + ASYMPTOTIC_TIME_RANGE_EXHAUSTED = 2 +} AsymptoticStatus; + +/* Unified canonical photon state in the asymptotic reference frame. `w` is + * the unit past-propagation direction: along the renderer's backward + * integration the spatial position moves as x(s) = x0 + s w, s = t0 - t. */ +typedef struct { + SpacetimeEndId end_id; + double t; + double x[3]; + double w[3]; + double log_alpha_p0; +} AsymptoticPhotonState; + +typedef enum { + ASYMPTOTIC_ROUTE_INSIDE, /* camera in a worldtube: activate at the camera */ + ASYMPTOTIC_ROUTE_ENTRY, /* camera outside, entry event produced */ + ASYMPTOTIC_ROUTE_ESCAPED, + ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED, + ASYMPTOTIC_ROUTE_INVALID +} AsymptoticRouteKind; + +typedef struct { + AsymptoticRouteKind kind; + SpacetimeEndId end_id; + /* Activation state, backend coordinates, for INSIDE and ENTRY. */ + double activate_t; + double x[3]; + double Pi[3]; + double log_alpha_p0; + /* Terminal infinity endpoint for ESCAPED. */ + double n_infinity[3]; + double frequency_ratio; +} AsymptoticRoute; + +/* Pre-route one camera ray against every declared end's worldtube. */ +AsymptoticStatus asymptotic_route_camera(const SpacetimeSource *source, + const ObserverState *observer, + const double direction[3], + AsymptoticRoute *route); + +/* Directed inside->outside crossing helper for the geodesic lifecycle. + * Returns ASYMPTOTIC_OK, ASYMPTOTIC_TIME_RANGE_EXHAUSTED (the backend cannot + * describe the worldtube at this time), or ASYMPTOTIC_INVALID. */ +int asymptotic_worldtube_value(const SpacetimeSource *source, + SpacetimeEndId end_id, double t, + const double x[3], double *value); +/* Finish an interior inside->outside crossing to an infinity endpoint. + * Returns ASYMPTOTIC_UNSUPPORTED for exterior models not implemented yet. */ +AsymptoticStatus asymptotic_finish_escape(const SpacetimeSource *source, + SpacetimeEndId end_id, double t, + const double x[3], const double Pi[3], + double log_alpha_p0, + RayEndpoint *endpoint); + +/* Canonical <-> backend bridge, valid only inside a supported exterior. */ +int asymptotic_canonical_from_backend(const SpacetimeSource *source, + SpacetimeEndId end_id, + const MetricData *metric, double t, + const double x[3], const double Pi[3], + double log_alpha_p0, + AsymptoticPhotonState *out); +int asymptotic_backend_from_canonical(const SpacetimeSource *source, + const MetricData *metric, + const AsymptoticPhotonState *canonical, + double x[3], double Pi[3], + double *log_alpha_p0); + +#endif diff --git a/src/asymptotic_gl48.h b/src/asymptotic_gl48.h new file mode 100644 index 0000000..b0b8464 --- /dev/null +++ b/src/asymptotic_gl48.h @@ -0,0 +1,42 @@ +#ifndef ASYMPTOTIC_GL48_H +#define ASYMPTOTIC_GL48_H + +/* 48-point Gauss-Legendre nodes and weights on [-1, 1], used only for the + * bounded residual of the Schwarzschild coordinate-time transfer. Generated + * with numpy.polynomial.legendre.leggauss(48); double precision. */ +static const double gl48_nodes[48] = { + -0.99877100725242607, -0.99353017226635076, -0.98412458372282685, + -0.97059159254624727, -0.9529877031604308, -0.93138669070655433, + -0.90587913671556963, -0.87657202027424785, -0.84358826162439349, + -0.80706620402944262, -0.76715903251574036, -0.72403413092381463, + -0.67787237963266389, -0.6288673967765136, -0.57722472608397268, + -0.523160974722233, -0.46690290475095841, -0.40868648199071672, + -0.34875588629216075, -0.28736248735545555, -0.22476379039468905, + -0.16122235606889174, -0.097004699209462697, -0.032380170962869367, + 0.032380170962869367, 0.097004699209462697, 0.16122235606889174, + 0.22476379039468905, 0.28736248735545555, 0.34875588629216075, + 0.40868648199071672, 0.46690290475095841, 0.523160974722233, + 0.57722472608397268, 0.6288673967765136, 0.67787237963266389, + 0.72403413092381463, 0.76715903251574036, 0.80706620402944262, + 0.84358826162439349, 0.87657202027424785, 0.90587913671556963, + 0.93138669070655433, 0.9529877031604308, 0.97059159254624727, + 0.98412458372282685, 0.99353017226635076, 0.99877100725242607}; +static const double gl48_weights[48] = { + 0.0031533460523098418, 0.0073275539012758505, 0.011477234579234699, + 0.015579315722943481, 0.019616160457356105, 0.023570760839324009, + 0.027426509708357052, 0.031167227832798117, 0.034777222564770421, + 0.038241351065830473, 0.041545082943464533, 0.044674560856694245, + 0.04761665849249027, 0.050359035553854216, 0.052890189485193424, + 0.05519950369998404, 0.057277292100402881, 0.059114839698395358, + 0.060704439165893562, 0.062039423159892415, 0.063114192286253756, + 0.06392423858464788, 0.06446616443594981, 0.064737696812683626, + 0.064737696812683626, 0.06446616443594981, 0.06392423858464788, + 0.063114192286253756, 0.062039423159892415, 0.060704439165893562, + 0.059114839698395358, 0.057277292100402881, 0.05519950369998404, + 0.052890189485193424, 0.050359035553854216, 0.04761665849249027, + 0.044674560856694245, 0.041545082943464533, 0.038241351065830473, + 0.034777222564770421, 0.031167227832798117, 0.027426509708357052, + 0.023570760839324009, 0.019616160457356105, 0.015579315722943481, + 0.011477234579234699, 0.0073275539012758505, 0.0031533460523098418}; + +#endif diff --git a/src/asymptotic_schwarzschild.c b/src/asymptotic_schwarzschild.c new file mode 100644 index 0000000..7f6a439 --- /dev/null +++ b/src/asymptotic_schwarzschild.c @@ -0,0 +1,528 @@ +#include "asymptotic_schwarzschild.h" + +#include +#include +#include + +static const double kPi = 3.14159265358979323846; + +static double dot3(const double a[3], const double b[3]) { + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; +} + +static void cross3(const double a[3], const double b[3], double out[3]) { + out[0] = a[1] * b[2] - a[2] * b[1]; + out[1] = a[2] * b[0] - a[0] * b[2]; + out[2] = a[0] * b[1] - a[1] * b[0]; +} + +static double normalize3(double v[3]) { + const double length = sqrt(dot3(v, v)); + if (length > 0.0) + for (int i = 0; i < 3; ++i) + v[i] /= length; + return length; +} + +/* ------------------------------------------------------------------------- * + * Complex arithmetic and Carlson R_F. + * ------------------------------------------------------------------------- */ + +typedef struct { + double re, im; +} cs; + +static cs cs_add(cs a, cs b) { return (cs){a.re + b.re, a.im + b.im}; } +static cs cs_sub(cs a, cs b) { return (cs){a.re - b.re, a.im - b.im}; } +static cs cs_mul(cs a, cs b) { + return (cs){a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re}; +} +static cs cs_scale(cs a, double s) { return (cs){a.re * s, a.im * s}; } +static double cs_abs(cs a) { return hypot(a.re, a.im); } + +static cs cs_inv(cs z) { + const double d = z.re * z.re + z.im * z.im; + return (cs){z.re / d, -z.im / d}; +} + +static cs cs_sqrt(cs z) { + const double r = hypot(z.re, z.im); + double re = sqrt(0.5 * (r + fabs(z.re))); + double im = sqrt(0.5 * (r - fabs(z.re))); + if (z.re < 0.0) { + const double t = re; + re = im; + im = t; + } + if (z.im < 0.0) + im = -im; + return (cs){re, im}; +} + +static cs cs_cbrt(cs z) { + /* Cardano needs the real cube root of real arguments (disc > 0); the + * principal complex root is correct for the conjugate pair (disc < 0). */ + if (z.im == 0.0) + return (cs){cbrt(z.re), 0.0}; + const double r = hypot(z.re, z.im); + const double theta = atan2(z.im, z.re); + const double cr = cbrt(r); + return (cs){cr * cos(theta / 3.0), cr * sin(theta / 3.0)}; +} + +static cs rf_naive(cs x, cs y, cs z) { + for (int iteration = 0; iteration < 80; ++iteration) { + const cs sx = cs_sqrt(x), sy = cs_sqrt(y), sz = cs_sqrt(z); + const cs lambda = + cs_add(cs_add(cs_mul(sx, sy), cs_mul(sy, sz)), cs_mul(sz, sx)); + x = cs_scale(cs_add(x, lambda), 0.25); + y = cs_scale(cs_add(y, lambda), 0.25); + z = cs_scale(cs_add(z, lambda), 0.25); + const cs a = cs_scale(cs_add(cs_add(x, y), z), 1.0 / 3.0); + const cs X = cs_sub((cs){1.0, 0.0}, cs_mul(x, cs_inv(a))); + const cs Y = cs_sub((cs){1.0, 0.0}, cs_mul(y, cs_inv(a))); + const cs Z = cs_sub((cs){1.0, 0.0}, cs_mul(z, cs_inv(a))); + if (fmax(fmax(cs_abs(X), cs_abs(Y)), cs_abs(Z)) < 1e-12) { + const cs e2 = cs_add(cs_add(cs_mul(X, Y), cs_mul(Y, Z)), + cs_mul(Z, X)); + const cs e3 = cs_mul(cs_mul(X, Y), Z); + const cs series = cs_add( + cs_add((cs){1.0, 0.0}, cs_scale(cs_mul(e2, e3), -3.0 / 44.0)), + cs_add(cs_scale(cs_mul(e2, e2), 1.0 / 24.0), + cs_add(cs_scale(e2, -1.0 / 10.0), cs_scale(e3, 1.0 / 14.0)))); + return cs_mul(series, cs_inv(cs_sqrt(a))); + } + } + return (cs){NAN, NAN}; +} + +/* R_F via Carlson duplication. The three-real-root branch is handled by the + * real Legendre form, so the only complex calls here come from the conjugate + * root pair, whose arguments are off the real axis and take the principal + * square-root branch consistently. */ +static cs rf(cs x, cs y, cs z) { return rf_naive(x, y, z); } + +/* Incomplete elliptic integral of the first kind with parameter m = k^2: + * F(phi,m) = sin(phi) R_F(cos^2 phi, 1 - m sin^2 phi, 1). */ +static double ellipf(double phi, double m) { + const double s = sin(phi), c = cos(phi); + const cs r = rf((cs){c * c, 0.0}, (cs){1.0 - m * s * s, 0.0}, + (cs){1.0, 0.0}); + return s * r.re; +} + +/* Leading-order estimate of phi for a candidate ordered real-root branch. */ +static double phi_three_real(double u0, double A, double B, double C) { + if (!(u0 > A) || !(u0 < B) || !(A < B) || !(B < C)) + return NAN; + const double sA = sqrt((0.0 - A) / (B - A)); + const double s0 = sqrt((u0 - A) / (B - A)); + const double m = (B - A) / (C - A); + return sqrt(2.0) / sqrt(C - A) * (ellipf(asin(s0), m) - ellipf(asin(sA), m)); +} + +/* Roots of 2 beta^2 u^3 - beta^2 u^2 + 1 = 0 through the depressed cubic + * w^3 + P w + Q = 0 with u = w + 1/6. */ +static void cubic_roots(double beta, cs e[3]) { + const double b2 = beta * beta; + const double c = 1.0 / (2.0 * b2); + const double P = -1.0 / 12.0; + const double Q = c - 1.0 / 108.0; + const double halfQ = 0.5 * Q; + const cs disc = (cs){halfQ * halfQ + (P * P * P) / 27.0, 0.0}; + const cs sq = cs_sqrt(disc); + const cs u1 = cs_cbrt(cs_add((cs){-halfQ, 0.0}, sq)); + const cs u2 = cs_cbrt(cs_add((cs){-halfQ, 0.0}, cs_scale(sq, -1.0))); + const cs omega = (cs){cos(2.0 * kPi / 3.0), sin(2.0 * kPi / 3.0)}; + const cs omega2 = cs_mul(omega, omega); + e[0] = cs_add(cs_add(u1, u2), (cs){1.0 / 6.0, 0.0}); + e[1] = cs_add(cs_add(cs_mul(omega, u1), cs_mul(omega2, u2)), + (cs){1.0 / 6.0, 0.0}); + e[2] = cs_add(cs_add(cs_mul(omega2, u1), cs_mul(omega, u2)), + (cs){1.0 / 6.0, 0.0}); + /* Cardano loses relative accuracy in the near-double-root regime. Polish + * the real roots with Newton so the grazing turning root keeps full + * relative precision. */ + for (int i = 0; i < 3; ++i) { + if (fabs(e[i].im) > 1e-9 * fmax(1.0, fabs(e[i].re))) + continue; + double u = e[i].re; + for (int step = 0; step < 20; ++step) { + const double p = 2.0 * b2 * u * u * u - b2 * u * u + 1.0; + const double dp = 6.0 * b2 * u * u - 2.0 * b2 * u; + if (dp == 0.0) + break; + const double du = p / dp; + u -= du; + if (fabs(du) <= 1e-18 * fmax(1.0, fabs(u))) + break; + } + e[i] = (cs){u, 0.0}; + } +} + +double asymptotic_schwarzschild_phi(double rho, double beta) { + if (!isfinite(rho) || rho <= 0.0 || !isfinite(beta) || beta < 0.0) + return NAN; + if (beta == 0.0) + return 0.0; + const double u0 = 1.0 / rho; + cs e[3]; + cubic_roots(beta, e); + const double imag_tol = 1e-11 * fmax(1.0, fabs(e[0].re)); + if (fabs(e[0].im) < imag_tol && fabs(e[1].im) < imag_tol && + fabs(e[2].im) < imag_tol) { + /* Three real roots: use the real Legendre form, which is accurate up to + * and through the grazing limit. */ + double r[3] = {e[0].re, e[1].re, e[2].re}; + for (int i = 0; i < 2; ++i) + for (int j = i + 1; j < 3; ++j) + if (r[j] < r[i]) { + const double t = r[i]; + r[i] = r[j]; + r[j] = t; + } + const double real_value = phi_three_real(u0, r[0], r[1], r[2]); + if (isfinite(real_value)) + return real_value; + } + const cs a = rf(cs_scale(e[0], -1.0), cs_scale(e[1], -1.0), + cs_scale(e[2], -1.0)); + const cs b = rf(cs_sub((cs){u0, 0.0}, e[0]), + cs_sub((cs){u0, 0.0}, e[1]), + cs_sub((cs){u0, 0.0}, e[2])); + const cs s = cs_scale(cs_sub(a, b), 2.0); + /* The real integral requires a real S; a non-negligible imaginary part + * means the principal branch failed. Report it instead of silently using + * a wrong angle. */ + if (!isfinite(s.re) || fabs(s.im) > 1e-6 * fmax(1.0, fabs(s.re))) + return NAN; + return fabs(s.re) / sqrt(2.0); +} + +double asymptotic_schwarzschild_turning_rho(double beta) { + if (!isfinite(beta) || beta <= 3.0 * sqrt(3.0)) + return INFINITY; + /* Larger positive root of f(rho) = rho^3 - beta^2 rho + 2 beta^2. + * f(3) = 27 - beta^2 < 0 and f(beta+2) > 0, and f is monotone on the + * bracket beyond its local minimum, so a bracketed bisection is safe. */ + const double b2 = beta * beta; + double lo = 3.0, hi = beta + 2.0; + for (int iteration = 0; iteration < 200; ++iteration) { + const double mid = 0.5 * (lo + hi); + const double f = mid * mid * mid - b2 * mid + 2.0 * b2; + if (f < 0.0) + lo = mid; + else + hi = mid; + if (hi - lo <= 4.0 * DBL_EPSILON * hi) + break; + } + double rho = 0.5 * (lo + hi); + for (int step = 0; step < 20; ++step) { + const double f = rho * rho * rho - b2 * rho + 2.0 * b2; + const double fp = 3.0 * rho * rho - b2; + if (fp == 0.0) + break; + const double next = rho - f / fp; + if (!(next > lo && next < hi)) + break; + rho = next; + } + return rho; +} + +/* ------------------------------------------------------------------------- * + * Coordinate-time transfer (ingoing Kerr-Schild time, M = 1 units). + * ------------------------------------------------------------------------- */ + +#include "asymptotic_gl48.h" + +static double sch_i2(double u_cam, double u_R, double beta) { + /* I2 = int_{u_cam}^{u_R} du / (1 + sqrt(P(u))). Substitute + * u = u_R - (u_R - u_cam) t^2 to remove the grazing branch point. */ + const double span = u_R - u_cam; + if (!(span > 0.0)) + return 0.0; + double sum = 0.0; + for (int i = 0; i < 48; ++i) { + const double t = 0.5 * (gl48_nodes[i] + 1.0); + const double u = u_R - span * t * t; + const double P = 1.0 - beta * beta * u * u + 2.0 * beta * beta * u * u * u; + const double f = 1.0 / (1.0 + sqrt(P)); + sum += gl48_weights[i] * f * 2.0 * span * t; + } + return 0.5 * sum; +} + +/* Positive coordinate time to travel outward from R to rho_cam. */ +static double sch_time_transfer(double rho_cam, double rho_R, double beta) { + const double u_cam = 1.0 / rho_cam, u_R = 1.0 / rho_R; + const double dphi = + asymptotic_schwarzschild_phi(rho_R, beta) - + asymptotic_schwarzschild_phi(rho_cam, beta); + const double elementary = + (rho_cam - rho_R) + 4.0 * log(u_R / u_cam) - + 4.0 * log((1.0 - 2.0 * u_R) / (1.0 - 2.0 * u_cam)); + return elementary + (beta * dphi - beta * beta * sch_i2(u_cam, u_R, beta)); +} + +/* ------------------------------------------------------------------------- * + * Rotation and canonical <-> backend bridging. + * ------------------------------------------------------------------------- */ + +static void rotate_axis(const double v[3], const double axis[3], double angle, + double out[3]) { + const double c = cos(angle), s = sin(angle); + double cross[3]; + cross3(axis, v, cross); + const double adotv = dot3(axis, v); + for (int i = 0; i < 3; ++i) + out[i] = v[i] * c + cross[i] * s + axis[i] * adotv * (1.0 - c); +} + +/* The algebraic monopole formulas below assume the asymptotic frame axes are + * the backend Cartesian axes; a rotated frame would require rotating the + * momentum and the worldtube. */ +static int sch_frame_is_aligned(const SpacetimeAsymptoticEnd *end) { + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) { + const double expected = i == j ? 1.0 : 0.0; + if (fabs(end->frame_axes[i][j] - expected) > 1e-12) + return 0; + } + return 1; +} + +int asymptotic_schwarzschild_canonical_from_state( + const SpacetimeAsymptoticEnd *end, const MetricData *metric, double t, + const double x[3], const double Pi[3], double log_alpha_p0, + SchwarzschildCanonical *out) { + if (end == NULL || metric == NULL || out == NULL || end->mass <= 0.0 || + !sch_frame_is_aligned(end)) + return -1; + double beta_dot_pi = 0.0; + for (int i = 0; i < 3; ++i) + beta_dot_pi += metric->beta[i] * Pi[i]; + const double energy = exp(log_alpha_p0) * (metric->alpha - beta_dot_pi); + if (!isfinite(energy) || energy <= 0.0) + return -1; + double rel[3]; + for (int i = 0; i < 3; ++i) + rel[i] = x[i] - end->frame_origin[i]; + const double radius = sqrt(dot3(rel, rel)); + if (!(radius > 0.0)) + return -1; + double Lvec[3]; + cross3(rel, Pi, Lvec); + const double Lmag = sqrt(dot3(Lvec, Lvec)); + const double denom = metric->alpha - beta_dot_pi; + out->end_id = end->end_id; + out->t = t; + out->rho = radius / end->mass; + out->energy = energy; + for (int i = 0; i < 3; ++i) + out->rhat[i] = rel[i] / radius; + if (Lmag > 0.0) { + out->beta = (Lmag / denom) / end->mass; + for (int i = 0; i < 3; ++i) + out->Lhat[i] = Lvec[i] / Lmag; + } else { + out->beta = 0.0; + out->Lhat[0] = out->Lhat[1] = out->Lhat[2] = 0.0; + } + double inv[3][3]; + const double det = + metric->gamma[0][0] * (metric->gamma[1][1] * metric->gamma[2][2] - + metric->gamma[1][2] * metric->gamma[2][1]) - + metric->gamma[0][1] * (metric->gamma[1][0] * metric->gamma[2][2] - + metric->gamma[1][2] * metric->gamma[2][0]) + + metric->gamma[0][2] * (metric->gamma[1][0] * metric->gamma[2][1] - + metric->gamma[1][1] * metric->gamma[2][0]); + inv[0][0] = (metric->gamma[1][1] * metric->gamma[2][2] - + metric->gamma[1][2] * metric->gamma[2][1]) / det; + inv[0][1] = (metric->gamma[0][2] * metric->gamma[2][1] - + metric->gamma[0][1] * metric->gamma[2][2]) / det; + inv[0][2] = (metric->gamma[0][1] * metric->gamma[1][2] - + metric->gamma[0][2] * metric->gamma[1][1]) / det; + inv[1][0] = (metric->gamma[1][2] * metric->gamma[2][0] - + metric->gamma[1][0] * metric->gamma[2][2]) / det; + inv[1][1] = (metric->gamma[0][0] * metric->gamma[2][2] - + metric->gamma[0][2] * metric->gamma[2][0]) / det; + inv[1][2] = (metric->gamma[0][2] * metric->gamma[1][0] - + metric->gamma[0][0] * metric->gamma[1][2]) / det; + inv[2][0] = (metric->gamma[1][0] * metric->gamma[2][1] - + metric->gamma[1][1] * metric->gamma[2][0]) / det; + inv[2][1] = (metric->gamma[0][1] * metric->gamma[2][0] - + metric->gamma[0][0] * metric->gamma[2][1]) / det; + inv[2][2] = (metric->gamma[0][0] * metric->gamma[1][1] - + metric->gamma[0][1] * metric->gamma[1][0]) / det; + double dxdt[3]; + for (int i = 0; i < 3; ++i) { + dxdt[i] = -metric->beta[i]; + for (int j = 0; j < 3; ++j) + dxdt[i] += metric->alpha * inv[i][j] * Pi[j]; + } + double radial = 0.0; + for (int i = 0; i < 3; ++i) + radial += -dxdt[i] * out->rhat[i]; + out->radial_sign = radial > 0.0 ? 1 : (radial < 0.0 ? -1 : 0); + return 0; +} + +int asymptotic_schwarzschild_state_from_canonical( + const SpacetimeAsymptoticEnd *end, const SchwarzschildCanonical *c, + double x[3], double Pi[3], double *log_alpha_p0) { + if (end == NULL || c == NULL || x == NULL || Pi == NULL || + !sch_frame_is_aligned(end)) + return -1; + const double rho = c->rho; + if (!(rho > 2.0)) + return -1; + const double Q = 1.0 - c->beta * c->beta * (1.0 - 2.0 / rho) / (rho * rho); + if (!(Q >= 0.0)) + return -1; + const double sqrtQ = sqrt(Q); + double e_phi[3] = {0.0, 0.0, 0.0}; + if (c->beta > 0.0) + cross3(c->Lhat, c->rhat, e_phi); + const double s_aff = -(double)c->radial_sign; /* physical (future) radial */ + const double kr = s_aff * c->energy * sqrtQ; + const double ktang = c->beta * c->energy / rho; + double kvec[3]; + for (int i = 0; i < 3; ++i) + kvec[i] = kr * c->rhat[i] + ktang * e_phi[i]; + const double kt_s = c->energy / (1.0 - 2.0 / rho); + const double kt_ks = kt_s + (2.0 / (rho - 2.0)) * kr; + const double alpha = 1.0 / sqrt(1.0 + 2.0 / rho); + const double ak0 = alpha * kt_ks; + if (!isfinite(ak0) || ak0 <= 0.0) + return -1; + for (int i = 0; i < 3; ++i) { + x[i] = end->frame_origin[i] + end->mass * rho * c->rhat[i]; + const double kcov = kvec[i] + (2.0 / rho) * c->rhat[i] * (kr + kt_ks); + Pi[i] = kcov / ak0; + } + if (log_alpha_p0 != NULL) + *log_alpha_p0 = log(ak0); + return 0; +} + +int asymptotic_schwarzschild_finish(const SpacetimeAsymptoticEnd *end, + const SchwarzschildCanonical *canonical, + double n_infinity[3], + double *frequency_ratio) { + if (end == NULL || canonical == NULL || n_infinity == NULL) + return -1; + const double phi = asymptotic_schwarzschild_phi(canonical->rho, + canonical->beta); + if (!isfinite(phi)) + return -1; + if (canonical->beta > 0.0) { + double e_phi[3]; + cross3(canonical->Lhat, canonical->rhat, e_phi); + for (int i = 0; i < 3; ++i) + n_infinity[i] = cos(phi) * canonical->rhat[i] - + sin(phi) * e_phi[i]; + } else { + const double s = canonical->radial_sign >= 0 ? 1.0 : -1.0; + for (int i = 0; i < 3; ++i) + n_infinity[i] = s * canonical->rhat[i]; + } + normalize3(n_infinity); + if (frequency_ratio != NULL) + *frequency_ratio = 1.0 / canonical->energy; + return 0; +} + +int asymptotic_schwarzschild_preroute( + const SpacetimeAsymptoticEnd *end, double worldtube_radius, + const SchwarzschildCanonical *camera, SchwarzschildRouteKind *kind, + double *activate_t, double x[3], double Pi[3], double *log_alpha_p0, + double n_infinity[3], double *frequency_ratio) { + if (end == NULL || camera == NULL || kind == NULL) + return -1; + const double R = worldtube_radius; + if (!(R > 2.0) || !(camera->rho >= R)) + return -1; + if (R / 1.0 < 64.0) { + *kind = SCH_ROUTE_UNSUPPORTED; + return 0; + } + const double beta_R = R / sqrt(1.0 - 2.0 / R); + + if (camera->radial_sign >= 0) { + /* Past propagation is outward or tangent: no entry, immediate infinity + * endpoint. (radial_sign == 0 means the camera is on the boundary with a + * tangent ray.) */ + const double phi = asymptotic_schwarzschild_phi(camera->rho, camera->beta); + if (!isfinite(phi)) { + *kind = SCH_ROUTE_UNSUPPORTED; + return 0; + } + if (camera->beta > 0.0) { + double e_phi[3]; + cross3(camera->Lhat, camera->rhat, e_phi); + for (int i = 0; i < 3; ++i) + n_infinity[i] = cos(phi) * camera->rhat[i] - sin(phi) * e_phi[i]; + } else { + for (int i = 0; i < 3; ++i) + n_infinity[i] = camera->rhat[i]; + } + normalize3(n_infinity); + *frequency_ratio = 1.0 / camera->energy; + *kind = SCH_ROUTE_ESCAPED; + return 0; + } + + if (camera->beta < beta_R) { + const double dphi = asymptotic_schwarzschild_phi(R, camera->beta) - + asymptotic_schwarzschild_phi(camera->rho, camera->beta); + double rhat_entry[3]; + if (camera->beta > 0.0) + rotate_axis(camera->rhat, camera->Lhat, -dphi, rhat_entry); + else + for (int i = 0; i < 3; ++i) + rhat_entry[i] = camera->rhat[i]; + SchwarzschildCanonical entry = *camera; + entry.rho = R; + for (int i = 0; i < 3; ++i) + entry.rhat[i] = rhat_entry[i]; + entry.radial_sign = -1; + if (asymptotic_schwarzschild_state_from_canonical(end, &entry, x, Pi, + log_alpha_p0)) + return -1; + const double T = + sch_time_transfer(camera->rho, R, camera->beta); + *activate_t = camera->t - end->mass * T; + *kind = SCH_ROUTE_ENTRY; + return 0; + } + + /* Inward but misses: turn before R and escape. */ + const double rho_turn = asymptotic_schwarzschild_turning_rho(camera->beta); + if (!isfinite(rho_turn) || rho_turn > camera->rho) { + *kind = SCH_ROUTE_UNSUPPORTED; + return 0; + } + const double total = + 2.0 * asymptotic_schwarzschild_phi(rho_turn, camera->beta) - + asymptotic_schwarzschild_phi(camera->rho, camera->beta); + if (!isfinite(total)) { + *kind = SCH_ROUTE_UNSUPPORTED; + return 0; + } + if (camera->beta > 0.0) { + double e_phi[3]; + cross3(camera->Lhat, camera->rhat, e_phi); + for (int i = 0; i < 3; ++i) + n_infinity[i] = cos(total) * camera->rhat[i] - sin(total) * e_phi[i]; + } else { + for (int i = 0; i < 3; ++i) + n_infinity[i] = camera->rhat[i]; + } + normalize3(n_infinity); + *frequency_ratio = 1.0 / camera->energy; + *kind = SCH_ROUTE_ESCAPED; + return 0; +} diff --git a/src/asymptotic_schwarzschild.h b/src/asymptotic_schwarzschild.h new file mode 100644 index 0000000..66232a0 --- /dev/null +++ b/src/asymptotic_schwarzschild.h @@ -0,0 +1,66 @@ +#ifndef ASYMPTOTIC_SCHWARZSCHILD_H +#define ASYMPTOTIC_SCHWARZSCHILD_H + +#include "geodesic.h" +#include "spacetime.h" + +/* Canonical photon state for a fixed, concentric Schwarzschild monopole + * exterior. All radial quantities are in units of the mass: rho = r / M. + * `Lhat` is the (unit) conserved angular-momentum direction = normalize(x x + * Pi); `beta` is the impact parameter b/M > 0. `radial_sign` is the sign of + * dr/ds along the renderer's past propagation (s = t_camera - t): +1 outward + * into the past, -1 inward into the past. `energy` is E = -p_t with the + * camera normalization E_camera = 1. */ +typedef struct { + SpacetimeEndId end_id; + double t; + double rho; + double rhat[3]; + double Lhat[3]; + double beta; + double energy; + int radial_sign; +} SchwarzschildCanonical; + +/* Angular primitive Phi(rho, beta): the azimuth swept on the outward branch + * from radius rho to infinity. Returns NAN outside the supported domain. */ +double asymptotic_schwarzschild_phi(double rho, double beta); + +/* Larger positive turning radius for the given impact parameter, or INFINITY + * when no turning point exists (beta <= 3 sqrt(3)). */ +double asymptotic_schwarzschild_turning_rho(double beta); + +/* Convert a backend state into the canonical form. `metric` must be the + * Schwarzschild Kerr-Schild metric at (t, x). */ +int asymptotic_schwarzschild_canonical_from_state( + const SpacetimeAsymptoticEnd *end, const MetricData *metric, double t, + const double x[3], const double Pi[3], double log_alpha_p0, + SchwarzschildCanonical *out); + +/* Rebuild the backend state at the stored radius / radial directions. */ +int asymptotic_schwarzschild_state_from_canonical( + const SpacetimeAsymptoticEnd *end, const SchwarzschildCanonical *canonical, + double x[3], double Pi[3], double *log_alpha_p0); + +/* Infinity endpoint for an outward crossing at the canonical radius. */ +int asymptotic_schwarzschild_finish(const SpacetimeAsymptoticEnd *end, + const SchwarzschildCanonical *canonical, + double n_infinity[3], + double *frequency_ratio); + +/* Pre-route a camera ray outside the worldtube. Fills one of the route + * kinds. `worldtube_radius` is R/M. */ +typedef enum { + SCH_ROUTE_ENTRY, + SCH_ROUTE_ESCAPED, + SCH_ROUTE_TIME_RANGE_EXHAUSTED, + SCH_ROUTE_UNSUPPORTED +} SchwarzschildRouteKind; + +int asymptotic_schwarzschild_preroute( + const SpacetimeAsymptoticEnd *end, double worldtube_radius, + const SchwarzschildCanonical *camera, SchwarzschildRouteKind *kind, + double *activate_t, double x[3], double Pi[3], double *log_alpha_p0, + double n_infinity[3], double *frequency_ratio); + +#endif diff --git a/src/frame.c b/src/frame.c index 4c18584..d1706a8 100644 --- a/src/frame.c +++ b/src/frame.c @@ -116,6 +116,7 @@ int frame_lens_mesh_trace(FrameLensMesh *mesh, const SpacetimeSource *spacetime, RayEndpoint endpoint = geodesic_trace_past(spacetime, observer, vertex->camera_direction, trace); vertex->status = endpoint.status; + vertex->end_id = endpoint.end_id; vertex->traced = 1; if (endpoint.status == RAY_ENDPOINT_ESCAPED) { for (int axis = 0; axis < 3; ++axis) @@ -437,11 +438,20 @@ static int all_vertices_traced(const FrameLensMesh *mesh) { static int terminal_mismatch(const LensVertex *a, const LensVertex *b, const LensVertex *c) { - int escaped = 0, captured = 0; + int escaped = 0, captured = 0, have_end = 0; + SpacetimeEndId end = SPACETIME_END_NONE; const LensVertex *vertices[] = {a, b, c}; for (size_t i = 0; i < 3; ++i) { escaped |= vertices[i]->status == RAY_ENDPOINT_ESCAPED; captured |= vertices[i]->status == RAY_ENDPOINT_CAPTURED; + if (vertices[i]->status == RAY_ENDPOINT_ESCAPED) { + if (!have_end) { + end = vertices[i]->end_id; + have_end = 1; + } else if (vertices[i]->end_id != end) { + return 1; /* two different infinity ends must not be interpolated */ + } + } } return escaped && captured; } @@ -576,6 +586,7 @@ int frame_lens_mesh_install_sample(FrameLensMesh *mesh, size_t sample_id, ? &mesh->vertices[sample->vertex_id] : &sample->vertex; vertex->status = endpoint->status; + vertex->end_id = endpoint->end_id; vertex->traced = 1; if (endpoint->status == RAY_ENDPOINT_ESCAPED) { for (int axis = 0; axis < 3; ++axis) @@ -610,6 +621,8 @@ static int discrete_jacobian(const FrameLensMesh *mesh, if (a->status != RAY_ENDPOINT_ESCAPED || b->status != RAY_ENDPOINT_ESCAPED || c->status != RAY_ENDPOINT_ESCAPED) return 0; + if (a->end_id != b->end_id || a->end_id != c->end_id) + return 0; const double image_area = spherical_signed_area( a->camera_direction, b->camera_direction, c->camera_direction); if (!isfinite(image_area) || fabs(image_area) <= 1e-15) diff --git a/src/frame.h b/src/frame.h index 61a6e15..2d0b21c 100644 --- a/src/frame.h +++ b/src/frame.h @@ -15,6 +15,9 @@ typedef struct { double n_infinity[3]; double log_frequency_ratio; RayEndpointStatus status; + /* Asymptotic end this escaped vertex belongs to; a triangle must not + * interpolate across two different ends. */ + SpacetimeEndId end_id; int traced; } LensVertex; diff --git a/src/geodesic.c b/src/geodesic.c index 3adeb89..4262693 100644 --- a/src/geodesic.c +++ b/src/geodesic.c @@ -1,4 +1,5 @@ #include "geodesic.h" +#include "asymptotic.h" #include #include @@ -105,15 +106,14 @@ static int rk4(const MetricSlab *slab, double t, double h, State *s) { return 0; } -int geodesic_initialize_past_ray(const MetricSlab *slab, - const ObserverState *o, const double n[3], - State *s) { - MetricData m; +int geodesic_initialize_past_ray_metric(const MetricData *metric, + const ObserverState *o, + const double n[3], State *s) { + const MetricData *m = metric; + if (m->alpha <= 0) + return -1; double k[4] = {o->tetrad[0][0], o->tetrad[0][1], o->tetrad[0][2], o->tetrad[0][3]}; - if (spacetime_slab_eval(slab, o->coordinate_time, o->coordinate_position, &m) || - m.alpha <= 0) - return -1; for (int a = 0; a < 3; a++) for (int mu = 0; mu < 4; mu++) k[mu] -= n[a] * o->tetrad[a + 1][mu]; @@ -123,15 +123,24 @@ int geodesic_initialize_past_ray(const MetricSlab *slab, s->x[i] = o->coordinate_position[i]; s->Pi[i] = 0; for (int j = 0; j < 3; j++) - s->Pi[i] += m.gamma[i][j] * (k[j + 1] + m.beta[j] * k[0]); - s->Pi[i] /= m.alpha * k[0]; + s->Pi[i] += m->gamma[i][j] * (k[j + 1] + m->beta[j] * k[0]); + s->Pi[i] /= m->alpha * k[0]; } - s->log_alpha_p0 = log(m.alpha * k[0]); + s->log_alpha_p0 = log(m->alpha * k[0]); s->coordinate_time = o->coordinate_time; s->steps = 0; return isfinite(s->log_alpha_p0) ? 0 : -1; } +int geodesic_initialize_past_ray(const MetricSlab *slab, + const ObserverState *o, const double n[3], + State *s) { + MetricData m; + if (spacetime_slab_eval(slab, o->coordinate_time, o->coordinate_position, &m)) + return -1; + return geodesic_initialize_past_ray_metric(&m, o, n, s); +} + static int escaped_direction(const MetricSlab *slab, double t, const State *s, double n[3]) { MetricData m; @@ -151,6 +160,82 @@ static int escaped_direction(const MetricSlab *slab, double t, return 0; } +/* Three-state lifecycle selection. Only a backend that declares no ends at + * all may use the legacy region test; a declared but inconsistent or + * unsupported end is an explicit protocol error, never a silent fallback. */ +typedef enum { + ASYM_LIFECYCLE_NONE, + ASYM_LIFECYCLE_READY, + ASYM_LIFECYCLE_PROTOCOL_ERROR +} AsymLifecycleMode; + +static AsymLifecycleMode asym_lifecycle_mode(const SpacetimeSource *source) { + const size_t count = spacetime_asymptotic_end_count(source); + if (count == 0) + return ASYM_LIFECYCLE_NONE; + for (size_t i = 0; i < count; ++i) { + SpacetimeAsymptoticEnd end; + if (spacetime_asymptotic_end(source, i, &end)) + return ASYM_LIFECYCLE_PROTOCOL_ERROR; + if (end.exterior_kind != ASYMPTOTIC_EXTERIOR_MINKOWSKI && + end.exterior_kind != ASYMPTOTIC_EXTERIOR_SCHWARZSCHILD_MONOPOLE) + return ASYM_LIFECYCLE_PROTOCOL_ERROR; + } + return ASYM_LIFECYCLE_READY; +} + +/* Bounded root localization of the first inside->outside worldtube crossing + * within one accepted step. Re-integrates from `before` with fractional step + * sizes; `after` lands on the outside end of the bracket. Returns + * ASYMPTOTIC_OK, ASYMPTOTIC_TIME_RANGE_EXHAUSTED (a midpoint fell into a + * history hole), or ASYMPTOTIC_INVALID. */ +static AsymptoticStatus localize_worldtube_crossing(const MetricSlab *slab, + SpacetimeEndId end_id, + const State *before, + double h, State *after) { + double f_lo = 0.0, f_hi = 1.0; + for (int iteration = 0; iteration < 64; ++iteration) { + const double f = 0.5 * (f_lo + f_hi); + State mid = *before; + if (rk4(slab, before->coordinate_time, h * f, &mid)) + return ASYMPTOTIC_INVALID; + mid.coordinate_time = before->coordinate_time + h * f; + double value; + const int status = asymptotic_worldtube_value( + slab->source, end_id, mid.coordinate_time, mid.x, &value); + if (status != ASYMPTOTIC_OK) + return (AsymptoticStatus)status; + if (value >= 0.0) + f_hi = f; + else + f_lo = f; + } + *after = *before; + if (rk4(slab, before->coordinate_time, h * f_hi, after)) + return ASYMPTOTIC_INVALID; + after->coordinate_time = before->coordinate_time + h * f_hi; + after->steps = before->steps + 1; + return ASYMPTOTIC_OK; +} + +static GeodesicAdvanceResult legacy_escape_or_capture(const MetricSlab *slab, + const State *s, + SpacetimeRayStatus status, + RayEndpoint *out) { + out->status = + status == SPACETIME_RAY_ESCAPED ? RAY_ENDPOINT_ESCAPED + : RAY_ENDPOINT_CAPTURED; + if (out->status == RAY_ENDPOINT_ESCAPED) { + if (escaped_direction(slab, s->coordinate_time, s, out->n_infinity) == 0) + out->frequency_ratio = exp(-s->log_alpha_p0); + else + out->status = RAY_ENDPOINT_INTEGRATION_FAILURE; + } + return out->status == RAY_ENDPOINT_INTEGRATION_FAILURE + ? GEODESIC_ADVANCE_FAILED + : GEODESIC_ADVANCE_TERMINATED; +} + GeodesicAdvanceResult geodesic_advance_past_ray( const MetricSlab *slab, State *s, double slab_left_time, const GeodesicTraceConfig *config, RayEndpoint *out) { @@ -158,36 +243,105 @@ GeodesicAdvanceResult geodesic_advance_past_ray( !config->max_steps || !isfinite(slab_left_time) || slab_left_time > s->coordinate_time) return GEODESIC_ADVANCE_FAILED; + const AsymLifecycleMode mode = asym_lifecycle_mode(slab->source); + if (mode == ASYM_LIFECYCLE_PROTOCOL_ERROR) { + out->status = RAY_ENDPOINT_INVALID; + return GEODESIC_ADVANCE_FAILED; + } + const int directed = mode == ASYM_LIFECYCLE_READY; + const size_t end_count = + directed ? spacetime_asymptotic_end_count(slab->source) : 0; while (s->coordinate_time > slab_left_time) { if (config->capture_log_alpha_p0 > 0.0 && s->log_alpha_p0 >= config->capture_log_alpha_p0) { out->status = RAY_ENDPOINT_CAPTURED; return GEODESIC_ADVANCE_TERMINATED; } - SpacetimeRayStatus status = + const SpacetimeRayStatus status = spacetime_slab_classify(slab, s->coordinate_time, s->x); - if (status != SPACETIME_RAY_ACTIVE) { - out->status = status == SPACETIME_RAY_ESCAPED ? RAY_ENDPOINT_ESCAPED - : RAY_ENDPOINT_CAPTURED; - if (out->status == RAY_ENDPOINT_ESCAPED && - escaped_direction(slab, s->coordinate_time, s, out->n_infinity) == 0) - out->frequency_ratio = exp(-s->log_alpha_p0); - else if (out->status == RAY_ENDPOINT_ESCAPED) - out->status = RAY_ENDPOINT_INTEGRATION_FAILURE; - return out->status == RAY_ENDPOINT_INTEGRATION_FAILURE - ? GEODESIC_ADVANCE_FAILED - : GEODESIC_ADVANCE_TERMINATED; + if (status == SPACETIME_RAY_CAPTURED) { + out->status = RAY_ENDPOINT_CAPTURED; + return GEODESIC_ADVANCE_TERMINATED; } + if (!directed && status != SPACETIME_RAY_ACTIVE) + return legacy_escape_or_capture(slab, s, status, out); if (s->steps >= config->max_steps) { out->status = RAY_ENDPOINT_MAX_STEPS; return GEODESIC_ADVANCE_TERMINATED; } const double h = -fmin(config->coordinate_time_step, s->coordinate_time - slab_left_time); + const State before = *s; if (rk4(slab, s->coordinate_time, h, s)) return GEODESIC_ADVANCE_FAILED; s->coordinate_time += h; ++s->steps; + if (!directed) + continue; + + if (spacetime_slab_classify(slab, s->coordinate_time, s->x) == + SPACETIME_RAY_CAPTURED) { + out->status = RAY_ENDPOINT_CAPTURED; + return GEODESIC_ADVANCE_TERMINATED; + } + for (size_t i = 0; i < end_count; ++i) { + SpacetimeAsymptoticEnd end; + if (spacetime_asymptotic_end(slab->source, i, &end)) { + out->status = RAY_ENDPOINT_INVALID; + return GEODESIC_ADVANCE_FAILED; + } + double f_before, f_after; + const int before_status = asymptotic_worldtube_value( + slab->source, end.end_id, before.coordinate_time, before.x, + &f_before); + const int after_status = asymptotic_worldtube_value( + slab->source, end.end_id, s->coordinate_time, s->x, &f_after); + if (before_status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED || + after_status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) { + /* A first-class terminal reason, matching pre-route exhaustion: + * preserve the end id and install it as terminated provenance. */ + out->status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED; + out->end_id = end.end_id; + return GEODESIC_ADVANCE_TERMINATED; + } + if (before_status != ASYMPTOTIC_OK || after_status != ASYMPTOTIC_OK) { + /* The backend cannot describe its own worldtube; this is an explicit + * failure, not a physical escape. */ + out->status = RAY_ENDPOINT_INVALID; + return GEODESIC_ADVANCE_FAILED; + } + /* Strict inside->outside: the step must end strictly outside, so a + * single touch at F == 0 (a tangent) is not accepted as a crossing. + * A crossing whose root lands exactly on a step boundary is picked up + * on the following step as f_before == 0, f_after > 0. */ + if (f_before > 0.0 || f_after <= 0.0) + continue; + State crossing; + const AsymptoticStatus localized = localize_worldtube_crossing( + slab, end.end_id, &before, h, &crossing); + if (localized == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) { + out->status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED; + out->end_id = end.end_id; + return GEODESIC_ADVANCE_TERMINATED; + } + if (localized != ASYMPTOTIC_OK) { + out->status = RAY_ENDPOINT_INVALID; + return GEODESIC_ADVANCE_FAILED; + } + const AsymptoticStatus transfer = asymptotic_finish_escape( + slab->source, end.end_id, crossing.coordinate_time, crossing.x, + crossing.Pi, crossing.log_alpha_p0, out); + if (transfer == ASYMPTOTIC_OK) + return GEODESIC_ADVANCE_TERMINATED; + out->status = transfer == ASYMPTOTIC_TIME_RANGE_EXHAUSTED + ? RAY_ENDPOINT_TIME_RANGE_EXHAUSTED + : RAY_ENDPOINT_INVALID; + if (transfer == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) { + out->end_id = end.end_id; + return GEODESIC_ADVANCE_TERMINATED; + } + return GEODESIC_ADVANCE_FAILED; + } } return GEODESIC_ADVANCE_ACTIVE; } @@ -198,22 +352,57 @@ RayEndpoint geodesic_trace_past(const SpacetimeSource *source, const GeodesicTraceConfig *config) { RayEndpoint out = {.frequency_ratio = 0, .magnification = 1, + .end_id = SPACETIME_END_NONE, .status = RAY_ENDPOINT_INTEGRATION_FAILURE}; - State state; - MetricSlab *slab = NULL; if (!source || !observer || !config || config->coordinate_time_step <= 0 || !config->max_steps || fabs(dot(n, n) - 1) > 1e-10) return out; - if (spacetime_load_slab(source, observer->coordinate_time, - observer->coordinate_time - - config->coordinate_time_step * config->max_steps - 1.0, - &slab) || - geodesic_initialize_past_ray(slab, observer, n, &state)) { - spacetime_free_slab(slab); + AsymptoticRoute route; + const AsymptoticStatus route_status = + asymptotic_route_camera(source, observer, n, &route); + if (route_status == ASYMPTOTIC_UNSUPPORTED) { + out.status = RAY_ENDPOINT_INVALID; return out; } + if (route_status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) { + out.status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED; + out.end_id = route.end_id; + return out; + } + if (route_status != ASYMPTOTIC_OK) { + out.status = RAY_ENDPOINT_INVALID; + return out; + } + if (route.kind == ASYMPTOTIC_ROUTE_ESCAPED) { + for (int i = 0; i < 3; ++i) + out.n_infinity[i] = route.n_infinity[i]; + out.frequency_ratio = route.frequency_ratio; + out.end_id = route.end_id; + out.status = RAY_ENDPOINT_ESCAPED; + return out; + } + if (route.kind == ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED) { + out.status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED; + out.end_id = route.end_id; + return out; + } + if (route.kind != ASYMPTOTIC_ROUTE_INSIDE && + route.kind != ASYMPTOTIC_ROUTE_ENTRY) { + out.status = RAY_ENDPOINT_INVALID; + return out; + } + State state = {.coordinate_time = route.activate_t, + .x = {route.x[0], route.x[1], route.x[2]}, + .Pi = {route.Pi[0], route.Pi[1], route.Pi[2]}, + .log_alpha_p0 = route.log_alpha_p0, + .steps = 0}; const double last_time = - observer->coordinate_time - config->coordinate_time_step * config->max_steps; + route.activate_t - config->coordinate_time_step * config->max_steps; + MetricSlab *slab = NULL; + if (spacetime_load_slab(source, route.activate_t, last_time - 1.0, &slab)) { + out.status = RAY_ENDPOINT_INVALID; + return out; + } if (geodesic_advance_past_ray(slab, &state, last_time, config, &out) == GEODESIC_ADVANCE_ACTIVE) out.status = RAY_ENDPOINT_MAX_STEPS; diff --git a/src/geodesic.h b/src/geodesic.h index c8da3bd..7d8bfb9 100644 --- a/src/geodesic.h +++ b/src/geodesic.h @@ -8,13 +8,18 @@ typedef enum { RAY_ENDPOINT_ESCAPED, RAY_ENDPOINT_CAPTURED, RAY_ENDPOINT_MAX_STEPS, - RAY_ENDPOINT_INTEGRATION_FAILURE + RAY_ENDPOINT_INTEGRATION_FAILURE, + RAY_ENDPOINT_TIME_RANGE_EXHAUSTED, + RAY_ENDPOINT_INVALID } RayEndpointStatus; typedef struct { double n_infinity[3]; double frequency_ratio; /* E_camera / E_infinity */ double magnification; /* Filled by the future local inverse lens map. */ + /* Meaningful for RAY_ENDPOINT_ESCAPED and for + * RAY_ENDPOINT_TIME_RANGE_EXHAUSTED; SPACETIME_END_NONE otherwise. */ + SpacetimeEndId end_id; RayEndpointStatus status; } RayEndpoint; @@ -52,6 +57,12 @@ int geodesic_initialize_past_ray(const MetricSlab *slab, const ObserverState *observer, const double camera_direction[3], GeodesicRayState *state); +/* Metric-based core of the initialization above; used by the asymptotic + * pre-route, which evaluates the metric at the camera event directly. */ +int geodesic_initialize_past_ray_metric(const MetricData *metric, + const ObserverState *observer, + const double camera_direction[3], + GeodesicRayState *state); GeodesicAdvanceResult geodesic_advance_past_ray( const MetricSlab *slab, GeodesicRayState *state, double slab_left_time, const GeodesicTraceConfig *config, diff --git a/src/main.c b/src/main.c index 8b56860..4cdae04 100644 --- a/src/main.c +++ b/src/main.c @@ -870,8 +870,12 @@ static double alcubierre_step_budget(const Settings *s) { static GeodesicTraceConfig trace_config(const Settings *s) { #ifdef SPACETIME_SCHWARZSCHILD (void)s; + /* The directed worldtube crossing makes an escaping ray traverse the + * interior as a round trip from the entry sphere (in, turn, back out), + * rather than the old one-way stop at the first outside sample. The step + * budget must cover roughly twice the escape sphere plus margin. */ return (GeodesicTraceConfig){.coordinate_time_step = 0.1, - .max_steps = 4096, + .max_steps = 65536, .capture_log_alpha_p0 = 8.0}; #elif defined(SPACETIME_ALCUBIERRE) const double step = alcubierre_time_step(s); @@ -949,11 +953,9 @@ static int build_observer(const Settings *s, const SpacetimeSource *spacetime, fputs("Camera position is inside the backend capture cutoff or invalid.\n", stderr); return -1; } - if (camera_status == SPACETIME_RAY_ESCAPED) { - fputs("Camera position is outside this backend's finite escape radius; " - "move the camera inward or enlarge the spacetime domain.\n", stderr); - return -1; - } + /* A camera outside the escape sphere is supported by the asymptotic + * exterior module for every declared end kind; unsupported exteriors are + * reported through the ray endpoints instead. */ MetricData metric; if (spacetime_eval(spacetime, camera.coordinate_time, camera.position, &metric)) { fputs("Could not evaluate metric at the camera event.\n", stderr); @@ -1133,6 +1135,7 @@ static int trace_movie_generation(Movie *movie, const Settings *s, return -1; } } + ray_pool_preroute(&rays, spacetime); double slab_hi = movie->frames[movie->frame_count - 1].coordinate_time; size_t slab_id = 0; while (ray_pool_has_live(&rays)) { diff --git a/src/ray.c b/src/ray.c index 8f167c2..68902ab 100644 --- a/src/ray.c +++ b/src/ray.c @@ -1,5 +1,8 @@ #include "ray.h" +#include "asymptotic.h" + +#include #include #include @@ -11,8 +14,9 @@ int ray_pool_init(RayPool *p, size_t capacity) { if (!(RAY_ALLOC(t) && RAY_ALLOC(x0) && RAY_ALLOC(x1) && RAY_ALLOC(x2) && RAY_ALLOC(p0) && RAY_ALLOC(p1) && RAY_ALLOC(p2) && RAY_ALLOC(observer) && RAY_ALLOC(direction0) && RAY_ALLOC(direction1) && RAY_ALLOC(direction2) && - RAY_ALLOC(log_alpha_p0) && RAY_ALLOC(steps) && RAY_ALLOC(frame_id) && - RAY_ALLOC(vertex_id) && RAY_ALLOC(status) && RAY_ALLOC(endpoint))) { + RAY_ALLOC(log_alpha_p0) && RAY_ALLOC(activate_t) && RAY_ALLOC(steps) && + RAY_ALLOC(frame_id) && RAY_ALLOC(vertex_id) && RAY_ALLOC(status) && + RAY_ALLOC(endpoint))) { ray_pool_destroy(p); return -1; } @@ -29,6 +33,7 @@ int ray_pool_append(RayPool *p, const ObserverState *observer, if (observer == NULL || direction == NULL) return -1; p->t[i] = observer->coordinate_time; + p->activate_t[i] = observer->coordinate_time; p->observer[i] = observer; p->direction0[i] = direction[0]; p->direction1[i] = direction[1]; @@ -37,28 +42,75 @@ int ray_pool_append(RayPool *p, const ObserverState *observer, p->vertex_id[i] = vertex_id; p->status[i] = RAY_POOL_PENDING; p->endpoint[i] = (RayEndpoint){.magnification = 1.0, + .end_id = SPACETIME_END_NONE, .status = RAY_ENDPOINT_INTEGRATION_FAILURE}; ++p->count; return 0; } -void ray_pool_activate_in_time_range(RayPool *p, const MetricSlab *slab) { +void ray_pool_preroute(RayPool *p, const SpacetimeSource *source) { + if (p == NULL || source == NULL) + return; +#pragma omp parallel for schedule(static) for (size_t i = 0; i < p->count; ++i) { - if (p->status[i] != RAY_POOL_PENDING || p->t[i] > slab->t_hi || - p->t[i] <= slab->t_lo) + if (p->status[i] != RAY_POOL_PENDING) continue; - GeodesicRayState state; - if (geodesic_initialize_past_ray( - slab, p->observer[i], - (double[]){p->direction0[i], p->direction1[i], p->direction2[i]}, - &state)) { + AsymptoticRoute route; + const AsymptoticStatus status = asymptotic_route_camera( + source, p->observer[i], + (double[]){p->direction0[i], p->direction1[i], p->direction2[i]}, + &route); + if (status == ASYMPTOTIC_UNSUPPORTED || status == ASYMPTOTIC_INVALID) { + p->endpoint[i].status = RAY_ENDPOINT_INVALID; + p->endpoint[i].end_id = route.end_id; p->status[i] = RAY_POOL_FAILED; continue; } - p->x0[i] = state.x[0]; p->x1[i] = state.x[1]; p->x2[i] = state.x[2]; - p->p0[i] = state.Pi[0]; p->p1[i] = state.Pi[1]; p->p2[i] = state.Pi[2]; - p->log_alpha_p0[i] = state.log_alpha_p0; - p->steps[i] = state.steps; + if (status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) { + p->endpoint[i].status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED; + p->endpoint[i].end_id = route.end_id; + p->status[i] = RAY_POOL_TERMINATED; + continue; + } + if (route.kind == ASYMPTOTIC_ROUTE_ESCAPED) { + for (int axis = 0; axis < 3; ++axis) + p->endpoint[i].n_infinity[axis] = route.n_infinity[axis]; + p->endpoint[i].frequency_ratio = route.frequency_ratio; + p->endpoint[i].end_id = route.end_id; + p->endpoint[i].status = RAY_ENDPOINT_ESCAPED; + p->status[i] = RAY_POOL_TERMINATED; + continue; + } + if (route.kind == ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED) { + p->endpoint[i].status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED; + p->endpoint[i].end_id = route.end_id; + p->status[i] = RAY_POOL_TERMINATED; + continue; + } + if (route.kind != ASYMPTOTIC_ROUTE_INSIDE && + route.kind != ASYMPTOTIC_ROUTE_ENTRY) { + p->endpoint[i].status = RAY_ENDPOINT_INVALID; + p->status[i] = RAY_POOL_FAILED; + continue; + } + p->activate_t[i] = route.activate_t; + p->x0[i] = route.x[0]; + p->x1[i] = route.x[1]; + p->x2[i] = route.x[2]; + p->p0[i] = route.Pi[0]; + p->p1[i] = route.Pi[1]; + p->p2[i] = route.Pi[2]; + p->log_alpha_p0[i] = route.log_alpha_p0; + } +} + +void ray_pool_activate_in_time_range(RayPool *p, const MetricSlab *slab) { + for (size_t i = 0; i < p->count; ++i) { + if (p->status[i] != RAY_POOL_PENDING || p->activate_t[i] > slab->t_hi || + p->activate_t[i] <= slab->t_lo) + continue; + p->t[i] = p->activate_t[i]; + p->steps[i] = 0; p->status[i] = RAY_POOL_ACTIVE; } } @@ -107,7 +159,8 @@ void ray_pool_destroy(RayPool *p) { free(p->t); free(p->x0); free(p->x1); free(p->x2); free(p->observer); free(p->direction0); free(p->direction1); free(p->direction2); free(p->p0); free(p->p1); free(p->p2); free(p->log_alpha_p0); - free(p->steps); free(p->frame_id); free(p->vertex_id); free(p->status); + free(p->activate_t); free(p->steps); free(p->frame_id); free(p->vertex_id); + free(p->status); free(p->endpoint); *p = (RayPool){0}; } diff --git a/src/ray.h b/src/ray.h index baf780a..eb38a70 100644 --- a/src/ray.h +++ b/src/ray.h @@ -15,6 +15,10 @@ typedef enum { typedef struct { double *t, *x0, *x1, *x2, *p0, *p1, *p2, *log_alpha_p0; + /* Coordinate time at which the pre-routed interior state becomes valid. + * For a camera inside a worldtube this equals the camera time; for an + * exterior hit it is the earlier entry time. */ + double *activate_t; const ObserverState **observer; double *direction0, *direction1, *direction2; unsigned int *steps; @@ -28,6 +32,8 @@ int ray_pool_init(RayPool *pool, size_t capacity); int ray_pool_append(RayPool *pool, const ObserverState *observer, const double direction[3], size_t frame_id, size_t vertex_id); +/* Pre-route every still-PENDING ray once, before the slab sweep. */ +void ray_pool_preroute(RayPool *pool, const SpacetimeSource *source); void ray_pool_activate_in_time_range(RayPool *pool, const MetricSlab *slab); void ray_pool_advance_active(RayPool *pool, const MetricSlab *slab, const GeodesicTraceConfig *config); diff --git a/src/spacetime.h b/src/spacetime.h index 96def36..c5e0228 100644 --- a/src/spacetime.h +++ b/src/spacetime.h @@ -1,6 +1,9 @@ #ifndef SPACETIME_H #define SPACETIME_H +#include +#include + typedef struct { double alpha; double beta[3]; @@ -17,6 +20,42 @@ typedef enum { SPACETIME_RAY_CAPTURED } SpacetimeRayStatus; +/* Stable identifier for one asymptotic end (infinity) of a backend. Backends + * may describe more than one; the current analytic backends expose one. */ +typedef uint32_t SpacetimeEndId; +#define SPACETIME_END_NONE ((SpacetimeEndId)0xffffffffu) + +typedef enum { + ASYMPTOTIC_EXTERIOR_MINKOWSKI, + ASYMPTOTIC_EXTERIOR_SCHWARZSCHILD_MONOPOLE +} AsymptoticExteriorKind; + +/* Declared asymptotic end. `frame_origin` and the columns of `frame_axes` + * express the asymptotic reference frame in backend coordinates; spatial + * `n_infinity` values use the same coordinate axes as the observer tetrad. */ +typedef struct { + SpacetimeEndId end_id; + AsymptoticExteriorKind exterior_kind; + double mass; + double frame_origin[3]; + double frame_axes[3][3]; +} SpacetimeAsymptoticEnd; + +/* Escape worldtube sample at one coordinate time. A zero `radius_rate` and a + * time-independent `velocity` describe the fixed/constant-velocity cases used + * in this phase. `valid == 0` means the backend cannot describe the worldtube + * at this time (history exhausted); callers must not treat that as a miss. */ +typedef struct { + double center[3]; + double velocity[3]; + double radius; + double radius_rate; + /* Nonzero when `velocity` and `radius_rate` are exact throughout the + * current motion segment, so the first entry has a closed quadratic form. */ + int velocity_constant; + int valid; +} SpacetimeEscapeWorldtubeSample; + typedef struct SpacetimeSource SpacetimeSource; typedef struct MetricSlab MetricSlab; @@ -38,6 +77,20 @@ typedef struct { MetricData *metric); SpacetimeRayStatus (*classify_slab)(const MetricSlab *slab, double t, const double x[3]); + /* Declared asymptotic ends and their moving escape worldtubes. Backends + * without an escape sphere may leave these NULL. */ + size_t (*asymptotic_end_count)(const SpacetimeSource *source); + int (*asymptotic_end)(const SpacetimeSource *source, size_t index, + SpacetimeAsymptoticEnd *out); + int (*escape_worldtube_sample)(const SpacetimeSource *source, + SpacetimeEndId end_id, double t, + SpacetimeEscapeWorldtubeSample *out); + /* Coordinate time of the next motion-segment boundary reached while + * integrating backward in time, i.e. the largest boundary strictly less + * than `t`. Return NAN when the worldtube description has a single open + * segment. */ + double (*escape_worldtube_next_segment)(const SpacetimeSource *source, + SpacetimeEndId end_id, double t); /* Analytic backends have negligible per-ray metric state. A numerical * backend must opt in once its metric slabs and evaluator workspaces need * to reserve memory alongside the private HDR render buffers. */ @@ -76,6 +129,19 @@ int spacetime_slab_eval(const MetricSlab *slab, double t, const double x[3], MetricData *metric); SpacetimeRayStatus spacetime_slab_classify(const MetricSlab *slab, double t, const double x[3]); +size_t spacetime_asymptotic_end_count(const SpacetimeSource *source); +int spacetime_asymptotic_end(const SpacetimeSource *source, size_t index, + SpacetimeAsymptoticEnd *out); +int spacetime_escape_worldtube_sample(const SpacetimeSource *source, + SpacetimeEndId end_id, double t, + SpacetimeEscapeWorldtubeSample *out); +double spacetime_escape_worldtube_next_segment(const SpacetimeSource *source, + SpacetimeEndId end_id, double t); +/* Common structural validation that every successful constructor must pass + * before returning. A source that passes is a promise that it can safely + * enter ray tracing; backend-specific history/segment validation stays in the + * backend constructor. On failure the constructor must destroy its context. */ +int spacetime_source_finalize(SpacetimeSource *source); int spacetime_limits_render_workers_by_memory(const SpacetimeSource *source); #endif diff --git a/src/spacetime_alcubierre.c b/src/spacetime_alcubierre.c index 8e739c0..0d7843a 100644 --- a/src/spacetime_alcubierre.c +++ b/src/spacetime_alcubierre.c @@ -127,9 +127,48 @@ static void alcubierre_destroy(SpacetimeSource *source) { source->ops = NULL; } +static size_t alcubierre_asymptotic_end_count(const SpacetimeSource *source) { + (void)source; + return 1; +} + +static int alcubierre_asymptotic_end(const SpacetimeSource *source, + size_t index, + SpacetimeAsymptoticEnd *out) { + (void)source; + if (index != 0) + return -1; + *out = (SpacetimeAsymptoticEnd){ + .end_id = 0, + .exterior_kind = ASYMPTOTIC_EXTERIOR_MINKOWSKI, + .mass = 0.0, + .frame_origin = {0.0, 0.0, 0.0}, + .frame_axes = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}}; + return 0; +} + +static int alcubierre_escape_worldtube_sample( + const SpacetimeSource *source, SpacetimeEndId end_id, double t, + SpacetimeEscapeWorldtubeSample *out) { + const AlcubierreContext *context = source->context; + if (end_id != 0) + return -1; + *out = (SpacetimeEscapeWorldtubeSample){ + .center = {context->vs * t, 0.0, 0.0}, + .velocity = {context->vs, 0.0, 0.0}, + .radius = context->escape_radius, + .radius_rate = 0.0, + .velocity_constant = 1, + .valid = 1}; + return 0; +} + static const SpacetimeOps alcubierre_ops = { .eval = alcubierre_eval, .classify = alcubierre_classify, + .asymptotic_end_count = alcubierre_asymptotic_end_count, + .asymptotic_end = alcubierre_asymptotic_end, + .escape_worldtube_sample = alcubierre_escape_worldtube_sample, .destroy = alcubierre_destroy, }; @@ -156,6 +195,10 @@ int spacetime_create_alcubierre(SpacetimeSource *source, double vs, context->escape_radius = escape_radius; source->ops = &alcubierre_ops; source->context = context; + if (spacetime_source_finalize(source)) { + alcubierre_destroy(source); + return -1; + } return 0; } diff --git a/src/spacetime_common.c b/src/spacetime_common.c index fc606fd..364e38a 100644 --- a/src/spacetime_common.c +++ b/src/spacetime_common.c @@ -1,5 +1,6 @@ #include "spacetime.h" +#include #include #include @@ -64,6 +65,79 @@ SpacetimeRayStatus spacetime_slab_classify(const MetricSlab *slab, double t, return spacetime_classify(slab->source, t, x); } +size_t spacetime_asymptotic_end_count(const SpacetimeSource *source) { + return source == NULL || source->ops == NULL || + source->ops->asymptotic_end_count == NULL + ? 0 + : source->ops->asymptotic_end_count(source); +} + +int spacetime_asymptotic_end(const SpacetimeSource *source, size_t index, + SpacetimeAsymptoticEnd *out) { + return source == NULL || source->ops == NULL || out == NULL || + source->ops->asymptotic_end == NULL + ? -1 + : source->ops->asymptotic_end(source, index, out); +} + +int spacetime_escape_worldtube_sample(const SpacetimeSource *source, + SpacetimeEndId end_id, double t, + SpacetimeEscapeWorldtubeSample *out) { + if (source == NULL || source->ops == NULL || out == NULL || + source->ops->escape_worldtube_sample == NULL) + return -1; + return source->ops->escape_worldtube_sample(source, end_id, t, out); +} + +double spacetime_escape_worldtube_next_segment(const SpacetimeSource *source, + SpacetimeEndId end_id, + double t) { + if (source == NULL || source->ops == NULL || + source->ops->escape_worldtube_next_segment == NULL) + return NAN; + return source->ops->escape_worldtube_next_segment(source, end_id, t); +} + +int spacetime_source_finalize(SpacetimeSource *source) { + if (source == NULL || source->ops == NULL || source->context == NULL) + return -1; + const SpacetimeOps *ops = source->ops; + if (ops->eval == NULL || ops->classify == NULL || ops->destroy == NULL) + return -1; + const size_t count = spacetime_asymptotic_end_count(source); + if (count == 0) + return 0; /* legacy backend without asymptotic ends */ + if (ops->asymptotic_end == NULL || ops->escape_worldtube_sample == NULL) + return -1; + if (count > 64) + return -1; + SpacetimeEndId ids[64]; + for (size_t i = 0; i < count; ++i) { + SpacetimeAsymptoticEnd end; + if (spacetime_asymptotic_end(source, i, &end)) + return -1; + if (end.end_id == SPACETIME_END_NONE) + return -1; + for (size_t j = 0; j < i; ++j) + if (ids[j] == end.end_id) + return -1; + ids[i] = end.end_id; + if (end.exterior_kind != ASYMPTOTIC_EXTERIOR_MINKOWSKI && + end.exterior_kind != ASYMPTOTIC_EXTERIOR_SCHWARZSCHILD_MONOPOLE) + return -1; + if (!isfinite(end.mass) || end.mass < 0.0) + return -1; + for (int k = 0; k < 3; ++k) { + if (!isfinite(end.frame_origin[k])) + return -1; + for (int l = 0; l < 3; ++l) + if (!isfinite(end.frame_axes[k][l])) + return -1; + } + } + return 0; +} + int spacetime_limits_render_workers_by_memory(const SpacetimeSource *source) { return source != NULL && source->ops != NULL && source->ops->limit_render_workers_by_memory; diff --git a/src/spacetime_minkowski.c b/src/spacetime_minkowski.c index 2058754..3b6007b 100644 --- a/src/spacetime_minkowski.c +++ b/src/spacetime_minkowski.c @@ -33,9 +33,47 @@ static void minkowski_destroy(SpacetimeSource *source) { source->ops = NULL; } +static size_t minkowski_asymptotic_end_count(const SpacetimeSource *source) { + (void)source; + return 1; +} + +static int minkowski_asymptotic_end(const SpacetimeSource *source, + size_t index, SpacetimeAsymptoticEnd *out) { + (void)source; + if (index != 0) + return -1; + *out = (SpacetimeAsymptoticEnd){ + .end_id = 0, + .exterior_kind = ASYMPTOTIC_EXTERIOR_MINKOWSKI, + .mass = 0.0, + .frame_origin = {0.0, 0.0, 0.0}, + .frame_axes = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}}; + return 0; +} + +static int minkowski_escape_worldtube_sample( + const SpacetimeSource *source, SpacetimeEndId end_id, double t, + SpacetimeEscapeWorldtubeSample *out) { + const MinkowskiContext *context = source->context; + if (end_id != 0) + return -1; + *out = (SpacetimeEscapeWorldtubeSample){.center = {0.0, 0.0, 0.0}, + .velocity = {0.0, 0.0, 0.0}, + .radius = context->escape_radius, + .radius_rate = 0.0, + .velocity_constant = 1, + .valid = 1}; + (void)t; + return 0; +} + static const SpacetimeOps minkowski_ops = { .eval = minkowski_eval, .classify = minkowski_classify, + .asymptotic_end_count = minkowski_asymptotic_end_count, + .asymptotic_end = minkowski_asymptotic_end, + .escape_worldtube_sample = minkowski_escape_worldtube_sample, .destroy = minkowski_destroy, }; @@ -48,6 +86,10 @@ int spacetime_create_minkowski(SpacetimeSource *source, double escape_radius) { context->escape_radius = escape_radius; source->ops = &minkowski_ops; source->context = context; + if (spacetime_source_finalize(source)) { + minkowski_destroy(source); + return -1; + } return 0; } diff --git a/src/spacetime_schwarzschild.c b/src/spacetime_schwarzschild.c index ffdc0ab..9137c25 100644 --- a/src/spacetime_schwarzschild.c +++ b/src/spacetime_schwarzschild.c @@ -104,9 +104,48 @@ static void schwarzschild_ks_destroy(SpacetimeSource *source) { source->ops = NULL; } +static size_t schwarzschild_ks_asymptotic_end_count( + const SpacetimeSource *source) { + (void)source; + return 1; +} + +static int schwarzschild_ks_asymptotic_end( + const SpacetimeSource *source, size_t index, SpacetimeAsymptoticEnd *out) { + const SchwarzschildKsContext *context = source->context; + if (index != 0) + return -1; + *out = (SpacetimeAsymptoticEnd){ + .end_id = 0, + .exterior_kind = ASYMPTOTIC_EXTERIOR_SCHWARZSCHILD_MONOPOLE, + .mass = context->mass, + .frame_origin = {0.0, 0.0, 0.0}, + .frame_axes = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}}; + return 0; +} + +static int schwarzschild_ks_escape_worldtube_sample( + const SpacetimeSource *source, SpacetimeEndId end_id, double t, + SpacetimeEscapeWorldtubeSample *out) { + const SchwarzschildKsContext *context = source->context; + if (end_id != 0) + return -1; + *out = (SpacetimeEscapeWorldtubeSample){.center = {0.0, 0.0, 0.0}, + .velocity = {0.0, 0.0, 0.0}, + .radius = context->escape_radius, + .radius_rate = 0.0, + .velocity_constant = 1, + .valid = 1}; + (void)t; + return 0; +} + static const SpacetimeOps schwarzschild_ks_ops = { .eval = schwarzschild_ks_eval, .classify = schwarzschild_ks_classify, + .asymptotic_end_count = schwarzschild_ks_asymptotic_end_count, + .asymptotic_end = schwarzschild_ks_asymptotic_end, + .escape_worldtube_sample = schwarzschild_ks_escape_worldtube_sample, .destroy = schwarzschild_ks_destroy, }; @@ -123,6 +162,10 @@ int spacetime_create_schwarzschild_ks(SpacetimeSource *source, double mass, *context = (SchwarzschildKsContext){mass, escape_radius, capture_radius}; source->ops = &schwarzschild_ks_ops; source->context = context; + if (spacetime_source_finalize(source)) { + schwarzschild_ks_destroy(source); + return -1; + } return 0; } diff --git a/tests/data/psf_event_sink_reference/schwarzschild_ra1_dec1_fov60_640x360.png b/tests/data/psf_event_sink_reference/schwarzschild_ra1_dec1_fov60_640x360.png index 83e4bc2..1477018 100644 Binary files a/tests/data/psf_event_sink_reference/schwarzschild_ra1_dec1_fov60_640x360.png and b/tests/data/psf_event_sink_reference/schwarzschild_ra1_dec1_fov60_640x360.png differ diff --git a/tests/data/psf_event_sink_reference/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits b/tests/data/psf_event_sink_reference/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits index 45ab3bd..7d4ec8a 100644 Binary files a/tests/data/psf_event_sink_reference/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits and b/tests/data/psf_event_sink_reference/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits differ diff --git a/tests/test_asymptotic.c b/tests/test_asymptotic.c new file mode 100644 index 0000000..b4b786f --- /dev/null +++ b/tests/test_asymptotic.c @@ -0,0 +1,747 @@ +#include "asymptotic.h" +#include "observer.h" +#include "ray.h" +#include "spacetime.h" + +#include +#include + +static int failures = 0; + +#define CHECK(condition, message) \ + do { \ + if (!(condition)) { \ + fprintf(stderr, "FAIL %s:%d: %s\n", __FILE__, __LINE__, message); \ + ++failures; \ + } \ + } while (0) + +static ObserverState flat_observer(double x, double y, double z) { + ObserverState o = {0}; + o.coordinate_position[0] = x; + o.coordinate_position[1] = y; + o.coordinate_position[2] = z; + o.tetrad[0][0] = 1.0; + o.tetrad[1][1] = 1.0; + o.tetrad[2][2] = 1.0; + o.tetrad[3][3] = 1.0; + return o; +} + +/* Synthetic flat exterior with a Minkowski end whose worldtube center follows + * x_c(t) = vx t + accel t^2 / 2. `constant` selects the closed quadratic path; + * otherwise the generic bracketed driver runs. */ +typedef struct { + double vx; + double accel; + double radius; + double radius_rate; + double valid_t_min; + int constant; + double segment_t; /* Motion-segment boundary for the cross-segment test. */ + int has_segment; + int end_descriptor_fails; /* Protocol-error injection. */ + int unsupported_kind; + double invalid_center, invalid_halfwidth; /* Isolated invalid time window. */ + int schwarzschild_kind; /* Declare a Schwarzschild monopole end. */ + int sample_callback_fails; /* make escape_worldtube_sample return -1 */ + int sample_invalid; /* valid = 0 */ + int sample_nan_radius; + int sample_nonpositive_radius; + int fail_on_sample_call; /* 1-based callback invocation to fail. */ + int sample_call_count; +} SyntheticContext; + +static int synthetic_eval(const SpacetimeSource *source, double t, + const double x[3], MetricData *metric) { + (void)source; + (void)t; + (void)x; + *metric = (MetricData){.alpha = 1.0, + .gamma = {{1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {0.0, 0.0, 1.0}}}; + return 0; +} + +static SpacetimeRayStatus synthetic_classify(const SpacetimeSource *source, + double t, const double x[3]) { + (void)source; + (void)t; + (void)x; + return SPACETIME_RAY_ACTIVE; +} + +static size_t synthetic_end_count(const SpacetimeSource *source) { + (void)source; + return 1; +} + +static int synthetic_end(const SpacetimeSource *source, size_t index, + SpacetimeAsymptoticEnd *out) { + const SyntheticContext *context = source->context; + if (index != 0 || context->end_descriptor_fails) + return -1; + AsymptoticExteriorKind kind = ASYMPTOTIC_EXTERIOR_MINKOWSKI; + double mass = 0.0; + if (context->unsupported_kind) { + kind = (AsymptoticExteriorKind)999; + } else if (context->schwarzschild_kind) { + kind = ASYMPTOTIC_EXTERIOR_SCHWARZSCHILD_MONOPOLE; + mass = 1.0; + } + *out = (SpacetimeAsymptoticEnd){ + .end_id = 0, + .exterior_kind = kind, + .mass = mass, + .frame_origin = {0.0, 0.0, 0.0}, + .frame_axes = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}}; + return 0; +} + +static int synthetic_worldtube(const SpacetimeSource *source, + SpacetimeEndId end_id, double t, + SpacetimeEscapeWorldtubeSample *out) { + SyntheticContext *mutable_context = source->context; + const SyntheticContext *context = mutable_context; + if (end_id != 0) + return -1; + ++mutable_context->sample_call_count; + if (context->fail_on_sample_call > 0 && + mutable_context->sample_call_count == context->fail_on_sample_call) + return -1; + if (context->sample_callback_fails) + return -1; + if (context->sample_invalid) { + *out = (SpacetimeEscapeWorldtubeSample){.valid = 0}; + return 0; + } + if (context->sample_nan_radius) { + *out = (SpacetimeEscapeWorldtubeSample){.radius = NAN, .valid = 1}; + return 0; + } + if (context->sample_nonpositive_radius) { + *out = (SpacetimeEscapeWorldtubeSample){.radius = 0.0, .valid = 1}; + return 0; + } + if (!isfinite(t) || t < context->valid_t_min) { + *out = (SpacetimeEscapeWorldtubeSample){.valid = 0}; + return 0; + } + if (context->invalid_halfwidth > 0.0 && + fabs(t - context->invalid_center) <= context->invalid_halfwidth) { + *out = (SpacetimeEscapeWorldtubeSample){.valid = 0}; + return 0; + } + if (context->has_segment && t < context->segment_t) { + /* Second segment: center moves toward +x as t decreases. */ + *out = (SpacetimeEscapeWorldtubeSample){ + .center = {context->segment_t - t, 0.0, 0.0}, + .velocity = {-1.0, 0.0, 0.0}, + .radius = context->radius, + .radius_rate = context->radius_rate, + .velocity_constant = context->constant, + .valid = 1}; + return 0; + } + *out = (SpacetimeEscapeWorldtubeSample){ + .center = {context->vx * t + 0.5 * context->accel * t * t, 0.0, 0.0}, + .velocity = {context->vx + context->accel * t, 0.0, 0.0}, + .radius = context->radius, + .radius_rate = context->radius_rate, + .velocity_constant = context->constant, + .valid = 1}; + return 0; +} + +static double synthetic_next_segment(const SpacetimeSource *source, + SpacetimeEndId end_id, double t) { + const SyntheticContext *context = source->context; + (void)end_id; + if (context->has_segment && t > context->segment_t) + return context->segment_t; + return NAN; +} + +static void synthetic_destroy(SpacetimeSource *source) { + /* The test context lives on the stack, so it is not freed; but match the + * real destroy postcondition. */ + source->context = NULL; + source->ops = NULL; +} + +static const SpacetimeOps synthetic_ops = { + .eval = synthetic_eval, + .classify = synthetic_classify, + .asymptotic_end_count = synthetic_end_count, + .asymptotic_end = synthetic_end, + .escape_worldtube_sample = synthetic_worldtube, + .escape_worldtube_next_segment = synthetic_next_segment, + .destroy = synthetic_destroy, +}; + +static void test_fixed_sphere(void) { + SpacetimeSource source = {0}; + CHECK(spacetime_create_minkowski(&source, 10.0) == 0, "create minkowski"); + AsymptoticRoute route; + + const ObserverState inside = flat_observer(0.0, 0.0, 0.0); + CHECK(asymptotic_route_camera(&source, &inside, (double[]){1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_INSIDE, + "origin camera is inside"); + + const ObserverState outside = flat_observer(50.0, 0.0, 0.0); + CHECK(asymptotic_route_camera(&source, &outside, (double[]){-1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ENTRY, + "outside ray toward sphere enters"); + CHECK(fabs(route.activate_t + 40.0) < 1e-9, "fixed-sphere entry time"); + CHECK(fabs(route.x[0] - 10.0) < 1e-9 && fabs(route.x[1]) < 1e-9 && + fabs(route.x[2]) < 1e-9, + "fixed-sphere entry position"); + + CHECK(asymptotic_route_camera(&source, &outside, (double[]){1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ESCAPED, + "outside ray away misses"); + CHECK(fabs(route.n_infinity[0] - 1.0) < 1e-12 && + fabs(route.n_infinity[1]) < 1e-12, + "miss direction"); + CHECK(fabs(route.frequency_ratio - 1.0) < 1e-12, "flat frequency ratio"); + + const ObserverState tangent = flat_observer(50.0, 10.0, 0.0); + CHECK(asymptotic_route_camera(&source, &tangent, (double[]){-1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ESCAPED, + "tangent ray is not a crossing"); + + const ObserverState near_miss = flat_observer(50.0, 10.0 + 1e-6, 0.0); + CHECK(asymptotic_route_camera(&source, &near_miss, + (double[]){-1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ESCAPED, + "near-tangent outside ray misses"); + const ObserverState near_hit = flat_observer(50.0, 10.0 - 1e-6, 0.0); + CHECK(asymptotic_route_camera(&source, &near_hit, + (double[]){-1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ENTRY, + "near-tangent inside ray enters"); + + RayEndpoint endpoint; + CHECK(asymptotic_finish_escape(&source, 0, 0.0, (double[]){10.0, 0.0, 0.0}, + (double[]){-1.0, 0.0, 0.0}, 0.0, + &endpoint) == ASYMPTOTIC_OK && + endpoint.status == RAY_ENDPOINT_ESCAPED && endpoint.end_id == 0, + "finish outward crossing"); + CHECK(fabs(endpoint.n_infinity[0] - 1.0) < 1e-12 && + fabs(endpoint.frequency_ratio - 1.0) < 1e-12, + "finish direction and frequency"); + + spacetime_destroy(&source); +} + +static void test_large_radius_quadratic(void) { + SpacetimeSource source = {0}; + CHECK(spacetime_create_minkowski(&source, 1.0e12) == 0, + "create huge minkowski sphere"); + AsymptoticRoute route; + const ObserverState hit = flat_observer(2.0e12, 5.0e11, 0.0); + CHECK(asymptotic_route_camera(&source, &hit, (double[]){-1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ENTRY, + "large-radius hit stays quadratic"); + const ObserverState miss = flat_observer(2.0e12, 2.0e12, 0.0); + CHECK(asymptotic_route_camera(&source, &miss, (double[]){-1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ESCAPED, + "large-radius miss stays quadratic"); + spacetime_destroy(&source); + + /* Small entry root a hair outside a large sphere: the cancellation-prone + * case for the naive formula. */ + SpacetimeSource big = {0}; + CHECK(spacetime_create_minkowski(&big, 1.0e9) == 0, "create 1e9 sphere"); + const ObserverState just_outside = flat_observer(1.0e9 + 1e-3, 0.0, 0.0); + CHECK(asymptotic_route_camera(&big, &just_outside, + (double[]){-1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ENTRY, + "just-outside hit"); + const double expected_delta = + just_outside.coordinate_position[0] - 1.0e9; + CHECK(fabs(-route.activate_t - expected_delta) < + 1e-7 + 1e-11 * fabs(expected_delta), + "just-outside entry time within budget"); + double residual; + CHECK(asymptotic_worldtube_value(&big, route.end_id, route.activate_t, + route.x, &residual) == 0 && + fabs(residual) <= 1e-12 * 1.0e9 * 1.0e9, + "just-outside entry on worldtube"); + spacetime_destroy(&big); +} + +static void test_boundary_semantics_minkowski(void) { + SpacetimeSource source = {0}; + CHECK(spacetime_create_minkowski(&source, 10.0) == 0, + "create minkowski"); + AsymptoticRoute route; + const ObserverState on_boundary = flat_observer(10.0, 0.0, 0.0); + CHECK(asymptotic_route_camera(&source, &on_boundary, + (double[]){-1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_INSIDE, + "on-boundary past-inward is inside"); + CHECK(asymptotic_route_camera(&source, &on_boundary, + (double[]){1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ESCAPED, + "on-boundary past-outward escapes"); + CHECK(asymptotic_route_camera(&source, &on_boundary, + (double[]){0.0, 1.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ESCAPED, + "on-boundary tangent escapes"); + spacetime_destroy(&source); +} + +static void test_boundary_semantics_generic(void) { + /* velocity_constant == 0 forces the generic bracketed driver. A finite + * history bounds the outward/tangent searches, which must not be reported + * as entries (they end as TIME_RANGE_EXHAUSTED instead). */ + SyntheticContext context = {.radius = 10.0, + .valid_t_min = -100.0, + .constant = 0}; + SpacetimeSource source = {.ops = &synthetic_ops, .context = &context}; + const ObserverState on_boundary = flat_observer(10.0, 0.0, 0.0); + AsymptoticRoute route; + const AsymptoticStatus inward = asymptotic_route_camera( + &source, &on_boundary, (double[]){-1.0, 0.0, 0.0}, &route); + CHECK(inward == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_INSIDE, + "generic on-boundary inward is inside"); + const AsymptoticStatus outward = asymptotic_route_camera( + &source, &on_boundary, (double[]){1.0, 0.0, 0.0}, &route); + CHECK(outward == ASYMPTOTIC_TIME_RANGE_EXHAUSTED, + "generic on-boundary outward is not an entry"); + const AsymptoticStatus tangent = asymptotic_route_camera( + &source, &on_boundary, (double[]){0.0, 1.0, 0.0}, &route); + CHECK(tangent == ASYMPTOTIC_TIME_RANGE_EXHAUSTED, + "generic on-boundary tangent is not an entry"); +} + +static void test_negative_radius_root_guard(void) { + /* Deliberately bypasses a constructor: every sampled radius is finite and + * positive, but the algebraic root sits where R < 0. The cheap root-level + * guard must reject it instead of fabricating a negative-radius entry. */ + SyntheticContext context = {.radius = 10.0, + .radius_rate = 2.0, + .valid_t_min = -1.0e30, + .constant = 1}; + SpacetimeSource source = {.ops = &synthetic_ops, .context = &context}; + const ObserverState on_boundary = flat_observer(10.0, 0.0, 0.0); + AsymptoticRoute route; + CHECK(asymptotic_route_camera(&source, &on_boundary, + (double[]){1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_INVALID, + "negative-radius algebraic root is rejected"); +} + +static void test_source_finalize(void) { + /* A real constructor already finalizes: finalize is idempotent. */ + SpacetimeSource good = {0}; + CHECK(spacetime_create_minkowski(&good, 10.0) == 0, "create minkowski"); + CHECK(spacetime_source_finalize(&good) == 0, "valid source finalizes"); + spacetime_destroy(&good); + + /* A structurally valid synthetic source must pass, so the failure cases + * below are attributable to their specific defect rather than to the test + * ops themselves. */ + SyntheticContext well_formed = {.radius = 10.0, + .valid_t_min = -1.0e30, + .constant = 1}; + SpacetimeSource valid_source = {.ops = &synthetic_ops, + .context = &well_formed}; + CHECK(spacetime_source_finalize(&valid_source) == 0, + "well-formed synthetic source finalizes"); + + /* Structural protocol errors must be rejected before any ray trace. */ + SyntheticContext bad_kind = {.radius = 10.0, + .valid_t_min = -1.0e30, + .constant = 1, + .unsupported_kind = 1}; + SpacetimeSource source = {.ops = &synthetic_ops, .context = &bad_kind}; + CHECK(spacetime_source_finalize(&source) != 0, + "unsupported exterior kind fails finalize"); + + SyntheticContext bad_desc = {.radius = 10.0, + .valid_t_min = -1.0e30, + .constant = 1, + .end_descriptor_fails = 1}; + source = (SpacetimeSource){.ops = &synthetic_ops, .context = &bad_desc}; + CHECK(spacetime_source_finalize(&source) != 0, + "broken end descriptor fails finalize"); +} + +static void test_motion_segment_domain(void) { + /* Segment 1 (t >= -50) is a static R=10 sphere; its quadratic root lies at + * s = 90, past the segment boundary. Segment 2 (t < -50) moves the center + * with velocity -1, so the true entry is at s = 70. The result must come + * from segment 2. */ + SyntheticContext context = {.radius = 10.0, + .valid_t_min = -1.0e30, + .constant = 1, + .segment_t = -50.0, + .has_segment = 1}; + SpacetimeSource source = {.ops = &synthetic_ops, .context = &context}; + const ObserverState camera = flat_observer(100.0, 0.0, 0.0); + AsymptoticRoute route; + CHECK(asymptotic_route_camera(&source, &camera, + (double[]){-1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ENTRY, + "entry found in the second motion segment"); + CHECK(fabs(route.activate_t + 70.0) < 1e-6, + "second-segment entry, not the stale first-segment root"); +} + +static void test_schwarzschild_sample_failures(void) { + const ObserverState camera = flat_observer(0.0, 0.0, 0.0); + AsymptoticRoute route; + SyntheticContext base = {.radius = 256.0, + .valid_t_min = -1.0e30, + .constant = 1, + .schwarzschild_kind = 1}; + + SyntheticContext callback = base; + callback.sample_callback_fails = 1; + SpacetimeSource s1 = {.ops = &synthetic_ops, .context = &callback}; + CHECK(asymptotic_route_camera(&s1, &camera, (double[]){1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_INVALID, + "schwarzschild callback failure is invalid"); + + SyntheticContext invalid = base; + invalid.sample_invalid = 1; + SpacetimeSource s2 = {.ops = &synthetic_ops, .context = &invalid}; + CHECK(asymptotic_route_camera(&s2, &camera, (double[]){1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_TIME_RANGE_EXHAUSTED, + "schwarzschild valid=0 is exhausted"); + + SyntheticContext nan = base; + nan.sample_nan_radius = 1; + SpacetimeSource s3 = {.ops = &synthetic_ops, .context = &nan}; + CHECK(asymptotic_route_camera(&s3, &camera, (double[]){1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_INVALID, + "schwarzschild NaN radius is invalid"); + + SyntheticContext zero = base; + zero.sample_nonpositive_radius = 1; + SpacetimeSource s4 = {.ops = &synthetic_ops, .context = &zero}; + CHECK(asymptotic_route_camera(&s4, &camera, (double[]){1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_INVALID, + "schwarzschild non-positive radius is invalid"); + + /* Second call fails: the containment sample (#1) succeeds with the camera + * outside, and schwarzschild_route's own sample (#2) is the one that fails. + * This locks the dedicated Schwarzschild sample handling. */ + SyntheticContext second = base; + second.fail_on_sample_call = 2; + SpacetimeSource s5 = {.ops = &synthetic_ops, .context = &second}; + const ObserverState outside = flat_observer(500.0, 0.0, 0.0); + CHECK(asymptotic_route_camera(&s5, &outside, (double[]){1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_INVALID, + "schwarzschild_route second-sample failure is invalid"); +} + +static void test_end_protocol_error(void) { + const ObserverState inside = flat_observer(0.0, 0.0, 0.0); + AsymptoticRoute route; + SyntheticContext bad = {.radius = 20.0, + .valid_t_min = -1.0e30, + .constant = 1, + .end_descriptor_fails = 1}; + SpacetimeSource bad_source = {.ops = &synthetic_ops, .context = &bad}; + CHECK(asymptotic_route_camera(&bad_source, &inside, + (double[]){1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_INVALID, + "bad end descriptor is an explicit protocol error"); + + SyntheticContext unsupported = {.radius = 20.0, + .valid_t_min = -1.0e30, + .constant = 1, + .unsupported_kind = 1}; + SpacetimeSource unsupported_source = {.ops = &synthetic_ops, + .context = &unsupported}; + CHECK(asymptotic_route_camera(&unsupported_source, &inside, + (double[]){1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_UNSUPPORTED, + "unsupported exterior with camera inside is not silently accepted"); + + /* The lifecycle layer, not just the pre-route, must refuse legacy fallback + * whenever ends are declared but broken. */ + MetricSlab *slab = NULL; + CHECK(spacetime_load_slab(&bad_source, 0.0, -10.0, &slab) == 0, + "bad-descriptor slab"); + GeodesicRayState state = {.coordinate_time = 0.0, + .x = {1.0, 0.0, 0.0}, + .Pi = {0.0, 0.0, 0.0}, + .log_alpha_p0 = 0.0, + .steps = 0}; + const GeodesicTraceConfig config = {.coordinate_time_step = 1.0, + .max_steps = 10}; + RayEndpoint endpoint = {.frequency_ratio = 0.0, + .magnification = 1.0, + .end_id = SPACETIME_END_NONE, + .status = RAY_ENDPOINT_INVALID}; + CHECK(geodesic_advance_past_ray(slab, &state, -10.0, &config, &endpoint) == + GEODESIC_ADVANCE_FAILED && + endpoint.status == RAY_ENDPOINT_INVALID, + "advance rejects a declared-but-broken end without legacy"); + spacetime_free_slab(slab); +} + +static void test_interior_crossing_bisection_failure(void) { + /* radius 20.3 makes the exit land strictly between steps: the accepted + * step goes from F < 0 (t = -119.7) to F > 0 (t = -120.7). The invalid + * window sits on the first bisection midpoint (t = -120.2), while both + * accepted-step endpoints stay valid. */ + SyntheticContext context = {.radius = 20.3, + .valid_t_min = -1.0e30, + .constant = 1, + .invalid_center = -120.2, + .invalid_halfwidth = 0.05}; + SpacetimeSource source = {.ops = &synthetic_ops, .context = &context}; + const ObserverState observer = flat_observer(100.0, 0.0, 0.0); + const GeodesicTraceConfig config = {.coordinate_time_step = 1.0, + .max_steps = 2048}; + const RayEndpoint endpoint = geodesic_trace_past( + &source, &observer, (double[]){-1.0, 0.0, 0.0}, &config); + CHECK(endpoint.status == RAY_ENDPOINT_TIME_RANGE_EXHAUSTED && + endpoint.end_id == 0, + "interior crossing bisection propagates history exhaustion"); +} + +static void test_generic_bisection_failure(void) { + /* The isolated invalid window lands on a bisection midpoint while the + * bracket endpoints stay valid, so only the bisection can see it. With the + * strict F < 0 entry test, the bracket is s = 80 (F == 0) to s = 90 + * (F < 0), so the first midpoint is t = -85. */ + SyntheticContext context = {.radius = 20.0, + .valid_t_min = -1.0e30, + .constant = 0, + .invalid_center = -85.0, + .invalid_halfwidth = 1.0}; + SpacetimeSource source = {.ops = &synthetic_ops, .context = &context}; + const ObserverState camera = flat_observer(100.0, 0.0, 0.0); + AsymptoticRoute route; + CHECK(asymptotic_route_camera(&source, &camera, (double[]){-1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_TIME_RANGE_EXHAUSTED, + "generic worldtube bisection propagates sample failure"); +} + +static void test_interior_history_exhaustion(void) { + SyntheticContext context = {.radius = 20.0, + .valid_t_min = -100.0, + .constant = 1}; + SpacetimeSource source = {.ops = &synthetic_ops, .context = &context}; + const ObserverState observer = flat_observer(100.0, 0.0, 0.0); + const GeodesicTraceConfig config = {.coordinate_time_step = 1.0, + .max_steps = 2048}; + const RayEndpoint endpoint = geodesic_trace_past( + &source, &observer, (double[]){-1.0, 0.0, 0.0}, &config); + CHECK(endpoint.status == RAY_ENDPOINT_TIME_RANGE_EXHAUSTED && + endpoint.end_id == 0, + "interior worldtube history exhaustion on a single trace"); + + RayPool pool; + CHECK(ray_pool_init(&pool, 1) == 0, "pool init"); + CHECK(ray_pool_append(&pool, &observer, (double[]){-1.0, 0.0, 0.0}, 0, 0) == + 0, + "append exhaustion ray"); + ray_pool_preroute(&pool, &source); + CHECK(pool.status[0] == RAY_POOL_PENDING, "exhaustion ray pends entry"); + MetricSlab *slab = NULL; + CHECK(spacetime_load_slab(&source, -80.0, -3000.0, &slab) == 0, + "exhaustion slab"); + ray_pool_activate_in_time_range(&pool, slab); + CHECK(pool.status[0] == RAY_POOL_ACTIVE, "exhaustion ray activates"); + ray_pool_advance_active(&pool, slab, &config); + CHECK(pool.endpoint[0].status == RAY_ENDPOINT_TIME_RANGE_EXHAUSTED && + pool.endpoint[0].end_id == 0 && + pool.status[0] == RAY_POOL_TERMINATED, + "interior worldtube history exhaustion on a RayPool"); + spacetime_free_slab(slab); + ray_pool_destroy(&pool); +} + +static void test_round_trip(void) { + SpacetimeSource source = {0}; + CHECK(spacetime_create_minkowski(&source, 10.0) == 0, "create minkowski"); + MetricData metric = {.alpha = 1.0, + .gamma = {{1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {0.0, 0.0, 1.0}}}; + AsymptoticPhotonState canonical; + CHECK(asymptotic_canonical_from_backend(&source, 0, &metric, 0.0, + (double[]){3.0, 4.0, 0.0}, + (double[]){-0.6, 0.8, 0.0}, 0.25, + &canonical) == 0, + "backend to canonical"); + double x[3], Pi[3], log_alpha_p0; + CHECK(asymptotic_backend_from_canonical(&source, &metric, &canonical, x, Pi, + &log_alpha_p0) == 0, + "canonical to backend"); + CHECK(fabs(x[0] - 3.0) < 1e-14 && fabs(x[1] - 4.0) < 1e-14 && + fabs(Pi[0] + 0.6) < 1e-14 && fabs(Pi[1] - 0.8) < 1e-14 && + fabs(log_alpha_p0 - 0.25) < 1e-14, + "round trip matches"); + spacetime_destroy(&source); +} + +static void test_moving_sphere(void) { + SyntheticContext context = {.vx = 0.5, .accel = 0.0, .radius = 25.0, + .valid_t_min = -1.0e30, .constant = 1}; + SpacetimeSource source = {.ops = &synthetic_ops, .context = &context}; + AsymptoticRoute route; + + const ObserverState head_on = flat_observer(100.0, 0.0, 0.0); + CHECK(asymptotic_route_camera(&source, &head_on, (double[]){-1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ENTRY, + "head-on moving-sphere entry"); + CHECK(fabs(route.activate_t + 150.0) < 1e-9, "head-on entry time"); + CHECK(fabs(route.x[0] + 50.0) < 1e-9, "head-on entry position"); + double value; + CHECK(asymptotic_worldtube_value(&source, route.end_id, route.activate_t, + route.x, &value) == 0 && + fabs(value) <= 1e-13 * 25.0 * 25.0, + "head-on entry lies on worldtube"); + + const ObserverState transverse = flat_observer(0.0, 40.0, 0.0); + CHECK(asymptotic_route_camera(&source, &transverse, + (double[]){0.0, -1.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ENTRY, + "transverse moving-sphere entry"); + CHECK(asymptotic_worldtube_value(&source, route.end_id, route.activate_t, + route.x, &value) == 0 && + fabs(value) <= 1e-13 * 25.0 * 25.0, + "transverse entry lies on worldtube"); + + const ObserverState away = flat_observer(100.0, 0.0, 0.0); + CHECK(asymptotic_route_camera(&source, &away, (double[]){1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ESCAPED, + "co-moving ray misses"); +} + +static void test_accelerated_worldtube(void) { + SyntheticContext context = {.vx = 0.0, .accel = 0.02, .radius = 20.0, + .valid_t_min = -1.0e30, .constant = 0}; + SpacetimeSource source = {.ops = &synthetic_ops, .context = &context}; + const ObserverState camera = flat_observer(100.0, 0.0, 0.0); + AsymptoticRoute route; + CHECK(asymptotic_route_camera(&source, &camera, (double[]){-1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ENTRY, + "accelerated worldtube entry"); + double value; + CHECK(asymptotic_worldtube_value(&source, route.end_id, route.activate_t, + route.x, &value) == 0 && + fabs(value) <= 1e-13 * 20.0 * 20.0, + "accelerated entry lies on worldtube"); + CHECK(route.activate_t < -40.0 && route.activate_t > -60.0, + "accelerated entry time in range"); + + context.valid_t_min = -30.0; + CHECK(asymptotic_route_camera(&source, &camera, (double[]){-1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_TIME_RANGE_EXHAUSTED, + "exhausted history is not a miss"); +} + +static void test_accelerated_segment(void) { + SyntheticContext context = {.vx = 0.0, + .accel = 0.0, + .radius = 20.0, + .valid_t_min = -1.0e30, + .constant = 0, + .segment_t = -50.0, + .has_segment = 1}; + SpacetimeSource source = {.ops = &synthetic_ops, .context = &context}; + const ObserverState camera = flat_observer(100.0, 0.0, 0.0); + AsymptoticRoute route; + CHECK(asymptotic_route_camera(&source, &camera, (double[]){-1.0, 0.0, 0.0}, + &route) == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ENTRY, + "cross-segment entry"); + CHECK(route.activate_t < context.segment_t, + "entry lies past the motion-segment boundary"); + CHECK(fabs(route.activate_t + 65.0) < 1e-6, "cross-segment entry time"); + double value; + CHECK(asymptotic_worldtube_value(&source, route.end_id, route.activate_t, + route.x, &value) == 0 && + fabs(value) <= 1e-13 * 20.0 * 20.0, + "cross-segment entry on worldtube"); +} + +static void test_ray_pool_lifecycle(void) { + SpacetimeSource source = {0}; + CHECK(spacetime_create_minkowski(&source, 10.0) == 0, "create minkowski"); + const ObserverState observer = flat_observer(50.0, 0.0, 0.0); + RayPool pool; + CHECK(ray_pool_init(&pool, 2) == 0, "pool init"); + CHECK(ray_pool_append(&pool, &observer, (double[]){-1.0, 0.0, 0.0}, 0, 0) == 0, + "append hit"); + CHECK(ray_pool_append(&pool, &observer, (double[]){1.0, 0.0, 0.0}, 0, 1) == 0, + "append miss"); + ray_pool_preroute(&pool, &source); + CHECK(pool.status[0] == RAY_POOL_PENDING && + pool.activate_t[0] < observer.coordinate_time - 1.0, + "entry ray stays pending until entry time"); + CHECK(pool.status[1] == RAY_POOL_TERMINATED && + pool.endpoint[1].status == RAY_ENDPOINT_ESCAPED, + "miss ray escapes during pre-route"); + + MetricSlab *early = NULL; + CHECK(spacetime_load_slab(&source, -20.0, -30.0, &early) == 0, "early slab"); + ray_pool_activate_in_time_range(&pool, early); + CHECK(pool.status[0] == RAY_POOL_PENDING, "entry ray not active early"); + spacetime_free_slab(early); + + MetricSlab *covering = NULL; + CHECK(spacetime_load_slab(&source, 0.0, -100.0, &covering) == 0, + "covering slab"); + ray_pool_activate_in_time_range(&pool, covering); + CHECK(pool.status[0] == RAY_POOL_ACTIVE && + fabs(pool.t[0] - pool.activate_t[0]) < 1e-30, + "entry ray activates at entry time"); + spacetime_free_slab(covering); + ray_pool_destroy(&pool); + spacetime_destroy(&source); +} + +int main(void) { + test_fixed_sphere(); + test_large_radius_quadratic(); + test_round_trip(); + test_moving_sphere(); + test_accelerated_worldtube(); + test_accelerated_segment(); + test_boundary_semantics_minkowski(); + test_boundary_semantics_generic(); + test_negative_radius_root_guard(); + test_source_finalize(); + test_motion_segment_domain(); + test_schwarzschild_sample_failures(); + test_end_protocol_error(); + test_generic_bisection_failure(); + test_interior_history_exhaustion(); + test_interior_crossing_bisection_failure(); + test_ray_pool_lifecycle(); + if (failures == 0) + puts("asymptotic regression passed"); + else + fprintf(stderr, "%d asymptotic regression failures\n", failures); + return failures == 0 ? 0 : 1; +} diff --git a/tests/test_asymptotic_schwarzschild.c b/tests/test_asymptotic_schwarzschild.c new file mode 100644 index 0000000..fd8bfb0 --- /dev/null +++ b/tests/test_asymptotic_schwarzschild.c @@ -0,0 +1,494 @@ +#include "asymptotic.h" +#include "asymptotic_schwarzschild.h" +#include "geodesic.h" +#include "observer.h" +#include "spacetime.h" + +#include +#include + +static int failures = 0; +#define CHECK(condition, message) \ + do { \ + if (!(condition)) { \ + fprintf(stderr, "FAIL %s:%d: %s\n", __FILE__, __LINE__, message); \ + ++failures; \ + } \ + } while (0) + +static double angle_between(const double a[3], const double b[3]) { + const double dot = a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; + const double cx = a[1] * b[2] - a[2] * b[1]; + const double cy = a[2] * b[0] - a[0] * b[2]; + const double cz = a[0] * b[1] - a[1] * b[0]; + return atan2(sqrt(cx * cx + cy * cy + cz * cz), dot); +} + +static void test_round_trip(void) { + SpacetimeSource source = {0}; + CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0, 1.5) == 0, + "create schwarzschild"); + SpacetimeAsymptoticEnd end; + CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end descriptor"); + SchwarzschildCanonical in = {.end_id = 0, + .t = 0.0, + .rho = 256.0, + .rhat = {1.0, 0.0, 0.0}, + .Lhat = {0.0, 1.0, 0.0}, + .beta = 5.0, + .energy = 1.0, + .radial_sign = 1}; + double x[3], Pi[3], log_alpha_p0; + CHECK(asymptotic_schwarzschild_state_from_canonical(&end, &in, x, Pi, + &log_alpha_p0) == 0, + "state from canonical"); + MetricData metric; + CHECK(spacetime_eval(&source, in.t, x, &metric) == 0, "metric"); + SchwarzschildCanonical out; + CHECK(asymptotic_schwarzschild_canonical_from_state( + &end, &metric, in.t, x, Pi, log_alpha_p0, &out) == 0, + "canonical from state"); + CHECK(fabs(out.beta - in.beta) < 1e-13, "beta round trip"); + CHECK(fabs(out.energy - in.energy) < 1e-13, "energy round trip"); + CHECK(out.radial_sign == in.radial_sign, "radial sign round trip"); + const double axis = angle_between(out.rhat, in.rhat); + CHECK(axis < 1e-13, "position direction round trip"); + spacetime_destroy(&source); +} + +static void test_finish_matches_integration(void) { + SpacetimeSource near = {0}, far = {0}; + CHECK(spacetime_create_schwarzschild_ks(&near, 1.0, 256.0, 1.5) == 0, + "create near"); + CHECK(spacetime_create_schwarzschild_ks(&far, 1.0, 1.0e5, 1.5) == 0, + "create far"); + SpacetimeAsymptoticEnd end; + CHECK(spacetime_asymptotic_end(&near, 0, &end) == 0, "near end"); + + const double betas[] = {0.0, 0.5, 4.0, 10.0, 30.0, 100.0, 250.0}; + const int beta_count = (int)(sizeof betas / sizeof betas[0]); + for (int k = 0; k < beta_count; ++k) { + SchwarzschildCanonical canonical = {.end_id = 0, + .t = 0.0, + .rho = 256.0, + .rhat = {0.8, 0.6, 0.0}, + .Lhat = {0.0, 0.0, 1.0}, + .beta = betas[k], + .energy = 1.0, + .radial_sign = 1}; + double x[3], Pi[3], log_alpha_p0; + CHECK(asymptotic_schwarzschild_state_from_canonical( + &end, &canonical, x, Pi, &log_alpha_p0) == 0, + "finish state build"); + + double n_analytic[3], freq_analytic; + CHECK(asymptotic_schwarzschild_finish(&end, &canonical, n_analytic, + &freq_analytic) == 0, + "analytic finish"); + + GeodesicRayState state = {.coordinate_time = 0.0, + .x = {x[0], x[1], x[2]}, + .Pi = {Pi[0], Pi[1], Pi[2]}, + .log_alpha_p0 = log_alpha_p0, + .steps = 0}; + const GeodesicTraceConfig config = {.coordinate_time_step = 5.0, + .max_steps = 100000}; + MetricSlab *slab = NULL; + CHECK(spacetime_load_slab(&far, 0.0, -1.0e6, &slab) == 0, "far slab"); + RayEndpoint endpoint = {.frequency_ratio = 0, .magnification = 1.0, + .end_id = SPACETIME_END_NONE, + .status = RAY_ENDPOINT_INVALID}; + const GeodesicAdvanceResult result = + geodesic_advance_past_ray(slab, &state, -1.0e6, &config, &endpoint); + spacetime_free_slab(slab); + CHECK(result == GEODESIC_ADVANCE_TERMINATED && + endpoint.status == RAY_ENDPOINT_ESCAPED, + "far integration escapes"); + /* Pipeline check only: the far integration at step 5 and escape radius + * 1e5 has its own O(1e-5..1e-3) error. Quantitative accuracy is checked + * against the high-precision reference constants below. */ + const double angle_error = + angle_between(n_analytic, endpoint.n_infinity); + CHECK(angle_error < 1e-2, "finish direction matches far integration"); + CHECK(fabs(freq_analytic - endpoint.frequency_ratio) / + freq_analytic < 1e-2, + "finish frequency matches far integration"); + (void)angle_error; + } + spacetime_destroy(&near); + spacetime_destroy(&far); +} + +/* Independent quadrature of the KS coordinate-time transfer for a camera + * outside the worldtube, used to check the analytic primitive. */ +static double simpson(const double a, const double b, int panels, + double (*f)(double, const void *), const void *ctx) { + if (panels < 2) + panels = 2; + if (panels % 2) + ++panels; + const double h = (b - a) / panels; + double sum = f(a, ctx) + f(b, ctx); + for (int i = 1; i < panels; ++i) + sum += (i % 2 ? 4.0 : 2.0) * f(a + i * h, ctx); + return sum * h / 3.0; +} + +typedef struct { + double beta; +} TransferContext; + +static double transfer_dt(double r, const void *context) { + const TransferContext *c = context; + const double Q = 1.0 - c->beta * c->beta * (1.0 - 2.0 / r) / (r * r); + return 1.0 / ((1.0 - 2.0 / r) * sqrt(Q)) + 2.0 / (r - 2.0); +} + +static double transfer_dphi(double r, const void *context) { + const TransferContext *c = context; + const double Q = 1.0 - c->beta * c->beta * (1.0 - 2.0 / r) / (r * r); + return c->beta / (r * r * sqrt(Q)); +} + +static void test_preroute_entry(void) { + SpacetimeSource source = {0}; + CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0, 1.5) == 0, + "create schwarzschild"); + const ObserverCamera camera = {.look_ra_deg = 0.0, .look_dec_deg = 0.0}; + ObserverCamera positioned = camera; + positioned.position[0] = 500.0; + positioned.look_ra_deg = 180.0; + positioned.look_dec_deg = 0.0; + const double direction[3] = {cos(0.3), sin(0.3), 0.0}; + MetricData metric; + CHECK(spacetime_eval(&source, 0.0, positioned.position, &metric) == 0, + "camera metric"); + ObserverState observer; + CHECK(observer_from_coordinate_camera(&metric, &positioned, &observer, + NULL) == OBSERVER_BUILD_OK, + "camera observer"); + + MetricSlab *camera_slab = NULL; + CHECK(spacetime_load_slab(&source, 0.0, -1.0, &camera_slab) == 0, + "camera slab"); + GeodesicRayState camera_state; + CHECK(geodesic_initialize_past_ray(camera_slab, &observer, direction, + &camera_state) == 0, + "camera state"); + SpacetimeAsymptoticEnd end; + CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end"); + SchwarzschildCanonical camera_can; + CHECK(asymptotic_schwarzschild_canonical_from_state( + &end, &metric, 0.0, camera_state.x, camera_state.Pi, + camera_state.log_alpha_p0, &camera_can) == 0, + "camera canonical"); + spacetime_free_slab(camera_slab); + + AsymptoticRoute route; + CHECK(asymptotic_route_camera(&source, &observer, direction, &route) == + ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ENTRY, + "outside camera enters"); + double value; + CHECK(asymptotic_worldtube_value(&source, route.end_id, route.activate_t, + route.x, &value) == 0 && + fabs(value) < 1e-3, + "entry on worldtube"); + + MetricData entry_metric; + CHECK(spacetime_eval(&source, route.activate_t, route.x, &entry_metric) == + 0, + "entry metric"); + SchwarzschildCanonical entry_can; + CHECK(asymptotic_schwarzschild_canonical_from_state( + &end, &entry_metric, route.activate_t, route.x, route.Pi, + route.log_alpha_p0, &entry_can) == 0, + "entry canonical"); + CHECK(fabs(entry_can.beta - camera_can.beta) < + 1e-12 * fmax(1.0, camera_can.beta), + "entry conserves impact parameter"); + CHECK(fabs(entry_can.energy - camera_can.energy) < 1e-12, + "entry conserves energy"); + CHECK(entry_can.radial_sign == -1, "entry is past-inward"); + CHECK(route.activate_t < 0.0, "entry time is in the past"); + + const TransferContext context = {.beta = camera_can.beta}; + const double t_analytic = -route.activate_t; + const double t_numeric = + simpson(256.0, 500.0, 20000, transfer_dt, &context); + CHECK(fabs(t_analytic - t_numeric) < 1e-9 * fmax(1.0, t_numeric), + "entry time matches quadrature"); + const double dphi_numeric = + simpson(256.0, 500.0, 20000, transfer_dphi, &context); + const double dphi_entry = angle_between(camera_can.rhat, entry_can.rhat); + CHECK(fabs(dphi_entry - dphi_numeric) < 1e-9, + "entry azimuth matches quadrature"); + if (fabs(t_analytic - t_numeric) >= 1e-9 * fmax(1.0, t_numeric) || + fabs(dphi_entry - dphi_numeric) >= 1e-9) + fprintf(stderr, " beta=%.6g t_an=%.12g t_num=%.12g dphi_an=%.12g " + "dphi_num=%.12g\n", + camera_can.beta, t_analytic, t_numeric, dphi_entry, + dphi_numeric); + spacetime_destroy(&source); +} + +/* High-precision (mpmath, 60 digits) reference values fixed into the ordinary + * C test: radial, complex-pair, three-real, grazing, and large-radius angle + * cases. */ +static void test_phi_reference_constants(void) { + static const struct { + double rho, beta, value; + } cases[] = { + {256.0, 0.0, 0.0}, + {256.0, 5.0, 0.019532484697919191145}, + {256.0, 60.0, 0.23656231243306290715}, + {64.0, 64.0, 1.4199914058161304301}, + {256.0, 255.0, 1.4527184167466732533}, + {1.0e6, 1.0, 1.0000000000001666664e-6}, + {300.0, 3.0, 0.010000165840750676787}, + {100.0, 5.3, 0.053024471018799209953}, + }; + for (size_t i = 0; i < sizeof cases / sizeof cases[0]; ++i) { + const double got = + asymptotic_schwarzschild_phi(cases[i].rho, cases[i].beta); + CHECK(fabs(got - cases[i].value) < 2e-13, "phi high-precision reference"); + } +} + +/* High-precision (mpmath, 60 digits) finish references covering radial, + * complex-pair, three-real, grazing, and large-radius scattering. The + * acceptance standard here is the error-budget-driven 1e-8 rad, not the + * measured ~1e-13. */ +static void test_finish_reference_constants(void) { + SpacetimeSource source = {0}; + CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0, 1.5) == 0, + "create schwarzschild"); + SpacetimeAsymptoticEnd end; + CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end"); + static const struct { + double rho, beta, n[3]; + } cases[] = { + {256.0, 0.0, {0.8, 0.6, 0.0}}, + {256.0, 3.0, {0.80697631554502468, 0.59058380112341107, 0.0}}, + {256.0, 60.0, {0.91833674808504193, 0.39579997109220491, 0.0}}, + {256.0, 255.0, {0.69006511550899122, -0.72374728763399356, 0.0}}, + {1.0e6, 1.0, {0.8000005999996, 0.5999991999997, 0.0}}, + }; + for (size_t i = 0; i < sizeof cases / sizeof cases[0]; ++i) { + SchwarzschildCanonical canonical = {.end_id = 0, + .t = 0.0, + .rho = cases[i].rho, + .rhat = {0.8, 0.6, 0.0}, + .Lhat = {0.0, 0.0, 1.0}, + .beta = cases[i].beta, + .energy = 2.5, + .radial_sign = 1}; + double n_inf[3], frequency = 0.0; + CHECK(asymptotic_schwarzschild_finish(&end, &canonical, n_inf, + &frequency) == 0, + "finish reference runs"); + CHECK(angle_between(n_inf, cases[i].n) < 1e-8, + "finish n_inf high-precision reference"); + CHECK(fabs(frequency - 0.4) < 1e-10 * 0.4, + "finish frequency high-precision reference"); + } + spacetime_destroy(&source); +} + +/* Turning equation residual |Q| at the computed turning radius. The final + * scattering direction is validated by test_grazing_reference(). */ +static void test_turning_reference(void) { + const double betas[] = {3.0 * sqrt(3.0) + 1e-9, 5.5, 6.0, 10.0, + 60.0, 255.0, 3890.44}; + for (size_t i = 0; i < sizeof betas / sizeof betas[0]; ++i) { + const double rho = asymptotic_schwarzschild_turning_rho(betas[i]); + CHECK(isfinite(rho) && rho > 3.0, "turning radius exists and is exterior"); + const double Q = + 1.0 - betas[i] * betas[i] * (1.0 - 2.0 / rho) / (rho * rho); + CHECK(fabs(Q) <= 1e-11, "turning equation residual"); + } +} + +/* High-precision entry coordinate-time and swept-azimuth references, checking + * both the KS time transfer and the entry direction construction. */ +static void test_time_reference(void) { + SpacetimeSource source = {0}; + CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0, 1.5) == 0, + "create schwarzschild"); + SpacetimeAsymptoticEnd end; + CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end"); + static const struct { + double rho_cam, beta, time, dphi; + } cases[] = { + {500.0, 10.0, 246.7884398934447041137, 0.01907105306677434549856}, + {500.0, 0.3, 246.693149021326385798, 0.0005718752307574524347376}, + {256.5, 10.0, 0.5082474340167056157528, 0.00007620281793853560952548}, + {256.5, 0.3, 0.5078666185420216617125, 0.000002284358278417004222584}, + {1000.0, 50.0, 753.1528987272233278083, 0.1465476883815797019938}, + {1.0e6, 10.0, 999777.308033789182372, + 0.03906238263856681534781}, + }; + for (size_t i = 0; i < sizeof cases / sizeof cases[0]; ++i) { + SchwarzschildCanonical camera = {.end_id = 0, + .t = 0.0, + .rho = cases[i].rho_cam, + .rhat = {1.0, 0.0, 0.0}, + .Lhat = {0.0, 0.0, 1.0}, + .beta = cases[i].beta, + .energy = 1.0, + .radial_sign = -1}; + SchwarzschildRouteKind kind = SCH_ROUTE_UNSUPPORTED; + double activate_t = 0.0, x[3], Pi[3], log_alpha_p0 = 0.0, n_inf[3], + frequency = 0.0; + CHECK(asymptotic_schwarzschild_preroute( + &end, 256.0, &camera, &kind, &activate_t, x, Pi, &log_alpha_p0, + n_inf, &frequency) == 0 && + kind == SCH_ROUTE_ENTRY, + "reference pre-route entry"); + /* Error-budget-driven mixed tolerance, well below one ODE step (0.1 M) + * and future metric cadence. */ + const double time_tol = 1e-7 + 1e-11 * fabs(cases[i].time); + CHECK(fabs(-activate_t - cases[i].time) < time_tol, + "entry time high-precision reference"); + const double radius = + sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]); + const double rhat[3] = {x[0] / radius, x[1] / radius, x[2] / radius}; + CHECK(fabs(angle_between(camera.rhat, rhat) - cases[i].dphi) < 2e-11, + "entry azimuth high-precision reference"); + } + spacetime_destroy(&source); +} + +/* Near-grazing references where the exterior integrals are most sensitive: + * the two sides of beta_R enter through different branches and the KS time + * integral has a near-singular endpoint. (A photon-sphere turning is not + * reachable from a camera outside R/M >= 64, so it is not tested here.) */ +static void test_grazing_reference(void) { + SpacetimeSource source = {0}; + CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0, 1.5) == 0, + "create schwarzschild"); + SpacetimeAsymptoticEnd end; + CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end"); + const double beta_R = 256.0 / sqrt(1.0 - 2.0 / 256.0); + SchwarzschildCanonical hit = {.end_id = 0, + .t = 0.0, + .rho = 500.0, + .rhat = {1.0, 0.0, 0.0}, + .Lhat = {0.0, 0.0, 1.0}, + .beta = beta_R * (1.0 - 1e-12), + .energy = 1.0, + .radial_sign = -1}; + SchwarzschildRouteKind kind; + double activate_t, x[3], Pi[3], log_alpha_p0, n_inf[3], frequency; + CHECK(asymptotic_schwarzschild_preroute(&end, 256.0, &hit, &kind, + &activate_t, x, Pi, &log_alpha_p0, + n_inf, &frequency) == 0 && + kind == SCH_ROUTE_ENTRY, + "near-grazing inside enters"); + const double dphi_ref = 1.0389037630217253661; + const double time_ref = 434.0116073725480308524; + const double radius = sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]); + const double rhat[3] = {x[0] / radius, x[1] / radius, x[2] / radius}; + CHECK(fabs(angle_between(hit.rhat, rhat) - dphi_ref) < 1e-8, + "near-grazing entry azimuth"); + CHECK(fabs(-activate_t - time_ref) < 1e-7 + 1e-11 * time_ref, + "near-grazing entry time"); + + SchwarzschildCanonical miss = hit; + miss.beta = beta_R * (1.0 + 1e-12); + CHECK(asymptotic_schwarzschild_preroute(&end, 256.0, &miss, &kind, + &activate_t, x, Pi, &log_alpha_p0, + n_inf, &frequency) == 0 && + kind == SCH_ROUTE_ESCAPED, + "near-grazing outside misses"); + const double n_ref[3] = {-0.86581533530640297059, + -0.50036367289028986187, 0.0}; + CHECK(angle_between(n_inf, n_ref) < 1e-8, "near-grazing miss n_inf"); + spacetime_destroy(&source); +} + +/* Deterministic coverage of the three pre-route branches: past-outward, + * past-inward hit, and past-inward miss (turn before the worldtube). */ +static void test_preroute_branches(void) { + SpacetimeSource source = {0}; + CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0, 1.5) == 0, + "create schwarzschild"); + SpacetimeAsymptoticEnd end; + CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end"); + const double beta_R = 256.0 / sqrt(1.0 - 2.0 / 256.0); + + SchwarzschildCanonical base = {.end_id = 0, + .t = 0.0, + .rho = 500.0, + .rhat = {1.0, 0.0, 0.0}, + .Lhat = {0.0, 0.0, 1.0}, + .beta = 10.0, + .energy = 1.0, + .radial_sign = -1}; + SchwarzschildRouteKind kind; + double activate_t, x[3], Pi[3], log_alpha_p0, n_inf[3], frequency; + const double outward_eps = 1e-12; + + SchwarzschildCanonical outward = base; + outward.radial_sign = 1; + CHECK(asymptotic_schwarzschild_preroute(&end, 256.0, &outward, &kind, + &activate_t, x, Pi, &log_alpha_p0, + n_inf, &frequency) == 0 && + kind == SCH_ROUTE_ESCAPED, + "past-outward branch escapes"); + CHECK(fabs(sqrt(n_inf[0]*n_inf[0]+n_inf[1]*n_inf[1]+n_inf[2]*n_inf[2]) - + 1.0) < outward_eps, + "outward n_inf is unit"); + CHECK(fabs(frequency - 1.0) < 1e-12, "outward frequency"); + + SchwarzschildCanonical hit = base; + CHECK(asymptotic_schwarzschild_preroute(&end, 256.0, &hit, &kind, + &activate_t, x, Pi, &log_alpha_p0, + n_inf, &frequency) == 0 && + kind == SCH_ROUTE_ENTRY, + "past-inward hit branch enters"); + + /* Genuine on-boundary tangent: rho = R, beta = beta_R (so Q = 0), zero + * radial past component. It must not enter. */ + SchwarzschildCanonical tangent = base; + tangent.rho = 256.0; + tangent.beta = beta_R; + tangent.radial_sign = 0; + CHECK(asymptotic_schwarzschild_preroute(&end, 256.0, &tangent, &kind, + &activate_t, x, Pi, &log_alpha_p0, + n_inf, &frequency) == 0 && + kind == SCH_ROUTE_ESCAPED, + "on-boundary tangent escapes"); + + SchwarzschildCanonical miss = base; + miss.beta = beta_R + 5.0; + CHECK(asymptotic_schwarzschild_preroute(&end, 256.0, &miss, &kind, + &activate_t, x, Pi, &log_alpha_p0, + n_inf, &frequency) == 0 && + kind == SCH_ROUTE_ESCAPED, + "past-inward miss branch escapes"); + CHECK(fabs(sqrt(n_inf[0]*n_inf[0]+n_inf[1]*n_inf[1]+n_inf[2]*n_inf[2]) - + 1.0) < outward_eps, + "miss n_inf is unit"); + /* A turning ray is deflected away from the radial direction. */ + CHECK(angle_between(n_inf, miss.rhat) > 1e-3, + "miss n_inf is deflected"); + spacetime_destroy(&source); +} + +int main(void) { + test_round_trip(); + test_finish_matches_integration(); + test_preroute_entry(); + test_phi_reference_constants(); + test_finish_reference_constants(); + test_turning_reference(); + test_time_reference(); + test_grazing_reference(); + test_preroute_branches(); + if (failures == 0) + puts("asymptotic schwarzschild regression passed"); + else + fprintf(stderr, "%d asymptotic schwarzschild failures\n", failures); + return failures == 0 ? 0 : 1; +} diff --git a/tests/test_observer.c b/tests/test_observer.c index f3f43cc..8813c4d 100644 --- a/tests/test_observer.c +++ b/tests/test_observer.c @@ -123,10 +123,11 @@ int main(int argc, char **argv) { const RayEndpoint ray = geodesic_trace_past(&source, &state, (double[]){1, 0, 0}, &trace); CHECK(ray.status == RAY_ENDPOINT_ESCAPED); CHECK(fabs(ray.n_infinity[0] - 1) < 1e-12); - /* Radial ingoing KS photon has k^r=-k^t and conserved E=k^t. - * Current escape convention measures Eulerian energy at finite R=256. */ + /* Radial ingoing KS photon has k^r=-k^t and conserved E=k^t. The + * asymptotic exterior transfers the photon to infinity, where + * g = E_camera / E_infinity = 1 / k^t. */ const double energy = state.tetrad[0][0] - state.tetrad[1][0]; - CHECK(fabs(ray.frequency_ratio - sqrt(1 + 2.0 / 256) / energy) < 2e-6); + CHECK(fabs(ray.frequency_ratio - 1.0 / energy) < 1e-10 * (1.0 / energy)); memset(camera.velocity, 0, sizeof camera.velocity); if (i > 0) CHECK(observer_from_coordinate_camera(&metric, &camera, &state, NULL) == OBSERVER_BUILD_NON_TIMELIKE);