Compare commits
3
Commits
94149d75e4
...
b69c9cfd45
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
b69c9cfd45 | ||
|
|
8d011b7ec5 | ||
|
|
6ae223a642 |
No files matched your search
@@ -85,6 +85,7 @@ catalog 内部数据保留 `(direction, temperature, amplitude)`,而非 RGB。
|
||||
|
||||
## Git 提交消息
|
||||
|
||||
- Git 提交消息必须使用英文,包括标题和正文。
|
||||
- 提交消息必须以 category 开头,格式为 `Category: 简洁说明`。
|
||||
- 常用 category 包括 `Doc:`、`Fix:`、`Feat:`、`Makefile:` 等。
|
||||
- 对不属于修复或功能的普通小更新,使用所涉及的模块作为 category,例如 `Catalog:`、`Geodesic:`、`Spacetime:` 或 `Output:`。
|
||||
@@ -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
|
||||
|
||||
@@ -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
|
||||
```
|
||||
|
||||
[](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
|
||||
[](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
@@ -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
|
||||
```
|
||||
|
||||
[](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
|
||||
[](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 |
@@ -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:
|
||||
|
||||
|
||||
@@ -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
|
||||
@@ -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 调度,不把本征时冒充坐标时间。
|
||||
@@ -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
@@ -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. */
|
||||
|
||||
+138
-38
@@ -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) ||
|
||||
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;
|
||||
@@ -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
@@ -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) {
|
||||
|
||||
@@ -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);
|
||||
|
||||
+85
-75
@@ -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;
|
||||
return 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;
|
||||
}
|
||||
|
||||
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);
|
||||
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;
|
||||
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
@@ -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
|
||||
@@ -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)
|
||||
@@ -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;
|
||||
}
|
||||
@@ -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
@@ -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);
|
||||
|
||||
@@ -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()
|
||||
@@ -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
|
||||
|
||||
|
||||
Reference in new issue
Block a user