GR-raytracing/tests/test_schwarzschild.c

#include "geodesic.h"
#include "frame.h"

#include <math.h>
#include <stdio.h>

static int camera_at(const SpacetimeSource *source, double radius,
                     double ra, double dec, ObserverState *out) {
  ObserverCamera camera = {.look_ra_deg = ra, .look_dec_deg = dec};
  const ObserverState pointing = observer_fixed_at_origin_look_at(ra, dec);
  for (int i = 0; i < 3; ++i)
    camera.position[i] = -radius * pointing.tetrad[1][i + 1];
  MetricData metric;
  return spacetime_eval(source, 0.0, camera.position, &metric) ||
         observer_from_coordinate_camera(&metric, &camera, out, NULL);
}

int main(void) {
  SpacetimeSource spacetime = {0};
  MetricData metric;
  ObserverState observer;
  ObserverState oriented_observer;
  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;
}