GR-raytracing/tests/test_termination_oracle.c

/*
 * Independent physics oracle for the ray-termination policy (plan P0).
 *
 * This test does not read production endpoints for its central assertions.
 * It builds Schwarzschild-KS states independently and checks:
 *   1. the two radial null branches dr/ds = 1 and dr/ds = (2M-r)/(2M+r);
 *   2. the critical impact parameter b = 3 sqrt(3) M and photon sphere r = 3M;
 *   3. camera energy normalization E_camera = 1 and tetrad orthonormality;
 *   4. the threshold proxy identity ln(p^0) = L - ln(alpha) with
 *      L = ln(alpha p^0).
 *
 * The radial-branch assertions are integrated with the production RK4 RHS in
 * src/geodesic.c so that an error in the 3+1 reduction is caught against a
 * closed-form invariant rather than against a second copy of the same algebra.
 */
#include "asymptotic_schwarzschild.h"
#include "geodesic.h"
#include "observer.h"
#include "spacetime.h"

#include <math.h>
#include <stdio.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 Eulerian orthonormal tetrad at x = (r0, 0, 0) for r0 > 0.  At this
 * point the KS spatial metric is diagonal, so the principal axes are already
 * orthonormal (up to the radial scale sqrt(gamma_rr)). */
static void radial_static_observer(const MetricData *metric, double r0,
                                   ObserverState *out) {
  *out = (ObserverState){.coordinate_time = 0.0,
                         .coordinate_position = {r0, 0.0, 0.0}};
  const double alpha = metric->alpha;
  out->tetrad[0][0] = 1.0 / alpha;
  for (int i = 0; i < 3; ++i)
    out->tetrad[0][i + 1] = -metric->beta[i] / alpha;
  const double radial_scale = sqrt(metric->gamma[0][0]);
  out->tetrad[1][1] = 1.0 / radial_scale;
  out->tetrad[2][2] = 1.0;
  out->tetrad[3][3] = 1.0;
}

/* Integrate a purely radial past ray with the production stepper and return
 * its final state.  Output endpoint is not inspected. */
static int trace_radial(const SpacetimeSource *source, const ObserverState *o,
                        double direction, GeodesicRayState *state) {
  MetricData metric;
  if (spacetime_eval(source, o->coordinate_time, o->coordinate_position,
                     &metric) != SPACETIME_POINT_OK)
    return -1;
  const double n[3] = {direction, 0.0, 0.0};
  if (geodesic_initialize_past_ray_metric(&metric, o, n, state))
    return -1;
  const GeodesicTraceConfig config = {.coordinate_time_step = 0.02,
                                      .max_steps = 400,
                                      .threshold = {.kind = THRESHOLD_DISABLED, .value = 0.0, .policy_version = 0}};
  MetricSlab *slab = NULL;
  if (spacetime_load_slab(source, 0.0, -1000.0, &slab))
    return -1;
  RayEndpoint endpoint = {.end_id = SPACETIME_END_NONE,
                          .outcome = RAY_OUTCOME_INCOMPLETE};
  const GeodesicAdvanceResult result =
      geodesic_advance_past_ray(slab, state, -1000.0, &config, &endpoint);
  spacetime_free_slab(slab);
  return result == GEODESIC_ADVANCE_FAILED ? -1 : 0;
}

static void test_radial_branches(void) {
  SpacetimeSource source = {0};
  CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0,
        "create schwarzschild");
  const double r0 = 10.0;
  MetricData metric;
  CHECK(spacetime_eval(&source, 0.0, (double[]){r0, 0.0, 0.0}, &metric) ==
            SPACETIME_POINT_OK,
        "metric at r0");
  ObserverState observer;
  radial_static_observer(&metric, r0, &observer);

  /* Branch dr/ds = +1: the closed-form solution is r = r0 + s. */
  GeodesicRayState outward;
  CHECK(trace_radial(&source, &observer, 1.0, &outward) == 0, "outward trace");
  const double s_out = -outward.coordinate_time;
  const double invariant_out = outward.x[0] - r0 - s_out;
  CHECK(fabs(invariant_out) < 1e-6, "outward branch r = r0 + s");

  /* Branch dr/ds = (2M-r)/(2M+r): the closed-form invariant is
   * (r-2M) + 4M ln(r-2M) + s = const. */
  GeodesicRayState inward;
  CHECK(trace_radial(&source, &observer, -1.0, &inward) == 0, "inward trace");
  const double s_in = -inward.coordinate_time;
  const double c0 = (r0 - 2.0) + 4.0 * log(r0 - 2.0);
  const double c1 = (inward.x[0] - 2.0) + 4.0 * log(inward.x[0] - 2.0) + s_in;
  CHECK(inward.x[0] > 2.0, "inward branch stays outside the horizon");
  CHECK(inward.x[0] < r0, "inward branch decreases r");
  CHECK(fabs(c1 - c0) < 1e-6, "inward branch closed-form invariant");

  /* Both branches are time-reversal partners: the outward and inward states
   * reach the same |dr/ds| magnitude in opposite senses at r0. */
  CHECK(outward.x[0] > r0, "outward branch increases r");
  spacetime_destroy(&source);
}

