#include "geodesic.h" #include "observer.h" #include #include #include #define CHECK(condition) do { if (!(condition)) { \ fprintf(stderr, "alcubierre regression failed at line %d: %s\n", \ __LINE__, #condition); \ return 1; } } while (0) static double shape(double r, double radius, double sigma) { const double sr = sigma * r; const double sR = sigma * radius; return (tanh(sr + sR) - tanh(sr - sR)) / (2.0 * tanh(sR)); } static double metric_g00(const MetricData *m) { double g = -m->alpha * m->alpha; for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) g += m->gamma[i][j] * m->beta[i] * m->beta[j]; return g; } static int eval(const SpacetimeSource *source, double t, const double x[3], MetricData *metric) { return spacetime_eval(source, t, x, metric); } int main(void) { const double vs = 0.5, radius = 5.0, sigma = 1.0; const double escape = spacetime_alcubierre_escape_radius(radius, sigma); SpacetimeSource source = {0}; MetricData metric; CHECK(spacetime_create_alcubierre(&source, vs, radius, sigma) == 0); /* Sub- and super-luminal velocities are both accepted. A successful * constructor installs a context, so each acceptance uses its own temporary * source that is destroyed immediately; the shared `source` above is never * overwritten with a second live context. */ { SpacetimeSource luminal = {0}; CHECK(spacetime_create_alcubierre(&luminal, 1.0, radius, sigma) == 0); spacetime_destroy(&luminal); } { SpacetimeSource superluminal = {0}; CHECK(spacetime_create_alcubierre(&superluminal, -1.5, radius, sigma) == 0); spacetime_destroy(&superluminal); } /* Non-finite velocities stay rejected. */ CHECK(spacetime_create_alcubierre(&source, NAN, radius, sigma) != 0); CHECK(spacetime_create_alcubierre(&source, INFINITY, radius, sigma) != 0); CHECK(spacetime_create_alcubierre(&source, -INFINITY, radius, sigma) != 0); CHECK(spacetime_create_alcubierre(&source, vs, 0.0, sigma) != 0); CHECK(spacetime_create_alcubierre(&source, vs, radius, 0.0) != 0); /* A derived escape radius that overflows or does not exceed R is rejected. */ CHECK(spacetime_create_alcubierre(&source, vs, 1.0, DBL_MIN) != 0); CHECK(spacetime_create_alcubierre(&source, vs, DBL_MAX, 1.0) != 0); /* At t = 0 the bubble is centered on the origin: f = 1, beta^x = -v_s, * flat spatial metric, K = 0. */ CHECK(eval(&source, 0.0, (double[]){0, 0, 0}, &metric) == 0); CHECK(metric.alpha == 1.0); CHECK(fabs(metric.beta[0] + vs) < 1e-15); CHECK(metric.beta[1] == 0.0 && metric.beta[2] == 0.0); for (int i = 0; i < 3; ++i) { CHECK(metric.d_alpha[i] == 0.0); for (int j = 0; j < 3; ++j) { CHECK(metric.gamma[i][j] == (i == j ? 1.0 : 0.0)); CHECK(metric.K[i][j] == 0.0); for (int k = 0; k < 3; ++k) CHECK(metric.d_gamma[i][j][k] == 0.0); } } /* Exact translation symmetry of the moving metric: * g(t, x, y, z) == g(0, x - v_s t, y, z) for the 3+1 data. */ { const double samples[3][4] = {{-2.0, 2.0, 3.0, -1.5}, {4.0, -6.5, 1.0, 2.0}, {-1.0, 0.5, -0.25, 0.75}}; for (int s = 0; s < 3; ++s) { const double t = samples[s][0]; double x[3] = {samples[s][1], samples[s][2], samples[s][3]}; double shifted[3] = {x[0] - vs * t, x[1], x[2]}; MetricData mt, m0; CHECK(eval(&source, t, x, &mt) == 0); CHECK(eval(&source, 0.0, shifted, &m0) == 0); CHECK(fabs(mt.beta[0] - m0.beta[0]) < 1e-14); for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) { CHECK(fabs(mt.d_beta[i][j] - m0.d_beta[i][j]) < 1e-13); CHECK(fabs(mt.K[i][j] - m0.K[i][j]) < 1e-13); } } } /* A generic off-axis point: 3+1 data must reconstruct the literal metric * ds^2 = -dt^2 + (dx - v_s f(r_s) dt)^2 + dy^2 + dz^2. */ const double t = -2.0; const double x[3] = {2.0, 3.0, -1.5}; const double dx = x[0] - vs * t; const double r = sqrt(dx * dx + x[1] * x[1] + x[2] * x[2]); const double f = shape(r, radius, sigma); CHECK(eval(&source, t, x, &metric) == 0); CHECK(fabs(metric_g00(&metric) - (-1.0 + vs * vs * f * f)) < 1e-14); for (int i = 0; i < 3; ++i) { double beta_lower = 0.0; for (int j = 0; j < 3; ++j) beta_lower += metric.gamma[i][j] * metric.beta[j]; const double target = (i == 0) ? -vs * f : 0.0; CHECK(fabs(metric.beta[i] - target) < 1e-14); CHECK(fabs(beta_lower - target) < 1e-14); CHECK(metric.d_alpha[i] == 0.0); } /* d_beta and K against central differences of beta at fixed t. The spatial * metric is flat and constant in time, so K_ij = * (d_i beta_j + d_j beta_i) / 2. */ { const double h = 1e-5; for (int direction = 0; direction < 3; ++direction) { double xp[3] = {x[0], x[1], x[2]}; double xm[3] = {x[0], x[1], x[2]}; MetricData mp, mm; xp[direction] += h; xm[direction] -= h; CHECK(eval(&source, t, xp, &mp) == 0); CHECK(eval(&source, t, xm, &mm) == 0); for (int j = 0; j < 3; ++j) { const double finite_difference = (mp.beta[j] - mm.beta[j]) / (2.0 * h); CHECK(fabs(metric.d_beta[direction][j] - finite_difference) < 1e-6); } } for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) { const double expected = 0.5 * (metric.d_beta[i][j] + metric.d_beta[j][i]); CHECK(fabs(metric.K[i][j] - expected) < 1e-14); CHECK(fabs(metric.K[i][j] - metric.K[j][i]) < 1e-15); } } /* Shape and derivative across the whole sigma*R domain, including the tiny * sigma*R regime where the direct tanh difference loses all its digits. A * long-double cosh form is cancellation-free and serves as the reference. */ { static const double cases[][2] = { {1e-20, 1.0}, {1e-8, 0.5}, {1e-3, 2.0}, {0.1, 0.3}, {0.24, 1.0}, {0.26, 1.0}, {0.5, 0.5}, {1.0, 1.0}, {5.0, 5.0}, {5.0, 8.0}}; for (size_t k = 0; k < sizeof cases / sizeof cases[0]; ++k) { const double sigma_r = cases[k][0]; const double sr = cases[k][1]; SpacetimeSource local = {0}; CHECK(spacetime_create_alcubierre(&local, vs, sigma_r, 1.0) == 0); const double r = sr; /* sigma = 1, so R = sigma_r and r = sigma_r_test */ MetricData m; CHECK(eval(&local, 0.0, (double[]){r, 0.0, 0.0}, &m) == 0); const double f = -m.beta[0] / vs; const double df = -m.d_beta[0][0] / vs; const long double C = coshl(2.0L * (long double)sigma_r); const long double fref = (C + 1.0L) / (coshl(2.0L * (long double)r) + C); const long double dfref = -(C + 1.0L) * 2.0L * sinhl(2.0L * (long double)r) / ((coshl(2.0L * (long double)r) + C) * (coshl(2.0L * (long double)r) + C)); CHECK(fabsl((long double)f - fref) < 1e-12L); CHECK(fabsl((long double)df - dfref) < 1e-9L); /* The escape sphere must be flat to below binary64 epsilon. */ const double local_escape = spacetime_alcubierre_escape_radius(sigma_r, 1.0); CHECK(eval(&local, 0.0, (double[]){local_escape, 0.0, 0.0}, &m) == 0); CHECK(fabs(m.beta[0] / vs) < 1e-15); spacetime_destroy(&local); } } /* Classification follows the bubble and is never CAPTURED. */ CHECK(spacetime_classify(&source, 0.0, (double[]){0, 0, 0}) == SPACETIME_RAY_ACTIVE); CHECK(spacetime_classify(&source, 0.0, (double[]){escape - 0.5, 0, 0}) == SPACETIME_RAY_ACTIVE); CHECK(spacetime_classify(&source, 0.0, (double[]){escape + 1.0, 0, 0}) == SPACETIME_RAY_ESCAPED); CHECK(spacetime_classify(&source, 0.0, (double[]){0, 0, 1000}) == SPACETIME_RAY_ESCAPED); /* At t = 3 the bubble center is at x_s = 1.5; the sphere moves with it. */ CHECK(spacetime_classify(&source, 3.0, (double[]){vs * 3.0, 0, 0}) == SPACETIME_RAY_ACTIVE); CHECK(spacetime_classify(&source, 3.0, (double[]){vs * 3.0 + escape + 1.0, 0, 0}) == SPACETIME_RAY_ESCAPED); /* Isometry check: the moving metric is invariant under the spacetime * translation (t, x) -> (t + T, x + v_s T). Two static observers related by * this isometry must therefore see identical escaping directions and * frequency ratios. This exercises the x_s(t) time dependence end to end. */ { const double T = 3.0; ObserverCamera camera0 = {.coordinate_time = 0.0, .position = {0.0, 0.0, 15.0}, .look_ra_deg = 90.0, .look_dec_deg = -90.0}; ObserverCamera camera1 = {.coordinate_time = T, .position = {vs * T, 0.0, 15.0}, .look_ra_deg = 90.0, .look_dec_deg = -90.0}; MetricData m0, m1; ObserverState o0, o1; CHECK(eval(&source, camera0.coordinate_time, camera0.position, &m0) == 0); CHECK(eval(&source, camera1.coordinate_time, camera1.position, &m1) == 0); CHECK(observer_from_coordinate_camera(&m0, &camera0, &o0, NULL) == OBSERVER_BUILD_OK); CHECK(observer_from_coordinate_camera(&m1, &camera1, &o1, NULL) == OBSERVER_BUILD_OK); const GeodesicTraceConfig trace = {.coordinate_time_step = 0.02, .max_steps = 1u << 20}; const double directions[3][3] = {{1, 0, 0}, {1, 0.25, 0}, {1, 0, 0.3}}; for (int i = 0; i < 3; ++i) { double n[3] = {directions[i][0], directions[i][1], directions[i][2]}; const double norm = sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]); for (int k = 0; k < 3; ++k) n[k] /= norm; const RayEndpoint r0 = geodesic_trace_past(&source, &o0, n, &trace); const RayEndpoint r1 = geodesic_trace_past(&source, &o1, n, &trace); CHECK(r0.outcome == RAY_OUTCOME_ESCAPED); CHECK(r1.outcome == RAY_OUTCOME_ESCAPED); for (int k = 0; k < 3; ++k) CHECK(fabs(r0.n_infinity[k] - r1.n_infinity[k]) < 1e-6); CHECK(fabs(r0.frequency_ratio - r1.frequency_ratio) < 1e-6); } } /* Flat limit v_s = 0 is exactly Minkowski. */ { SpacetimeSource flat = {0}; CHECK(spacetime_create_alcubierre(&flat, 0.0, radius, sigma) == 0); MetricData flat_metric; CHECK(eval(&flat, 0.0, (double[]){2, 3, 4}, &flat_metric) == 0); CHECK(flat_metric.alpha == 1.0); CHECK(flat_metric.beta[0] == 0.0 && flat_metric.beta[1] == 0.0 && flat_metric.beta[2] == 0.0); const GeodesicTraceConfig trace = {.coordinate_time_step = 0.25, .max_steps = 200}; const ObserverState observer = observer_fixed_at_origin(); const RayEndpoint ray = geodesic_trace_past( &flat, &observer, (double[]){1, 0, 0}, &trace); CHECK(ray.outcome == RAY_OUTCOME_ESCAPED); CHECK(fabs(ray.n_infinity[0]) < 1e-12); CHECK(fabs(ray.n_infinity[1]) < 1e-12); CHECK(fabs(ray.n_infinity[2] + 1.0) < 1e-12); CHECK(fabs(ray.frequency_ratio - 1.0) < 1e-12); spacetime_destroy(&flat); } /* Reflection symmetry at fixed t: invariant under y -> -y, so transverse * beta derivatives and K components flip sign. */ { MetricData mirrored; CHECK(eval(&source, t, (double[]){x[0], -x[1], x[2]}, &mirrored) == 0); CHECK(fabs(metric.beta[0] - mirrored.beta[0]) < 1e-15); CHECK(fabs(metric.d_beta[0][0] - mirrored.d_beta[0][0]) < 1e-14); CHECK(fabs(metric.d_beta[1][0] + mirrored.d_beta[1][0]) < 1e-14); CHECK(fabs(metric.K[0][1] + mirrored.K[0][1]) < 1e-14); CHECK(fabs(metric.K[0][0] - mirrored.K[0][0]) < 1e-14); } /* Near-luminal bubble: a photon that propagates along +x with the bubble * separates from its center at only 1 - |v_s| and, traced backwards, meets * the bubble again near t ~ -15/(1-v_s) = -15000. It must still reach the * escape sphere; with a fixed 2^18 budget it would end in MAX_STEPS. The * budget below is the one main.c derives: 1.25 * 4*escape/((1-|v_s|)*step) * = 1.25 * 4*25/(0.001*0.05) = 2.5e6. */ { SpacetimeSource fast = {0}; CHECK(spacetime_create_alcubierre(&fast, 0.999, radius, sigma) == 0); ObserverCamera camera = {.position = {15.0, 0.0, 0.0}, .look_ra_deg = 0.0, .look_dec_deg = 0.0}; MetricData camera_metric; ObserverState observer; CHECK(eval(&fast, 0.0, camera.position, &camera_metric) == 0); CHECK(observer_from_coordinate_camera(&camera_metric, &camera, &observer, NULL) == OBSERVER_BUILD_OK); const GeodesicTraceConfig trace = {.coordinate_time_step = 0.05, .max_steps = 2500000u}; const RayEndpoint ray = geodesic_trace_past( &fast, &observer, (double[]){-1, 0, 0}, &trace); CHECK(ray.outcome == RAY_OUTCOME_ESCAPED); spacetime_destroy(&fast); } /* Refinement convergence: a ray grazing the bubble wall must converge in * n_infinity as the coordinate step is halved. */ { ObserverCamera camera = {.position = {-15.0, 0.0, 0.0}, .look_ra_deg = 0.0, .look_dec_deg = 0.0}; MetricData camera_metric; ObserverState observer; CHECK(eval(&source, 0.0, camera.position, &camera_metric) == 0); CHECK(observer_from_coordinate_camera(&camera_metric, &camera, &observer, NULL) == OBSERVER_BUILD_OK); double n[3] = {0.9995, 0.0316, 0.0}; { const double norm = sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]); for (int k = 0; k < 3; ++k) n[k] /= norm; } RayEndpoint previous = {0}; double previous_error = INFINITY; for (int level = 0; level < 3; ++level) { const GeodesicTraceConfig trace = { .coordinate_time_step = 0.08 / (1 << level), .max_steps = 1u << 20}; const RayEndpoint ray = geodesic_trace_past(&source, &observer, n, &trace); CHECK(ray.outcome == RAY_OUTCOME_ESCAPED); if (level > 0) { double error = 0.0; for (int k = 0; k < 3; ++k) { const double difference = ray.n_infinity[k] - previous.n_infinity[k]; error += difference * difference; } error = sqrt(error); CHECK(error <= previous_error); previous_error = error; } if (level == 2) CHECK(previous_error < 1e-5); previous = ray; } } /* Production DP54 axial superluminal check. A comoving bubble-center camera * at R = 1, sigma = 1 (escape radius 21) sees Pi_x = +-1 along the bubble * axis. The shared camera-relative dark policy must fire at threshold 8 for * the direction along the bubble motion and the opposite direction must * escape. The vs = +-2 constants were computed independently with mpmath; * this test has no dependency on any local experiment fixture. */ { const double R1 = 1.0, sig1 = 1.0; static const double vlist[] = {1.0, -1.0, 2.0, -2.0, 0.9999}; for (size_t k = 0; k < sizeof vlist / sizeof vlist[0]; ++k) { const double v = vlist[k]; const double sgn = v > 0.0 ? 1.0 : -1.0; SpacetimeSource fast = {0}; CHECK(spacetime_create_alcubierre(&fast, v, R1, sig1) == 0); ObserverCamera cam = {.coordinate_time = 0.0, .position = {0.0, 0.0, 0.0}, .velocity = {v, 0.0, 0.0}, .look_ra_deg = 0.0, .look_dec_deg = 0.0, .roll_deg = 0.0}; MetricData m; CHECK(eval(&fast, 0.0, cam.position, &m) == 0); ObserverState o; CHECK(observer_from_coordinate_camera(&m, &cam, &o, NULL) == OBSERVER_BUILD_OK); /* A coordinate-static center camera is timelike only for |v| < 1. */ ObserverCamera stat = cam; stat.velocity[0] = stat.velocity[1] = stat.velocity[2] = 0.0; ObserverState so; const int static_ok = observer_from_coordinate_camera(&m, &stat, &so, NULL) == OBSERVER_BUILD_OK; CHECK(static_ok == (fabs(v) < 1.0)); /* An outer static camera in the flat exterior is always legal. */ { const double pos[3] = {26.0, 0.0, 0.0}; MetricData om; ObserverCamera oc = {.coordinate_time = 0.0, .position = {26.0, 0.0, 0.0}, .look_ra_deg = 0.0, .look_dec_deg = 0.0, .roll_deg = 0.0}; CHECK(eval(&fast, 0.0, pos, &om) == 0); ObserverState oo; CHECK(observer_from_coordinate_camera(&om, &oc, &oo, NULL) == OBSERVER_BUILD_OK); } const GeodesicTraceConfig trace = { .coordinate_time_step = 0.05, .max_steps = 100000u, .threshold = {.kind = THRESHOLD_LOG_ENERGY_GROWTH, .value = 8.0, .policy_version = 1}, .stepper = GEODESIC_STEPPER_DP54, .atol_x = 1e-9, .atol_Pi = 1e-9, .atol_L = 1e-9, .rtol = 1e-9, .min_step = 1e-12, .max_step = 0.4, .consecutive_rejection_limit = 32, .max_lookback_time = 30000.0}; const double n_dark[3] = {-sgn, 0.0, 0.0}; const double n_esc[3] = {sgn, 0.0, 0.0}; GeodesicRayState idark, iesc; CHECK(geodesic_initialize_past_ray_metric(&m, &o, n_dark, &idark) == 0); CHECK(geodesic_initialize_past_ray_metric(&m, &o, n_esc, &iesc) == 0); printf("alcubierre vs=%.6g dark Pi_x=%.17g escape Pi_x=%.17g\n", v, idark.Pi[0], iesc.Pi[0]); CHECK(sgn * idark.Pi[0] > 0.999 && sgn * idark.Pi[0] < 1.000000001); CHECK(sgn * iesc.Pi[0] < -0.999 && sgn * iesc.Pi[0] > -1.000000001); const RayEndpoint dark = geodesic_trace_past(&fast, &o, n_dark, &trace); CHECK(dark.outcome == RAY_OUTCOME_DARK); CHECK(isfinite(dark.stop_coordinate_time)); CHECK(fabs(dark.threshold_value - 8.0) < 1e-6); CHECK(fabs((dark.final_log_alpha_p0 - dark.final_log_alpha_p0_0) - 8.0) < 1e-6); if (fabs(v) == 2.0) { const double q = dark.final_x[0] - v * dark.stop_coordinate_time; CHECK(fabs(fabs(q) - 1.2181434100155241) < 1e-6); CHECK(fabs(dark.stop_coordinate_time + 7.09915163394274) < 1e-6); } const RayEndpoint esc = geodesic_trace_past(&fast, &o, n_esc, &trace); CHECK(esc.outcome == RAY_OUTCOME_ESCAPED); CHECK(isfinite(esc.frequency_ratio) && esc.frequency_ratio > 0.0); if (fabs(v) == 2.0) CHECK(fabs(esc.frequency_ratio - 3.0) < 1e-6); for (int i = 0; i < 3; ++i) CHECK(isfinite(esc.n_infinity[i])); spacetime_destroy(&fast); } } spacetime_destroy(&source); puts("alcubierre regression passed"); return 0; }