diff --git a/Makefile b/Makefile index 853e7fe..3fa7e4f 100644 --- a/Makefile +++ b/Makefile @@ -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 diff --git a/README.md b/README.md index 4631adc..1f3b7ed 100644 --- a/README.md +++ b/README.md @@ -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 diff --git a/README.zh-CN.md b/README.zh-CN.md index 9cdafc2..83161b2 100644 --- a/README.zh-CN.md +++ b/README.zh-CN.md @@ -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 的测试恒星都被推向蓝白色。 ### 示例：叠加网格的合成测试星表 diff --git a/assets/images/schwarzschild_near_horizon_outward_R2.1_0.01_refine3.png b/assets/images/schwarzschild_near_horizon_outward_R2.1_0.01_refine3.png new file mode 100644 index 0000000..de7d6af Binary files /dev/null and b/assets/images/schwarzschild_near_horizon_outward_R2.1_0.01_refine3.png differ diff --git a/build.md b/build.md index 7e1d69d..153dbcb 100644 --- a/build.md +++ b/build.md @@ -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: diff --git a/mk/reference_images.mk b/mk/reference_images.mk index ee6be30..9aa408d 100644 --- a/mk/reference_images.mk +++ b/mk/reference_images.mk @@ -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 diff --git a/nr_spacetime_movie_renderer_design.md b/nr_spacetime_movie_renderer_design.md index f9b315c..561b562 100644 --- a/nr_spacetime_movie_renderer_design.md +++ b/nr_spacetime_movie_renderer_design.md @@ -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 初始化 diff --git a/src/frame.c b/src/frame.c index 3ceb8a1..1612d71 100644 --- a/src/frame.c +++ b/src/frame.c @@ -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. */ diff --git a/src/main.c b/src/main.c index aa35b57..f0172b4 100644 --- a/src/main.c +++ b/src/main.c @@ -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; @@ -83,14 +85,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 +109,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 +134,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 +258,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 && @@ -339,8 +353,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" @@ -483,15 +507,78 @@ static GeodesicTraceConfig trace_config(void) { #endif } -static int default_observer(const Settings *s, ObserverState *observer) { +static int resolve_camera(Settings *s) { + 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 +686,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 @@ -895,7 +974,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] " @@ -946,16 +1026,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 +1066,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); diff --git a/src/observer.c b/src/observer.c index c4cbd11..c6e03cd 100644 --- a/src/observer.c +++ b/src/observer.c @@ -1,5 +1,6 @@ #include "observer.h" +#include #include #include @@ -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; } diff --git a/src/observer.h b/src/observer.h index 6b7cb15..ea26598 100644 --- a/src/observer.h +++ b/src/observer.h @@ -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 diff --git a/tests/test_camera_cli.py b/tests/test_camera_cli.py new file mode 100644 index 0000000..102fd0e --- /dev/null +++ b/tests/test_camera_cli.py @@ -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(' +#include +#include + +#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; +} diff --git a/tests/test_schwarzschild.c b/tests/test_schwarzschild.c index 3a7109d..c8dea39 100644 --- a/tests/test_schwarzschild.c +++ b/tests/test_schwarzschild.c @@ -4,12 +4,22 @@ #include #include +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); diff --git a/usage.md b/usage.md index cf52e23..49ef529 100644 --- a/usage.md +++ b/usage.md @@ -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