diff --git a/README.md b/README.md index 1f3b7ed..3e40a7b 100644 --- a/README.md +++ b/README.md @@ -185,3 +185,49 @@ mkdir -p output/imgs [![Synthetic stellar grid lensed by a Schwarzschild black hole, with adaptive mesh overlay](assets/images/schwarzschild_test_grid.png)](assets/images/schwarzschild_test_grid.png) *4K test-grid reference image. Click to view at full resolution.* + +### Freely falling Schwarzschild movie camera + +`scripts/schwarzschild_camera_track.py` generates the canonical 21-column observer +CSV in ingoing Cartesian Kerr–Schild coordinates (`G=c=M=1`, matching the renderer). +It requires Python 3, NumPy and SciPy. Position is `(x,y,z)`; velocity is coordinate +`(dx/dt,dy/dt,dz/dt)`. Look RA/Dec and roll construct the initial rest-frame +forward/up/right legs with the same convention as the single-image camera. Without +look angles, the camera initially points toward the origin. The orientation then +follows Fermi–Walker transport; for free fall this is parallel transport, so it does +not keep pointing at the black hole. `--tetrad` alternatively accepts 16 row-major +components `(e0t,e0x,...,e3z)` of an orthonormal tetrad, with `e0` matching velocity. + +```sh +python3 scripts/schwarzschild_camera_track.py \ + --position 8 0 0 --velocity 0 0 0 \ + --look-ra-deg 0 --look-dec-deg 0 --roll-deg 0 \ + --fps 30 --duration 10 --output /tmp/freefall_camera.csv +make SPACETIME=schwarzschild all +mkdir -p output/freefall_frames +./build/Release/schwarzschild_sky \ + --observer-track /tmp/freefall_camera.csv --movie-track-samples \ + --frames-dir output/freefall_frames --catalog assets/sky_grid_5deg.csv \ + --width 640 --height 360 --fov-deg 60 --exposure 0.2 +ffmpeg -framerate 30 -i output/freefall_frames/frame_%06d.png \ + -c:v libx264 -pix_fmt yuv420p output/freefall.mp4 +``` + +Here `--duration` is elapsed **proper time**, and `--fps` is samples per unit proper +time. CSV rows occur at `tau=k/fps <= duration`, including the initial event and an +endpoint only if it lies on that cadence (10 at 30 fps gives 301 rows). `tau` starts +at zero; `--t0` sets the initial coordinate time. `--movie-track-samples` uses each +row exactly once and ignores renderer `--start-time`, `--duration`, and `--fps`; +encode the PNG sequence at the generator's fps. Without this flag the existing +movie mode resamples at uniform coordinate time, which changes the proper-time cadence. + +DOP853 jointly integrates the geodesic and all tetrad legs using analytic metric +derivatives. Defaults: `--rtol 1e-10 --atol 1e-12 --stop-radius 0.001`. +Integration crosses the horizon and stops at this numerical guard before the +singularity, reporting its proper time and retaining only regular cadence samples. +The guard is not the exact singularity; reduce it and tolerances to check convergence. +The renderer's independent ray capture cutoff remains `r=1.5M`: rows inside it are +valid trajectory data but the current renderer captures those rays immediately. +The script reports maximum tetrad drift and rejects errors above `1e-6` rather +than silently repairing the transported frame. Run the orbit, transport and CSV +render regressions after building with `python3 tests/test_schwarzschild_camera_track.py`. diff --git a/README.zh-CN.md b/README.zh-CN.md index 83161b2..5986406 100644 --- a/README.zh-CN.md +++ b/README.zh-CN.md @@ -132,3 +132,44 @@ mkdir -p output/imgs [![Schwarzschild 黑洞对合成恒星网格的透镜效果，叠加自适应网格](assets/images/schwarzschild_test_grid.png)](assets/images/schwarzschild_test_grid.png) *4K 测试网格参考图像。点击查看完整分辨率。* + +### Schwarzschild 自由落体相机轨迹 + +`scripts/schwarzschild_camera_track.py` 生成 movie 使用的 21 列相机 CSV，采用 +入射 Cartesian Kerr–Schild 坐标及 `G=c=M=1`，依赖 Python 3、NumPy、SciPy。 +`--position` 是坐标位置，`--velocity` 是 `dx/dt,dy/dt,dz/dt`，不是局域三速度。 +初始标架由 RA/Dec 与 roll 构造，约定与单张相机一致；省略指向时初始朝向原点。 +也可用 `--tetrad` 提供 16 个按行排列的四标架分量，顺序为 `e0t,e0x,...,e3z`， +要求正交归一、定向正确，且 `e0` 与指定速度一致。自由落体中费米–沃克输运等同于 +平行输运；初始之后不再强制朝向黑洞。 + +```sh +python3 scripts/schwarzschild_camera_track.py \ + --position 8 0 0 --velocity 0 0 0 \ + --look-ra-deg 0 --look-dec-deg 0 --roll-deg 0 \ + --fps 30 --duration 10 --output /tmp/freefall_camera.csv +make SPACETIME=schwarzschild all +mkdir -p output/freefall_frames +./build/Release/schwarzschild_sky \ + --observer-track /tmp/freefall_camera.csv --movie-track-samples \ + --frames-dir output/freefall_frames --catalog assets/sky_grid_5deg.csv \ + --width 640 --height 360 --fov-deg 60 --exposure 0.2 +ffmpeg -framerate 30 -i output/freefall_frames/frame_%06d.png \ + -c:v libx264 -pix_fmt yuv420p output/freefall.mp4 +``` + +脚本的 `--duration` 是持续本征时（单位 M），`--fps` 是每单位本征时的采样数。 +采样为 `tau=k/fps <= duration`，包含初始帧，只在终点恰好满足采样节奏时包含终点； +10 M、30 fps 共 301 行。`tau` 从零开始，`--t0` 指定初始坐标时间。 +新增 `--movie-track-samples` 让每行恰好对应一帧，保留真实坐标时间，并忽略 renderer +的 `--start-time/--duration/--fps`；编码时使用生成脚本的 fps。 +不加此选项时原有 movie 仍按坐标时间等间隔采样，不能保持这里的本征时节奏。 + +脚本以解析 metric 导数及 DOP853 联合积分测地线与四标架，默认 +`--rtol 1e-10 --atol 1e-12 --stop-radius 0.001`。轨迹穿过视界继续积分，在该小半径 +数值保护边界停止，报告终止本征时，只保留此前的规则采样；不声称到达精确的 `r=0`。 +可减小半径及误差容限检查收敛。现有 renderer 的光线捕获边界仍为 `r=1.5M`： +CSV 可以记录更深处的相机，但目前 renderer 会把这些相机发出的光线立即判为捕获。 +脚本报告标架正交归一误差，超过 `1e-6` 时直接报错，不自动修正输运后的标架。 +构建后运行 `python3 tests/test_schwarzschild_camera_track.py`，检查解析径向落体、 +圆轨道及标架输运收敛、非法初值和实际 CSV 渲染。 diff --git a/nr_spacetime_movie_renderer_design.md b/nr_spacetime_movie_renderer_design.md index 561b562..f8383d9 100644 --- a/nr_spacetime_movie_renderer_design.md +++ b/nr_spacetime_movie_renderer_design.md @@ -1359,3 +1359,26 @@ renderer 顶层架构原则上不应为 BBH 重新设计。 13. **所有高开销 mutable cache 都 thread-local。** 14. **第一版 CPU-only，先把物理与数据流做正确，再谈 GPU。** 15. **从 Minkowski → analytic Schwarzschild → numerical Schwarzschild → BBH 逐级验证。** + +## Schwarzschild 自由落体轨迹辅助脚本 + +`scripts/schwarzschild_camera_track.py` 在入射 Cartesian Kerr–Schild 坐标中使用 +$g_{\mu\nu}=\eta_{\mu\nu}+(2/r)\ell_\mu\ell_\nu$、$\ell_\mu=(1,x_i/r)$， +固定 $G=c=M=1$ 与解析 backend 一致。解析微分 metric 构造四维 Christoffel， +以本征时联合积分 +$dz^\mu/d\tau=u^\mu$ 与 +$de_{(a)}^\mu/d\tau=-\Gamma^\mu_{\alpha\beta}u^\alpha e_{(a)}^\beta$。 +$e_{(0)}=u$ 同时满足自由落体方程；四加速度为零时费米–沃克输运就是平行输运。 +初始坐标速度及指向构造沿用 10.1 的约定，也接受经验证的完整标架。 +不在积分期间重新正交化以隐藏误差；输出正交归一误差超过 $10^{-6}$ 时失败。 + +DOP853 默认 rtol=$10^{-10}$、atol=$10^{-12}$，可配置并通过解析径向自由落体、 +圆轨道和圆轨道平行输运的收敛回归验证。视界不终止相机；默认 $r=10^{-3}M$ +只是可配置的奇点数值保护边界，不等于精确撞击奇点。它独立于光线的 $1.5M$ +捕获 cutoff；该 cutoff 内的轨迹可输出，但当前 renderer 的光线会立即被捕获。 + +CSV 保持 21 列不变，记录 $\tau=k/\mathrm{fps}$ 与积分得到的真实坐标时间。 +只保留不超过请求持续本征时或提前终止时刻的规则采样，包含 $\tau=0$。 +`--movie-track-samples` 选择每个 CSV 样本直接生成一帧，绕过坐标时间均匀插值， +忽略 renderer 的 start-time/duration/fps。默认 movie 路径不变；两条路径随后共用 +按真实 coordinate time 的 time-slab 调度，不把本征时冒充坐标时间。 diff --git a/scripts/schwarzschild_camera_track.py b/scripts/schwarzschild_camera_track.py new file mode 100644 index 0000000..51a7c05 --- /dev/null +++ b/scripts/schwarzschild_camera_track.py @@ -0,0 +1,165 @@ +#!/usr/bin/env python3 +"""Free-fall camera with parallel (= geodesic Fermi-Walker) transport. + +Ingoing Cartesian Kerr-Schild, signature -+++, G=c=M=1. Requires numpy/scipy. +The renderer's Schwarzschild backend fixes M=1 as well. +""" +import argparse +import csv +import math +from pathlib import Path +import sys + +import numpy as np +from scipy.integrate import solve_ivp + +ETA = np.diag([-1., 1., 1., 1.]) +HEADER = ['t', 'tau', 'x', 'y', 'z'] + [f'e{a}{c}' for a in range(4) for c in 'txyz'] + + +def metric_connection(x): + """Analytic g and Christoffels; dg[k,mu,nu] = partial_k g_mu_nu.""" + r = np.linalg.norm(x) + if not np.isfinite(r) or r <= 0: + raise ValueError('metric undefined at r=0 or nonfinite position') + n = x / r + ell = np.r_[1., n] + f = 2 / r + g = ETA + f * np.outer(ell, ell) + raised = ETA @ ell + inverse = ETA - f * np.outer(raised, raised) + dg = np.zeros((4, 4, 4)) + for k in range(3): + dl = np.r_[0., (np.eye(3)[k] - n[k] * n) / r] + dg[k+1] = f * (np.outer(dl, ell) + np.outer(ell, dl) + - n[k] / r * np.outer(ell, ell)) + connection = .5 * np.einsum('ml,alb->mab', inverse, + dg + dg.transpose(2, 1, 0) - dg.transpose(1, 0, 2)) + return g, connection + + +def initial_state(position, velocity, ra=None, dec=None, roll=0., tetrad=None, t0=0.): + g, _ = metric_connection(position) + u = np.r_[1., velocity] + q = u @ g @ u + if not np.isfinite(q) or q >= 0: + raise ValueError('coordinate velocity must be future timelike: g(1,v;1,v) < 0') + u /= math.sqrt(-q) + if tetrad is not None: + e = np.asarray(tetrad, dtype=float) + if e.shape != (4, 4) or not np.all(np.isfinite(e)): + raise ValueError('initial tetrad must contain four rows of four finite components') + if not np.allclose(e[0], u, rtol=1e-9, atol=1e-9): + raise ValueError('initial tetrad e0 must agree with the specified coordinate velocity') + if np.max(np.abs(e @ g @ e.T - ETA)) > 1e-8: + raise ValueError('initial tetrad must be Lorentz orthonormal') + if np.linalg.det(e) <= 0: + raise ValueError('initial tetrad must have forward cross up = right orientation') + else: + if ra is None: + direction = -np.asarray(position) / np.linalg.norm(position) + ra = math.degrees(math.atan2(direction[1], direction[0])) % 360 + dec = math.degrees(math.asin(direction[2])) + a, d = np.deg2rad([ra, dec]) + e = np.array([u, [0, math.cos(d)*math.cos(a), math.cos(d)*math.sin(a), math.sin(d)], + [0, -math.sin(d)*math.cos(a), -math.sin(d)*math.sin(a), math.cos(d)], + [0, math.sin(a), -math.cos(a), 0]]) + for i in range(1, 4): + for _ in range(2): + for j in range(i): + e[i] -= (-1 if j == 0 else 1) * (e[i] @ g @ e[j]) * e[j] + e[i] /= math.sqrt(e[i] @ g @ e[i]) + angle = math.radians(roll) + up, right = e[2].copy(), e[3].copy() + e[2] = math.cos(angle)*up + math.sin(angle)*right + e[3] = -math.sin(angle)*up + math.cos(angle)*right + return np.r_[t0, position, e.ravel()] + + +def rhs(tau, state): + _, connection = metric_connection(state[1:4]) + e = state[4:].reshape(4, 4) + # e0 is u: transporting all four legs also integrates the timelike geodesic. + de = -np.einsum('mab,a,ib->im', connection, e[0], e) + return np.r_[e[0], de.ravel()] + + +def integrate(state, duration, fps, stop_radius=1e-3, rtol=1e-10, atol=1e-12): + if np.linalg.norm(state[1:4]) <= stop_radius: + raise ValueError('initial radius must exceed stop-radius') + if duration == 0: + return np.array([0.]), state[None, :], None + + def stop(tau, y): + return np.linalg.norm(y[1:4]) - stop_radius + stop.terminal = True + stop.direction = -1 + solution = solve_ivp(rhs, (0., duration), state, method='DOP853', + rtol=rtol, atol=atol, events=stop, dense_output=True) + if not solution.success: + raise ValueError(f'integration failed at tau={solution.t[-1]:.17g}: {solution.message}') + end = solution.t[-1] + # Include tau=0 and every full 1/fps interval; never add an off-cadence endpoint. + count = int(math.floor(np.nextafter(end * fps, np.inf))) + 1 + tau = np.arange(count, dtype=float) / fps + tau = tau[tau <= end] + states = solution.sol(tau).T + if not np.all(np.isfinite(states)) or np.any(np.diff(states[:, 0]) <= 0): + raise ValueError('nonfinite trajectory or coordinate time lost monotonicity') + return tau, states, end if solution.t_events[0].size else None + + +def main(): + p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.ArgumentDefaultsHelpFormatter) + p.add_argument('--output', type=Path, required=True, help='21-column movie CSV') + p.add_argument('--position', nargs=3, type=float, required=True, metavar=('X', 'Y', 'Z'), help='initial KS Cartesian position in M') + p.add_argument('--velocity', nargs=3, type=float, default=[0., 0., 0.], help='coordinate dx/dt, dy/dt, dz/dt (not local 3-speed)') + p.add_argument('--look-ra-deg', type=float, help='initial forward RA; default points toward origin') + p.add_argument('--look-dec-deg', type=float, help='initial forward Dec; specify together with RA') + p.add_argument('--roll-deg', type=float, default=0., help='initial roll, same sign as renderer') + p.add_argument('--tetrad', type=float, nargs=16, help='explicit row-major e0,e1,e2,e3 in t,x,y,z; replaces look/roll') + p.add_argument('--t0', type=float, default=0., help='initial KS coordinate time') + p.add_argument('--fps', type=float, default=30., help='samples per unit proper time M') + p.add_argument('--duration', type=float, required=True, help='requested elapsed proper time in M') + p.add_argument('--stop-radius', type=float, default=1e-3, help='numerical singularity guard in M, strictly between 0 and 2') + p.add_argument('--rtol', type=float, default=1e-10, help='DOP853 relative error tolerance') + p.add_argument('--atol', type=float, default=1e-12, help='DOP853 absolute error tolerance') + args = p.parse_args() + try: + numbers = [*args.position, *args.velocity, args.roll_deg, args.t0, args.fps, + args.duration, args.stop_radius, args.rtol, args.atol] + numbers += [v for v in (args.look_ra_deg, args.look_dec_deg) if v is not None] + if not all(math.isfinite(v) for v in numbers): + raise ValueError('all numeric arguments must be finite') + if args.fps <= 0 or args.duration < 0 or not 0 < args.stop_radius < 2 or min(args.rtol, args.atol) <= 0: + raise ValueError('require fps,rtol,atol > 0, duration >= 0 and 0 < stop-radius < 2') + if (args.look_ra_deg is None) != (args.look_dec_deg is None): + raise ValueError('specify both look-ra-deg and look-dec-deg') + if args.look_ra_deg is not None and not (0 <= args.look_ra_deg < 360 and abs(args.look_dec_deg) <= 90): + raise ValueError('require 0 <= RA < 360 and -90 <= Dec <= 90') + if args.tetrad is not None and (args.look_ra_deg is not None or args.roll_deg != 0): + raise ValueError('tetrad conflicts with look/roll') + e = None if args.tetrad is None else np.array(args.tetrad).reshape(4, 4) + state = initial_state(args.position, args.velocity, args.look_ra_deg, + args.look_dec_deg, args.roll_deg, e, args.t0) + tau, states, stopped = integrate(state, args.duration, args.fps, + args.stop_radius, args.rtol, args.atol) + error = max(np.max(np.abs(s[4:].reshape(4, 4) @ metric_connection(s[1:4])[0] + @ s[4:].reshape(4, 4).T - ETA)) for s in states) + if error > 1e-6: + raise ValueError(f'tetrad norm drift {error:.3g} exceeds 1e-6; tighten tolerances') + with args.output.open('w', newline='') as stream: + writer = csv.writer(stream) + writer.writerow(HEADER) + for t, s in zip(tau, states): + writer.writerow(format(v, '.17g') for v in np.r_[s[0], t, s[1:]]) + print(f'Wrote {len(tau)} frames; tau=0..{tau[-1]:.17g}, t={states[0,0]:.17g}..{states[-1,0]:.17g}; max tetrad error={error:.3g}', file=sys.stderr) + if stopped is not None: + print(f'Stopped early at r={args.stop_radius:g} M, tau={stopped:.17g} (numerical guard before r=0).', file=sys.stderr) + print('Render with --observer-track --movie-track-samples --frames-dir ; encode at the chosen fps.', file=sys.stderr) + except (ValueError, OSError, OverflowError) as exc: + p.error(str(exc)) + + +if __name__ == '__main__': + main() diff --git a/src/main.c b/src/main.c index f0172b4..210e2ac 100644 --- a/src/main.c +++ b/src/main.c @@ -44,6 +44,7 @@ typedef struct { char hdr_output_path[PATH_MAX]; #endif const char *observer_track_path; + int movie_track_samples; const char *frames_dir; const char *frames_prefix; const char *write_minkowski_accel_track_path; @@ -299,6 +300,8 @@ static int parse_args(int argc, char **argv, Settings *s, *write_path = argv[++i]; else if (!strcmp(argv[i], "--observer-track") && i + 1 < argc) s->observer_track_path = argv[++i]; + else if (!strcmp(argv[i], "--movie-track-samples")) + s->movie_track_samples = 1; else if (!strcmp(argv[i], "--frames-dir") && i + 1 < argc) s->frames_dir = argv[++i]; else if (!strcmp(argv[i], "--frames-prefix") && i + 1 < argc) @@ -385,6 +388,7 @@ static void print_help(const char *program) { " --draw-mesh Draw the final lens mesh overlay (default: disabled)\n" "\nMovie and observer track:\n" " --observer-track PATH Observer worldline/tetrad CSV for movie rendering (default: disabled)\n" + " --movie-track-samples One frame per CSV row; ignores start-time/duration/fps (default: disabled)\n" " --frames-dir DIR Write a movie image sequence to this directory (default: disabled)\n" " --frames-prefix NAME Movie frame filename prefix (default: frame)\n" " --start-time T Movie start coordinate time (default: 0)\n" @@ -508,6 +512,11 @@ static GeodesicTraceConfig trace_config(void) { } static int resolve_camera(Settings *s) { + if (s->movie_track_samples && (!s->observer_track_path || !s->frames_dir || + s->lens_map_input_path)) { + fputs("--movie-track-samples requires --observer-track and --frames-dir, without --lens-map-input.\n", stderr); + return -1; + } const int camera_specified = s->position_specified || s->look_specified || s->radius_specified || s->velocity_specified || s->roll_specified; if (s->position_specified && s->radius_specified) { @@ -798,8 +807,9 @@ static int render_movie(const Settings *s, StarCatalog *catalog, int result = -1; if (s->observer_track_path == NULL || observer_track_load_csv(&track, s->observer_track_path) || - movie_init(&movie, &track, s->movie_start_time, s->movie_duration, - s->movie_fps) || + (s->movie_track_samples ? movie_init_track_samples(&movie, &track) : + movie_init(&movie, &track, s->movie_start_time, s->movie_duration, + s->movie_fps)) || movie_build_coarse_meshes(&movie, s->width, s->height, s->coarse_cell_pixels, s->horizontal_fov_deg)) goto done; @@ -990,7 +1000,7 @@ int main(int argc, char **argv) { "[--draw-mesh] [--write-catalog PATH] " "[--catalog-load-workers N] " "[--observer-track PATH --frames-dir DIR --frames-prefix NAME " - "--start-time T --duration T --fps N] " + "--start-time T --duration T --fps N | --movie-track-samples] " "[--proper-acceleration A --write-minkowski-accel-track PATH]\n", argv[0]); return 2; diff --git a/src/movie.c b/src/movie.c index eeb2c43..4c7774f 100644 --- a/src/movie.c +++ b/src/movie.c @@ -3,6 +3,28 @@ #include #include #include +#include + +int movie_init_track_samples(Movie *movie, const ObserverTrack *track) { + if (!movie || !track || !track->samples || !track->count || + track->count > SIZE_MAX / sizeof *movie->frames) + return -1; + *movie = (Movie){0}; + movie->frames = calloc(track->count, sizeof *movie->frames); + if (!movie->frames) return -1; + movie->frame_count = track->count; + for (size_t i = 0; i < track->count; ++i) { + const ObserverSample *s = &track->samples[i]; + MovieFrame *f = &movie->frames[i]; + f->frame_id = i; + f->coordinate_time = f->observer.coordinate_time = s->coordinate_time; + f->proper_time = s->proper_time; + memcpy(f->observer.coordinate_position, s->coordinate_position, + sizeof s->coordinate_position); + memcpy(f->observer.tetrad, s->tetrad, sizeof s->tetrad); + } + return 0; +} int movie_init(Movie *movie, const ObserverTrack *track, double start_time, double duration, double frames_per_second) { diff --git a/src/movie.h b/src/movie.h index 4cc17d8..d09fc1a 100644 --- a/src/movie.h +++ b/src/movie.h @@ -21,6 +21,8 @@ typedef struct { int movie_init(Movie *movie, const ObserverTrack *track, double start_time, double duration, double frames_per_second); +/* Preserve the supplied event and tetrad exactly: one CSV sample per frame. */ +int movie_init_track_samples(Movie *movie, const ObserverTrack *track); int movie_build_coarse_meshes(Movie *movie, int width, int height, int cell_pixels, double horizontal_fov_deg); void movie_destroy(Movie *movie); diff --git a/tests/test_observer_track.c b/tests/test_observer_track.c index 5dc911e..474723e 100644 --- a/tests/test_observer_track.c +++ b/tests/test_observer_track.c @@ -30,6 +30,23 @@ int main(void) { -(sqrt(1.0 + 1.52 * 1.52) - 1.0) / 1.52) || !nearly_equal(proper_time, asinh(1.52) / 1.52)) goto done; + movie_destroy(&movie); + /* Nonuniform coordinate times must survive the row-per-frame path exactly. */ + for (size_t i = 0; i < loaded.count; ++i) + loaded.samples[i].coordinate_time += 0.001 * i * i; + if (movie_init_track_samples(&movie, &loaded) || + movie.frame_count != loaded.count) + goto done; + for (size_t i = 0; i < loaded.count; ++i) { + if (movie.frames[i].coordinate_time != loaded.samples[i].coordinate_time || + movie.frames[i].observer.coordinate_time != loaded.samples[i].coordinate_time || + movie.frames[i].proper_time != loaded.samples[i].proper_time) + goto done; + for (int a = 0; a < 4; ++a) + for (int mu = 0; mu < 4; ++mu) + if (movie.frames[i].observer.tetrad[a][mu] != loaded.samples[i].tetrad[a][mu]) + goto done; + } { FILE *bad = fopen(path, "w"); ObserverTrack invalid = {0}; diff --git a/tests/test_schwarzschild_camera_track.py b/tests/test_schwarzschild_camera_track.py new file mode 100644 index 0000000..4361dda --- /dev/null +++ b/tests/test_schwarzschild_camera_track.py @@ -0,0 +1,117 @@ +#!/usr/bin/env python3 +"""Independent orbit/transport invariants and actual renderer CSV consumption.""" +import importlib.util +import os +from pathlib import Path +import subprocess +import sys +import tempfile +import unittest + +import numpy as np + +ROOT = Path(__file__).resolve().parents[1] +spec = importlib.util.spec_from_file_location('track', ROOT / 'scripts/schwarzschild_camera_track.py') +track = importlib.util.module_from_spec(spec) +spec.loader.exec_module(track) + + +class TrackTests(unittest.TestCase): + def test_sampling_endpoints(self): + state = track.initial_state([8, 0, 0], [0, 0, 0]) + for duration, count in ((0., 1), (.29, 30), (.295, 30)): + tau, states, stopped = track.integrate(state, duration, 100) + self.assertEqual(len(tau), count) + self.assertLessEqual(tau[-1], duration) + self.assertIsNone(stopped) + + def test_radial_infall_through_horizon(self): + # E=1 infall from rest at infinity: dr/dtau=-sqrt(2/r). + r = 8. + w = np.sqrt(2/r) + ut = (1 + w + w*w) / (1+w) + state = track.initial_state([r, 0, 0], [-w/ut, 0, 0]) + tau, states, stopped = track.integrate(state, 20, 20, stop_radius=.01) + expected_stop = 2/(3*np.sqrt(2)) * (r**1.5 - .01**1.5) + self.assertAlmostEqual(stopped, expected_stop, delta=2e-8) + radii = np.linalg.norm(states[:, 1:4], axis=1) + np.testing.assert_allclose(radii, (r**1.5 - 1.5*np.sqrt(2)*tau)**(2/3), atol=2e-8, rtol=2e-8) + self.assertLess(radii[-1], 2) + for s in states: + g, _ = track.metric_connection(s[1:4]) + e = s[4:].reshape(4, 4) + np.testing.assert_allclose(e @ g @ e.T, track.ETA, atol=2e-8) + self.assertAlmostEqual(-(g @ e[0])[0], 1, delta=2e-8) + + def test_circular_orbit_and_transport_convergence(self): + r = 8. + omega = r**-1.5 + state = track.initial_state([r, 0, 0], [0, r*omega, 0]) + # Choose a pure Schwarzschild radial leg, transformed to KS time. + # Projecting a zero-KS-time radial seed would produce a different leg. + e = state[4:].reshape(4, 4) + root = np.sqrt(1-2/r) + e[1] = [-2/r/root, -root, 0, 0] + e[2] = [0, 0, 0, 1] + e[3] = [e[0, 2]/root, 0, e[0, 0]*root, 0] + duration = 2*np.pi/omega*np.sqrt(1-3/r) + results = [] + for tol in (1e-6, 1e-10): + tau, states, stopped = track.integrate(state, duration, 2, rtol=tol, atol=tol*.01) + self.assertIsNone(stopped) + angle = omega*tau/np.sqrt(1-3/r) + expected = r*np.column_stack([np.cos(angle), np.sin(angle), np.zeros_like(angle)]) + results.append(np.max(np.abs(states[:, 1:4]-expected))) + self.assertLess(results[1], 2e-7) + self.assertLess(results[1], results[0]/100) + s = states[-1] + e = s[4:].reshape(4, 4) + g, _ = track.metric_connection(s[1:4]) + np.testing.assert_allclose(e @ g @ e.T, track.ETA, atol=2e-8) + # Analytic parallel transport of initially inward radial e1 on circular orbit. + # In Schwarzschild components e1^r=-sqrt(1-2/r) cos(omega*tau). + n = s[1:4]/np.linalg.norm(s[1:4]) + self.assertAlmostEqual(n @ e[1, 1:], -np.sqrt(1-2/r)*np.cos(omega*tau[-1]), delta=2e-8) + + def test_initial_tetrad_and_invalid_velocity(self): + state = track.initial_state([2, 0, 0], [-.5, .1, 0], 40, 30, 17) + explicit = track.initial_state([2, 0, 0], [-.5, .1, 0], tetrad=state[4:].reshape(4, 4)) + np.testing.assert_array_equal(state, explicit) + with self.assertRaises(ValueError): + track.initial_state([2, 0, 0], [0, 0, 0]) + with self.assertRaises(ValueError): + track.initial_state([8, 0, 0], [0, 0, 0], tetrad=np.eye(4)) + + def test_csv_movie(self): + binary = ROOT / 'build/Release/schwarzschild_sky' + if not binary.exists(): + self.fail('build Schwarzschild Release renderer before running this test') + with tempfile.TemporaryDirectory() as directory: + directory = Path(directory) + csv = directory / 'camera.csv' + command = [sys.executable, str(ROOT/'scripts/schwarzschild_camera_track.py'), + '--output', str(csv), '--position', '8', '0', '0', + '--look-ra-deg', '0', '--look-dec-deg', '0', + '--fps', '10', '--duration', '.21', '--t0', '7'] + subprocess.run(command, check=True, capture_output=True, text=True) + rows = np.loadtxt(csv, delimiter=',', skiprows=1) + self.assertEqual(rows.shape, (3, 21)) + np.testing.assert_allclose(rows[:, 1], [0, .1, .2]) + result = subprocess.run([str(binary), '--observer-track', str(csv), + '--movie-track-samples', '--frames-dir', str(directory), + '--catalog', str(ROOT/'assets/sky_grid_5deg.csv'), '--width', '16', + '--height', '16', '--coarse-cell-pixels', '8', '--refine-max-level', '0'], + env=dict(os.environ, OMP_NUM_THREADS='2'), capture_output=True, text=True) + self.assertEqual(result.returncode, 0, result.stderr) + self.assertEqual(len(list(directory.glob('frame_*.png'))), 3) + for png in directory.glob('frame_*.png'): + self.assertEqual(png.read_bytes()[:8], b'\x89PNG\r\n\x1a\n') + # Error must happen before writing an output file. + csv.unlink() + bad = subprocess.run(command + ['--velocity', '2', '0', '0'], capture_output=True) + self.assertNotEqual(bad.returncode, 0) + self.assertFalse(csv.exists()) + + +if __name__ == '__main__': + unittest.main()