Compare commits

...
3 Commits
Author SHA1 Message Date
wyj b69c9cfd45 Feat: generate freely falling Schwarzschild camera tracks
Integrate timelike geodesics and Fermi-Walker tetrads in ingoing Kerr-Schild coordinates, sampled at a configurable proper-time cadence.

Add movie-track-samples to preserve CSV events as frames. Document usage and singularity guards, and cover analytic orbits, transport convergence, sampling, and CSV rendering.
2026-09-06 05:56:51 -04:00
wyj 8d011b7ec5 Doc: require English commit messages 2026-09-06 04:01:38 -04:00
wyj 6ae223a642 Feat: support general single-frame cameras and add a near-horizon example 2026-09-06 04:00:21 -04:00
21 changed files with 1218 additions and 169 deletions

No files matched your search

+1
View File
@@ -85,6 +85,7 @@ catalog 内部数据保留 `(direction, temperature, amplitude)`,而非 RGB。
## Git 提交消息
- Git 提交消息必须使用英文,包括标题和正文。
- 提交消息必须以 category 开头,格式为 `Category: 简洁说明`。
- 常用 category 包括 `Doc:`、`Fix:`、`Feat:`、`Makefile:` 等。
- 对不属于修复或功能的普通小更新,使用所涉及的模块作为 category,例如 `Catalog:`、`Geodesic:`、`Spacetime:` 或 `Output:`。
+11 -1
View File
@@ -41,6 +41,7 @@ SCHWARZSCHILD_TEST_TARGET := $(BUILD_DIR)/test_schwarzschild
OBSERVER_TRACK_TEST_TARGET := $(BUILD_DIR)/test_observer_track
CATALOG_PREFETCH_TEST_TARGET := $(BUILD_DIR)/test_catalog_prefetch
HIP_PSF_TEST_TARGET := $(BUILD_DIR)/test_hip_psf
CAMERA_TEST_TARGETS := $(BUILD_DIR)/test_observer_minkowski $(BUILD_DIR)/test_observer_schwarzschild
.PHONY: all backend clean run test hip-psf-test minkowski schwarzschild FORCE
@@ -155,12 +156,21 @@ $(OBSERVER_TRACK_TEST_TARGET): tests/test_observer_track.c $(COMMON_SOURCES) | $
$(CATALOG_PREFETCH_TEST_TARGET): tests/test_catalog_prefetch.c $(CORE_MINKOWSKI_SOURCES) | $(BUILD_DIR)
$(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@
test: $(TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET)
$(BUILD_DIR)/test_observer_minkowski: tests/test_observer.c $(CORE_MINKOWSKI_SOURCES) | $(BUILD_DIR)
$(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@
$(BUILD_DIR)/test_observer_schwarzschild: tests/test_observer.c $(COMMON_SOURCES) src/spacetime_schwarzschild.c | $(BUILD_DIR)
$(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -DSPACETIME_SCHWARZSCHILD -Isrc $^ $(LDLIBS) -o $@
test: $(CAMERA_TEST_TARGETS) $(TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET)
./$(BUILD_DIR)/test_observer_minkowski
./$(BUILD_DIR)/test_observer_schwarzschild
./$(TEST_TARGET)
./$(FRAME_TEST_TARGET)
./$(SCHWARZSCHILD_TEST_TARGET)
./$(OBSERVER_TRACK_TEST_TARGET)
./$(CATALOG_PREFETCH_TEST_TARGET)
python3 tests/test_camera_cli.py $(BUILD_DIR)
clean:
rm -rf build
+85 -2
View File
@@ -17,7 +17,13 @@ Stars are individual catalog point sources with direction, temperature, and
amplitude. The renderer maps them into the camera image, accounting for
multiple images, gravitational lensing magnification, and frequency shifts,
then accumulates their sub-pixel point-spread functions (PSFs) into an HDR
image. Moving cameras are described by worldline and tetrad tracks.
image. Movie cameras are described by worldline and tetrad tracks.
Single images support independent coordinate position and look direction,
coordinate velocity, and camera roll in both backends. For example,
`--observer-position 1.75 0 0 --observer-velocity -0.5 0 0 --look-ra-deg 0 --look-dec-deg 0`
specifies an inward-moving Schwarzschild camera looking outward from inside
the horizon. See [single-frame camera parameters and complete commands](usage.md#single-frame-camera).
## Current status
@@ -119,12 +125,43 @@ to the cache radius. The [original benchmark record](benchmarks/2mass_galactic_c
preserves the command and terminal output; the command above omits the optional
HDR export.
The Schwarzschild camera points toward the hole. Adaptive refinement should
This example infers a position facing the hole from its look direction and radius. Adaptive refinement should
be configured for the desired image accuracy; it is disabled by default.
Camera controls, movie sequences, lens-map reuse, PSF settings, and HDR output
are described in [usage.md](usage.md). Both binaries provide a complete option
list with `--help`.
### Example: Looking outward just above the Schwarzschild horizon
This camera is static at `(2.1, 0, 0)` in Cartesian Kerr–Schild coordinates,
just outside the horizon at `r=2M` (`M=1`). RA=0°, Dec=0° points along +X,
radially outward. Its coordinate velocity is zero: this is a static
observer, not a freely falling one. The scene uses the bundled synthetic
stellar grid, so no survey download is needed.
```sh
mkdir -p output/imgs
./build/Release/schwarzschild_sky \
--catalog assets/sky_grid_5deg.csv \
--observer-position 2.1 0 0 \
--observer-velocity 0 0 0 \
--look-ra-deg 0 --look-dec-deg 0 \
--camera-roll-deg 0 \
--width 3840 --height 2160 --fov-deg 90 \
--coarse-cell-pixels 16 --refine-max-level 3 \
--exposure 0.01 \
--output output/imgs/schwarzschild_near_horizon_outward_R2.1_0.01_refine3.png
```
[![Synthetic stellar sky seen outward by a static observer at r=2.1M](assets/images/schwarzschild_near_horizon_outward_R2.1_0.01_refine3.png)](assets/images/schwarzschild_near_horizon_outward_R2.1_0.01_refine3.png)
*4K render with a 90° horizontal field of view, exposure `0.01`, and maximum
refinement level 3. Click to view at full resolution.*
The distant sky occupies a bounded angular region around the outward
direction, with repeated images crowded near its edge. Strong gravitational
blueshift pushes both the 3000 K and 12000 K test stars toward blue-white.
### Example: Synthetic test grid with mesh overlay
This example uses `assets/sky_grid_5deg.csv` to inspect lensing and adaptive
@@ -148,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`.
+75 -2
View File
@@ -8,7 +8,12 @@
*从半径 100 M 处的相机观察 Schwarzschild 时空中的 2MASS 银心方向星场。点击图片查看完整 4K 图像;[渲染命令见下文](#示例schwarzschild-时空中的银心方向星场)。*
恒星以独立的星表点源表示,保留方向、温度和振幅。渲染器将它们映射到相机图像中,处理多像、引力透镜放大和频移,再将保留亚像素位置的点扩散函数(PSF)累积到 HDR 图像。运动相机由世界线与四标架(tetrad)轨迹描述。
恒星以独立的星表点源表示,保留方向、温度和振幅。渲染器将它们映射到相机图像中,处理多像、引力透镜放大和频移,再将保留亚像素位置的点扩散函数(PSF)累积到 HDR 图像。电影中的运动相机由世界线与四标架(tetrad)轨迹描述。
单张模式在两个 backend 中都支持独立的位置、指向、坐标速度和 roll。例如
`--observer-position 1.75 0 0 --observer-velocity -0.5 0 0 --look-ra-deg 0 --look-dec-deg 0`
表示位于史瓦西视界内、向内运动且朝外看的相机,无需借用电影模式。
参数语义与完整命令见 [单张相机指南](usage.md#single-frame-camera)。
## 当前进展
@@ -77,7 +82,34 @@ mkdir -p output/imgs
这里使用参考图像的渲染设置,包括 `--max-cache-psf-flux 1e8` 这一预览近似:它将明亮 PSF 的翼部截断在缓存半径内。[原始 benchmark 记录](benchmarks/2mass_galactic_center_blackhole.md)保留了命令与终端输出;上述命令省略了可选的 HDR 导出。
Schwarzschild 相机朝向黑洞。自适应细分默认关闭,应根据所需图像精度配置。相机控制、图像序列、透镜映射复用、PSF 设置及 HDR 输出参见 [usage.md](usage.md)(英文)。两个可执行文件都可通过 `--help` 查看完整选项。
上述 Schwarzschild 示例省略位置,因此按指向与半径推导出朝向黑洞的相机。自适应细分默认关闭,应根据所需图像精度配置。相机控制、图像序列、透镜映射复用、PSF 设置及 HDR 输出参见 [usage.md](usage.md)(英文)。两个可执行文件都可通过 `--help` 查看完整选项。
### 示例:紧贴史瓦西视界向外看
相机位于 Cartesian Kerr–Schild 坐标 `(2.1, 0, 0)`,紧贴 `r=2M` 的视界外侧
(`M=1`)。RA=0°、Dec=0° 对应沿 +X 径向向外看;坐标速度为零,表示静态观者,
而非自由落体相机。这里使用仓库自带的合成测试星表,无需下载巡天数据。
```sh
mkdir -p output/imgs
./build/Release/schwarzschild_sky \
--catalog assets/sky_grid_5deg.csv \
--observer-position 2.1 0 0 \
--observer-velocity 0 0 0 \
--look-ra-deg 0 --look-dec-deg 0 \
--camera-roll-deg 0 \
--width 3840 --height 2160 --fov-deg 90 \
--coarse-cell-pixels 16 --refine-max-level 3 \
--exposure 0.01 \
--output output/imgs/schwarzschild_near_horizon_outward_R2.1_0.01_refine3.png
```
[![r=2.1M 处的静态观者向外观察合成恒星天空](assets/images/schwarzschild_near_horizon_outward_R2.1_0.01_refine3.png)](assets/images/schwarzschild_near_horizon_outward_R2.1_0.01_refine3.png)
*水平视场角 90°、曝光 `0.01`、最大细分级别 3 的 4K 图像。点击查看完整分辨率。*
远方天空集中在朝外方向的有限角域中,多级成像在边缘附近密集堆叠。
强烈的引力蓝移使 3000 K 和 12000 K 的测试恒星都被推向蓝白色。
### 示例:叠加网格的合成测试星表
@@ -100,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 渲染。
Binary file not shown.

After

Width:  |  Height:  |  Size: 5.1 MiB

+8
View File
@@ -110,6 +110,14 @@ Run CPU regression checks:
make test
```
The CPU suite includes both coordinate-camera builders, CLI validation,
small PNG renders, and single-frame/movie lens-map agreement. The reference
image checks require Python 3 and CFITSIO. The original Schwarzschild HDR
fixture is retained with a float32 relative tolerance of `2^-23` (zero absolute
tolerance): replacing the specialized static tetrad with metric-based
orthonormalization changes four components by one ULP. The Minkowski fixture
still requires exact equality.
The HIP regression requires an actual compatible GPU. The Make target builds
the test executable; run it separately:
+6 -3
View File
@@ -9,8 +9,11 @@ test: test-reference-images
test-reference-images: $(FLOATDIFF_SCRIPT) $(REFERENCE_DIR)/minkowski_ra1_dec1_640x360_HDR.fits $(REFERENCE_DIR)/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits
$(MAKE) SPACETIME=minkowski ENABLE_HDR=1 backend
mkdir -p $(REFERENCE_TMP_DIR)
OMP_NUM_THREADS=16 $(BUILD_DIR)/minkowski_sky --catalog assets/sky_grid_5deg.csv --output $(REFERENCE_TMP_DIR)/minkowski_ra1_dec1_640x360.png --hdr-output --width 640 --height 360 --fov-deg 30 --look-ra-deg 1 --look-dec-deg 1 --exposure 0.1 --observer-radius 30 --observer-inward-speed 0 --psf-fwhm-pixels 2.7 --psf-moffat-beta 4.5 --max-magnification 1e300 --max-cache-psf-flux 1 --psf-relative-tail 1e-8 --psf-min-y 0 --coarse-cell-pixels 16 --refine-max-level 0 --refine-angle-abs-deg 0.001 --refine-angle-rel 0.1 --refine-jacobian-min 1e-3 --refine-min-edge-pixels 0.5 --refine-min-area-pixels2 0.25 --catalog-load-workers 4
OMP_NUM_THREADS=16 $(BUILD_DIR)/minkowski_sky --catalog assets/sky_grid_5deg.csv --output $(REFERENCE_TMP_DIR)/minkowski_ra1_dec1_640x360.png --hdr-output --width 640 --height 360 --fov-deg 30 --look-ra-deg 1 --look-dec-deg 1 --exposure 0.1 --observer-radius 30 --observer-velocity 0 0 0 --psf-fwhm-pixels 2.7 --psf-moffat-beta 4.5 --max-magnification 1e300 --max-cache-psf-flux 1 --psf-relative-tail 1e-8 --psf-min-y 0 --coarse-cell-pixels 16 --refine-max-level 0 --refine-angle-abs-deg 0.001 --refine-angle-rel 0.1 --refine-jacobian-min 1e-3 --refine-min-edge-pixels 0.5 --refine-min-area-pixels2 0.25 --catalog-load-workers 4
python3 $(FLOATDIFF_SCRIPT) $(REFERENCE_DIR)/minkowski_ra1_dec1_640x360_HDR.fits $(REFERENCE_TMP_DIR)/minkowski_ra1_dec1_640x360_HDR.fits
$(MAKE) SPACETIME=schwarzschild ENABLE_HDR=1 backend
OMP_NUM_THREADS=16 $(BUILD_DIR)/schwarzschild_sky --catalog assets/sky_grid_5deg.csv --output $(REFERENCE_TMP_DIR)/schwarzschild_ra1_dec1_fov60_640x360.png --hdr-output --width 640 --height 360 --fov-deg 60 --look-ra-deg 1 --look-dec-deg 1 --exposure 0.1 --observer-radius 30 --observer-inward-speed 0 --psf-fwhm-pixels 2.7 --psf-moffat-beta 4.5 --max-magnification 1e300 --max-cache-psf-flux 1 --psf-relative-tail 1e-8 --psf-min-y 0 --coarse-cell-pixels 16 --refine-max-level 0 --refine-angle-abs-deg 0.001 --refine-angle-rel 0.1 --refine-jacobian-min 1e-3 --refine-min-edge-pixels 0.5 --refine-min-area-pixels2 0.25 --catalog-load-workers 4
python3 $(FLOATDIFF_SCRIPT) $(REFERENCE_DIR)/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits $(REFERENCE_TMP_DIR)/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits
OMP_NUM_THREADS=16 $(BUILD_DIR)/schwarzschild_sky --catalog assets/sky_grid_5deg.csv --output $(REFERENCE_TMP_DIR)/schwarzschild_ra1_dec1_fov60_640x360.png --hdr-output --width 640 --height 360 --fov-deg 60 --look-ra-deg 1 --look-dec-deg 1 --exposure 0.1 --observer-radius 30 --observer-velocity 0 0 0 --psf-fwhm-pixels 2.7 --psf-moffat-beta 4.5 --max-magnification 1e300 --max-cache-psf-flux 1 --psf-relative-tail 1e-8 --psf-min-y 0 --coarse-cell-pixels 16 --refine-max-level 0 --refine-angle-abs-deg 0.001 --refine-angle-rel 0.1 --refine-jacobian-min 1e-3 --refine-min-edge-pixels 0.5 --refine-min-area-pixels2 0.25 --catalog-load-workers 4
# The general metric-based tetrad changes four float32 samples by one ULP
# versus the legacy analytic static tetrad. Keep the original fixture and
# allow at most float32 relative rounding (2^-23), with zero abs tolerance.
python3 $(FLOATDIFF_SCRIPT) $(REFERENCE_DIR)/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits $(REFERENCE_TMP_DIR)/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits --rel-tolerance 1.1920928955078125e-7
+73 -1
View File
@@ -340,7 +340,7 @@ trajectory generator 可以采用不同物理规则:
- 人为给定加速度;
- 人为叠加 camera pointing/rotation。
renderer 本身只读取轨迹。
电影 renderer 本身只读取轨迹;单张入口可以用下述参数直接构造一个事件处的完整 tetrad。
示意:
@@ -361,6 +361,55 @@ typedef struct {
每帧根据 `t_camera` 插值得到完整 observer state。
## 10.1 单张的瞬时相机
单张与电影共享 `ObserverState` 和 ray 初始化。单张不是只指定三维位置:
由事件、坐标速度、指向与 roll 生成完整四速度和 tetrad。
目前解析单张事件取 $t=0$;此构造不要求静态时空或静止观测者。
`--observer-position X Y Z` 与 `--look-ra-deg` / `--look-dec-deg` 独立。
只有位置时取朝原点的坐标方向;只有指向时令 $\mathbf x=-R\mathbf d$。
`--observer-radius R` 仅用于缺省位置的推导(默认 30),与显式位置互斥。
几何参数全部省略时保留 Minkowski 原点、Schwarzschild 半径 30 的默认相机。
原点无法推导指向,要求显式 look;由位置推导方向且在极轴上时固定 RA=0。
RA 或 Dec 只给一个时,另一个取默认值(90°、−90°)。
`--observer-velocity VX VY VZ` 表示 $v^i=dx^i/dt$,默认零。
由事件处的 $\alpha,\beta^i,\gamma_{ij}$ 计算
$$
Q=-\alpha^2+\gamma_{ij}(v^i+\beta^i)(v^j+\beta^j),\qquad
u^\mu=\frac{(1,v^i)}{\sqrt{-Q}}.
$$
要求 $Q<0$,未来指向选择 $u^t>0$。不能用欧氏 $|\mathbf v|<1$ 或
$r>2M$ 代替类时性检查;删除旧的 `--observer-inward-speed` 参数。
构造失败直接报错,不修改速度或退回静止相机。
采用标准 ICRS 坐标种子
$$
\mathbf d=(\cos\delta\cos\alpha_{\rm RA},\cos\delta\sin\alpha_{\rm RA},\sin\delta),
$$
将 $D^\mu=(0,\mathbf d)$ 投影为
$F^\mu=D^\mu+u^\mu(u_\nu D^\nu)$,归一化得到 forward。
随后把天球北向和西向的坐标种子投影、正交化为 up、right,保留标准画面定向。
因此 look 表示坐标方向在相机静止空间内的投影,不保证特定远方物体位于画面中心。
`--camera-roll-deg` 默认零,正角按
$e'_2=\cos\rho\,e_2+\sin\rho\,e_3$、
$e'_3=-\sin\rho\,e_2+\cos\rho\,e_3$ 定义。
observer 构造只接收当地 metric 和已补全的参数,不加载 slab、不分类 ray。
调用方在昂贵的 catalog/PSF 初始化前验证相机及 backend 数据域。
现有 Schwarzschild cutoff 为 $r=1.5M$;相机必须在 cutoff 外,但允许在视界内。
此功能不改变捕获 cutoff 或向过去追踪的高红移终止条件。
单张相机参数与轨迹输入、lens-map 导入互斥;导入仍跳过 metric 与 observer 初始化。
验证包括 tetrad 正交归一和 null 初始化、平直时空平移不变性与解析光行差/多普勒、
KS 视界处和视界内的向过去径向逃逸及频移,以及相同 observer 的单张/电影单帧一致性。
---
# 11. Ray 初始化
@@ -1310,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 调度,不把本征时冒充坐标时间。
+165
View File
@@ -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 <CSV> --movie-track-samples --frames-dir <DIR>; encode at the chosen fps.', file=sys.stderr)
except (ValueError, OSError, OverflowError) as exc:
p.error(str(exc))
if __name__ == '__main__':
main()
+2 -1
View File
@@ -413,7 +413,8 @@ int frame_lens_mesh_prepare_generation(FrameLensMesh *mesh,
return -1;
}
if (config->max_level == 0)
return 0;
/* Coarse vertices still need tracing when refinement is disabled. */
return (int)mesh->sample_count;
/* Every generation may batch newly inserted vertices with probes for its
* new leaves: probe positions depend only on image-plane geometry. Their
* endpoints are considered only after this complete generation finishes. */
+137 -37
View File
@@ -28,7 +28,9 @@ typedef struct {
double psf_relative_tail;
double psf_min_y;
double observer_radius;
double observer_inward_speed;
double observer_position[3], observer_velocity[3], camera_roll_deg;
int position_specified, velocity_specified, look_specified;
int radius_specified, roll_specified;
PointSpreadFunction psf;
PsfKernelCache psf_cache;
const char *catalog_path;
@@ -42,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;
@@ -83,14 +86,14 @@ static int parse_ra_deg(const char *text, double *value) {
char *end;
errno = 0;
*value = strtod(text, &end);
return errno || *end || *value < 0.0 || *value >= 360.0 ? -1 : 0;
return errno || end == text || *end || !isfinite(*value) || *value < 0.0 || *value >= 360.0 ? -1 : 0;
}
static int parse_dec_deg(const char *text, double *value) {
char *end;
errno = 0;
*value = strtod(text, &end);
return errno || *end || *value < -90.0 || *value > 90.0 ? -1 : 0;
return errno || end == text || *end || !isfinite(*value) || *value < -90.0 || *value > 90.0 ? -1 : 0;
}
static int parse_positive(const char *text, double *value) {
@@ -107,11 +110,11 @@ static int parse_nonnegative(const char *text, double *value) {
return errno || *end || !isfinite(*value) || *value < 0.0 ? -1 : 0;
}
static int parse_speed(const char *text, double *value) {
static int parse_finite(const char *text, double *value) {
char *end;
errno = 0;
*value = strtod(text, &end);
return errno || *end || *value < 0.0 || *value >= 1.0 ? -1 : 0;
return errno || end == text || *end || !isfinite(*value) ? -1 : 0;
}
static int parse_moffat_beta(const char *text, double *value) {
@@ -132,7 +135,7 @@ static int parse_finite_positive(const char *text, double *value) {
char *end;
errno = 0;
*value = strtod(text, &end);
return errno || *end || !isfinite(*value) || *value <= 0.0 ? -1 : 0;
return errno || end == text || *end || !isfinite(*value) || *value <= 0.0 ? -1 : 0;
}
static int validate_tonemapped_output_path(const char *path) {
@@ -256,14 +259,26 @@ static int parse_args(int argc, char **argv, Settings *s,
s->fov_specified = 1;
} else if (!strcmp(argv[i], "--look-ra-deg") && i + 1 < argc &&
!parse_ra_deg(argv[++i], &s->look_ra_deg)) {
s->look_specified = 1;
} else if (!strcmp(argv[i], "--look-dec-deg") && i + 1 < argc &&
!parse_dec_deg(argv[++i], &s->look_dec_deg)) {
s->look_specified = 1;
} else if (!strcmp(argv[i], "--exposure") && i + 1 < argc &&
!parse_positive(argv[++i], &s->exposure)) {
} else if (!strcmp(argv[i], "--observer-inward-speed") && i + 1 < argc &&
!parse_speed(argv[++i], &s->observer_inward_speed)) {
} else if ((!strcmp(argv[i], "--observer-position") ||
!strcmp(argv[i], "--observer-velocity")) && i + 3 < argc) {
const int position = !strcmp(argv[i], "--observer-position");
double *v = position ? s->observer_position : s->observer_velocity;
for (int component = 0; component < 3; ++component)
if (parse_finite(argv[++i], &v[component])) return -1;
if (position) s->position_specified = 1;
else s->velocity_specified = 1;
} else if (!strcmp(argv[i], "--camera-roll-deg") && i + 1 < argc &&
!parse_finite(argv[++i], &s->camera_roll_deg)) {
s->roll_specified = 1;
} else if (!strcmp(argv[i], "--observer-radius") && i + 1 < argc &&
!parse_positive(argv[++i], &s->observer_radius)) {
!parse_finite_positive(argv[++i], &s->observer_radius)) {
s->radius_specified = 1;
} else if (!strcmp(argv[i], "--psf-fwhm-pixels") && i + 1 < argc &&
!parse_positive(argv[++i], &s->psf.fwhm_pixels)) {
} else if (!strcmp(argv[i], "--psf-moffat-beta") && i + 1 < argc &&
@@ -285,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)
@@ -339,8 +356,18 @@ static void print_help(const char *program) {
" --look-ra-deg D ICRS look direction right ascension in degrees (default: 90)\n"
" --look-dec-deg D ICRS look direction declination in degrees (default: -90)\n"
" --exposure E Linear exposure multiplier (default: 1e-3)\n"
" --observer-radius R Observer radius in Schwarzschild units (default: 30)\n"
" --observer-inward-speed V Inward observer speed as a fraction of c (default: 0)\n"
" --observer-position X Y Z Coordinate position; alone implies looking at the origin\n"
" --observer-radius R Infer position = -R * look direction (default R: 30); conflicts with position\n"
" --observer-velocity VX VY VZ Coordinate dx/dt, dy/dt, dz/dt (default: 0 0 0); must be timelike\n"
" --camera-roll-deg ANGLE Rotate up toward right about forward (default: 0)\n"
" Explicit look alone implies position = -R * look direction.\n"
" Look is projected into the moving camera rest space.\n"
#ifdef SPACETIME_SCHWARZSCHILD
" Default position: (0,0,30); look RA=90, Dec=-90.\n"
#else
" Default position: (0,0,0); look RA=90, Dec=-90.\n"
#endif
" Single-frame camera options conflict with track/map input.\n"
"\nPSF and catalog splatting:\n"
" --psf-fwhm-pixels N Moffat PSF FWHM in pixels (default: 2.7)\n"
" --psf-moffat-beta N Moffat PSF beta, greater than 1 (default: 4.5)\n"
@@ -361,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"
@@ -483,15 +511,83 @@ static GeodesicTraceConfig trace_config(void) {
#endif
}
static int default_observer(const Settings *s, ObserverState *observer) {
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) {
fputs("--observer-position and --observer-radius are mutually exclusive.\n", stderr);
return -1;
}
if (camera_specified && (s->observer_track_path || s->frames_dir ||
s->lens_map_input_path)) {
fputs("Single-frame camera options cannot be combined with movie/observer-track or --lens-map-input.\n", stderr);
return -1;
}
if (s->position_specified && !s->look_specified) {
const double *x = s->observer_position;
const double radius = hypot(hypot(x[0], x[1]), x[2]);
if (!isfinite(radius) || radius == 0.0) {
fputs("Cannot infer a look direction from this position; specify --look-ra-deg and/or --look-dec-deg.\n", stderr);
return -1;
}
const double degrees = 180.0 / 3.14159265358979323846;
s->look_ra_deg = (x[0] == 0.0 && x[1] == 0.0)
? 0.0 : atan2(-x[1], -x[0]) * degrees;
if (s->look_ra_deg < 0.0) s->look_ra_deg += 360.0;
if (s->look_ra_deg >= 360.0) s->look_ra_deg = 0.0;
s->look_dec_deg = atan2(-x[2], hypot(x[0], x[1])) * degrees;
}
int infer_position = s->look_specified || s->radius_specified;
#ifdef SPACETIME_SCHWARZSCHILD
return observer_inward_schwarzschild_ks_look_at(
1.0, s->observer_radius, s->look_ra_deg, s->look_dec_deg,
s->observer_inward_speed, observer);
#else
*observer = observer_fixed_at_origin_look_at(s->look_ra_deg, s->look_dec_deg);
return 0;
infer_position = 1;
#endif
if (!s->position_specified && infer_position) {
const ObserverState pointing = observer_fixed_at_origin_look_at(
s->look_ra_deg, s->look_dec_deg);
for (int i = 0; i < 3; ++i)
s->observer_position[i] = -s->observer_radius * pointing.tetrad[1][i + 1];
}
return 0;
}
static int build_observer(const Settings *s, const SpacetimeSource *spacetime,
ObserverState *observer) {
ObserverCamera camera = {.look_ra_deg = s->look_ra_deg,
.look_dec_deg = s->look_dec_deg,
.roll_deg = s->camera_roll_deg};
for (int i = 0; i < 3; ++i) {
camera.position[i] = s->observer_position[i];
camera.velocity[i] = s->observer_velocity[i];
}
if (spacetime_classify(spacetime, camera.coordinate_time, camera.position) ==
SPACETIME_RAY_CAPTURED) {
fputs("Camera position is inside the backend capture cutoff or invalid.\n", stderr);
return -1;
}
MetricData metric;
if (spacetime_eval(spacetime, camera.coordinate_time, camera.position, &metric)) {
fputs("Could not evaluate metric at the camera event.\n", stderr);
return -1;
}
double q;
const ObserverBuildResult result = observer_from_coordinate_camera(
&metric, &camera, observer, &q);
if (result != OBSERVER_BUILD_OK) {
fprintf(stderr, "Camera construction failed (%s): position=(%.17g, %.17g, %.17g), "
"coordinate velocity=(%.17g, %.17g, %.17g), Q=%.17g.\n",
result == OBSERVER_BUILD_NON_TIMELIKE ? "velocity is not timelike" :
result == OBSERVER_BUILD_INVALID_INPUT ? "invalid input/metric" :
"invalid tetrad",
camera.position[0], camera.position[1], camera.position[2],
camera.velocity[0], camera.velocity[1], camera.velocity[2], q);
return -1;
}
return 0;
}
static int render_observer_frame(const Settings *s, StarCatalog *catalog,
@@ -599,14 +695,6 @@ static int render_observer_frame(const Settings *s, StarCatalog *catalog,
return result;
}
static int render_frame(const Settings *s, StarCatalog *catalog,
const SpacetimeSource *spacetime) {
ObserverState observer;
return default_observer(s, &observer) ? -1
: render_observer_frame(s, catalog, spacetime,
&observer, s->output_path);
}
static int frame_output_path(char path[PATH_MAX], const Settings *s,
size_t frame_id) {
#ifdef ENABLE_PNG
@@ -719,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) ||
(s->movie_track_samples ? movie_init_track_samples(&movie, &track) :
movie_init(&movie, &track, s->movie_start_time, s->movie_duration,
s->movie_fps) ||
s->movie_fps)) ||
movie_build_coarse_meshes(&movie, s->width, s->height,
s->coarse_cell_pixels, s->horizontal_fov_deg))
goto done;
@@ -895,7 +984,8 @@ int main(int argc, char **argv) {
"Usage: %s [--catalog PATH | --all-sky-catalog DIR] [--output PATH] [--width N] [--height "
"N] [--fov-deg D] [--look-ra-deg D] [--look-dec-deg D] "
"[--lens-map-input FILE | --lens-map-output FILE] "
"[--exposure E] [--observer-radius R] [--observer-inward-speed V] "
"[--exposure E] [--observer-radius R | --observer-position X Y Z] "
"[--observer-velocity VX VY VZ] [--camera-roll-deg ANGLE] "
"[--psf-fwhm-pixels N] [--psf-moffat-beta N] "
"[--max-magnification M] [--max-cache-psf-flux F] "
"[--psf-relative-tail R] [--psf-min-y Y] "
@@ -910,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;
@@ -946,16 +1036,32 @@ int main(int argc, char **argv) {
return 2;
}
#endif
if (resolve_camera(&settings)) return 2;
SpacetimeSource spacetime = {0};
ObserverState observer;
if (settings.lens_map_input_path == NULL) {
if (spacetime_create_default(&spacetime)) {
fputs("Could not create spacetime source\n", stderr);
return 1;
}
if (settings.frames_dir == NULL &&
build_observer(&settings, &spacetime, &observer)) {
spacetime_destroy(&spacetime);
return 2;
}
}
StarCatalog catalog = {0};
if (settings.all_sky_catalog_path != NULL) {
if (catalog_load_all_sky(&catalog, settings.all_sky_catalog_path)) {
perror(settings.all_sky_catalog_path);
spacetime_destroy(&spacetime);
return 1;
}
} else if (catalog_load_csv(&catalog, settings.catalog_path)) {
if (catalog_write_octant_grid(settings.catalog_path) ||
catalog_load_csv(&catalog, settings.catalog_path)) {
perror(settings.catalog_path);
spacetime_destroy(&spacetime);
return 1;
}
fprintf(stderr, "Created test catalog: %s\n", settings.catalog_path);
@@ -970,16 +1076,10 @@ int main(int argc, char **argv) {
psf_kernel_cache_destroy(&settings.psf_cache);
return result == 0 ? 0 : 1;
}
SpacetimeSource spacetime = {0};
if (spacetime_create_default(&spacetime)) {
fputs("Could not create spacetime source\n", stderr);
catalog_destroy(&catalog);
psf_kernel_cache_destroy(&settings.psf_cache);
return 1;
}
int result = settings.frames_dir != NULL
? render_movie(&settings, &catalog, &spacetime)
: render_frame(&settings, &catalog, &spacetime);
: render_observer_frame(&settings, &catalog, &spacetime,
&observer, settings.output_path);
spacetime_destroy(&spacetime);
catalog_destroy(&catalog);
psf_kernel_cache_destroy(&settings.psf_cache);
+22
View File
@@ -3,6 +3,28 @@
#include <math.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
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) {
+2
View File
@@ -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);
+86 -76
View File
@@ -1,5 +1,6 @@
#include "observer.h"
#include <float.h>
#include <math.h>
#include <stddef.h>
@@ -32,84 +33,93 @@ ObserverState observer_fixed_at_origin_look_at(double ra_deg, double dec_deg) {
{0.0, right[0], right[1], right[2]}}};
}
static int schwarzschild_look_direction(double ra_deg, double dec_deg,
double direction[3], double up[3],
double right[3]) {
if (!isfinite(ra_deg) || !isfinite(dec_deg) || ra_deg < 0.0 ||
ra_deg >= 360.0 || dec_deg < -90.0 || dec_deg > 90.0)
return -1;
const double ra = ra_deg * pi / 180.0;
const double dec = dec_deg * pi / 180.0;
const double cos_ra = cos(ra), sin_ra = sin(ra);
const double cos_dec = cos(dec), sin_dec = sin(dec);
direction[0] = cos_dec * cos_ra;
direction[1] = cos_dec * sin_ra;
direction[2] = sin_dec;
up[0] = -sin_dec * cos_ra;
up[1] = -sin_dec * sin_ra;
up[2] = cos_dec;
right[0] = sin_ra;
right[1] = -cos_ra;
right[2] = 0.0;
/* Use the 3+1 form directly, including the shift in every four-vector. */
static double inner(const MetricData *m, const double a[4], const double b[4]) {
double value = -m->alpha * m->alpha * a[0] * b[0];
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j)
value += m->gamma[i][j] * (a[i + 1] + m->beta[i] * a[0]) *
(b[j + 1] + m->beta[j] * b[0]);
return value;
}
static int valid_metric(const MetricData *m) {
if (!isfinite(m->alpha) || m->alpha <= 0.0) return 0;
double l[3][3] = {{0}};
for (int i = 0; i < 3; ++i) {
if (!isfinite(m->beta[i])) return 0;
for (int j = 0; j <= i; ++j) {
double value = m->gamma[i][j];
const double transposed = m->gamma[j][i];
if (!isfinite(value) || !isfinite(transposed) ||
fabs(value - transposed) > 32 * DBL_EPSILON *
fmax(fabs(value), fabs(transposed)))
return 0;
}
int observer_static_schwarzschild_ks_look_at(double mass, double radius,
double look_ra_deg,
double look_dec_deg,
ObserverState *out) {
double direction[3], up[3], right[3];
if (out == NULL || mass <= 0.0 || radius <= 2.0 * mass ||
schwarzschild_look_direction(look_ra_deg, look_dec_deg, direction, up,
right))
return -1;
const double f = 2.0 * mass / radius;
const double normalization = sqrt(1.0 - f);
*out = (ObserverState){
.coordinate_time = 0.0,
.coordinate_position = {-radius * direction[0], -radius * direction[1],
-radius * direction[2]},
/* e_(0) is the static four-velocity. e_(1) points inward, toward the
* origin; its time component makes the tetrad orthonormal in the KS
* metric. */
.tetrad = {{1.0 / normalization, 0.0, 0.0, 0.0},
{-f / normalization, normalization * direction[0],
normalization * direction[1], normalization * direction[2]},
{0.0, up[0], up[1], up[2]},
{0.0, right[0], right[1], right[2]}}};
return 0;
}
int observer_static_schwarzschild_ks(double mass, double radius,
ObserverState *out) {
return observer_static_schwarzschild_ks_look_at(mass, radius, 180.0, 0.0,
out);
}
int observer_inward_schwarzschild_ks_look_at(double mass, double radius,
double look_ra_deg,
double look_dec_deg,
double inward_speed,
ObserverState *out) {
ObserverState static_observer;
if (inward_speed < 0.0 || inward_speed >= 1.0 ||
observer_static_schwarzschild_ks_look_at(
mass, radius, look_ra_deg, look_dec_deg, &static_observer))
return -1;
const double gamma = 1.0 / sqrt(1.0 - inward_speed * inward_speed);
*out = static_observer;
for (int mu = 0; mu < 4; ++mu) {
const double e0 = static_observer.tetrad[0][mu];
const double forward = static_observer.tetrad[1][mu];
out->tetrad[0][mu] = gamma * (e0 + inward_speed * forward);
out->tetrad[1][mu] = gamma * (inward_speed * e0 + forward);
for (int k = 0; k < j; ++k) value -= l[i][k] * l[j][k];
if (i == j) {
if (!isfinite(value) || value <= 0.0) return 0;
l[i][j] = sqrt(value);
} else l[i][j] = value / l[j][j];
}
return 0;
}
return 1;
}
int observer_inward_schwarzschild_ks(double mass, double radius,
double inward_speed, ObserverState *out) {
return observer_inward_schwarzschild_ks_look_at(mass, radius, 180.0, 0.0,
inward_speed, out);
ObserverBuildResult observer_from_coordinate_camera(
const MetricData *metric, const ObserverCamera *camera,
ObserverState *out, double *q) {
if (q) *q = NAN;
if (!metric || !camera || !out || !valid_metric(metric) ||
!isfinite(camera->coordinate_time) ||
!isfinite(camera->look_ra_deg) || camera->look_ra_deg < 0.0 ||
camera->look_ra_deg >= 360.0 || !isfinite(camera->look_dec_deg) ||
fabs(camera->look_dec_deg) > 90.0 || !isfinite(camera->roll_deg))
return OBSERVER_BUILD_INVALID_INPUT;
ObserverState state = observer_fixed_at_origin_look_at(
camera->look_ra_deg, camera->look_dec_deg);
state.coordinate_time = camera->coordinate_time;
for (int i = 0; i < 3; ++i) {
if (!isfinite(camera->position[i]) || !isfinite(camera->velocity[i]))
return OBSERVER_BUILD_INVALID_INPUT;
state.coordinate_position[i] = camera->position[i];
state.tetrad[0][i + 1] = camera->velocity[i];
}
const double norm = inner(metric, state.tetrad[0], state.tetrad[0]);
if (q) *q = norm;
if (!isfinite(norm) || norm >= 0.0) return OBSERVER_BUILD_NON_TIMELIKE;
for (int mu = 0; mu < 4; ++mu) state.tetrad[0][mu] /= sqrt(-norm);
for (int a = 1; a < 4; ++a) {
/* Modified Gram-Schmidt with reorthogonalization in the observer rest
* space; the coordinate forward seed is given first priority. */
for (int pass = 0; pass < 2; ++pass)
for (int b = 0; b < a; ++b) {
const double projection = inner(metric, state.tetrad[a], state.tetrad[b]);
for (int mu = 0; mu < 4; ++mu)
state.tetrad[a][mu] -= (b == 0 ? -projection : projection) *
state.tetrad[b][mu];
}
const double length2 = inner(metric, state.tetrad[a], state.tetrad[a]);
if (!isfinite(length2) || length2 <= 0.0)
return OBSERVER_BUILD_INVALID_TETRAD;
for (int mu = 0; mu < 4; ++mu) state.tetrad[a][mu] /= sqrt(length2);
}
const double roll = remainder(camera->roll_deg, 360.0) * pi / 180.0;
for (int mu = 0; mu < 4; ++mu) {
const double up = state.tetrad[2][mu], right = state.tetrad[3][mu];
state.tetrad[2][mu] = cos(roll) * up + sin(roll) * right;
state.tetrad[3][mu] = -sin(roll) * up + cos(roll) * right;
}
for (int a = 0; a < 4; ++a) {
for (int mu = 0; mu < 4; ++mu)
if (!isfinite(state.tetrad[a][mu])) return OBSERVER_BUILD_INVALID_TETRAD;
for (int b = 0; b <= a; ++b) {
const double error = inner(metric, state.tetrad[a], state.tetrad[b]) -
(a == b ? (a == 0 ? -1.0 : 1.0) : 0.0);
if (!isfinite(error) || fabs(error) > 1e-8)
return OBSERVER_BUILD_INVALID_TETRAD;
}
}
if (state.tetrad[0][0] <= 0.0) return OBSERVER_BUILD_INVALID_TETRAD;
*out = state;
return OBSERVER_BUILD_OK;
}
+25 -19
View File
@@ -1,6 +1,8 @@
#ifndef OBSERVER_H
#define OBSERVER_H
#include "spacetime.h"
typedef struct {
double coordinate_time;
double coordinate_position[3];
@@ -13,24 +15,28 @@ ObserverState observer_fixed_at_origin(void);
* direction. The local spatial axes are (forward, celestial north,
* celestial west), so a north-up image has decreasing RA to the right. */
ObserverState observer_fixed_at_origin_look_at(double ra_deg, double dec_deg);
/* Static camera at Cartesian Kerr--Schild position -radius * look_direction,
* directed toward the Schwarzschild black hole at the origin. look_direction
* uses the same ICRS-style RA/Dec convention as the flat-space camera. */
int observer_static_schwarzschild_ks_look_at(double mass, double radius,
double look_ra_deg,
double look_dec_deg,
ObserverState *out);
/* Compatibility shortcut for the +X camera directed toward the origin. */
int observer_static_schwarzschild_ks(double mass, double radius,
ObserverState *out);
/* As above, with a local radial inward boost relative to the static observer. */
int observer_inward_schwarzschild_ks_look_at(double mass, double radius,
double look_ra_deg,
double look_dec_deg,
double inward_speed,
ObserverState *out);
/* Compatibility shortcut for the +X camera directed toward the origin. */
int observer_inward_schwarzschild_ks(double mass, double radius,
double inward_speed, ObserverState *out);
/* Fully resolved instantaneous camera; velocity is dx^i/dt, not a local boost.
* Look angles specify a coordinate direction projected into the rest space of
* the resulting four-velocity. No static reference observer is required. */
typedef struct {
double coordinate_time;
double position[3], velocity[3];
double look_ra_deg, look_dec_deg, roll_deg;
} ObserverCamera;
typedef enum {
OBSERVER_BUILD_OK = 0,
OBSERVER_BUILD_INVALID_INPUT,
OBSERVER_BUILD_NON_TIMELIKE,
OBSERVER_BUILD_INVALID_TETRAD
} ObserverBuildResult;
/* Metric belongs to the camera event. q, if non-NULL, receives g((1,v),(1,v)).
* Positive roll: up' = cos(roll) up + sin(roll) right,
* right' = -sin(roll) up + cos(roll) right.
* out is only written on success. */
ObserverBuildResult observer_from_coordinate_camera(
const MetricData *metric, const ObserverCamera *camera,
ObserverState *out, double *q);
#endif
+140
View File
@@ -0,0 +1,140 @@
#!/usr/bin/env python3
"""Exercise camera defaults/errors and single-frame/movie agreement (CPU PNG builds)."""
import os
from pathlib import Path
import struct
import subprocess
import sys
import tempfile
import zlib
BUILD = Path(sys.argv[1] if len(sys.argv) > 1 else 'build/Release').resolve()
ENV = dict(os.environ, OMP_NUM_THREADS='4')
def run(binary, *args, ok=True):
result = subprocess.run([str(binary), *map(str, args)], env=ENV,
capture_output=True, text=True)
if (result.returncode == 0) != ok:
raise AssertionError(f'{binary.name} {args}: {result.returncode}\n{result.stderr}')
return result
def image_payload(path):
data = path.read_bytes()
assert data[:8] == b'\x89PNG\r\n\x1a\n'
offset, compressed = 8, bytearray()
while offset < len(data):
count, kind = struct.unpack_from('>I4s', data, offset)
payload = data[offset + 8:offset + 8 + count]
if kind == b'IHDR':
assert struct.unpack_from('>II', payload) == (64, 48)
if kind == b'IDAT':
compressed.extend(payload)
offset += count + 12
raw = zlib.decompress(compressed)
assert any(raw), f'empty image: {path}'
return raw
def map_vertices(path):
data = path.read_bytes()
assert data[:8] == b'GRLENS\x01\x00'
assert struct.unpack_from('<Q', data, 32)[0] == 1
vertices, triangles = struct.unpack_from('<QQ', data, 64)
offset = 80
values = []
for _ in range(vertices):
values.append(struct.unpack_from('<9dI', data, offset))
offset += 76
return values, data[offset:offset + triangles * 28]
with tempfile.TemporaryDirectory(prefix='gr-camera-cli-') as directory:
tmp = Path(directory)
for backend in ('minkowski', 'schwarzschild'):
binary = BUILD / f'{backend}_sky'
help_text = run(binary, '--help').stdout
for option in ('--observer-position', '--observer-velocity', '--camera-roll-deg'):
assert option in help_text
assert '--observer-inward-speed' not in help_text
common = ['--catalog', 'assets/sky_grid_5deg.csv', '--width', 64,
'--height', 48, '--fov-deg', 80, '--exposure', 0.1,
'--coarse-cell-pixels', 8, '--refine-max-level', 0, '--psf-relative-tail', 1e-4]
def render(name, *options):
path = tmp / f'{backend}_{name}.png'
run(binary, *common, '--output', path, *options)
return image_payload(path)
# Equivalent independently specified and inferred camera geometry.
inferred = render('position', '--observer-position', -30, 0, 0)
explicit = render('explicit', '--observer-position', -30, 0, 0,
'--look-ra-deg', 0, '--look-dec-deg', 0)
angled = render('angle', '--look-ra-deg', 0, '--look-dec-deg', 0)
assert inferred == explicit == angled
pole = render('pole', '--observer-position', 0, 0, 30)
assert pole == render('pole_explicit', '--observer-position', 0, 0, 30,
'--look-ra-deg', 0, '--look-dec-deg', -90)
default = render('default')
pos = (0, 0, 0) if backend == 'minkowski' else (0, 0, 30)
assert default == render('default_explicit', '--observer-position', *pos,
'--look-ra-deg', 90, '--look-dec-deg', -90)
assert render('radius', '--observer-radius', 40) == render(
'radius_explicit', '--observer-radius', 40, '--look-ra-deg', 90,
'--look-dec-deg', -90)
for partial, value, ra, dec in [('--look-ra-deg', 37, 37, -90),
('--look-dec-deg', -23, 90, -23)]:
assert render('partial', partial, value) == render(
'complete', '--look-ra-deg', ra, '--look-dec-deg', dec)
errors = [
(['--observer-position', 1, 2], None),
(['--observer-position', 1, 2, 'nan'], None),
(['--observer-velocity', 0, 0, 'inf'], None),
(['--look-ra-deg', 'nan'], None),
(['--look-dec-deg', 'inf'], None),
(['--observer-radius', 'nan'], None),
(['--observer-radius', 0], None),
(['--camera-roll-deg', 'nan'], None),
(['--observer-position', 0, 0, 0], 'Cannot infer'),
(['--observer-position', 3, 4, 5, '--observer-radius', 30], 'mutually exclusive'),
(['--observer-velocity', 10, 0, 0], 'not timelike'),
(['--observer-inward-speed', 0], None),
(['--observer-track', 'missing.csv', '--observer-velocity', 0, 0, 0], 'cannot be combined'),
(['--frames-dir', tmp, '--look-ra-deg', 0], 'cannot be combined'),
(['--lens-map-input', 'missing.grlens', '--camera-roll-deg', 0], 'cannot be combined'),
]
if backend == 'schwarzschild':
errors += [(['--observer-position', 1.5, 0, 0, '--observer-velocity', -0.5, 0, 0], 'capture cutoff'),
(['--observer-position', 1.75, 0, 0], 'not timelike')]
render('inside', '--observer-position', 1.75, 0, 0,
'--observer-velocity', -0.5, 0, 0, '--look-ra-deg', 0, '--look-dec-deg', 0)
for options, message in errors:
missing_catalog = tmp / 'should_not_be_created.csv'
result = run(binary, '--catalog', missing_catalog, *options, ok=False)
if message:
assert message in result.stderr, result.stderr
assert not missing_catalog.exists(), result.stderr
assert 'PSF cache ready' not in result.stderr
track = tmp / f'{backend}.csv'
run(BUILD / f'test_observer_{backend}', track)
single_map, movie_map = tmp / 'single.grlens', tmp / 'movie.grlens'
single = render('moving', '--observer-position', 3, -4, 5,
'--observer-velocity', 0.2, -0.1, 0.3,
'--look-ra-deg', 37, '--look-dec-deg', -23,
'--camera-roll-deg', 19, '--lens-map-output', single_map)
run(binary, *common, '--observer-track', track, '--frames-dir', tmp,
'--frames-prefix', backend, '--duration', 0, '--fps', 1,
'--lens-map-output', movie_map)
movie = image_payload(tmp / f'{backend}_000000.png')
assert single == movie, f'{backend}: single/movie PNG mismatch'
a, ta = map_vertices(single_map)
b, tb = map_vertices(movie_map)
assert len(a) == len(b) and ta == tb
max_error = 0
for x, y in zip(a, b):
assert x[-1] == y[-1], 'ray classification mismatch'
max_error = max(max_error, *(abs(v - w) for v, w in zip(x[:-1], y[:-1])))
assert max_error < 1e-9, max_error
# A map import must still work without evaluating a camera/metric.
assert single == render('import', '--lens-map-input', single_map)
print(f'{backend}: CLI checks passed; single/movie PNG identical, map max error {max_error:.3g}', flush=True)
+157
View File
@@ -0,0 +1,157 @@
#include "geodesic.h"
#include "observer_track.h"
#include <math.h>
#include <stdio.h>
#include <string.h>
#define CHECK(condition) do { if (!(condition)) { \
fprintf(stderr, "observer regression failed at line %d: %s\n", __LINE__, #condition); \
return 1; } } while (0)
/* Independent covariant four-metric contraction in extended precision. */
static long double dot(const MetricData *m, const double a[4], const double b[4]) {
long double g[4][4] = {{0}};
g[0][0] = -(long double)m->alpha * m->alpha;
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j) {
g[i + 1][j + 1] = m->gamma[i][j];
g[0][0] += (long double)m->gamma[i][j] * m->beta[i] * m->beta[j];
g[0][i + 1] += (long double)m->gamma[i][j] * m->beta[j];
g[i + 1][0] = g[0][i + 1];
}
long double result = 0;
for (int i = 0; i < 4; ++i)
for (int j = 0; j < 4; ++j) result += g[i][j] * a[i] * b[j];
return result;
}
static int check_state(const MetricData *m, const ObserverCamera *c,
const ObserverState *o) {
CHECK(o->coordinate_time == c->coordinate_time);
CHECK(o->tetrad[0][0] > 0);
for (int i = 0; i < 3; ++i) {
CHECK(o->coordinate_position[i] == c->position[i]);
CHECK(fabs(o->tetrad[0][i + 1] / o->tetrad[0][0] - c->velocity[i]) < 1e-12);
}
for (int a = 0; a < 4; ++a)
for (int b = 0; b < 4; ++b)
CHECK(fabsl(dot(m, o->tetrad[a], o->tetrad[b]) -
(a == b ? (a == 0 ? -1 : 1) : 0)) < 1e-11L);
const double n[3] = {0.36, 0.48, 0.8};
double k[4];
for (int mu = 0; mu < 4; ++mu) {
k[mu] = o->tetrad[0][mu];
for (int a = 0; a < 3; ++a) k[mu] -= n[a] * o->tetrad[a + 1][mu];
}
CHECK(k[0] > 0 && fabsl(dot(m, k, k)) < 1e-11L);
CHECK(fabsl(dot(m, k, o->tetrad[0]) + 1) < 1e-11L);
return 0;
}
int main(int argc, char **argv) {
SpacetimeSource source = {0};
CHECK(spacetime_create_default(&source) == 0);
ObserverCamera camera = {.position = {3, -4, 5}, .velocity = {0.2, -0.1, 0.3},
.look_ra_deg = 37, .look_dec_deg = -23, .roll_deg = 19};
MetricData metric;
ObserverState state;
CHECK(spacetime_eval(&source, 0, camera.position, &metric) == 0);
CHECK(observer_from_coordinate_camera(&metric, &camera, &state, NULL) == OBSERVER_BUILD_OK);
CHECK(check_state(&metric, &camera, &state) == 0);
if (argc == 2) {
ObserverSample sample = {.coordinate_time = 0, .proper_time = 0};
memcpy(sample.coordinate_position, state.coordinate_position, sizeof sample.coordinate_position);
memcpy(sample.tetrad, state.tetrad, sizeof sample.tetrad);
ObserverTrack track = {.samples = &sample, .count = 1};
CHECK(observer_track_write_csv(&track, argv[1]) == 0);
}
/* A rescaled/shifted coordinate system can have timelike |dx/dt| > 1.
* The builder must use the metric, not impose a Euclidean speed limit. */
const MetricData shifted_metric = {.alpha = 2, .beta = {0.1, -0.2, 0.3},
.gamma = {{1, 0.1, 0}, {0.1, 1.2, 0.1}, {0, 0.1, 0.9}}};
const ObserverCamera fast_coordinate = {.velocity = {1.2, 0, 0},
.look_ra_deg = 123, .look_dec_deg = 45, .roll_deg = -31};
CHECK(observer_from_coordinate_camera(&shifted_metric, &fast_coordinate,
&state, NULL) == OBSERVER_BUILD_OK);
CHECK(check_state(&shifted_metric, &fast_coordinate, &state) == 0);
MetricData rounded_metric = shifted_metric;
rounded_metric.gamma[1][0] = nextafter(rounded_metric.gamma[1][0], INFINITY);
CHECK(observer_from_coordinate_camera(&rounded_metric, &fast_coordinate,
&state, NULL) == OBSERVER_BUILD_OK);
CHECK(check_state(&rounded_metric, &fast_coordinate, &state) == 0);
MetricData invalid_metric = shifted_metric;
invalid_metric.gamma[2][2] = -1;
CHECK(observer_from_coordinate_camera(&invalid_metric, &fast_coordinate,
&state, NULL) == OBSERVER_BUILD_INVALID_INPUT);
const double ras[] = {0, 37, 90, 180, 359.9};
const double decs[] = {-90, -23, 0, 45, 90};
for (int a = 0; a < 5; ++a)
for (int b = 0; b < 5; ++b) {
camera.look_ra_deg = ras[a]; camera.look_dec_deg = decs[b];
CHECK(observer_from_coordinate_camera(&metric, &camera, &state, NULL) == OBSERVER_BUILD_OK);
CHECK(check_state(&metric, &camera, &state) == 0);
}
camera.look_ra_deg = 0; camera.look_dec_deg = 0; camera.roll_deg = 0;
ObserverState unrolled;
CHECK(observer_from_coordinate_camera(&metric, &camera, &unrolled, NULL) == OBSERVER_BUILD_OK);
camera.roll_deg = 90;
CHECK(observer_from_coordinate_camera(&metric, &camera, &state, NULL) == OBSERVER_BUILD_OK);
for (int mu = 0; mu < 4; ++mu) {
CHECK(fabs(state.tetrad[2][mu] - unrolled.tetrad[3][mu]) < 1e-12);
CHECK(fabs(state.tetrad[3][mu] + unrolled.tetrad[2][mu]) < 1e-12);
}
camera.velocity[0] = NAN;
CHECK(observer_from_coordinate_camera(&metric, &camera, &state, NULL) == OBSERVER_BUILD_INVALID_INPUT);
camera.velocity[0] = 10;
CHECK(observer_from_coordinate_camera(&metric, &camera, &state, NULL) == OBSERVER_BUILD_NON_TIMELIKE);
camera.velocity[0] = 0;
camera.roll_deg = INFINITY;
CHECK(observer_from_coordinate_camera(&metric, &camera, &state, NULL) == OBSERVER_BUILD_INVALID_INPUT);
#ifdef SPACETIME_SCHWARZSCHILD
/* Ingoing radial light seen from the horizon and its interior must still
* trace backwards to the external sky, rather than be classified captured. */
const GeodesicTraceConfig trace = {.coordinate_time_step = 0.05,
.max_steps = 8192, .capture_log_alpha_p0 = 8};
for (int i = 0; i < 3; ++i) {
camera = (ObserverCamera){.position = {2.25 - 0.25 * i, 0, 0},
.velocity = {-0.5, 0, 0}};
CHECK(spacetime_eval(&source, 0, camera.position, &metric) == 0);
CHECK(observer_from_coordinate_camera(&metric, &camera, &state, NULL) == OBSERVER_BUILD_OK);
CHECK(check_state(&metric, &camera, &state) == 0);
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. */
const double energy = state.tetrad[0][0] - state.tetrad[1][0];
CHECK(fabs(ray.frequency_ratio - sqrt(1 + 2.0 / 256) / energy) < 2e-6);
memset(camera.velocity, 0, sizeof camera.velocity);
if (i > 0)
CHECK(observer_from_coordinate_camera(&metric, &camera, &state, NULL) == OBSERVER_BUILD_NON_TIMELIKE);
}
#else
/* For transverse velocity along Y, projected +X remains F=(0,1,0,0).
* Independently, k=(gamma,-1,gamma*v,0) gives aberration and Doppler. */
camera = (ObserverCamera){.velocity = {0, 0.6, 0}};
CHECK(spacetime_eval(&source, 0, camera.position, &metric) == 0);
CHECK(observer_from_coordinate_camera(&metric, &camera, &state, NULL) == OBSERVER_BUILD_OK);
const GeodesicTraceConfig trace = {.coordinate_time_step = 1, .max_steps = 2048};
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] - 0.8) < 1e-12);
CHECK(fabs(ray.n_infinity[1] + 0.6) < 1e-12);
CHECK(fabs(ray.frequency_ratio - 0.8) < 1e-12);
camera.position[0] = 25; camera.position[1] = -30; camera.position[2] = 10;
CHECK(observer_from_coordinate_camera(&metric, &camera, &state, NULL) == OBSERVER_BUILD_OK);
const RayEndpoint shifted = geodesic_trace_past(&source, &state, (double[]){1, 0, 0}, &trace);
CHECK(shifted.status == ray.status && fabs(shifted.frequency_ratio - ray.frequency_ratio) < 1e-12);
for (int i = 0; i < 3; ++i) CHECK(fabs(shifted.n_infinity[i] - ray.n_infinity[i]) < 1e-12);
camera.velocity[1] = 1;
CHECK(observer_from_coordinate_camera(&metric, &camera, &state, NULL) == OBSERVER_BUILD_NON_TIMELIKE);
#endif
spacetime_destroy(&source);
puts("coordinate-camera regression passed");
return 0;
}
+17
View File
@@ -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};
+14 -16
View File
@@ -4,12 +4,22 @@
#include <math.h>
#include <stdio.h>
static int camera_at(const SpacetimeSource *source, double radius,
double ra, double dec, ObserverState *out) {
ObserverCamera camera = {.look_ra_deg = ra, .look_dec_deg = dec};
const ObserverState pointing = observer_fixed_at_origin_look_at(ra, dec);
for (int i = 0; i < 3; ++i)
camera.position[i] = -radius * pointing.tetrad[1][i + 1];
MetricData metric;
return spacetime_eval(source, 0.0, camera.position, &metric) ||
observer_from_coordinate_camera(&metric, &camera, out, NULL);
}
int main(void) {
SpacetimeSource spacetime = {0};
MetricData metric;
ObserverState observer;
ObserverState oriented_observer;
ObserverState inward_observer;
const GeodesicTraceConfig trace = {.coordinate_time_step = 0.1,
.max_steps = 4096,
.capture_log_alpha_p0 = 8.0};
@@ -18,8 +28,8 @@ int main(void) {
spacetime_eval(&spacetime, 0.0, (double[]){2.0, 0.0, 0.0}, &metric) ||
!isfinite(metric.alpha) || !isfinite(metric.gamma[0][0]) ||
!isfinite(metric.K[0][0]) ||
observer_static_schwarzschild_ks(1.0, 30.0, &observer) ||
observer_static_schwarzschild_ks_look_at(1.0, 40.0, 270.0, 30.0,
camera_at(&spacetime, 30.0, 180.0, 0.0, &observer) ||
camera_at(&spacetime, 40.0, 270.0, 30.0,
&oriented_observer) ||
fabs(oriented_observer.coordinate_position[0]) > 1e-12 ||
fabs(oriented_observer.coordinate_position[1] - 20.0 * sqrt(3.0)) >
@@ -30,19 +40,7 @@ int main(void) {
fabs(oriented_observer.tetrad[1][2] + sqrt(0.95) * sqrt(3.0) / 2.0) >
1e-12 ||
fabs(oriented_observer.tetrad[1][3] - 0.5 * sqrt(0.95)) >
1e-12 ||
observer_inward_schwarzschild_ks(1.0, 30.0, 0.5,
&inward_observer) ||
fabs(inward_observer.tetrad[0][0] -
(2.0 / sqrt(3.0)) * (observer.tetrad[0][0] +
0.5 * observer.tetrad[1][0])) >
1e-12 ||
fabs(inward_observer.tetrad[1][1] -
(2.0 / sqrt(3.0)) * (0.5 * observer.tetrad[0][1] +
observer.tetrad[1][1])) >
1e-12 ||
!observer_inward_schwarzschild_ks(1.0, 30.0, 1.0,
&inward_observer))
1e-12)
goto done;
const RayEndpoint central = geodesic_trace_past(
&spacetime, &observer, (double[]){1.0, 0.0, 0.0}, &trace);
+117
View File
@@ -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()
+75 -11
View File
@@ -14,19 +14,83 @@ processing are documented in [assets/2mass/README.md](assets/2mass/README.md).
synthetic test catalog; `--exposure 1e15` is a starting point, to be adjusted
for the field and camera.
## Schwarzschild camera
## Single-frame camera
`schwarzschild_sky` uses an analytic Schwarzschild metric in Cartesian ingoing
Kerr–Schild coordinates, with mass `M=1`. The default camera is static at
coordinate radius `30`. `--look-ra-deg` and `--look-dec-deg` set the direction
from the camera to the hole; the camera is placed on the opposite side of the
origin and points radially inward. `--observer-radius R` requires `R > 2`.
`--observer-inward-speed V` applies a local inward boost relative to the static
observer, with `0 <= V < 1` and default `0`.
Both backends accept the same instantaneous camera parameters. Position and
velocity use the backend's coordinates; velocity means `dx/dt, dy/dt, dz/dt`,
not a local physical speed. The analytic single-frame event is at `t=0`.
The analytic setup uses escape radius `256` and capture radius `1.5`, inside
the horizon at `r=2`. These are current demonstration settings; they do not
establish termination criteria for future numerical-relativity data.
| Option | Meaning / default |
| --- | --- |
| `--observer-position X Y Z` | Coordinate position; if look is omitted, point toward the origin |
| `--look-ra-deg RA`, `--look-dec-deg DEC` | Coordinate look direction; missing angle defaults to RA=90°, Dec=-90° |
| `--observer-radius R` | Positive radius used only to infer position, default 30; conflicts with explicit position |
| `--observer-velocity VX VY VZ` | Coordinate velocity, default `0 0 0`; must be future-timelike in the local metric |
| `--camera-roll-deg ANGLE` | Rotate the camera up axis toward its right axis, default 0° |
When position and look are both supplied, they are independent. Look alone
implies position `-R * look_direction`, in either backend. Radius alone uses
the default look. With no position/look/radius options, Minkowski keeps the
origin camera and Schwarzschild keeps its radius-30 camera on +Z; both look
along -Z. Velocity and roll alone do not alter these defaults.
An explicit origin position requires a look angle. Position-derived directions
on the Z axis use RA=0° to fix the pole orientation deterministically.
The axes are right-handed ICRS: X=RA 0°/Dec 0°, Y=RA 90°/Dec 0°,
Z=Dec +90°. The coordinate look vector `(cos(dec) cos(ra), cos(dec) sin(ra),
sin(dec))` is projected into the moving camera's rest space. The coordinate
north and west seeds are then orthogonalized there to form up and right.
Consequently the look angles specify a projected coordinate orientation;
they do not lock a lensed sky object to the image center. Positive roll uses
`up' = cos(roll) up + sin(roll) right`,
`right' = -sin(roll) up + cos(roll) right`.
No static observer is needed. The metric determines whether `(1,VX,VY,VZ)`
is timelike; there is no Euclidean `|V| < 1` check in curved coordinates.
Invalid velocities are rejected before loading catalogs or building the PSF
cache, with the velocity and its metric norm in the diagnostic.
`--observer-inward-speed` has been removed.
Single-frame camera options cannot be combined with `--observer-track`,
`--frames-dir`, or `--lens-map-input`.
Schwarzschild uses Cartesian ingoing Kerr–Schild coordinates with `M=1`.
Cameras at and inside the horizon `r=2` are allowed with a valid timelike
coordinate velocity. The current backend excludes camera positions at or
inside its capture cutoff `r=1.5`; its finite escape radius is `256`.
These remain analytic demonstration settings, not criteria for future NR data.
Zero coordinate velocity at or inside the horizon is not timelike and is rejected.
The following complete examples use the bundled synthetic catalog:
```sh
make -j PSF_BACKEND=cpu SPACETIME=minkowski backend
mkdir -p output/imgs
./build/Release/minkowski_sky --catalog assets/sky_grid_5deg.csv \
--observer-position 3 -4 5 --observer-velocity 0.2 -0.1 0.3 \
--look-ra-deg 37 --look-dec-deg -23 --camera-roll-deg 19 \
--width 64 --height 48 --fov-deg 80 --exposure 0.1 \
--coarse-cell-pixels 8 --refine-max-level 0 --psf-direct \
--output output/imgs/minkowski_moving.png
make -j PSF_BACKEND=cpu SPACETIME=schwarzschild backend
./build/Release/schwarzschild_sky --catalog assets/sky_grid_5deg.csv \
--observer-position 3 -4 5 --observer-velocity 0.2 -0.1 0.3 \
--look-ra-deg 37 --look-dec-deg -23 --camera-roll-deg 19 \
--width 64 --height 48 --fov-deg 80 --exposure 0.1 \
--coarse-cell-pixels 8 --refine-max-level 0 --psf-direct \
--output output/imgs/schwarzschild_off_axis.png
./build/Release/schwarzschild_sky --catalog assets/sky_grid_5deg.csv \
--observer-position 1.75 0 0 --observer-velocity -0.5 0 0 \
--look-ra-deg 0 --look-dec-deg 0 \
--width 64 --height 48 --fov-deg 80 --exposure 0.1 \
--coarse-cell-pixels 8 --refine-max-level 0 --psf-direct \
--output output/imgs/schwarzschild_inside_horizon.png
```
The last example points outward from a camera moving inward inside the horizon.
It needs no CSV trajectory or movie wrapper. These small images are camera
checks; increase resolution and refinement for production renders.
## Movie image sequences