#include "spacetime.h"

#include <math.h>
#include <stddef.h>
#include <stdlib.h>

/* Escape sphere lies this many wall thicknesses 1/sigma beyond R.  At
 * r = R + span/sigma the shape has decayed to ~2*exp(-2*span) for a thin wall
 * and ~4*exp(-2*span) for a broad bump, i.e. below binary64 epsilon for
 * span = 20, so the escape sphere is Minkowski to machine accuracy. */
#define ALCUBIERRE_ESCAPE_SPAN 20.0

/* For 2 sigma R below this threshold the direct difference of two nearby
 * tanh values loses about 1/(2 sigma R) digits and can round the shape to
 * zero while it is still O(1).  Switch to an algebraically equivalent form
 * that is free of cancellation in that regime. */
#define ALCUBIERRE_SMALL_WALL 0.5

typedef struct {
  double vs;
  double radius;
  double sigma;
  double escape_radius;
} AlcubierreContext;

/* Alcubierre shape function
 *   f(r) = (tanh(sigma (r + R)) - tanh(sigma (r - R))) / (2 tanh(sigma R)),
 * positive, equal to 1 at r = 0 for any sigma R > 0, and decaying to zero
 * past r = R over a transition width ~1/sigma.
 *
 * The identity f = (1 - s^2) / (1 - s^2 t^2) with s = tanh(sigma r),
 * t = tanh(sigma R) is exact and has no cancellation when sigma R is small,
 * where the shape tends to sech^2(sigma r). */
static double alcubierre_shape(double r, double radius, double sigma) {
  const double sR = sigma * radius;
  if (2.0 * sR < ALCUBIERRE_SMALL_WALL) {
    const double s = tanh(sigma * r);
    const double t = tanh(sR);
    return (1.0 - s * s) / (1.0 - s * s * t * t);
  }
  return (tanh(sigma * (r + radius)) - tanh(sigma * (r - radius))) /
         (2.0 * tanh(sR));
}

/* d f / d r.  The thin-wall branch uses sech^2(x) = 1 - tanh(x)^2; the
 * broad-bump branch uses the cancellation-free derivative of the identity
 * above.  Both underflow to zero far outside the bubble, which is the
 * intended exactly-flat limit. */
static double alcubierre_shape_derivative(double r, double radius,
                                          double sigma) {
  const double sR = sigma * radius;
  if (2.0 * sR < ALCUBIERRE_SMALL_WALL) {
    const double s = tanh(sigma * r);
    const double c = cosh(2.0 * sR);
    const double denom = 1.0 + s * s + c * (1.0 - s * s);
    return -(c + 1.0) * 4.0 * sigma * s * (1.0 - s * s) / (denom * denom);
  }
  const double tanh_plus = tanh(sigma * (r + radius));
  const double tanh_minus = tanh(sigma * (r - radius));
  const double sech2_plus = 1.0 - tanh_plus * tanh_plus;
  const double sech2_minus = 1.0 - tanh_minus * tanh_minus;
  return sigma * (sech2_plus - sech2_minus) / (2.0 * tanh(sR));
}

/* Moving Alcubierre bubble in the lab coordinates
 *   ds^2 = -dt^2 + (dx - v_s f(r_s) dt)^2 + dy^2 + dz^2,
 *   r_s = sqrt((x - x_s)^2 + y^2 + z^2),   x_s(t) = v_s t,
 * with x_s(0) = 0.  This is not a comoving (x_s = 0) slicing: the bubble
 * propagates through the coordinates.  The spatial slices stay flat, so
 * alpha = 1, gamma_ij = delta_ij, beta^x = -v_s f(r_s), and
 *   K_ij = (D_i beta_j + D_j beta_i) / (2 alpha)
 *        = -v_s (delta_jx d_i f + delta_ix d_j f) / 2,
 * where d_i differentiates at fixed t (only the spatial argument of f moves
 * with t).  K encodes the time dependence required by the 3+1 null-ray RHS. */