/* The production dark policy is the camera-relative growth A_0 = L - L_0,
 * independent of the backend.  This oracle retains the stationary-KS
 * conserved Killing energy A_K = L - ln|E_K| as an independent cross-check of
 * the same ray: it verifies E_K conservation and the identity
 * A_K - A_0 = -ln|alpha_0 - beta_0.Pi_0|.  It is not the production
 * criterion. */
static void test_killing_energy_reference(void) {
  SpacetimeSource source = {0};
  CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0,
        "create schwarzschild for killing reference");
  const double r0 = 10.0;
  MetricData start_metric;
  CHECK(spacetime_eval(&source, 0.0, (double[]){r0, 0.0, 0.0}, &start_metric) ==
            SPACETIME_POINT_OK,
        "metric for killing reference");
  ObserverState observer;
  radial_static_observer(&start_metric, r0, &observer);
  GeodesicRayState start;
  CHECK(geodesic_initialize_past_ray_metric(&start_metric, &observer,
                                            (double[]){1.0, 0.0, 0.0},
                                            &start) == 0,
        "initialize killing ray");
  double beta0 = 0.0;
  for (int i = 0; i < 3; ++i)
    beta0 += start_metric.beta[i] * start.Pi[i];
  const double ek0 = exp(start.log_alpha_p0) * (start_metric.alpha - beta0);
  CHECK(isfinite(ek0) && fabs(ek0) > 0.0, "nonzero Killing energy");

  GeodesicRayState end;
  CHECK(trace_radial(&source, &observer, 1.0, &end) == 0,
        "trace killing reference ray");
  MetricData end_metric;
  CHECK(spacetime_eval(&source, end.coordinate_time, end.x, &end_metric) ==
            SPACETIME_POINT_OK,
        "metric at killing reference end");
  double beta1 = 0.0;
  for (int i = 0; i < 3; ++i)
    beta1 += end_metric.beta[i] * end.Pi[i];
  const double ek1 = exp(end.log_alpha_p0) * (end_metric.alpha - beta1);
  CHECK(fabs(ek1 / ek0 - 1.0) < 1e-6,
        "Killing energy conserved along the geodesic");
  const double a0 = end.log_alpha_p0 - start.log_alpha_p0;
  const double ak = end.log_alpha_p0 - log(fabs(ek1));
  const double predicted = -log(fabs(start_metric.alpha - beta0));
  CHECK(fabs((ak - a0) - predicted) < 1e-9,
        "A_K - A_0 equals the initial boost factor");
  spacetime_destroy(&source);
}

static void test_critical_parameters(void) {
  const double b_crit = 3.0 * sqrt(3.0);
  CHECK(!isfinite(asymptotic_schwarzschild_turning_rho(b_crit - 1e-6)),
        "no turning point below b_crit");
  CHECK(!isfinite(asymptotic_schwarzschild_turning_rho(3.0)),
        "no turning point for a deeply plunging ray");
  const double just_above = asymptotic_schwarzschild_turning_rho(b_crit + 1e-6);
  CHECK(isfinite(just_above) && just_above > 3.0 && just_above < 3.01,
        "turning radius approaches the photon sphere at b_crit");
  const double b6 = asymptotic_schwarzschild_turning_rho(6.0);
  CHECK(isfinite(b6) && b6 > 3.0, "turning radius above the photon sphere");
  /* Verify the turning radius is an independent root of
   * f(rho) = rho^3 - b^2 rho + 2 b^2. */
  const double residual = b6 * b6 * b6 - 36.0 * b6 + 72.0;
  CHECK(fabs(residual) < 1e-9, "turning radius satisfies the radial equation");
}

