wyj
789549ed2f
Use scaled long-double quadratic arithmetic without explicit FMA. Validate entry candidates against backend geometry and localize uncertain entries along the original exterior trajectory. Preserve conservative miss semantics and propagate concrete entry failures. Add production-sample and numerical regression coverage.
661 lines
28 KiB
C
661 lines
28 KiB
C
#include "asymptotic.h"
|
|
#include "asymptotic_schwarzschild.h"
|
|
#include "geodesic.h"
|
|
#include "observer.h"
|
|
#include "spacetime.h"
|
|
|
|
#include <math.h>
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
|
|
static int failures = 0;
|
|
#define CHECK(condition, message) \
|
|
do { \
|
|
if (!(condition)) { \
|
|
fprintf(stderr, "FAIL %s:%d: %s\n", __FILE__, __LINE__, message); \
|
|
++failures; \
|
|
} \
|
|
} while (0)
|
|
|
|
static double angle_between(const double a[3], const double b[3]) {
|
|
const double dot = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
|
|
const double cx = a[1] * b[2] - a[2] * b[1];
|
|
const double cy = a[2] * b[0] - a[0] * b[2];
|
|
const double cz = a[0] * b[1] - a[1] * b[0];
|
|
return atan2(sqrt(cx * cx + cy * cy + cz * cz), dot);
|
|
}
|
|
|
|
static void test_round_trip(void) {
|
|
SpacetimeSource source = {0};
|
|
CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0,
|
|
"create schwarzschild");
|
|
SpacetimeAsymptoticEnd end;
|
|
CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end descriptor");
|
|
SchwarzschildCanonical in = {.end_id = 0,
|
|
.t = 0.0,
|
|
.rho = 256.0,
|
|
.rhat = {1.0, 0.0, 0.0},
|
|
.Lhat = {0.0, 1.0, 0.0},
|
|
.beta = 5.0,
|
|
.energy = 1.0,
|
|
.radial_sign = 1};
|
|
double x[3], Pi[3], log_alpha_p0;
|
|
CHECK(asymptotic_schwarzschild_state_from_canonical(&end, &in, x, Pi,
|
|
&log_alpha_p0) == 0,
|
|
"state from canonical");
|
|
MetricData metric;
|
|
CHECK(spacetime_eval(&source, in.t, x, &metric) == 0, "metric");
|
|
SchwarzschildCanonical out;
|
|
CHECK(asymptotic_schwarzschild_canonical_from_state(
|
|
&end, &metric, in.t, x, Pi, log_alpha_p0, &out) == 0,
|
|
"canonical from state");
|
|
CHECK(fabs(out.beta - in.beta) < 1e-13, "beta round trip");
|
|
CHECK(fabs(out.energy - in.energy) < 1e-13, "energy round trip");
|
|
CHECK(out.radial_sign == in.radial_sign, "radial sign round trip");
|
|
const double axis = angle_between(out.rhat, in.rhat);
|
|
CHECK(axis < 1e-13, "position direction round trip");
|
|
spacetime_destroy(&source);
|
|
}
|
|
|
|
static void test_finish_matches_integration(void) {
|
|
SpacetimeSource near = {0}, far = {0};
|
|
CHECK(spacetime_create_schwarzschild_ks(&near, 1.0, 256.0) == 0,
|
|
"create near");
|
|
CHECK(spacetime_create_schwarzschild_ks(&far, 1.0, 1.0e5) == 0,
|
|
"create far");
|
|
SpacetimeAsymptoticEnd end;
|
|
CHECK(spacetime_asymptotic_end(&near, 0, &end) == 0, "near end");
|
|
|
|
const double betas[] = {0.0, 0.5, 4.0, 10.0, 30.0, 100.0, 250.0};
|
|
const int beta_count = (int)(sizeof betas / sizeof betas[0]);
|
|
for (int k = 0; k < beta_count; ++k) {
|
|
SchwarzschildCanonical canonical = {.end_id = 0,
|
|
.t = 0.0,
|
|
.rho = 256.0,
|
|
.rhat = {0.8, 0.6, 0.0},
|
|
.Lhat = {0.0, 0.0, 1.0},
|
|
.beta = betas[k],
|
|
.energy = 1.0,
|
|
.radial_sign = 1};
|
|
double x[3], Pi[3], log_alpha_p0;
|
|
CHECK(asymptotic_schwarzschild_state_from_canonical(
|
|
&end, &canonical, x, Pi, &log_alpha_p0) == 0,
|
|
"finish state build");
|
|
|
|
double n_analytic[3], freq_analytic;
|
|
CHECK(asymptotic_schwarzschild_finish(&end, &canonical, n_analytic,
|
|
&freq_analytic) == 0,
|
|
"analytic finish");
|
|
|
|
GeodesicRayState state = {.coordinate_time = 0.0,
|
|
.x = {x[0], x[1], x[2]},
|
|
.Pi = {Pi[0], Pi[1], Pi[2]},
|
|
.log_alpha_p0 = log_alpha_p0,
|
|
.steps = 0};
|
|
const GeodesicTraceConfig config = {.coordinate_time_step = 5.0,
|
|
.max_steps = 100000};
|
|
MetricSlab *slab = NULL;
|
|
CHECK(spacetime_load_slab(&far, 0.0, -1.0e6, &slab) == 0, "far slab");
|
|
RayEndpoint endpoint = {.frequency_ratio = 0, .magnification = 1.0,
|
|
.end_id = SPACETIME_END_NONE,
|
|
.outcome = RAY_OUTCOME_INCOMPLETE};
|
|
const GeodesicAdvanceResult result =
|
|
geodesic_advance_past_ray(slab, &state, -1.0e6, &config, &endpoint);
|
|
spacetime_free_slab(slab);
|
|
CHECK(result == GEODESIC_ADVANCE_TERMINATED &&
|
|
endpoint.outcome == RAY_OUTCOME_ESCAPED,
|
|
"far integration escapes");
|
|
/* Pipeline check only: the far integration at step 5 and escape radius
|
|
* 1e5 has its own O(1e-5..1e-3) error. Quantitative accuracy is checked
|
|
* against the high-precision reference constants below. */
|
|
const double angle_error =
|
|
angle_between(n_analytic, endpoint.n_infinity);
|
|
CHECK(angle_error < 1e-2, "finish direction matches far integration");
|
|
CHECK(fabs(freq_analytic - endpoint.frequency_ratio) /
|
|
freq_analytic < 1e-2,
|
|
"finish frequency matches far integration");
|
|
(void)angle_error;
|
|
}
|
|
spacetime_destroy(&near);
|
|
spacetime_destroy(&far);
|
|
}
|
|
|
|
/* Independent quadrature of the KS coordinate-time transfer for a camera
|
|
* outside the worldtube, used to check the analytic primitive. */
|
|
static double simpson(const double a, const double b, int panels,
|
|
double (*f)(double, const void *), const void *ctx) {
|
|
if (panels < 2)
|
|
panels = 2;
|
|
if (panels % 2)
|
|
++panels;
|
|
const double h = (b - a) / panels;
|
|
double sum = f(a, ctx) + f(b, ctx);
|
|
for (int i = 1; i < panels; ++i)
|
|
sum += (i % 2 ? 4.0 : 2.0) * f(a + i * h, ctx);
|
|
return sum * h / 3.0;
|
|
}
|
|
|
|
typedef struct {
|
|
double beta;
|
|
} TransferContext;
|
|
|
|
static double transfer_dt(double r, const void *context) {
|
|
const TransferContext *c = context;
|
|
const double Q = 1.0 - c->beta * c->beta * (1.0 - 2.0 / r) / (r * r);
|
|
return 1.0 / ((1.0 - 2.0 / r) * sqrt(Q)) + 2.0 / (r - 2.0);
|
|
}
|
|
|
|
static double transfer_dphi(double r, const void *context) {
|
|
const TransferContext *c = context;
|
|
const double Q = 1.0 - c->beta * c->beta * (1.0 - 2.0 / r) / (r * r);
|
|
return c->beta / (r * r * sqrt(Q));
|
|
}
|
|
|
|
static void test_preroute_entry(void) {
|
|
SpacetimeSource source = {0};
|
|
CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0,
|
|
"create schwarzschild");
|
|
const ObserverCamera camera = {.look_ra_deg = 0.0, .look_dec_deg = 0.0};
|
|
ObserverCamera positioned = camera;
|
|
positioned.position[0] = 500.0;
|
|
positioned.look_ra_deg = 180.0;
|
|
positioned.look_dec_deg = 0.0;
|
|
const double direction[3] = {cos(0.3), sin(0.3), 0.0};
|
|
MetricData metric;
|
|
CHECK(spacetime_eval(&source, 0.0, positioned.position, &metric) == 0,
|
|
"camera metric");
|
|
ObserverState observer;
|
|
CHECK(observer_from_coordinate_camera(&metric, &positioned, &observer,
|
|
NULL) == OBSERVER_BUILD_OK,
|
|
"camera observer");
|
|
|
|
MetricSlab *camera_slab = NULL;
|
|
CHECK(spacetime_load_slab(&source, 0.0, -1.0, &camera_slab) == 0,
|
|
"camera slab");
|
|
GeodesicRayState camera_state;
|
|
CHECK(geodesic_initialize_past_ray(camera_slab, &observer, direction,
|
|
&camera_state) == 0,
|
|
"camera state");
|
|
SpacetimeAsymptoticEnd end;
|
|
CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end");
|
|
SchwarzschildCanonical camera_can;
|
|
CHECK(asymptotic_schwarzschild_canonical_from_state(
|
|
&end, &metric, 0.0, camera_state.x, camera_state.Pi,
|
|
camera_state.log_alpha_p0, &camera_can) == 0,
|
|
"camera canonical");
|
|
spacetime_free_slab(camera_slab);
|
|
|
|
AsymptoticRoute route;
|
|
CHECK(asymptotic_route_camera(&source, &observer, direction, &route) ==
|
|
ASYMPTOTIC_OK &&
|
|
route.kind == ASYMPTOTIC_ROUTE_ENTRY,
|
|
"outside camera enters");
|
|
CHECK(route.entry_fallback_evaluations == 0,
|
|
"analytic entry stays on the fast path");
|
|
double value;
|
|
CHECK(asymptotic_worldtube_value(&source, route.end_id, route.activate_t,
|
|
route.x, &value) == 0 &&
|
|
fabs(value) < 1e-3,
|
|
"entry on worldtube");
|
|
|
|
MetricData entry_metric;
|
|
CHECK(spacetime_eval(&source, route.activate_t, route.x, &entry_metric) ==
|
|
0,
|
|
"entry metric");
|
|
SchwarzschildCanonical entry_can;
|
|
CHECK(asymptotic_schwarzschild_canonical_from_state(
|
|
&end, &entry_metric, route.activate_t, route.x, route.Pi,
|
|
route.log_alpha_p0, &entry_can) == 0,
|
|
"entry canonical");
|
|
CHECK(fabs(entry_can.beta - camera_can.beta) <
|
|
1e-12 * fmax(1.0, camera_can.beta),
|
|
"entry conserves impact parameter");
|
|
CHECK(fabs(entry_can.energy - camera_can.energy) < 1e-12,
|
|
"entry conserves energy");
|
|
CHECK(entry_can.radial_sign == -1, "entry is past-inward");
|
|
CHECK(route.activate_t < 0.0, "entry time is in the past");
|
|
|
|
const TransferContext context = {.beta = camera_can.beta};
|
|
const double t_analytic = -route.activate_t;
|
|
const double t_numeric =
|
|
simpson(256.0, 500.0, 20000, transfer_dt, &context);
|
|
CHECK(fabs(t_analytic - t_numeric) < 1e-9 * fmax(1.0, t_numeric),
|
|
"entry time matches quadrature");
|
|
const double dphi_numeric =
|
|
simpson(256.0, 500.0, 20000, transfer_dphi, &context);
|
|
const double dphi_entry = angle_between(camera_can.rhat, entry_can.rhat);
|
|
CHECK(fabs(dphi_entry - dphi_numeric) < 1e-9,
|
|
"entry azimuth matches quadrature");
|
|
if (fabs(t_analytic - t_numeric) >= 1e-9 * fmax(1.0, t_numeric) ||
|
|
fabs(dphi_entry - dphi_numeric) >= 1e-9)
|
|
fprintf(stderr, " beta=%.6g t_an=%.12g t_num=%.12g dphi_an=%.12g "
|
|
"dphi_num=%.12g\n",
|
|
camera_can.beta, t_analytic, t_numeric, dphi_entry,
|
|
dphi_numeric);
|
|
spacetime_destroy(&source);
|
|
}
|
|
|
|
/* An external camera's dark-threshold reference must be the L at the camera
|
|
* event, not the worldtube entry energy. Changing only the worldtube radius
|
|
* must not change the reference but may change the entry L. */
|
|
static void test_camera_reference_radius_independent(void) {
|
|
const double radii[2] = {128.0, 256.0};
|
|
double reference[2], entry[2];
|
|
for (int k = 0; k < 2; ++k) {
|
|
SpacetimeSource source = {0};
|
|
CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, radii[k]) == 0,
|
|
"create reference source");
|
|
ObserverCamera positioned = {.look_ra_deg = 180.0, .look_dec_deg = 0.0};
|
|
positioned.position[0] = 500.0;
|
|
MetricData metric;
|
|
CHECK(spacetime_eval(&source, 0.0, positioned.position, &metric) == 0,
|
|
"reference camera metric");
|
|
ObserverState observer;
|
|
CHECK(observer_from_coordinate_camera(&metric, &positioned, &observer,
|
|
NULL) == OBSERVER_BUILD_OK,
|
|
"reference camera observer");
|
|
const double direction[3] = {cos(0.05), sin(0.05), 0.0};
|
|
MetricSlab *slab = NULL;
|
|
GeodesicRayState camera_state;
|
|
CHECK(spacetime_load_slab(&source, 0.0, -1.0, &slab) == 0,
|
|
"reference camera slab");
|
|
CHECK(geodesic_initialize_past_ray(slab, &observer, direction,
|
|
&camera_state) == 0,
|
|
"reference camera state");
|
|
spacetime_free_slab(slab);
|
|
AsymptoticRoute route;
|
|
CHECK(asymptotic_route_camera(&source, &observer, direction, &route) ==
|
|
ASYMPTOTIC_OK &&
|
|
route.kind == ASYMPTOTIC_ROUTE_ENTRY,
|
|
"reference ray enters the worldtube");
|
|
CHECK(fabs(route.log_alpha_p0_camera - camera_state.log_alpha_p0) < 1e-12,
|
|
"route reference is the camera-event L");
|
|
reference[k] = route.log_alpha_p0_camera;
|
|
entry[k] = route.log_alpha_p0;
|
|
spacetime_destroy(&source);
|
|
}
|
|
CHECK(fabs(reference[0] - reference[1]) < 1e-12,
|
|
"camera reference is worldtube-radius independent");
|
|
CHECK(fabs(entry[0] - entry[1]) > 1e-6,
|
|
"entry energy depends on the worldtube radius");
|
|
}
|
|
|
|
/* High-precision (mpmath, 60 digits) reference values fixed into the ordinary
|
|
* C test: radial, complex-pair, three-real, grazing, and large-radius angle
|
|
* cases. */
|
|
static void test_phi_reference_constants(void) {
|
|
static const struct {
|
|
double rho, beta, value;
|
|
} cases[] = {
|
|
{256.0, 0.0, 0.0},
|
|
{256.0, 5.0, 0.019532484697919191145},
|
|
{256.0, 60.0, 0.23656231243306290715},
|
|
{64.0, 64.0, 1.4199914058161304301},
|
|
{256.0, 255.0, 1.4527184167466732533},
|
|
{1.0e6, 1.0, 1.0000000000001666664e-6},
|
|
{300.0, 3.0, 0.010000165840750676787},
|
|
{100.0, 5.3, 0.053024471018799209953},
|
|
};
|
|
for (size_t i = 0; i < sizeof cases / sizeof cases[0]; ++i) {
|
|
const double got =
|
|
asymptotic_schwarzschild_phi(cases[i].rho, cases[i].beta);
|
|
CHECK(fabs(got - cases[i].value) < 2e-13, "phi high-precision reference");
|
|
}
|
|
}
|
|
|
|
/* High-precision (mpmath, 60 digits) finish references covering radial,
|
|
* complex-pair, three-real, grazing, and large-radius scattering. The
|
|
* acceptance standard here is the error-budget-driven 1e-8 rad, not the
|
|
* measured ~1e-13. */
|
|
static void test_finish_reference_constants(void) {
|
|
SpacetimeSource source = {0};
|
|
CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0,
|
|
"create schwarzschild");
|
|
SpacetimeAsymptoticEnd end;
|
|
CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end");
|
|
static const struct {
|
|
double rho, beta, n[3];
|
|
} cases[] = {
|
|
{256.0, 0.0, {0.8, 0.6, 0.0}},
|
|
{256.0, 3.0, {0.80697631554502468, 0.59058380112341107, 0.0}},
|
|
{256.0, 60.0, {0.91833674808504193, 0.39579997109220491, 0.0}},
|
|
{256.0, 255.0, {0.69006511550899122, -0.72374728763399356, 0.0}},
|
|
{1.0e6, 1.0, {0.8000005999996, 0.5999991999997, 0.0}},
|
|
};
|
|
for (size_t i = 0; i < sizeof cases / sizeof cases[0]; ++i) {
|
|
SchwarzschildCanonical canonical = {.end_id = 0,
|
|
.t = 0.0,
|
|
.rho = cases[i].rho,
|
|
.rhat = {0.8, 0.6, 0.0},
|
|
.Lhat = {0.0, 0.0, 1.0},
|
|
.beta = cases[i].beta,
|
|
.energy = 2.5,
|
|
.radial_sign = 1};
|
|
double n_inf[3], frequency = 0.0;
|
|
CHECK(asymptotic_schwarzschild_finish(&end, &canonical, n_inf,
|
|
&frequency) == 0,
|
|
"finish reference runs");
|
|
CHECK(angle_between(n_inf, cases[i].n) < 1e-8,
|
|
"finish n_inf high-precision reference");
|
|
CHECK(fabs(frequency - 0.4) < 1e-10 * 0.4,
|
|
"finish frequency high-precision reference");
|
|
}
|
|
spacetime_destroy(&source);
|
|
}
|
|
|
|
/* Turning equation residual |Q| at the computed turning radius. The final
|
|
* scattering direction is validated by test_grazing_reference(). */
|
|
static void test_turning_reference(void) {
|
|
const double betas[] = {3.0 * sqrt(3.0) + 1e-9, 5.5, 6.0, 10.0,
|
|
60.0, 255.0, 3890.44};
|
|
for (size_t i = 0; i < sizeof betas / sizeof betas[0]; ++i) {
|
|
const double rho = asymptotic_schwarzschild_turning_rho(betas[i]);
|
|
CHECK(isfinite(rho) && rho > 3.0, "turning radius exists and is exterior");
|
|
const double Q =
|
|
1.0 - betas[i] * betas[i] * (1.0 - 2.0 / rho) / (rho * rho);
|
|
CHECK(fabs(Q) <= 1e-11, "turning equation residual");
|
|
}
|
|
}
|
|
|
|
/* High-precision entry coordinate-time and swept-azimuth references, checking
|
|
* both the KS time transfer and the entry direction construction. */
|
|
static void test_time_reference(void) {
|
|
SpacetimeSource source = {0};
|
|
CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0,
|
|
"create schwarzschild");
|
|
SpacetimeAsymptoticEnd end;
|
|
CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end");
|
|
static const struct {
|
|
double rho_cam, beta, time, dphi;
|
|
} cases[] = {
|
|
{500.0, 10.0, 246.7884398934447041137, 0.01907105306677434549856},
|
|
{500.0, 0.3, 246.693149021326385798, 0.0005718752307574524347376},
|
|
{256.5, 10.0, 0.5082474340167056157528, 0.00007620281793853560952548},
|
|
{256.5, 0.3, 0.5078666185420216617125, 0.000002284358278417004222584},
|
|
{1000.0, 50.0, 753.1528987272233278083, 0.1465476883815797019938},
|
|
{1.0e6, 10.0, 999777.308033789182372,
|
|
0.03906238263856681534781},
|
|
};
|
|
for (size_t i = 0; i < sizeof cases / sizeof cases[0]; ++i) {
|
|
SchwarzschildCanonical camera = {.end_id = 0,
|
|
.t = 0.0,
|
|
.rho = cases[i].rho_cam,
|
|
.rhat = {1.0, 0.0, 0.0},
|
|
.Lhat = {0.0, 0.0, 1.0},
|
|
.beta = cases[i].beta,
|
|
.energy = 1.0,
|
|
.radial_sign = -1};
|
|
SchwarzschildRouteKind kind = SCH_ROUTE_UNSUPPORTED;
|
|
double activate_t = 0.0, x[3], Pi[3], log_alpha_p0 = 0.0, n_inf[3],
|
|
frequency = 0.0;
|
|
CHECK(asymptotic_schwarzschild_preroute(
|
|
&end, 256.0, &camera, &kind, &activate_t, x, Pi, &log_alpha_p0,
|
|
n_inf, &frequency) == 0 &&
|
|
kind == SCH_ROUTE_ENTRY,
|
|
"reference pre-route entry");
|
|
/* Error-budget-driven mixed tolerance, well below one ODE step (0.1 M)
|
|
* and future metric cadence. */
|
|
const double time_tol = 1e-7 + 1e-11 * fabs(cases[i].time);
|
|
CHECK(fabs(-activate_t - cases[i].time) < time_tol,
|
|
"entry time high-precision reference");
|
|
const double radius =
|
|
sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
|
|
const double rhat[3] = {x[0] / radius, x[1] / radius, x[2] / radius};
|
|
CHECK(fabs(angle_between(camera.rhat, rhat) - cases[i].dphi) < 2e-11,
|
|
"entry azimuth high-precision reference");
|
|
}
|
|
spacetime_destroy(&source);
|
|
}
|
|
|
|
/* Near-grazing references where the exterior integrals are most sensitive:
|
|
* the two sides of beta_R enter through different branches and the KS time
|
|
* integral has a near-singular endpoint. (A photon-sphere turning is not
|
|
* reachable from a camera outside R/M >= 64, so it is not tested here.) */
|
|
static void test_grazing_reference(void) {
|
|
SpacetimeSource source = {0};
|
|
CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0,
|
|
"create schwarzschild");
|
|
SpacetimeAsymptoticEnd end;
|
|
CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end");
|
|
const double beta_R = 256.0 / sqrt(1.0 - 2.0 / 256.0);
|
|
SchwarzschildCanonical hit = {.end_id = 0,
|
|
.t = 0.0,
|
|
.rho = 500.0,
|
|
.rhat = {1.0, 0.0, 0.0},
|
|
.Lhat = {0.0, 0.0, 1.0},
|
|
.beta = beta_R * (1.0 - 1e-12),
|
|
.energy = 1.0,
|
|
.radial_sign = -1};
|
|
SchwarzschildRouteKind kind;
|
|
double activate_t, x[3], Pi[3], log_alpha_p0, n_inf[3], frequency;
|
|
CHECK(asymptotic_schwarzschild_preroute(&end, 256.0, &hit, &kind,
|
|
&activate_t, x, Pi, &log_alpha_p0,
|
|
n_inf, &frequency) == 0 &&
|
|
kind == SCH_ROUTE_ENTRY,
|
|
"near-grazing inside enters");
|
|
const double dphi_ref = 1.0389037630217253661;
|
|
const double time_ref = 434.0116073725480308524;
|
|
const double radius = sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
|
|
const double rhat[3] = {x[0] / radius, x[1] / radius, x[2] / radius};
|
|
CHECK(fabs(angle_between(hit.rhat, rhat) - dphi_ref) < 1e-8,
|
|
"near-grazing entry azimuth");
|
|
CHECK(fabs(-activate_t - time_ref) < 1e-7 + 1e-11 * time_ref,
|
|
"near-grazing entry time");
|
|
|
|
SchwarzschildCanonical miss = hit;
|
|
miss.beta = beta_R * (1.0 + 1e-12);
|
|
CHECK(asymptotic_schwarzschild_preroute(&end, 256.0, &miss, &kind,
|
|
&activate_t, x, Pi, &log_alpha_p0,
|
|
n_inf, &frequency) == 0 &&
|
|
kind == SCH_ROUTE_ESCAPED,
|
|
"near-grazing outside misses");
|
|
const double n_ref[3] = {-0.86581533530640297059,
|
|
-0.50036367289028986187, 0.0};
|
|
CHECK(angle_between(n_inf, n_ref) < 1e-8, "near-grazing miss n_inf");
|
|
spacetime_destroy(&source);
|
|
}
|
|
|
|
/* Deterministic coverage of the three pre-route branches: past-outward,
|
|
* past-inward hit, and past-inward miss (turn before the worldtube). */
|
|
static void test_preroute_branches(void) {
|
|
SpacetimeSource source = {0};
|
|
CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0,
|
|
"create schwarzschild");
|
|
SpacetimeAsymptoticEnd end;
|
|
CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end");
|
|
const double beta_R = 256.0 / sqrt(1.0 - 2.0 / 256.0);
|
|
|
|
SchwarzschildCanonical base = {.end_id = 0,
|
|
.t = 0.0,
|
|
.rho = 500.0,
|
|
.rhat = {1.0, 0.0, 0.0},
|
|
.Lhat = {0.0, 0.0, 1.0},
|
|
.beta = 10.0,
|
|
.energy = 1.0,
|
|
.radial_sign = -1};
|
|
SchwarzschildRouteKind kind;
|
|
double activate_t, x[3], Pi[3], log_alpha_p0, n_inf[3], frequency;
|
|
const double outward_eps = 1e-12;
|
|
|
|
SchwarzschildCanonical outward = base;
|
|
outward.radial_sign = 1;
|
|
CHECK(asymptotic_schwarzschild_preroute(&end, 256.0, &outward, &kind,
|
|
&activate_t, x, Pi, &log_alpha_p0,
|
|
n_inf, &frequency) == 0 &&
|
|
kind == SCH_ROUTE_ESCAPED,
|
|
"past-outward branch escapes");
|
|
CHECK(fabs(sqrt(n_inf[0]*n_inf[0]+n_inf[1]*n_inf[1]+n_inf[2]*n_inf[2]) -
|
|
1.0) < outward_eps,
|
|
"outward n_inf is unit");
|
|
CHECK(fabs(frequency - 1.0) < 1e-12, "outward frequency");
|
|
|
|
SchwarzschildCanonical hit = base;
|
|
CHECK(asymptotic_schwarzschild_preroute(&end, 256.0, &hit, &kind,
|
|
&activate_t, x, Pi, &log_alpha_p0,
|
|
n_inf, &frequency) == 0 &&
|
|
kind == SCH_ROUTE_ENTRY,
|
|
"past-inward hit branch enters");
|
|
|
|
/* Genuine on-boundary tangent: rho = R, beta = beta_R (so Q = 0), zero
|
|
* radial past component. It must not enter. */
|
|
SchwarzschildCanonical tangent = base;
|
|
tangent.rho = 256.0;
|
|
tangent.beta = beta_R;
|
|
tangent.radial_sign = 0;
|
|
CHECK(asymptotic_schwarzschild_preroute(&end, 256.0, &tangent, &kind,
|
|
&activate_t, x, Pi, &log_alpha_p0,
|
|
n_inf, &frequency) == 0 &&
|
|
kind == SCH_ROUTE_ESCAPED,
|
|
"on-boundary tangent escapes");
|
|
|
|
SchwarzschildCanonical miss = base;
|
|
miss.beta = beta_R + 5.0;
|
|
CHECK(asymptotic_schwarzschild_preroute(&end, 256.0, &miss, &kind,
|
|
&activate_t, x, Pi, &log_alpha_p0,
|
|
n_inf, &frequency) == 0 &&
|
|
kind == SCH_ROUTE_ESCAPED,
|
|
"past-inward miss branch escapes");
|
|
CHECK(fabs(sqrt(n_inf[0]*n_inf[0]+n_inf[1]*n_inf[1]+n_inf[2]*n_inf[2]) -
|
|
1.0) < outward_eps,
|
|
"miss n_inf is unit");
|
|
/* A turning ray is deflected away from the radial direction. */
|
|
CHECK(angle_between(n_inf, miss.rhat) > 1e-3,
|
|
"miss n_inf is deflected");
|
|
spacetime_destroy(&source);
|
|
}
|
|
|
|
/* Translated-origin wrapper around the analytic Schwarzschild KS source: the
|
|
* inner metric is evaluated at x - origin and the worldtube/end are shifted by
|
|
* the same origin. This models a black hole at a large coordinate offset; the
|
|
* double reconstruction of the boundary state loses the sub-ULP offset and
|
|
* trips the common entry fallback, while the exact inward orbit transfer stays
|
|
* valid. It exercises the production fallback path, not a synthetic
|
|
* nonlinearity. */
|
|
typedef struct {
|
|
SpacetimeSource inner;
|
|
double origin[3];
|
|
} ShiftedOriginContext;
|
|
|
|
static SpacetimePointStatus shifted_origin_eval(const SpacetimeSource *source,
|
|
double t, const double x[3],
|
|
MetricData *metric) {
|
|
const ShiftedOriginContext *ctx = source->context;
|
|
const double local[3] = {x[0] - ctx->origin[0], x[1] - ctx->origin[1],
|
|
x[2] - ctx->origin[2]};
|
|
return spacetime_eval(&ctx->inner, t, local, metric);
|
|
}
|
|
|
|
static SpacetimeRayStatus shifted_origin_classify(const SpacetimeSource *source,
|
|
double t,
|
|
const double x[3]) {
|
|
const ShiftedOriginContext *ctx = source->context;
|
|
const double local[3] = {x[0] - ctx->origin[0], x[1] - ctx->origin[1],
|
|
x[2] - ctx->origin[2]};
|
|
return spacetime_classify(&ctx->inner, t, local);
|
|
}
|
|
|
|
static size_t shifted_origin_end_count(const SpacetimeSource *source) {
|
|
const ShiftedOriginContext *ctx = source->context;
|
|
return spacetime_asymptotic_end_count(&ctx->inner);
|
|
}
|
|
|
|
static int shifted_origin_end(const SpacetimeSource *source, size_t index,
|
|
SpacetimeAsymptoticEnd *out) {
|
|
const ShiftedOriginContext *ctx = source->context;
|
|
if (spacetime_asymptotic_end(&ctx->inner, index, out))
|
|
return -1;
|
|
for (int i = 0; i < 3; ++i)
|
|
out->frame_origin[i] = ctx->origin[i];
|
|
return 0;
|
|
}
|
|
|
|
static int shifted_origin_worldtube(const SpacetimeSource *source,
|
|
SpacetimeEndId end_id, double t,
|
|
SpacetimeEscapeWorldtubeSample *out) {
|
|
const ShiftedOriginContext *ctx = source->context;
|
|
if (spacetime_escape_worldtube_sample(&ctx->inner, end_id, t, out))
|
|
return -1;
|
|
for (int i = 0; i < 3; ++i)
|
|
out->center[i] += ctx->origin[i];
|
|
return 0;
|
|
}
|
|
|
|
static void shifted_origin_destroy(SpacetimeSource *source) {
|
|
ShiftedOriginContext *ctx = source->context;
|
|
if (ctx != NULL) {
|
|
spacetime_destroy(&ctx->inner);
|
|
free(ctx);
|
|
}
|
|
source->context = NULL;
|
|
source->ops = NULL;
|
|
}
|
|
|
|
static const SpacetimeOps shifted_origin_ops = {
|
|
.eval = shifted_origin_eval,
|
|
.classify = shifted_origin_classify,
|
|
.asymptotic_end_count = shifted_origin_end_count,
|
|
.asymptotic_end = shifted_origin_end,
|
|
.escape_worldtube_sample = shifted_origin_worldtube,
|
|
.destroy = shifted_origin_destroy};
|
|
|
|
static void test_translated_origin_fallback(void) {
|
|
ShiftedOriginContext *ctx = malloc(sizeof *ctx);
|
|
CHECK(ctx != NULL, "shifted-origin context");
|
|
if (ctx == NULL)
|
|
return;
|
|
ctx->origin[0] = 1.0e6;
|
|
ctx->origin[1] = 2.0e6;
|
|
ctx->origin[2] = -3.0e6;
|
|
CHECK(spacetime_create_schwarzschild_ks(&ctx->inner, 1.0, 256.0) == 0,
|
|
"shifted-origin inner source");
|
|
SpacetimeSource source = {.ops = &shifted_origin_ops, .context = ctx};
|
|
|
|
ObserverCamera cam = {.look_ra_deg = 180.0, .look_dec_deg = 0.0};
|
|
for (int i = 0; i < 3; ++i)
|
|
cam.position[i] = ctx->origin[i];
|
|
cam.position[0] += 500.0;
|
|
MetricData metric;
|
|
CHECK(spacetime_eval(&source, 0.0, cam.position, &metric) == 0,
|
|
"shifted-origin camera metric");
|
|
ObserverState observer;
|
|
CHECK(observer_from_coordinate_camera(&metric, &cam, &observer, NULL) ==
|
|
OBSERVER_BUILD_OK,
|
|
"shifted-origin camera observer");
|
|
const double direction[3] = {cos(0.3), sin(0.3), 0.0};
|
|
AsymptoticRoute route;
|
|
CHECK(asymptotic_route_camera(&source, &observer, direction, &route) ==
|
|
ASYMPTOTIC_OK &&
|
|
route.kind == ASYMPTOTIC_ROUTE_ENTRY,
|
|
"shifted-origin entry found");
|
|
CHECK(route.entry_fallback_evaluations > 0,
|
|
"shifted-origin entry used the common fallback");
|
|
CHECK(route.failure_reason == RAY_REASON_NONE,
|
|
"shifted-origin fallback has no failure reason");
|
|
double value;
|
|
CHECK(asymptotic_worldtube_value(&source, route.end_id, route.activate_t,
|
|
route.x, &value) == 0 &&
|
|
value <= 0.0,
|
|
"shifted-origin fallback state is inside the worldtube");
|
|
CHECK(route.activate_t < 0.0, "shifted-origin entry is in the past");
|
|
spacetime_destroy(&source);
|
|
}
|
|
|
|
int main(void) {
|
|
test_round_trip();
|
|
test_finish_matches_integration();
|
|
test_preroute_entry();
|
|
test_camera_reference_radius_independent();
|
|
test_phi_reference_constants();
|
|
test_finish_reference_constants();
|
|
test_turning_reference();
|
|
test_time_reference();
|
|
test_grazing_reference();
|
|
test_preroute_branches();
|
|
test_translated_origin_fallback();
|
|
if (failures == 0)
|
|
puts("asymptotic schwarzschild regression passed");
|
|
else
|
|
fprintf(stderr, "%d asymptotic schwarzschild failures\n", failures);
|
|
return failures == 0 ? 0 : 1;
|
|
}
|