static SpacetimePointStatus alcubierre_eval(const SpacetimeSource *source,
                                            double t, const double x[3],
                                            MetricData *metric) {
  const AlcubierreContext *context = source->context;
  const double vs = context->vs;
  const double dx = x[0] - vs * t;
  const double r2 = dx * dx + x[1] * x[1] + x[2] * x[2];
  double df[3] = {0.0, 0.0, 0.0};
  double f;
  if (!isfinite(r2))
    return SPACETIME_POINT_INVALID_METRIC;
  const double r = sqrt(r2);
  *metric = (MetricData){
      .alpha = 1.0,
      .gamma = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}};
  if (r > 0.0) {
    f = alcubierre_shape(r, context->radius, context->sigma);
    const double radial_scale =
        alcubierre_shape_derivative(r, context->radius, context->sigma) / r;
    df[0] = radial_scale * dx;
    df[1] = radial_scale * x[1];
    df[2] = radial_scale * x[2];
  } else {
    f = alcubierre_shape(0.0, context->radius, context->sigma);
  }
  metric->beta[0] = -vs * f;
  for (int i = 0; i < 3; ++i) {
    metric->d_beta[i][0] = -vs * df[i];
    for (int j = 0; j < 3; ++j)
      metric->K[i][j] =
          -0.5 * vs * ((j == 0 ? df[i] : 0.0) + (i == 0 ? df[j] : 0.0));
  }
  return SPACETIME_POINT_OK;
}

/* The exotic matter that would source the bubble is treated as optically
 * transparent.  The escape sphere follows the moving bubble. */
static SpacetimeRayStatus alcubierre_classify(const SpacetimeSource *source,
                                              double t, const double x[3]) {
  const AlcubierreContext *context = source->context;
  const double dx = x[0] - context->vs * t;
  const double r2 = dx * dx + x[1] * x[1] + x[2] * x[2];
  return r2 >= context->escape_radius * context->escape_radius
             ? SPACETIME_RAY_ESCAPED
             : SPACETIME_RAY_ACTIVE;
}

static void alcubierre_destroy(SpacetimeSource *source) {
  free(source->context);
  source->context = NULL;
  source->ops = NULL;
}

static size_t alcubierre_asymptotic_end_count(const SpacetimeSource *source) {
  (void)source;
  return 1;
}

static int alcubierre_asymptotic_end(const SpacetimeSource *source,
                                     size_t index,
                                     SpacetimeAsymptoticEnd *out) {
  (void)source;
  if (index != 0)
    return -1;
  *out = (SpacetimeAsymptoticEnd){
      .end_id = 0,
      .exterior_kind = ASYMPTOTIC_EXTERIOR_MINKOWSKI,
      .mass = 0.0,
      .frame_origin = {0.0, 0.0, 0.0},
      .frame_axes = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}};
  return 0;
}

static int alcubierre_escape_worldtube_sample(
    const SpacetimeSource *source, SpacetimeEndId end_id, double t,
    SpacetimeEscapeWorldtubeSample *out) {
  const AlcubierreContext *context = source->context;
  if (end_id != 0)
    return -1;
  *out = (SpacetimeEscapeWorldtubeSample){
      .center = {context->vs * t, 0.0, 0.0},
      .velocity = {context->vs, 0.0, 0.0},
      .radius = context->escape_radius,
      .radius_rate = 0.0,
      .velocity_constant = 1,
      .valid = 1};
  return 0;
}

static const SpacetimeOps alcubierre_ops = {
    .eval = alcubierre_eval,
    .classify = alcubierre_classify,
    .asymptotic_end_count = alcubierre_asymptotic_end_count,
    .asymptotic_end = alcubierre_asymptotic_end,
    .escape_worldtube_sample = alcubierre_escape_worldtube_sample,
    .destroy = alcubierre_destroy,
};

double spacetime_alcubierre_escape_radius(double radius, double sigma) {
  return radius + ALCUBIERRE_ESCAPE_SPAN / sigma;
}

int spacetime_create_alcubierre(SpacetimeSource *source, double vs,
                                double radius, double sigma) {
  if (source == NULL || !isfinite(vs) ||
      !isfinite(radius) || radius <= 0.0 || !isfinite(sigma) || sigma <= 0.0)
    return -1;
  /* Reject parameter combinations whose derived domain overflows or does not
   * actually extend beyond the bubble. */
  const double escape_radius = spacetime_alcubierre_escape_radius(radius, sigma);
  if (!isfinite(escape_radius) || escape_radius <= radius)
    return -1;
  AlcubierreContext *context = malloc(sizeof *context);
  if (context == NULL)
    return -1;
  context->vs = vs;
  context->radius = radius;
  context->sigma = sigma;
  context->escape_radius = escape_radius;
  source->ops = &alcubierre_ops;
  source->context = context;
  if (spacetime_source_finalize(source)) {
    alcubierre_destroy(source);
    return -1;
  }
  return 0;
}

int spacetime_create_default(SpacetimeSource *source) {
  return spacetime_create_alcubierre(source, 0.5, 5.0, 1.0);
}