static void test_observer_normalization(void) {
  SpacetimeSource source = {0};
  CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0,
        "create schwarzschild for observer");
  ObserverCamera camera = {.coordinate_time = 0.0,
                           .position = {30.0, 0.0, 0.0},
                           .velocity = {0.0, 0.0, 0.0},
                           .look_ra_deg = 0.0,
                           .look_dec_deg = 0.0};
  MetricData metric;
  ObserverState observer;
  CHECK(spacetime_eval(&source, 0.0, camera.position, &metric) ==
            SPACETIME_POINT_OK,
        "metric at camera");
  CHECK(observer_from_coordinate_camera(&metric, &camera, &observer, NULL) ==
            OBSERVER_BUILD_OK,
        "build observer");
  /* Orthonormality of the production tetrad, independently of the geodesic
   * layer: g(e_a, e_b) = diag(-1, 1, 1, 1). */
  for (int a = 0; a < 4; ++a) {
    for (int b = 0; b < 4; ++b) {
      const double *ea = observer.tetrad[a];
      const double *eb = observer.tetrad[b];
      double inner = -metric.alpha * metric.alpha * ea[0] * eb[0];
      for (int i = 0; i < 3; ++i)
        for (int j = 0; j < 3; ++j)
          inner += metric.gamma[i][j] * (ea[i + 1] + metric.beta[i] * ea[0]) *
                   (eb[j + 1] + metric.beta[j] * eb[0]);
      const double expected = a == b ? (a == 0 ? -1.0 : 1.0) : 0.0;
      CHECK(fabs(inner - expected) < 1e-10, "tetrad orthonormal");
    }
  }
  const double local[3] = {0.3, 0.5, 0.9};
  const double norm = sqrt(local[0] * local[0] + local[1] * local[1] +
                           local[2] * local[2]);
  const double direction[3] = {local[0] / norm, local[1] / norm,
                               local[2] / norm};
  GeodesicRayState state;
  CHECK(geodesic_initialize_past_ray_metric(&metric, &observer, direction,
                                            &state) == 0,
        "initialize past ray");
  /* gamma is diagonal at (30, 0, 0): gamma_xx = 1 + 2/r. */
  double gamma_inv[3][3];
  for (int i = 0; i < 3; ++i)
    for (int j = 0; j < 3; ++j)
      gamma_inv[i][j] = (i == j) ? 1.0 / metric.gamma[i][j] : 0.0;
  /* Null constraint gamma^{ij} Pi_i Pi_j = 1. */
  double null_residual = 0.0;
  for (int i = 0; i < 3; ++i)
    for (int j = 0; j < 3; ++j)
      null_residual += gamma_inv[i][j] * state.Pi[i] * state.Pi[j];
  CHECK(fabs(null_residual - 1.0) < 1e-10, "null constraint preserved");
  /* Observed energy -g(k, e0) = 1 for the unit observer four-velocity.  The
   * photon four-momentum is reconstructed from the stored state:
   * p^0 = exp(L)/alpha and p^i = alpha p^0 gamma^{ij} Pi_j - beta^i p^0. */
  const double k0 = exp(state.log_alpha_p0) / metric.alpha;
  double k[4] = {k0, 0.0, 0.0, 0.0};
  for (int i = 0; i < 3; ++i) {
    double covariant = 0.0;
    for (int j = 0; j < 3; ++j)
      covariant += gamma_inv[i][j] * state.Pi[j];
    k[i + 1] = metric.alpha * k0 * covariant - metric.beta[i] * k0;
  }
  const double *e0 = observer.tetrad[0];
  double inner = -metric.alpha * metric.alpha * k[0] * e0[0];
  for (int i = 0; i < 3; ++i)
    for (int j = 0; j < 3; ++j)
      inner += metric.gamma[i][j] * (k[i + 1] + metric.beta[i] * k[0]) *
               (e0[j + 1] + metric.beta[j] * e0[0]);
  CHECK(fabs(inner + 1.0) < 1e-10, "camera energy normalized to one");
  spacetime_destroy(&source);
}

static void test_threshold_proxies(void) {
  SpacetimeSource source = {0};
  CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0,
        "create schwarzschild for proxies");
  const double r0 = 30.0;
  MetricData metric;
  CHECK(spacetime_eval(&source, 0.0, (double[]){r0, 0.0, 0.0}, &metric) ==
            SPACETIME_POINT_OK,
        "metric for proxies");
  ObserverState observer;
  radial_static_observer(&metric, r0, &observer);
  const double direction[3] = {1.0, 0.0, 0.0};
  GeodesicRayState state;
  CHECK(geodesic_initialize_past_ray_metric(&metric, &observer, direction,
                                            &state) == 0,
        "initialize proxy ray");
  /* L = ln(alpha p^0) is stored; ln(p^0) = L - ln(alpha).  Recompute p^0 from
   * the tetrad and direction independently. */
  const double k0 = observer.tetrad[0][0] - direction[0] * observer.tetrad[1][0] -
                    direction[1] * observer.tetrad[2][0] -
                    direction[2] * observer.tetrad[3][0];
  const double log_p0 = state.log_alpha_p0 - log(metric.alpha);
  CHECK(fabs(log_p0 - log(k0)) < 1e-12,
        "ln(p^0) = L - ln(alpha) with L = ln(alpha p^0)");
  /* For a static observer far outside, alpha -> 1 and the two proxies agree
   * to O(M/r); this documents why a fixed L threshold is not a fixed p^0
   * threshold. */
  CHECK(fabs(state.log_alpha_p0 - log_p0) > 1e-3,
        "L and ln(p^0) differ near the hole");
  spacetime_destroy(&source);
}

int main(void) {
  test_radial_branches();
  test_killing_energy_reference();
  test_critical_parameters();
  test_observer_normalization();
  test_threshold_proxies();
  if (failures != 0) {
    fprintf(stderr, "termination oracle: %d failure(s)\n", failures);
    return 1;
  }
  return 0;
}