#include "geodesic.h"
#include "observer_track.h"

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

static int nearly_equal(double a, double b) { return fabs(a - b) < 1e-12; }

static int check_ray(const SpacetimeSource *source,
                     const ObserverState *observer,
                     const double local_direction[3],
                     const double expected[3]) {
  const GeodesicTraceConfig config = {.coordinate_time_step = 0.25,
                                      .max_steps = 100};
  RayEndpoint ray =
      geodesic_trace_past(source, observer, local_direction, &config);
  if (ray.status != RAY_ENDPOINT_ESCAPED ||
      !nearly_equal(ray.frequency_ratio, 1.0) ||
      !nearly_equal(ray.n_infinity[0], expected[0]) ||
      !nearly_equal(ray.n_infinity[1], expected[1]) ||
      !nearly_equal(ray.n_infinity[2], expected[2])) {
    fprintf(stderr, "flat-space geodesic regression failed\n");
    return 1;
  }
  return 0;
}

int main(void) {
  SpacetimeSource source = {0};
  const ObserverState observer = observer_fixed_at_origin();
  if (spacetime_create_minkowski(&source, 10.0))
    return 1;
  int result = check_ray(&source, &observer, (double[]){1.0, 0.0, 0.0},
                         (double[]){0.0, 0.0, -1.0}) ||
               check_ray(&source, &observer, (double[]){0.0, 0.0, 1.0},
                         (double[]){1.0, 0.0, 0.0});
  const ObserverState look_at_ra_zero =
      observer_fixed_at_origin_look_at(0.0, 0.0);
  result = result || check_ray(&source, &look_at_ra_zero,
                               (double[]){1.0, 0.0, 0.0},
                               (double[]){1.0, 0.0, 0.0});
  ObserverTrack accelerated = {0};
  ObserverState final_observer;
  if (observer_track_generate_minkowski_acceleration(&accelerated, 1.52, 2.0,
                                                      1.0 / 30.0) ||
      observer_track_interpolate(&accelerated, 2.0, &final_observer, NULL)) {
    result = 1;
  } else {
    const RayEndpoint forward = geodesic_trace_past(
        &source, &final_observer, (double[]){1.0, 0.0, 0.0},
        &(GeodesicTraceConfig){.coordinate_time_step = 0.25, .max_steps = 100});
    const double expected_g = sqrt(1.0 + 3.04 * 3.04) + 3.04;
    if (forward.status != RAY_ENDPOINT_ESCAPED ||
        !nearly_equal(forward.frequency_ratio, expected_g)) {
      fputs("accelerated-observer Doppler regression failed\n", stderr);
      result = 1;
    }
  }
  observer_track_destroy(&accelerated);
  spacetime_destroy(&source);
  return result;
}