diff --git a/Makefile b/Makefile index cd00a00..0dbd096 100644 --- a/Makefile +++ b/Makefile @@ -29,14 +29,14 @@ else IMAGE_EXT := ppm endif -COMMON_SOURCES := $(filter-out src/main.c src/dummy_psf.c src/fast_psf_fftw.c src/spacetime_minkowski.c src/spacetime_schwarzschild.c,$(wildcard src/*.c)) +COMMON_SOURCES := $(filter-out src/main.c src/dummy_psf.c src/fast_psf_fftw.c src/spacetime_minkowski.c src/spacetime_schwarzschild.c src/spacetime_alcubierre.c,$(wildcard src/*.c)) PROVIDER_SOURCE := src/spacetime_$(SPACETIME).c BUILD_DIR := build/$(BUILD_TYPE) TARGET_BASENAME := $(SPACETIME)_sky OBJECT_DIR := $(BUILD_DIR)/obj/$(SPACETIME) CORE_MINKOWSKI_SOURCES := $(COMMON_SOURCES) src/spacetime_minkowski.c -.PHONY: all backend clean run test tone-map-test sensor-bloom-test sensor-bloom-bench hip-psf-test hip-psf-bench fast-psf-fftw-bench minkowski schwarzschild FORCE +.PHONY: all backend clean run test tone-map-test sensor-bloom-test sensor-bloom-bench hip-psf-test hip-psf-bench fast-psf-fftw-bench minkowski schwarzschild alcubierre FORCE ifneq ($(filter 0 1,$(PSF_EVENT_SINK)),$(PSF_EVENT_SINK)) $(error Unknown PSF_EVENT_SINK '$(PSF_EVENT_SINK)'; choose 0 or 1) @@ -85,8 +85,10 @@ ifeq ($(SPACETIME),minkowski) BACKEND_CPPFLAGS := -DSPACETIME_MINKOWSKI else ifeq ($(SPACETIME),schwarzschild) BACKEND_CPPFLAGS := -DSPACETIME_SCHWARZSCHILD +else ifeq ($(SPACETIME),alcubierre) +BACKEND_CPPFLAGS := -DSPACETIME_ALCUBIERRE else -$(error Unknown SPACETIME '$(SPACETIME)'; choose minkowski or schwarzschild) +$(error Unknown SPACETIME '$(SPACETIME)'; choose minkowski, schwarzschild, or alcubierre) endif ifeq ($(ENABLE_HDR),1) @@ -106,6 +108,7 @@ TEST_OUT_DIR := $(OBJECT_DIR)/$(HDR_BUILD_TAG) TEST_TARGET := $(TEST_OUT_DIR)/test_geodesic FRAME_TEST_TARGET := $(TEST_OUT_DIR)/test_frame SCHWARZSCHILD_TEST_TARGET := $(TEST_OUT_DIR)/test_schwarzschild +ALCUBIERRE_TEST_TARGET := $(TEST_OUT_DIR)/test_alcubierre OBSERVER_TRACK_TEST_TARGET := $(TEST_OUT_DIR)/test_observer_track CATALOG_PREFETCH_TEST_TARGET := $(TEST_OUT_DIR)/test_catalog_prefetch HIP_PSF_TEST_TARGET := $(TEST_OUT_DIR)/test_hip_psf @@ -142,7 +145,7 @@ RENDER_DEPS := $(RENDER_OBJECTS:.o=.d) ifneq ($(filter command\ line environment environment\ override,$(origin SPACETIME)),) all: backend else -all: minkowski schwarzschild +all: minkowski schwarzschild alcubierre endif minkowski: @@ -151,6 +154,9 @@ minkowski: schwarzschild: $(MAKE) SPACETIME=schwarzschild ENABLE_HDR=$(ENABLE_HDR) backend +alcubierre: + $(MAKE) SPACETIME=alcubierre ENABLE_HDR=$(ENABLE_HDR) backend + # Build exactly the selected backend/configuration, e.g. # make SPACETIME=schwarzschild ENABLE_HDR=1 backend backend: $(TARGET) @@ -203,6 +209,9 @@ $(FRAME_TEST_TARGET): tests/test_frame.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SO $(SCHWARZSCHILD_TEST_TARGET): tests/test_schwarzschild.c $(COMMON_SOURCES) src/spacetime_schwarzschild.c $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ +$(ALCUBIERRE_TEST_TARGET): tests/test_alcubierre.c $(COMMON_SOURCES) src/spacetime_alcubierre.c $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) + $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ + $(OBSERVER_TRACK_TEST_TARGET): tests/test_observer_track.c $(COMMON_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ @@ -244,12 +253,13 @@ FAST_PSF_FFTW_TEST_DEP := FAST_PSF_FFTW_TEST_RUN := endif -test: $(CAMERA_TEST_TARGETS) $(TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET) $(FAST_PSF_FFTW_TEST_DEP) $(TONE_MAP_TEST_TARGET) $(SENSOR_BLOOM_TEST_TARGET) +test: $(CAMERA_TEST_TARGETS) $(TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(ALCUBIERRE_TEST_TARGET) $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET) $(FAST_PSF_FFTW_TEST_DEP) $(TONE_MAP_TEST_TARGET) $(SENSOR_BLOOM_TEST_TARGET) $(TEST_OUT_DIR)/test_observer_minkowski $(TEST_OUT_DIR)/test_observer_schwarzschild $(TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) + $(ALCUBIERRE_TEST_TARGET) $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET) $(FAST_PSF_FFTW_TEST_RUN) diff --git a/README.md b/README.md index e976ded..263eb9d 100644 --- a/README.md +++ b/README.md @@ -28,8 +28,9 @@ the horizon. See [single-frame camera parameters and complete commands](usage.md ## Current status -The current implementation supports analytic **Minkowski** and -**Schwarzschild** spacetimes, single images and observer-track image sequences, +The current implementation supports analytic **Minkowski**, +**Schwarzschild**, and moving **Alcubierre** warp-bubble spacetimes, +single images and observer-track image sequences, adaptive lens meshes, and reusable lens-map files. It is written primarily in C with OpenMP CPU parallelism; an optional HIP backend accelerates PSF accumulation. @@ -52,8 +53,8 @@ libpng development files. From the repository root: make -j ``` -This builds both `build/Release/minkowski_sky` and -`build/Release/schwarzschild_sky`. For individual backends, Debug builds, +This builds `build/Release/minkowski_sky`, `build/Release/schwarzschild_sky`, +and `build/Release/alcubierre_sky`. For individual backends, Debug builds, optional HDR/FITS output, HIP support, and regression checks, see [build.md](build.md). diff --git a/build.md b/build.md index 1b42fcb..e27a4e6 100644 --- a/build.md +++ b/build.md @@ -24,18 +24,20 @@ Optional dependencies are CFITSIO for HDR/FITS output, and HIP/ROCm with make -j ``` -With no explicit `SPACETIME` setting, this builds both supported spacetimes: +With no explicit `SPACETIME` setting, this builds all supported spacetimes: | Executable | Spacetime | | --- | --- | | `build/Release/minkowski_sky` | Flat Minkowski spacetime | | `build/Release/schwarzschild_sky` | Analytic Schwarzschild in ingoing Kerr–Schild coordinates | +| `build/Release/alcubierre_sky` | Analytic moving Alcubierre warp bubble, `x_s(t)=v_s t` (no capture) | To build only one: ```sh make -j SPACETIME=minkowski backend make -j SPACETIME=schwarzschild backend +make -j SPACETIME=alcubierre backend ``` Each executable contains one metric provider, selected at compile time. diff --git a/src/main.c b/src/main.c index 253f95a..4e253e2 100644 --- a/src/main.c +++ b/src/main.c @@ -57,6 +57,7 @@ typedef struct { double movie_start_time, movie_duration, movie_fps; double slab_duration; double minkowski_proper_acceleration; + double alcubierre_vs, alcubierre_radius, alcubierre_sigma; int catalog_load_workers; const char *blackbody_table_path; RefinementConfig refinement; @@ -340,6 +341,9 @@ static int parse_args(int argc, char **argv, Settings *s, .movie_fps = 30.0, .slab_duration = 64.0, .minkowski_proper_acceleration = 1.52, + .alcubierre_vs = 0.5, + .alcubierre_radius = 5.0, + .alcubierre_sigma = 1.0, .catalog_load_workers = 4, .refinement = {.angle_absolute_rad = 1e-3 * 3.14159265358979323846 / 180.0, @@ -348,6 +352,11 @@ static int parse_args(int argc, char **argv, Settings *s, .min_edge_pixels = 0.5, .min_area_pixels2 = 0.25}, .tone_map = {.op = TONE_MAP_SOFTCLIP, .p = 2.0}}; +#ifdef SPACETIME_ALCUBIERRE + /* The default bubble (R=5, sigma=1) has escape radius 25, so the generic + * radius-30 camera would sit outside the active domain. */ + s->observer_radius = 15.0; +#endif *write_path = NULL; int tone_map_p_specified = 0; for (int i = 1; i < argc; ++i) { @@ -481,6 +490,14 @@ static int parse_args(int argc, char **argv, Settings *s, } else if (!strcmp(argv[i], "--proper-acceleration") && i + 1 < argc && !parse_nonnegative(argv[++i], &s->minkowski_proper_acceleration)) { +#ifdef SPACETIME_ALCUBIERRE + } else if (!strcmp(argv[i], "--alcubierre-vs") && i + 1 < argc && + !parse_finite(argv[++i], &s->alcubierre_vs)) { + } else if (!strcmp(argv[i], "--alcubierre-radius") && i + 1 < argc && + !parse_positive(argv[++i], &s->alcubierre_radius)) { + } else if (!strcmp(argv[i], "--alcubierre-sigma") && i + 1 < argc && + !parse_positive(argv[++i], &s->alcubierre_sigma)) { +#endif } else if (!strcmp(argv[i], "--catalog-load-workers") && i + 1 < argc && !parse_int(argv[++i], &s->catalog_load_workers)) { } else if (!strcmp(argv[i], "--blackbody-table") && i + 1 < argc) { @@ -555,6 +572,8 @@ static void print_help(const char *program) { " Look is projected into the moving camera rest space.\n" #ifdef SPACETIME_SCHWARZSCHILD " Default position: (0,0,30); look RA=90, Dec=-90.\n" +#elif defined(SPACETIME_ALCUBIERRE) + " Default position: (0,0,15); look RA=90, Dec=-90.\n" #else " Default position: (0,0,0); look RA=90, Dec=-90.\n" #endif @@ -592,6 +611,16 @@ static void print_help(const char *program) { #else fputs(" --draw-mesh Also write the final lens-mesh overlay as _mesh.ppm\n", stdout); +#endif +#ifdef SPACETIME_ALCUBIERRE + fputs( + "\nAlcubierre warp bubble (moving x_s(t)=v_s*t; no capture):\n" + " --alcubierre-vs V Constant bubble velocity v_s, |v_s| < 1 (default: 0.5)\n" + " --alcubierre-radius R Bubble radius R > 0 (default: 5)\n" + " --alcubierre-sigma S Wall sharpness sigma > 0 (default: 1)\n" + " The escape radius R + 20/sigma is derived internally; the\n" + " camera must lie inside it.\n", + stdout); #endif fputs( "\nMovie and observer track:\n" @@ -713,12 +742,55 @@ static void report_frame_refinement(void *context, size_t generation, generation, added_vertices, vertex_count, triangle_count); } -static GeodesicTraceConfig trace_config(void) { +#ifdef SPACETIME_ALCUBIERRE +/* Upper bound on the per-ray step budget. Legal parameters whose worst-case + * near-comoving ray could need more than this are rejected at startup rather + * than silently terminating as RAY_ENDPOINT_MAX_STEPS. */ +#define ALCUBIERRE_MAX_TRACE_STEPS (1u << 24) + +/* Safety margin over the straight-line worst case: wall-region deflection can + * make a ray linger, and 1 - |v_s| is only the asymptotic separation rate. */ +#define ALCUBIERRE_BUDGET_MARGIN 1.25 + +static double alcubierre_time_step(const Settings *s) { + /* Resolve the wall transition ~1/sigma. */ + return fmin(0.1, 0.05 / s->alcubierre_sigma); +} + +/* Worst-case per-ray step budget, including ALCUBIERRE_BUDGET_MARGIN. A ray + * that is nearly comoving with the bubble separates from its center in the + * propagation direction at only ~1 - |v_s|, so crossing the ~4*escape domain + * can take ~4*escape/(1-|v_s|) in coordinate time. trace_config() and the + * startup rejection share this single value so the configured limit always + * carries the full margin when it is accepted. */ +static double alcubierre_step_budget(const Settings *s) { + const double escape = + spacetime_alcubierre_escape_radius(s->alcubierre_radius, + s->alcubierre_sigma); + const double separation = 1.0 - fabs(s->alcubierre_vs); + return ALCUBIERRE_BUDGET_MARGIN * 4.0 * escape / + (separation * alcubierre_time_step(s)); +} +#endif + +static GeodesicTraceConfig trace_config(const Settings *s) { #ifdef SPACETIME_SCHWARZSCHILD + (void)s; return (GeodesicTraceConfig){.coordinate_time_step = 0.1, .max_steps = 4096, .capture_log_alpha_p0 = 8.0}; +#elif defined(SPACETIME_ALCUBIERRE) + const double step = alcubierre_time_step(s); + const double budget = alcubierre_step_budget(s); + unsigned max_steps = ALCUBIERRE_MAX_TRACE_STEPS; + if (budget < (double)max_steps && isfinite(budget)) + max_steps = (unsigned)ceil(budget); + if (max_steps < 1024u) + max_steps = 1024u; + return (GeodesicTraceConfig){.coordinate_time_step = step, + .max_steps = max_steps}; #else + (void)s; return (GeodesicTraceConfig){.coordinate_time_step = 1.0, .max_steps = 2048}; #endif @@ -756,7 +828,7 @@ static int resolve_camera(Settings *s) { s->look_dec_deg = atan2(-x[2], hypot(x[0], x[1])) * degrees; } int infer_position = s->look_specified || s->radius_specified; -#ifdef SPACETIME_SCHWARZSCHILD +#if defined(SPACETIME_SCHWARZSCHILD) || defined(SPACETIME_ALCUBIERRE) infer_position = 1; #endif if (!s->position_specified && infer_position) { @@ -777,11 +849,17 @@ static int build_observer(const Settings *s, const SpacetimeSource *spacetime, camera.position[i] = s->observer_position[i]; camera.velocity[i] = s->observer_velocity[i]; } - if (spacetime_classify(spacetime, camera.coordinate_time, camera.position) == - SPACETIME_RAY_CAPTURED) { + const SpacetimeRayStatus camera_status = + spacetime_classify(spacetime, camera.coordinate_time, camera.position); + if (camera_status == SPACETIME_RAY_CAPTURED) { fputs("Camera position is inside the backend capture cutoff or invalid.\n", stderr); return -1; } + if (camera_status == SPACETIME_RAY_ESCAPED) { + fputs("Camera position is outside this backend's finite escape radius; " + "move the camera inward or enlarge the spacetime domain.\n", stderr); + return -1; + } MetricData metric; if (spacetime_eval(spacetime, camera.coordinate_time, camera.position, &metric)) { fputs("Could not evaluate metric at the camera event.\n", stderr); @@ -819,7 +897,7 @@ static int render_observer_frame(const Settings *s, StarCatalog *catalog, const SpacetimeSource *spacetime, const ObserverState *observer, const char *output_path) { - const GeodesicTraceConfig trace = trace_config(); + const GeodesicTraceConfig trace = trace_config(s); FrameLensMesh mesh = {0}; double *hdr = calloc((size_t)s->width * s->height * 3, sizeof *hdr); if (hdr == NULL || frame_lens_mesh_build_coarse(&mesh, s->width, s->height, @@ -1028,7 +1106,7 @@ static int render_movie(const Settings *s, StarCatalog *catalog, const SpacetimeSource *spacetime) { ObserverTrack track = {0}; Movie movie = {0}; - const GeodesicTraceConfig trace = trace_config(); + const GeodesicTraceConfig trace = trace_config(s); int result = -1; if (s->observer_track_path == NULL || observer_track_load_csv(&track, s->observer_track_path) || @@ -1341,10 +1419,39 @@ int main(int argc, char **argv) { SpacetimeSource spacetime = {0}; ObserverState observer; if (settings.lens_map_input_path == NULL) { +#ifdef SPACETIME_ALCUBIERRE + if (spacetime_create_alcubierre(&spacetime, settings.alcubierre_vs, + settings.alcubierre_radius, + settings.alcubierre_sigma)) { + fputs("Could not create Alcubierre spacetime source; require |v_s| < 1, " + "R > 0, sigma > 0.\n", stderr); + return 1; + } + if (alcubierre_step_budget(&settings) > + (double)ALCUBIERRE_MAX_TRACE_STEPS) { + /* The budget has a V-shaped minimum at sigma = 0.5, where the step + * stops being capped: below it the 20/sigma term dominates (increase + * sigma helps), above it the step scales as 1/sigma (decrease sigma + * helps), and at exactly 0.5 neither direction improves anything. */ + const char *sigma_advice = ""; + if (settings.alcubierre_sigma > 0.5) + sigma_advice = "decrease --alcubierre-sigma, "; + else if (settings.alcubierre_sigma < 0.5) + sigma_advice = "increase --alcubierre-sigma, "; + fprintf(stderr, + "Alcubierre trace budget exceeds the %u-step cap; decrease " + "--alcubierre-radius, %sor move --alcubierre-vs away from " + "+/-1.\n", + ALCUBIERRE_MAX_TRACE_STEPS, sigma_advice); + spacetime_destroy(&spacetime); + return 2; + } +#else if (spacetime_create_default(&spacetime)) { fputs("Could not create spacetime source\n", stderr); return 1; } +#endif if (settings.frames_dir == NULL && build_observer(&settings, &spacetime, &observer)) { spacetime_destroy(&spacetime); diff --git a/src/spacetime.h b/src/spacetime.h index 7fccebd..96def36 100644 --- a/src/spacetime.h +++ b/src/spacetime.h @@ -57,6 +57,13 @@ int spacetime_create_minkowski(SpacetimeSource *source, double escape_radius); int spacetime_create_schwarzschild_ks(SpacetimeSource *source, double mass, double escape_radius, double capture_radius); +/* Moving Alcubierre bubble with x_s(t) = vs*t and x_s(0) = 0. Requires + * |vs| < 1, R > 0, and sigma > 0. */ +int spacetime_create_alcubierre(SpacetimeSource *source, double vs, + double radius, double sigma); +/* Bubble-centered escape radius used by the Alcubierre backend; also lets + * callers size their integration step budget. */ +double spacetime_alcubierre_escape_radius(double radius, double sigma); void spacetime_destroy(SpacetimeSource *source); int spacetime_eval(const SpacetimeSource *source, double t, const double x[3], MetricData *metric); diff --git a/src/spacetime_alcubierre.c b/src/spacetime_alcubierre.c new file mode 100644 index 0000000..8e739c0 --- /dev/null +++ b/src/spacetime_alcubierre.c @@ -0,0 +1,164 @@ +#include "spacetime.h" + +#include +#include +#include + +/* 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 int 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 -1; + 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 0; +} + +/* A warp bubble has no curvature singularity or horizon for |v_s| < 1, so + * rays are only ever ACTIVE or ESCAPED; the exotic matter that would source + * the bubble is treated as optically transparent. The escape sphere follows + * the bubble, so rays terminate only once the metric is flat to machine + * precision at their current location. */ +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 const SpacetimeOps alcubierre_ops = { + .eval = alcubierre_eval, + .classify = alcubierre_classify, + .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) || fabs(vs) >= 1.0 || + !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; + return 0; +} + +int spacetime_create_default(SpacetimeSource *source) { + return spacetime_create_alcubierre(source, 0.5, 5.0, 1.0); +} diff --git a/tests/test_alcubierre.c b/tests/test_alcubierre.c new file mode 100644 index 0000000..8bf1716 --- /dev/null +++ b/tests/test_alcubierre.c @@ -0,0 +1,328 @@ +#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); + CHECK(spacetime_create_alcubierre(&source, 1.0, radius, sigma) != 0); + CHECK(spacetime_create_alcubierre(&source, -1.5, 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.status == RAY_ENDPOINT_ESCAPED); + CHECK(r1.status == RAY_ENDPOINT_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.status == RAY_ENDPOINT_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.status == RAY_ENDPOINT_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.status == RAY_ENDPOINT_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; + } + } + + spacetime_destroy(&source); + puts("alcubierre regression passed"); + return 0; +} diff --git a/usage.md b/usage.md index e5d9086..a82acfd 100644 --- a/usage.md +++ b/usage.md @@ -92,6 +92,59 @@ The last example points outward from a camera moving inward inside the horizon. It needs no CSV trajectory or movie wrapper. These small images are camera checks; increase resolution and refinement for production renders. +## Alcubierre warp-bubble spacetime + +`alcubierre_sky` renders the moving Alcubierre line element + +$$ds^2 = -dt^2 + \left[dx - v_s f(r_s)\,dt\right]^2 + dy^2 + dz^2,$$ + +with the bubble center following the constant-velocity worldline +`x_s(t) = v_s t` (a lab/non-comoving slicing fixed by `x_s(0) = 0`) and + +$$f(r) = \frac{\tanh(\sigma(r+R)) - \tanh(\sigma(r-R))}{2\tanh(\sigma R)},\qquad +r_s = \sqrt{(x-x_s)^2 + y^2 + z^2}.$$ + +The bubble therefore propagates through the coordinates, and the metric is +time-dependent: the renderer evaluates `f(r_s)` and its spatial derivatives at +each coordinate time, while the extrinsic curvature supplies the required +`d_t gamma` information to the 3+1 null-ray equations. The exotic matter that +would source the bubble is treated as optically transparent, so there is +**no capture**: rays are only active or escaped. This is why the backend +requires a sub-luminal `|v_s| < 1`; at or above `1` the metric develops an +ergoregion/event horizon and static observers cease to exist, which is outside +the current no-capture scope. + +| Option | Meaning / default | +| --- | --- | +| `--alcubierre-vs V` | Constant shift parameter, `|V| < 1` (default 0.5) | +| `--alcubierre-radius R` | Bubble radius `R > 0` (default 5) | +| `--alcubierre-sigma S` | Wall sharpness `S > 0` (default 1) | + +`f` decays to zero past `r_s = R` over a transition width `~1/sigma`, so the +finite escape sphere is bubble-centered with radius `R + 20/sigma` and needs no +CLI option; it follows the moving bubble, so rays terminate only once the local +metric is flat to below double precision. The single-frame camera default is +`(0,0,15)` at `t = 0`, when the bubble is still at the origin; it must lie +inside the escape sphere, or the observer build fails with an explicit error. +The per-ray step budget scales with the escape radius and `1/(1-|v_s|)`, so +near-luminal `v_s` still lets grazing rays escape; combinations whose +worst-case budget would exceed the internal cap are rejected at startup. + +Lensing and frequency shifts come from the bubble wall. The configuration is +invariant under the isometry `(t, x) -> (t + T, x + v_s T)`, so observers +related by it see identical escaping directions and frequency ratios. For +example: + +```sh +make -j PSF_BACKEND=cpu SPACETIME=alcubierre backend +./build/Release/alcubierre_sky --catalog assets/sky_grid_5deg.csv \ + --observer-radius 15 --look-ra-deg 90 --look-dec-deg -90 \ + --alcubierre-vs 0.5 --alcubierre-radius 5 --alcubierre-sigma 1 \ + --width 640 --height 360 --fov-deg 60 --exposure 1 \ + --coarse-cell-pixels 16 --refine-max-level 2 --psf-direct \ + --output output/imgs/alcubierre_wall.png +``` + ## Movie image sequences Movie mode reads an observer-track CSV containing coordinate time, proper