#include "geodesic.h" #include "frame.h" #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; const GeodesicTraceConfig trace = {.coordinate_time_step = 0.1, .max_steps = 4096, .threshold = {.kind = THRESHOLD_LOG_ENERGY_GROWTH, .value = 8.0, .policy_version = 3}}; int result = 1; if (spacetime_create_schwarzschild_ks(&spacetime, 1.0, 256.0) || 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]) || 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)) > 1e-12 || fabs(oriented_observer.coordinate_position[2] + 20.0) > 1e-12 || fabs(oriented_observer.tetrad[1][1]) > 1e-12 || 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) goto done; const RayEndpoint central = geodesic_trace_past( &spacetime, &observer, (double[]){1.0, 0.0, 0.0}, &trace); const RayEndpoint inside_shadow = geodesic_trace_past( &spacetime, &observer, (double[]){cos(0.10), sin(0.10), 0.0}, &trace); const RayEndpoint outside_shadow = geodesic_trace_past( &spacetime, &observer, (double[]){cos(0.30), sin(0.30), 0.0}, &trace); if (central.outcome != RAY_OUTCOME_DARK || inside_shadow.outcome != RAY_OUTCOME_DARK || outside_shadow.outcome != RAY_OUTCOME_ESCAPED) { fprintf(stderr, "Schwarzschild KS shadow regression failed (center=%d, inside=%d, " "outside=%d)\n", central.outcome, inside_shadow.outcome, outside_shadow.outcome); goto done; } /* The dark threshold must also be checked on the final accepted step when * that step lands exactly on the slab's left boundary. */ { ObserverState inner; if (camera_at(&spacetime, 3.0, 180.0, 0.0, &inner)) goto done; const GeodesicTraceConfig last_step = { .coordinate_time_step = 0.125, .max_steps = 1, .threshold = {.kind = THRESHOLD_LOG_ENERGY_GROWTH, .value = 0.01, .policy_version = 3}}; const RayEndpoint endpoint = geodesic_trace_past( &spacetime, &inner, (double[]){1.0, 0.0, 0.0}, &last_step); if (endpoint.outcome != RAY_OUTCOME_DARK || endpoint.reason != RAY_REASON_REDSHIFT_LIMIT || !(endpoint.threshold_value >= 0.01)) { fprintf(stderr, "last-step dark threshold regression failed (outcome=%d reason=%d " "value=%.12g)\n", endpoint.outcome, endpoint.reason, endpoint.threshold_value); goto done; } } /* A budget-exhausted ray is UNRESOLVED (retryable), keeps its last trusted * state, and resolves when resumed from that state. */ { ObserverCamera camera = {.position = {30,0,0}, .velocity = {-0.99999999,0,0}, .look_ra_deg = 0}; ObserverState boosted; MetricData m; if (spacetime_eval(&spacetime, 0, camera.position, &m) || observer_from_coordinate_camera(&m, &camera, &boosted, NULL)) goto done; GeodesicRayState initial; if (geodesic_initialize_past_ray_metric(&m, &boosted, (double[]){1,0,0}, &initial) || initial.log_alpha_p0 <= 8) goto done; GeodesicTraceConfig disabled = trace; disabled.threshold.kind = THRESHOLD_DISABLED; RayEndpoint enabled = geodesic_trace_past(&spacetime, &boosted, (double[]){1,0,0}, &trace); RayEndpoint reference = geodesic_trace_past(&spacetime, &boosted, (double[]){1,0,0}, &disabled); if (enabled.outcome != RAY_OUTCOME_ESCAPED || reference.outcome != RAY_OUTCOME_ESCAPED || fabs(enabled.frequency_ratio/reference.frequency_ratio-1) > 1e-10) { fputs("initial high-energy false-dark regression failed\n", stderr); goto done; } } const GeodesicTraceConfig tiny = { .coordinate_time_step = 0.1, .max_steps = 30, .threshold = {.kind = THRESHOLD_LOG_ENERGY_GROWTH, .value = 8.0, .policy_version = 3}}; const RayEndpoint unresolved = geodesic_trace_past( &spacetime, &observer, (double[]){cos(0.30), sin(0.30), 0.0}, &tiny); if (unresolved.outcome != RAY_OUTCOME_UNRESOLVED || unresolved.reason != RAY_REASON_BUDGET_EXHAUSTED || unresolved.end_id != SPACETIME_END_NONE) { fputs("budget-exhausted ray classification regression failed\n", stderr); goto done; } const GeodesicRayState continuation = { .coordinate_time = unresolved.stop_coordinate_time, .x = {unresolved.final_x[0], unresolved.final_x[1], unresolved.final_x[2]}, .Pi = {unresolved.final_Pi[0], unresolved.final_Pi[1], unresolved.final_Pi[2]}, .log_alpha_p0 = unresolved.final_log_alpha_p0, .log_alpha_p0_0 = unresolved.final_log_alpha_p0_0, .steps = unresolved.accepted_steps}; GeodesicTraceConfig more = tiny; more.max_steps = 8192; const RayEndpoint resumed = geodesic_trace_past_from_state(&spacetime, &continuation, &more); if (resumed.outcome != RAY_OUTCOME_ESCAPED) { fprintf(stderr, "resumed ray classification regression failed (outcome=%d)\n", (int)resumed.outcome); goto done; } /* A coarse field covering the shadow must genuinely refine: its initial * capture/escape-discontinuous triangles are a separate trigger from the * smooth direction-error criterion. */ FrameLensMesh mesh = {0}; const RefinementConfig refinement = {.max_level = 1, .angle_absolute_rad = 1e-5, .angle_relative = 1e-5, .jacobian_minimum = 1e-3, .min_edge_pixels = 1.0, .min_area_pixels2 = 1.0}; if (frame_lens_mesh_build_coarse(&mesh, 48, 48, 24, 40.0) || frame_lens_mesh_trace(&mesh, &spacetime, &observer, &trace) || frame_lens_mesh_refine(&mesh, &spacetime, &observer, &trace, &refinement) || mesh.vertex_count <= 9 || mesh.triangle_count <= 8) { fputs("Schwarzschild adaptive-refinement regression failed\n", stderr); frame_lens_mesh_destroy(&mesh); goto done; } frame_lens_mesh_destroy(&mesh); result = 0; done: spacetime_destroy(&spacetime); return result; }