wyj/GR-raytracing
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Feat: Add directed asymptotic escape and analytic Schwarzschild exterior

Browse Source 
Replace radius-only escape termination with a common asymptotic exterior protocol: declared ends, moving escape worldtubes, directed inside->outside crossings, and a PENDING_ENTRY lifecycle shared by single-frame and movie tracing.

Add an analytic Carlson-integral Schwarzschild monopole exterior (angle primitive, bracketed turning radius, ingoing Kerr-Schild coordinate-time transfer, conserved-energy frequency) so a camera outside the escape sphere is traced through an entry event.

Make the lifecycle tri-state (no ends / ready / protocol error), carry end_id through the endpoint and lens mesh, validate sources in constructors via spacetime_source_finalize(), and refresh the Schwarzschild reference images for the corrected finish.
This commit is contained in:
wyj committed 2026-10-04 03:32:28 -04:00
1 parent 04611e3e5a
commit 09a7417961
24 files changed
+3566 -58
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

src/asymptotic.c
+710
View File
@@ -0,0 +1,710 @@
#include "asymptotic.h"
#include "asymptotic_schwarzschild.h"
#include <float.h>
#include <math.h>
#include <stddef.h>
static double dot3(const double a[3], const double b[3]) {
return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
static double normalize3(double v[3]) {
const double length = sqrt(dot3(v, v));
if (length > 0.0)
for (int i = 0; i < 3; ++i)
v[i] /= length;
return length;
}
static int invert3(double a[3][3], double inv[3][3]) {
const double det =
a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1]) -
a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0]) +
a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);
if (!isfinite(det) || fabs(det) < 1e-300)
return -1;
inv[0][0] = (a[1][1] * a[2][2] - a[1][2] * a[2][1]) / det;
inv[0][1] = (a[0][2] * a[2][1] - a[0][1] * a[2][2]) / det;
inv[0][2] = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) / det;
inv[1][0] = (a[1][2] * a[2][0] - a[1][0] * a[2][2]) / det;
inv[1][1] = (a[0][0] * a[2][2] - a[0][2] * a[2][0]) / det;
inv[1][2] = (a[0][2] * a[1][0] - a[0][0] * a[1][2]) / det;
inv[2][0] = (a[1][0] * a[2][1] - a[1][1] * a[2][0]) / det;
inv[2][1] = (a[0][1] * a[2][0] - a[0][0] * a[2][1]) / det;
inv[2][2] = (a[0][0] * a[1][1] - a[0][1] * a[1][0]) / det;
return 0;
}
static int find_end(const SpacetimeSource *source, SpacetimeEndId end_id,
SpacetimeAsymptoticEnd *out) {
const size_t count = spacetime_asymptotic_end_count(source);
for (size_t i = 0; i < count; ++i) {
SpacetimeAsymptoticEnd end;
if (spacetime_asymptotic_end(source, i, &end))
continue;
if (end.end_id == end_id) {
if (out != NULL)
*out = end;
return 0;
}
}
return -1;
}
/* Backend vector -> asymptotic-frame vector, where the frame axes are the
* columns of end->frame_axes expressed in backend coordinates. */
static void backend_vector_to_frame(const SpacetimeAsymptoticEnd *end,
const double a[3], double out[3]) {
for (int j = 0; j < 3; ++j)
out[j] = end->frame_axes[0][j] * a[0] + end->frame_axes[1][j] * a[1] +
end->frame_axes[2][j] * a[2];
}
static void frame_vector_to_backend(const SpacetimeAsymptoticEnd *end,
const double a[3], double out[3]) {
for (int i = 0; i < 3; ++i)
out[i] = end->frame_axes[i][0] * a[0] + end->frame_axes[i][1] * a[1] +
end->frame_axes[i][2] * a[2];
}
static void backend_position_to_frame(const SpacetimeAsymptoticEnd *end,
const double x[3], double out[3]) {
const double shifted[3] = {x[0] - end->frame_origin[0],
x[1] - end->frame_origin[1],
x[2] - end->frame_origin[2]};
backend_vector_to_frame(end, shifted, out);
}
static void frame_position_to_backend(const SpacetimeAsymptoticEnd *end,
const double x[3], double out[3]) {
double rotated[3];
frame_vector_to_backend(end, x, rotated);
for (int i = 0; i < 3; ++i)
out[i] = rotated[i] + end->frame_origin[i];
}
static AsymptoticStatus worldtube_sample(
const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end, double t,
SpacetimeEscapeWorldtubeSample *out) {
if (spacetime_escape_worldtube_sample(source, end->end_id, t, out))
return ASYMPTOTIC_INVALID;
if (!out->valid)
return ASYMPTOTIC_TIME_RANGE_EXHAUSTED;
/* A declared worldtube must be a finite, positive-radius sphere. */
if (!(out->radius > 0.0) || !isfinite(out->radius) ||
!isfinite(out->radius_rate))
return ASYMPTOTIC_INVALID;
for (int i = 0; i < 3; ++i)
if (!isfinite(out->center[i]) || !isfinite(out->velocity[i]))
return ASYMPTOTIC_INVALID;
return ASYMPTOTIC_OK;
}
/* Worldtube value F and its derivative dF/ds along the past direction `w`
* (unit past spatial velocity, s = t0 - t). On the boundary F == 0 the sign
* of the slope decides inside vs outside. */
static AsymptoticStatus worldtube_value_and_slope(
const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end, double t,
const double x[3], const double w[3], double *value, double *slope) {
SpacetimeEscapeWorldtubeSample sample;
const AsymptoticStatus status = worldtube_sample(source, end, t, &sample);
if (status != ASYMPTOTIC_OK)
return status;
double d[3], q[3];
for (int i = 0; i < 3; ++i) {
d[i] = x[i] - sample.center[i];
q[i] = w[i] + sample.velocity[i];
}
*value = dot3(d, d) - sample.radius * sample.radius;
*slope = 2.0 * dot3(d, q) + 2.0 * sample.radius * sample.radius_rate;
return ASYMPTOTIC_OK;
}
int asymptotic_worldtube_value(const SpacetimeSource *source,
SpacetimeEndId end_id, double t,
const double x[3], double *value) {
SpacetimeAsymptoticEnd end;
if (value == NULL || find_end(source, end_id, &end))
return ASYMPTOTIC_INVALID;
double slope;
return worldtube_value_and_slope(source, &end, t, x,
(const double[3]){0.0, 0.0, 0.0}, value,
&slope);
}
int asymptotic_canonical_from_backend(const SpacetimeSource *source,
SpacetimeEndId end_id,
const MetricData *metric, double t,
const double x[3], const double Pi[3],
double log_alpha_p0,
AsymptoticPhotonState *out) {
SpacetimeAsymptoticEnd end;
double inv[3][3], gamma[3][3];
if (out == NULL || metric == NULL || find_end(source, end_id, &end))
return -1;
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j)
gamma[i][j] = metric->gamma[i][j];
if (invert3(gamma, inv))
return -1;
if (end.exterior_kind != ASYMPTOTIC_EXTERIOR_MINKOWSKI)
return -1;
double w_backend[3] = {0.0, 0.0, 0.0};
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j)
w_backend[i] -= inv[i][j] * Pi[j];
if (normalize3(w_backend) <= 0.0)
return -1;
out->end_id = end_id;
out->t = t;
out->log_alpha_p0 = log_alpha_p0;
backend_position_to_frame(&end, x, out->x);
backend_vector_to_frame(&end, w_backend, out->w);
return 0;
}
int asymptotic_backend_from_canonical(const SpacetimeSource *source,
const MetricData *metric,
const AsymptoticPhotonState *canonical,
double x[3], double Pi[3],
double *log_alpha_p0) {
SpacetimeAsymptoticEnd end;
(void)metric;
if (canonical == NULL || x == NULL || Pi == NULL ||
find_end(source, canonical->end_id, &end))
return -1;
if (end.exterior_kind != ASYMPTOTIC_EXTERIOR_MINKOWSKI)
return -1;
double w_backend[3];
frame_position_to_backend(&end, canonical->x, x);
frame_vector_to_backend(&end, canonical->w, w_backend);
/* Flat exterior: the covariant momentum is the unit past direction negated. */
for (int i = 0; i < 3; ++i)
Pi[i] = -w_backend[i];
if (log_alpha_p0 != NULL)
*log_alpha_p0 = canonical->log_alpha_p0;
return 0;
}
/* Solve |d + q s|^2 = (R0 - rr s)^2 for the smallest s >= 0 with outside ->
* inside crossing. Returns 1 on entry (sets s), 0 on miss, -1 on error. */
static int solve_entry_quadratic(const double d[3], const double q[3],
double R0, double rr, double *s_out) {
const double qq = dot3(q, q);
const double a = qq - rr * rr;
const double b = 2.0 * (dot3(d, q) + R0 * rr);
const double c = dot3(d, d) - R0 * R0;
if (c < 0.0) {
/* Strictly inside; the lifecycle normally handles this as INSIDE. */
*s_out = 0.0;
return 1;
}
if (c == 0.0) {
/* On the boundary: classify by dF/ds = b. Past-inward enters at once;
* outward/tangent rays may still re-enter later when the sphere shrinks
* (a < 0), so do not declare a permanent miss on the zero root. */
if (b < 0.0) {
*s_out = 0.0;
return 1;
}
if (b == 0.0) {
if (a < 0.0) {
*s_out = 0.0;
return 1;
}
return 0;
}
if (a < 0.0) {
*s_out = -b / a;
return 1;
}
return 0;
}
/* Compare the quadratic coefficient against the velocity-squared scale it
* is built from; mixing in R0^2 would let a large radius misclassify a
* genuinely quadratic entry as linear. */
const double scale = qq + rr * rr;
if (fabs(a) <= 32.0 * DBL_EPSILON * scale) {
if (!isfinite(b) || b >= 0.0)
return 0;
const double s = -c / b;
if (s <= 0.0)
return 0;
*s_out = s;
return 1;
}
const double disc = b * b - 4.0 * a * c;
if (!isfinite(disc) || disc <= 0.0)
return 0;
/* Numerically stable quadratic roots: q avoids cancellation in the root
* with the same sign as b, which is exactly the small entry root when the
* camera sits just outside a large sphere. */
const double root = sqrt(disc);
const double qq2 = -0.5 * (b + copysign(root, b));
const double r1 = qq2 / a;
const double r2 = c / qq2;
/* The first outside->inside crossing is the smallest positive root. */
double s = INFINITY;
if (r1 > 0.0)
s = r1;
if (r2 > 0.0 && r2 < s)
s = r2;
if (!(s < INFINITY))
return 0;
*s_out = s;
return 1;
}
static double worldtube_F_frame(const double c_frame[3], double radius,
const double x_frame[3]) {
const double d[3] = {x_frame[0] - c_frame[0], x_frame[1] - c_frame[1],
x_frame[2] - c_frame[2]};
return dot3(d, d) - radius * radius;
}
/* Bracketed first-entry search for a worldtube whose motion is not constant.
* This is the future interface for accelerated worldtubes; the current
* backends always take the closed quadratic path above. */
static AsymptoticStatus minkowski_generic_first_entry(
const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end,
double t0, const double x_frame[3], const double w_frame[3],
double *s_out) {
double s = 0.0;
SpacetimeEscapeWorldtubeSample sample;
AsymptoticStatus status = worldtube_sample(source, end, t0, &sample);
if (status != ASYMPTOTIC_OK)
return status;
double c_frame[3];
backend_position_to_frame(end, sample.center, c_frame);
const double f_start = worldtube_F_frame(c_frame, sample.radius, x_frame);
if (f_start < 0.0) {
*s_out = 0.0;
return ASYMPTOTIC_OK;
}
if (f_start == 0.0) {
/* On the boundary only an inward slope is an entry; outward and tangent
* rays keep searching. */
double v_frame[3], d0[3], q0[3];
backend_vector_to_frame(end, sample.velocity, v_frame);
for (int i = 0; i < 3; ++i) {
d0[i] = x_frame[i] - c_frame[i];
q0[i] = w_frame[i] + v_frame[i];
}
const double slope0 =
2.0 * dot3(d0, q0) + 2.0 * sample.radius * sample.radius_rate;
if (slope0 < 0.0) {
*s_out = 0.0;
return ASYMPTOTIC_OK;
}
}
const double speed = fabs(sample.velocity[0]) + fabs(sample.velocity[1]) +
fabs(sample.velocity[2]) + fabs(sample.radius_rate) + 1.0;
const double base_step = 0.5 * fmax(sample.radius, 1.0) / speed;
for (int iteration = 0; iteration < 1000000; ++iteration) {
double step = base_step;
const double boundary = spacetime_escape_worldtube_next_segment(
source, end->end_id, t0 - s);
if (isfinite(boundary)) {
/* The backward-integration distance to a past segment boundary. */
const double to_boundary = (t0 - boundary) - s;
if (to_boundary > 0.0)
step = fmin(step, to_boundary);
}
const double s_next = s + step;
const double t_next = t0 - s_next;
SpacetimeEscapeWorldtubeSample next;
status = worldtube_sample(source, end, t_next, &next);
if (status != ASYMPTOTIC_OK)
return status;
double c_next[3], ray_next[3];
backend_position_to_frame(end, next.center, c_next);
for (int i = 0; i < 3; ++i)
ray_next[i] = x_frame[i] + w_frame[i] * s_next;
const double f_next = worldtube_F_frame(c_next, next.radius, ray_next);
/* Only a strictly negative sample is an entry; a single touch at F == 0
* (tangent) is not. */
if (f_next < 0.0) {
double lo = s, hi = s_next;
AsymptoticStatus bisect_status = ASYMPTOTIC_OK;
for (int bisect = 0; bisect < 80; ++bisect) {
const double mid = 0.5 * (lo + hi);
SpacetimeEscapeWorldtubeSample mid_sample;
bisect_status = worldtube_sample(source, end, t0 - mid, &mid_sample);
if (bisect_status != ASYMPTOTIC_OK)
break;
double c_mid[3], ray_mid[3];
backend_position_to_frame(end, mid_sample.center, c_mid);
for (int i = 0; i < 3; ++i)
ray_mid[i] = x_frame[i] + w_frame[i] * mid;
const double f_mid =
worldtube_F_frame(c_mid, mid_sample.radius, ray_mid);
if (f_mid <= 0.0)
hi = mid;
else
lo = mid;
}
/* A sample failure inside the bracket must not be disguised as a
* normal entry. */
if (bisect_status != ASYMPTOTIC_OK)
return bisect_status;
*s_out = hi;
return ASYMPTOTIC_OK;
}
s = s_next;
}
return ASYMPTOTIC_INVALID;
}
static void minkowski_route_escaped(const SpacetimeAsymptoticEnd *end,
const double w_frame[3],
SpacetimeEndId end_id,
AsymptoticRoute *route) {
route->kind = ASYMPTOTIC_ROUTE_ESCAPED;
route->end_id = end_id;
double w_backend[3];
frame_vector_to_backend(end, w_frame, w_backend);
for (int i = 0; i < 3; ++i)
route->n_infinity[i] = w_backend[i];
}
static void minkowski_route_entry(const SpacetimeAsymptoticEnd *end, double t0,
const double x_frame[3],
const double w_frame[3], double s_entry,
SpacetimeEndId end_id,
AsymptoticRoute *route) {
double x_entry_frame[3];
for (int i = 0; i < 3; ++i)
x_entry_frame[i] = x_frame[i] + w_frame[i] * s_entry;
route->kind = ASYMPTOTIC_ROUTE_ENTRY;
route->end_id = end_id;
route->activate_t = t0 - s_entry;
frame_position_to_backend(end, x_entry_frame, route->x);
double w_backend[3];
frame_vector_to_backend(end, w_frame, w_backend);
for (int i = 0; i < 3; ++i)
route->Pi[i] = -w_backend[i];
}
static AsymptoticStatus minkowski_preroute(
const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end, double t0,
const double x_frame[3], const double w_frame[3], SpacetimeEndId end_id,
AsymptoticRoute *route) {
/* Walk constant-velocity motion segments. A quadratic root is only valid
* inside the current segment and while the radius stays positive; otherwise
* advance to the next segment boundary and re-sample. */
double s = 0.0;
for (int segment = 0; segment < 1000000; ++segment) {
const double t = t0 - s;
SpacetimeEscapeWorldtubeSample sample;
AsymptoticStatus status = worldtube_sample(source, end, t, &sample);
if (status != ASYMPTOTIC_OK)
return status;
double x_cur[3];
for (int i = 0; i < 3; ++i)
x_cur[i] = x_frame[i] + w_frame[i] * s;
if (!sample.velocity_constant) {
double s_rel;
status = minkowski_generic_first_entry(source, end, t, x_cur, w_frame,
&s_rel);
if (status != ASYMPTOTIC_OK)
return status;
minkowski_route_entry(end, t0, x_frame, w_frame, s + s_rel, end_id,
route);
return ASYMPTOTIC_OK;
}
double c_frame[3], v_frame[3], d[3], q[3];
backend_position_to_frame(end, sample.center, c_frame);
backend_vector_to_frame(end, sample.velocity, v_frame);
for (int i = 0; i < 3; ++i) {
d[i] = x_cur[i] - c_frame[i];
q[i] = w_frame[i] + v_frame[i];
}
const double boundary = spacetime_escape_worldtube_next_segment(
source, end->end_id, t);
const double s_segment = isfinite(boundary) ? (t - boundary) : INFINITY;
if (!(s_segment >= 0.0))
return ASYMPTOTIC_INVALID;
double sigma;
const int hit = solve_entry_quadratic(d, q, sample.radius,
sample.radius_rate, &sigma);
/* The backend constructor guarantees R > 0 throughout every segment, so
* a root inside the segment is a real entry. A root past the segment
* boundary is not adopted here; the next segment is sampled instead.
* The cheap R > 0 test at the root guards against a backend that
* bypasses its constructor. */
if (hit && sigma >= 0.0 && sigma <= s_segment) {
if (sample.radius - sample.radius_rate * sigma <= 0.0)
return ASYMPTOTIC_INVALID;
minkowski_route_entry(end, t0, x_frame, w_frame, s + sigma, end_id,
route);
return ASYMPTOTIC_OK;
}
if (!isfinite(s_segment)) {
/* Open final segment with no entry: a genuine miss. */
minkowski_route_escaped(end, w_frame, end_id, route);
return ASYMPTOTIC_OK;
}
const double next_s = t0 - nextafter(boundary, -INFINITY);
if (!(next_s > s))
return ASYMPTOTIC_INVALID;
s = next_s;
}
/* Segment budget exhausted without a conclusion: never disguise this as an
* escape. */
return ASYMPTOTIC_INVALID;
}
static AsymptoticStatus schwarzschild_route(
const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end,
const MetricData *metric, const GeodesicRayState *state,
AsymptoticRoute *route) {
SpacetimeEscapeWorldtubeSample sample;
const AsymptoticStatus sample_status = worldtube_sample(
source, end, state->coordinate_time, &sample);
if (sample_status != ASYMPTOTIC_OK)
return sample_status;
/* The analytic monopole exterior only covers a fixed, concentric sphere. */
const double center_tol = 1e-12 * fmax(1.0, sample.radius);
for (int i = 0; i < 3; ++i) {
if (fabs(sample.center[i] - end->frame_origin[i]) > center_tol ||
fabs(sample.velocity[i]) > 1e-12 ||
!(sample.radius > 2.0 * end->mass))
return ASYMPTOTIC_UNSUPPORTED;
}
if (fabs(sample.radius_rate) > 1e-12)
return ASYMPTOTIC_UNSUPPORTED;
if (sample.radius / end->mass < 64.0)
return ASYMPTOTIC_UNSUPPORTED;
SchwarzschildCanonical camera;
if (asymptotic_schwarzschild_canonical_from_state(
end, metric, state->coordinate_time, state->x, state->Pi,
state->log_alpha_p0, &camera))
return ASYMPTOTIC_INVALID;
SchwarzschildRouteKind kind = SCH_ROUTE_UNSUPPORTED;
double activate_t = 0.0, x[3] = {0.0, 0.0, 0.0}, Pi[3] = {0.0, 0.0, 0.0};
double log_alpha_p0 = 0.0, n_inf[3] = {0.0, 0.0, 0.0}, frequency = 0.0;
if (asymptotic_schwarzschild_preroute(
end, sample.radius / end->mass, &camera, &kind, &activate_t, x, Pi,
&log_alpha_p0, n_inf, &frequency))
return ASYMPTOTIC_INVALID;
if (kind == SCH_ROUTE_UNSUPPORTED)
return ASYMPTOTIC_UNSUPPORTED;
if (kind == SCH_ROUTE_TIME_RANGE_EXHAUSTED) {
route->kind = ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED;
route->end_id = end->end_id;
return ASYMPTOTIC_TIME_RANGE_EXHAUSTED;
}
route->end_id = end->end_id;
if (kind == SCH_ROUTE_ENTRY) {
route->kind = ASYMPTOTIC_ROUTE_ENTRY;
route->activate_t = activate_t;
for (int i = 0; i < 3; ++i) {
route->x[i] = x[i];
route->Pi[i] = Pi[i];
}
route->log_alpha_p0 = log_alpha_p0;
return ASYMPTOTIC_OK;
}
route->kind = ASYMPTOTIC_ROUTE_ESCAPED;
for (int i = 0; i < 3; ++i)
route->n_infinity[i] = n_inf[i];
route->frequency_ratio = frequency;
return ASYMPTOTIC_OK;
}
AsymptoticStatus asymptotic_route_camera(const SpacetimeSource *source,
const ObserverState *observer,
const double direction[3],
AsymptoticRoute *route) {
if (source == NULL || observer == NULL || direction == NULL || route == NULL)
return ASYMPTOTIC_INVALID;
*route = (AsymptoticRoute){.kind = ASYMPTOTIC_ROUTE_INVALID,
.end_id = SPACETIME_END_NONE};
MetricData metric;
if (spacetime_eval(source, observer->coordinate_time,
observer->coordinate_position, &metric))
return ASYMPTOTIC_INVALID;
GeodesicRayState state;
if (geodesic_initialize_past_ray_metric(&metric, observer, direction, &state))
return ASYMPTOTIC_INVALID;
const size_t count = spacetime_asymptotic_end_count(source);
if (count == 0) {
route->kind = ASYMPTOTIC_ROUTE_INSIDE;
route->end_id = SPACETIME_END_NONE;
route->activate_t = state.coordinate_time;
for (int i = 0; i < 3; ++i) {
route->x[i] = state.x[i];
route->Pi[i] = state.Pi[i];
}
route->log_alpha_p0 = state.log_alpha_p0;
return ASYMPTOTIC_OK;
}
/* A backend that declares ends must describe them consistently and use a
* supported exterior; otherwise the protocol is broken and no route may be
* fabricated (not even an "inside" one). */
for (size_t i = 0; i < count; ++i) {
SpacetimeAsymptoticEnd end;
if (spacetime_asymptotic_end(source, i, &end))
return ASYMPTOTIC_INVALID;
if (end.exterior_kind != ASYMPTOTIC_EXTERIOR_MINKOWSKI &&
end.exterior_kind != ASYMPTOTIC_EXTERIOR_SCHWARZSCHILD_MONOPOLE)
return ASYMPTOTIC_UNSUPPORTED;
}
double past_w[3];
{
double inv[3][3];
if (invert3(metric.gamma, inv))
return ASYMPTOTIC_INVALID;
for (int i = 0; i < 3; ++i) {
double dxdt = -metric.beta[i];
for (int j = 0; j < 3; ++j)
dxdt += metric.alpha * inv[i][j] * state.Pi[j];
past_w[i] = -dxdt;
}
}
for (size_t i = 0; i < count; ++i) {
SpacetimeAsymptoticEnd end;
if (spacetime_asymptotic_end(source, i, &end))
return ASYMPTOTIC_INVALID;
double value, slope;
const AsymptoticStatus ws = worldtube_value_and_slope(
source, &end, state.coordinate_time, state.x, past_w, &value, &slope);
if (ws == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) {
route->kind = ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED;
route->end_id = end.end_id;
return ASYMPTOTIC_TIME_RANGE_EXHAUSTED;
}
if (ws != ASYMPTOTIC_OK)
return ASYMPTOTIC_INVALID;
if (value < 0.0 || (value == 0.0 && slope < 0.0)) {
route->kind = ASYMPTOTIC_ROUTE_INSIDE;
route->end_id = end.end_id;
route->activate_t = state.coordinate_time;
for (int k = 0; k < 3; ++k) {
route->x[k] = state.x[k];
route->Pi[k] = state.Pi[k];
}
route->log_alpha_p0 = state.log_alpha_p0;
return ASYMPTOTIC_OK;
}
}
AsymptoticPhotonState canonical;
int have_entry = 0, have_miss = 0;
double best_s = INFINITY;
AsymptoticRoute best = {.kind = ASYMPTOTIC_ROUTE_INVALID};
SpacetimeEndId first_end = SPACETIME_END_NONE;
double miss_n_inf[3] = {0.0, 0.0, 0.0};
for (size_t i = 0; i < count; ++i) {
SpacetimeAsymptoticEnd end;
if (spacetime_asymptotic_end(source, i, &end))
return ASYMPTOTIC_INVALID;
if (first_end == SPACETIME_END_NONE)
first_end = end.end_id;
if (end.exterior_kind == ASYMPTOTIC_EXTERIOR_SCHWARZSCHILD_MONOPOLE)
return schwarzschild_route(source, &end, &metric, &state, route);
if (end.exterior_kind != ASYMPTOTIC_EXTERIOR_MINKOWSKI)
return ASYMPTOTIC_UNSUPPORTED;
if (asymptotic_canonical_from_backend(source, end.end_id, &metric,
state.coordinate_time, state.x,
state.Pi, state.log_alpha_p0,
&canonical))
return ASYMPTOTIC_INVALID;
AsymptoticRoute candidate = {.kind = ASYMPTOTIC_ROUTE_INVALID};
const AsymptoticStatus status = minkowski_preroute(
source, &end, state.coordinate_time, canonical.x, canonical.w,
end.end_id, &candidate);
if (status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) {
route->kind = ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED;
route->end_id = end.end_id;
return ASYMPTOTIC_TIME_RANGE_EXHAUSTED;
}
if (status != ASYMPTOTIC_OK)
return status;
if (candidate.kind == ASYMPTOTIC_ROUTE_ENTRY) {
const double s = state.coordinate_time - candidate.activate_t;
if (!have_entry || s < best_s) {
have_entry = 1;
best_s = s;
best = candidate;
}
} else if (!have_miss) {
have_miss = 1;
for (int k = 0; k < 3; ++k)
miss_n_inf[k] = candidate.n_infinity[k];
}
}
if (have_entry) {
*route = best;
route->log_alpha_p0 = state.log_alpha_p0;
return ASYMPTOTIC_OK;
}
if (!have_miss) {
/* Declared ends exist but none produced a route: broken protocol. */
return ASYMPTOTIC_INVALID;
}
route->kind = ASYMPTOTIC_ROUTE_ESCAPED;
route->end_id = first_end;
for (int k = 0; k < 3; ++k)
route->n_infinity[k] = miss_n_inf[k];
route->frequency_ratio = exp(-state.log_alpha_p0);
return ASYMPTOTIC_OK;
}
AsymptoticStatus asymptotic_finish_escape(const SpacetimeSource *source,
SpacetimeEndId end_id, double t,
const double x[3], const double Pi[3],
double log_alpha_p0,
RayEndpoint *endpoint) {
if (source == NULL || x == NULL || Pi == NULL || endpoint == NULL)
return ASYMPTOTIC_INVALID;
SpacetimeAsymptoticEnd end;
if (find_end(source, end_id, &end))
return ASYMPTOTIC_INVALID;
if (end.exterior_kind == ASYMPTOTIC_EXTERIOR_SCHWARZSCHILD_MONOPOLE) {
MetricData sch_metric;
if (spacetime_eval(source, t, x, &sch_metric))
return ASYMPTOTIC_INVALID;
SchwarzschildCanonical canonical;
if (asymptotic_schwarzschild_canonical_from_state(
&end, &sch_metric, t, x, Pi, log_alpha_p0, &canonical))
return ASYMPTOTIC_INVALID;
double n_inf[3], frequency;
if (asymptotic_schwarzschild_finish(&end, &canonical, n_inf, &frequency))
return ASYMPTOTIC_INVALID;
for (int i = 0; i < 3; ++i)
endpoint->n_infinity[i] = n_inf[i];
endpoint->frequency_ratio = frequency;
endpoint->end_id = end_id;
endpoint->status = RAY_ENDPOINT_ESCAPED;
endpoint->magnification = 1.0;
return ASYMPTOTIC_OK;
}
if (end.exterior_kind != ASYMPTOTIC_EXTERIOR_MINKOWSKI)
return ASYMPTOTIC_UNSUPPORTED;
MetricData metric;
double inv[3][3];
if (spacetime_eval(source, t, x, &metric) || invert3(metric.gamma, inv))
return ASYMPTOTIC_INVALID;
double n[3] = {0.0, 0.0, 0.0};
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j)
n[i] -= inv[i][j] * Pi[j];
if (normalize3(n) <= 0.0)
return ASYMPTOTIC_INVALID;
double beta_dot_pi = 0.0;
for (int i = 0; i < 3; ++i)
beta_dot_pi += metric.beta[i] * Pi[i];
const double energy = exp(log_alpha_p0) * (metric.alpha - beta_dot_pi);
if (!isfinite(energy) || energy <= 0.0)
return ASYMPTOTIC_INVALID;
for (int i = 0; i < 3; ++i)
endpoint->n_infinity[i] = n[i];
endpoint->frequency_ratio = 1.0 / energy;
endpoint->end_id = end_id;
endpoint->status = RAY_ENDPOINT_ESCAPED;
endpoint->magnification = 1.0;
return ASYMPTOTIC_OK;
}
src/asymptotic.h
+82
View File
@@ -0,0 +1,82 @@
#ifndef ASYMPTOTIC_H
#define ASYMPTOTIC_H
#include "geodesic.h"
#include "spacetime.h"
/* Status codes for the common asymptotic-exterior module. Unsupported and
* exhausted are reported explicitly; callers must not turn them into a
* plausible-looking escape. */
typedef enum {
ASYMPTOTIC_INVALID = -1,
ASYMPTOTIC_OK = 0,
ASYMPTOTIC_UNSUPPORTED = 1,
ASYMPTOTIC_TIME_RANGE_EXHAUSTED = 2
} AsymptoticStatus;
/* Unified canonical photon state in the asymptotic reference frame. `w` is
* the unit past-propagation direction: along the renderer's backward
* integration the spatial position moves as x(s) = x0 + s w, s = t0 - t. */
typedef struct {
SpacetimeEndId end_id;
double t;
double x[3];
double w[3];
double log_alpha_p0;
} AsymptoticPhotonState;
typedef enum {
ASYMPTOTIC_ROUTE_INSIDE, /* camera in a worldtube: activate at the camera */
ASYMPTOTIC_ROUTE_ENTRY, /* camera outside, entry event produced */
ASYMPTOTIC_ROUTE_ESCAPED,
ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED,
ASYMPTOTIC_ROUTE_INVALID
} AsymptoticRouteKind;
typedef struct {
AsymptoticRouteKind kind;
SpacetimeEndId end_id;
/* Activation state, backend coordinates, for INSIDE and ENTRY. */
double activate_t;
double x[3];
double Pi[3];
double log_alpha_p0;
/* Terminal infinity endpoint for ESCAPED. */
double n_infinity[3];
double frequency_ratio;
} AsymptoticRoute;
/* Pre-route one camera ray against every declared end's worldtube. */
AsymptoticStatus asymptotic_route_camera(const SpacetimeSource *source,
const ObserverState *observer,
const double direction[3],
AsymptoticRoute *route);
/* Directed inside->outside crossing helper for the geodesic lifecycle.
* Returns ASYMPTOTIC_OK, ASYMPTOTIC_TIME_RANGE_EXHAUSTED (the backend cannot
* describe the worldtube at this time), or ASYMPTOTIC_INVALID. */
int asymptotic_worldtube_value(const SpacetimeSource *source,
SpacetimeEndId end_id, double t,
const double x[3], double *value);
/* Finish an interior inside->outside crossing to an infinity endpoint.
* Returns ASYMPTOTIC_UNSUPPORTED for exterior models not implemented yet. */
AsymptoticStatus asymptotic_finish_escape(const SpacetimeSource *source,
SpacetimeEndId end_id, double t,
const double x[3], const double Pi[3],
double log_alpha_p0,
RayEndpoint *endpoint);
/* Canonical <-> backend bridge, valid only inside a supported exterior. */
int asymptotic_canonical_from_backend(const SpacetimeSource *source,
SpacetimeEndId end_id,
const MetricData *metric, double t,
const double x[3], const double Pi[3],
double log_alpha_p0,
AsymptoticPhotonState *out);
int asymptotic_backend_from_canonical(const SpacetimeSource *source,
const MetricData *metric,
const AsymptoticPhotonState *canonical,
double x[3], double Pi[3],
double *log_alpha_p0);
#endif
src/asymptotic_gl48.h
+42
View File
@@ -0,0 +1,42 @@
#ifndef ASYMPTOTIC_GL48_H
#define ASYMPTOTIC_GL48_H
/* 48-point Gauss-Legendre nodes and weights on [-1, 1], used only for the
* bounded residual of the Schwarzschild coordinate-time transfer. Generated
* with numpy.polynomial.legendre.leggauss(48); double precision. */
static const double gl48_nodes[48] = {
-0.99877100725242607, -0.99353017226635076, -0.98412458372282685,
-0.97059159254624727, -0.9529877031604308, -0.93138669070655433,
-0.90587913671556963, -0.87657202027424785, -0.84358826162439349,
-0.80706620402944262, -0.76715903251574036, -0.72403413092381463,
-0.67787237963266389, -0.6288673967765136, -0.57722472608397268,
-0.523160974722233, -0.46690290475095841, -0.40868648199071672,
-0.34875588629216075, -0.28736248735545555, -0.22476379039468905,
-0.16122235606889174, -0.097004699209462697, -0.032380170962869367,
0.032380170962869367, 0.097004699209462697, 0.16122235606889174,
0.22476379039468905, 0.28736248735545555, 0.34875588629216075,
0.40868648199071672, 0.46690290475095841, 0.523160974722233,
0.57722472608397268, 0.6288673967765136, 0.67787237963266389,
0.72403413092381463, 0.76715903251574036, 0.80706620402944262,
0.84358826162439349, 0.87657202027424785, 0.90587913671556963,
0.93138669070655433, 0.9529877031604308, 0.97059159254624727,
0.98412458372282685, 0.99353017226635076, 0.99877100725242607};
static const double gl48_weights[48] = {
0.0031533460523098418, 0.0073275539012758505, 0.011477234579234699,
0.015579315722943481, 0.019616160457356105, 0.023570760839324009,
0.027426509708357052, 0.031167227832798117, 0.034777222564770421,
0.038241351065830473, 0.041545082943464533, 0.044674560856694245,
0.04761665849249027, 0.050359035553854216, 0.052890189485193424,
0.05519950369998404, 0.057277292100402881, 0.059114839698395358,
0.060704439165893562, 0.062039423159892415, 0.063114192286253756,
0.06392423858464788, 0.06446616443594981, 0.064737696812683626,
0.064737696812683626, 0.06446616443594981, 0.06392423858464788,
0.063114192286253756, 0.062039423159892415, 0.060704439165893562,
0.059114839698395358, 0.057277292100402881, 0.05519950369998404,
0.052890189485193424, 0.050359035553854216, 0.04761665849249027,
0.044674560856694245, 0.041545082943464533, 0.038241351065830473,
0.034777222564770421, 0.031167227832798117, 0.027426509708357052,
0.023570760839324009, 0.019616160457356105, 0.015579315722943481,
0.011477234579234699, 0.0073275539012758505, 0.0031533460523098418};
#endif
src/asymptotic_schwarzschild.c
+528
View File
@@ -0,0 +1,528 @@
#include "asymptotic_schwarzschild.h"
#include <float.h>
#include <math.h>
#include <stddef.h>
static const double kPi = 3.14159265358979323846;
static double dot3(const double a[3], const double b[3]) {
return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
static void cross3(const double a[3], const double b[3], double out[3]) {
out[0] = a[1] * b[2] - a[2] * b[1];
out[1] = a[2] * b[0] - a[0] * b[2];
out[2] = a[0] * b[1] - a[1] * b[0];
}
static double normalize3(double v[3]) {
const double length = sqrt(dot3(v, v));
if (length > 0.0)
for (int i = 0; i < 3; ++i)
v[i] /= length;
return length;
}
/* ------------------------------------------------------------------------- *
* Complex arithmetic and Carlson R_F.
* ------------------------------------------------------------------------- */
typedef struct {
double re, im;
} cs;
static cs cs_add(cs a, cs b) { return (cs){a.re + b.re, a.im + b.im}; }
static cs cs_sub(cs a, cs b) { return (cs){a.re - b.re, a.im - b.im}; }
static cs cs_mul(cs a, cs b) {
return (cs){a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re};
}
static cs cs_scale(cs a, double s) { return (cs){a.re * s, a.im * s}; }
static double cs_abs(cs a) { return hypot(a.re, a.im); }
static cs cs_inv(cs z) {
const double d = z.re * z.re + z.im * z.im;
return (cs){z.re / d, -z.im / d};
}
static cs cs_sqrt(cs z) {
const double r = hypot(z.re, z.im);
double re = sqrt(0.5 * (r + fabs(z.re)));
double im = sqrt(0.5 * (r - fabs(z.re)));
if (z.re < 0.0) {
const double t = re;
re = im;
im = t;
}
if (z.im < 0.0)
im = -im;
return (cs){re, im};
}
static cs cs_cbrt(cs z) {
/* Cardano needs the real cube root of real arguments (disc > 0); the
* principal complex root is correct for the conjugate pair (disc < 0). */
if (z.im == 0.0)
return (cs){cbrt(z.re), 0.0};
const double r = hypot(z.re, z.im);
const double theta = atan2(z.im, z.re);
const double cr = cbrt(r);
return (cs){cr * cos(theta / 3.0), cr * sin(theta / 3.0)};
}
static cs rf_naive(cs x, cs y, cs z) {
for (int iteration = 0; iteration < 80; ++iteration) {
const cs sx = cs_sqrt(x), sy = cs_sqrt(y), sz = cs_sqrt(z);
const cs lambda =
cs_add(cs_add(cs_mul(sx, sy), cs_mul(sy, sz)), cs_mul(sz, sx));
x = cs_scale(cs_add(x, lambda), 0.25);
y = cs_scale(cs_add(y, lambda), 0.25);
z = cs_scale(cs_add(z, lambda), 0.25);
const cs a = cs_scale(cs_add(cs_add(x, y), z), 1.0 / 3.0);
const cs X = cs_sub((cs){1.0, 0.0}, cs_mul(x, cs_inv(a)));
const cs Y = cs_sub((cs){1.0, 0.0}, cs_mul(y, cs_inv(a)));
const cs Z = cs_sub((cs){1.0, 0.0}, cs_mul(z, cs_inv(a)));
if (fmax(fmax(cs_abs(X), cs_abs(Y)), cs_abs(Z)) < 1e-12) {
const cs e2 = cs_add(cs_add(cs_mul(X, Y), cs_mul(Y, Z)),
cs_mul(Z, X));
const cs e3 = cs_mul(cs_mul(X, Y), Z);
const cs series = cs_add(
cs_add((cs){1.0, 0.0}, cs_scale(cs_mul(e2, e3), -3.0 / 44.0)),
cs_add(cs_scale(cs_mul(e2, e2), 1.0 / 24.0),
cs_add(cs_scale(e2, -1.0 / 10.0), cs_scale(e3, 1.0 / 14.0))));
return cs_mul(series, cs_inv(cs_sqrt(a)));
}
}
return (cs){NAN, NAN};
}
/* R_F via Carlson duplication. The three-real-root branch is handled by the
* real Legendre form, so the only complex calls here come from the conjugate
* root pair, whose arguments are off the real axis and take the principal
* square-root branch consistently. */
static cs rf(cs x, cs y, cs z) { return rf_naive(x, y, z); }
/* Incomplete elliptic integral of the first kind with parameter m = k^2:
* F(phi,m) = sin(phi) R_F(cos^2 phi, 1 - m sin^2 phi, 1). */
static double ellipf(double phi, double m) {
const double s = sin(phi), c = cos(phi);
const cs r = rf((cs){c * c, 0.0}, (cs){1.0 - m * s * s, 0.0},
(cs){1.0, 0.0});
return s * r.re;
}
/* Leading-order estimate of phi for a candidate ordered real-root branch. */
static double phi_three_real(double u0, double A, double B, double C) {
if (!(u0 > A) || !(u0 < B) || !(A < B) || !(B < C))
return NAN;
const double sA = sqrt((0.0 - A) / (B - A));
const double s0 = sqrt((u0 - A) / (B - A));
const double m = (B - A) / (C - A);
return sqrt(2.0) / sqrt(C - A) * (ellipf(asin(s0), m) - ellipf(asin(sA), m));
}
/* Roots of 2 beta^2 u^3 - beta^2 u^2 + 1 = 0 through the depressed cubic
* w^3 + P w + Q = 0 with u = w + 1/6. */
static void cubic_roots(double beta, cs e[3]) {
const double b2 = beta * beta;
const double c = 1.0 / (2.0 * b2);
const double P = -1.0 / 12.0;
const double Q = c - 1.0 / 108.0;
const double halfQ = 0.5 * Q;
const cs disc = (cs){halfQ * halfQ + (P * P * P) / 27.0, 0.0};
const cs sq = cs_sqrt(disc);
const cs u1 = cs_cbrt(cs_add((cs){-halfQ, 0.0}, sq));
const cs u2 = cs_cbrt(cs_add((cs){-halfQ, 0.0}, cs_scale(sq, -1.0)));
const cs omega = (cs){cos(2.0 * kPi / 3.0), sin(2.0 * kPi / 3.0)};
const cs omega2 = cs_mul(omega, omega);
e[0] = cs_add(cs_add(u1, u2), (cs){1.0 / 6.0, 0.0});
e[1] = cs_add(cs_add(cs_mul(omega, u1), cs_mul(omega2, u2)),
(cs){1.0 / 6.0, 0.0});
e[2] = cs_add(cs_add(cs_mul(omega2, u1), cs_mul(omega, u2)),
(cs){1.0 / 6.0, 0.0});
/* Cardano loses relative accuracy in the near-double-root regime. Polish
* the real roots with Newton so the grazing turning root keeps full
* relative precision. */
for (int i = 0; i < 3; ++i) {
if (fabs(e[i].im) > 1e-9 * fmax(1.0, fabs(e[i].re)))
continue;
double u = e[i].re;
for (int step = 0; step < 20; ++step) {
const double p = 2.0 * b2 * u * u * u - b2 * u * u + 1.0;
const double dp = 6.0 * b2 * u * u - 2.0 * b2 * u;
if (dp == 0.0)
break;
const double du = p / dp;
u -= du;
if (fabs(du) <= 1e-18 * fmax(1.0, fabs(u)))
break;
}
e[i] = (cs){u, 0.0};
}
}
double asymptotic_schwarzschild_phi(double rho, double beta) {
if (!isfinite(rho) || rho <= 0.0 || !isfinite(beta) || beta < 0.0)
return NAN;
if (beta == 0.0)
return 0.0;
const double u0 = 1.0 / rho;
cs e[3];
cubic_roots(beta, e);
const double imag_tol = 1e-11 * fmax(1.0, fabs(e[0].re));
if (fabs(e[0].im) < imag_tol && fabs(e[1].im) < imag_tol &&
fabs(e[2].im) < imag_tol) {
/* Three real roots: use the real Legendre form, which is accurate up to
* and through the grazing limit. */
double r[3] = {e[0].re, e[1].re, e[2].re};
for (int i = 0; i < 2; ++i)
for (int j = i + 1; j < 3; ++j)
if (r[j] < r[i]) {
const double t = r[i];
r[i] = r[j];
r[j] = t;
}
const double real_value = phi_three_real(u0, r[0], r[1], r[2]);
if (isfinite(real_value))
return real_value;
}
const cs a = rf(cs_scale(e[0], -1.0), cs_scale(e[1], -1.0),
cs_scale(e[2], -1.0));
const cs b = rf(cs_sub((cs){u0, 0.0}, e[0]),
cs_sub((cs){u0, 0.0}, e[1]),
cs_sub((cs){u0, 0.0}, e[2]));
const cs s = cs_scale(cs_sub(a, b), 2.0);
/* The real integral requires a real S; a non-negligible imaginary part
* means the principal branch failed. Report it instead of silently using
* a wrong angle. */
if (!isfinite(s.re) || fabs(s.im) > 1e-6 * fmax(1.0, fabs(s.re)))
return NAN;
return fabs(s.re) / sqrt(2.0);
}
double asymptotic_schwarzschild_turning_rho(double beta) {
if (!isfinite(beta) || beta <= 3.0 * sqrt(3.0))
return INFINITY;
/* Larger positive root of f(rho) = rho^3 - beta^2 rho + 2 beta^2.
* f(3) = 27 - beta^2 < 0 and f(beta+2) > 0, and f is monotone on the
* bracket beyond its local minimum, so a bracketed bisection is safe. */
const double b2 = beta * beta;
double lo = 3.0, hi = beta + 2.0;
for (int iteration = 0; iteration < 200; ++iteration) {
const double mid = 0.5 * (lo + hi);
const double f = mid * mid * mid - b2 * mid + 2.0 * b2;
if (f < 0.0)
lo = mid;
else
hi = mid;
if (hi - lo <= 4.0 * DBL_EPSILON * hi)
break;
}
double rho = 0.5 * (lo + hi);
for (int step = 0; step < 20; ++step) {
const double f = rho * rho * rho - b2 * rho + 2.0 * b2;
const double fp = 3.0 * rho * rho - b2;
if (fp == 0.0)
break;
const double next = rho - f / fp;
if (!(next > lo && next < hi))
break;
rho = next;
}
return rho;
}
/* ------------------------------------------------------------------------- *
* Coordinate-time transfer (ingoing Kerr-Schild time, M = 1 units).
* ------------------------------------------------------------------------- */
#include "asymptotic_gl48.h"
static double sch_i2(double u_cam, double u_R, double beta) {
/* I2 = int_{u_cam}^{u_R} du / (1 + sqrt(P(u))). Substitute
* u = u_R - (u_R - u_cam) t^2 to remove the grazing branch point. */
const double span = u_R - u_cam;
if (!(span > 0.0))
return 0.0;
double sum = 0.0;
for (int i = 0; i < 48; ++i) {
const double t = 0.5 * (gl48_nodes[i] + 1.0);
const double u = u_R - span * t * t;
const double P = 1.0 - beta * beta * u * u + 2.0 * beta * beta * u * u * u;
const double f = 1.0 / (1.0 + sqrt(P));
sum += gl48_weights[i] * f * 2.0 * span * t;
}
return 0.5 * sum;
}
/* Positive coordinate time to travel outward from R to rho_cam. */
static double sch_time_transfer(double rho_cam, double rho_R, double beta) {
const double u_cam = 1.0 / rho_cam, u_R = 1.0 / rho_R;
const double dphi =
asymptotic_schwarzschild_phi(rho_R, beta) -
asymptotic_schwarzschild_phi(rho_cam, beta);
const double elementary =
(rho_cam - rho_R) + 4.0 * log(u_R / u_cam) -
4.0 * log((1.0 - 2.0 * u_R) / (1.0 - 2.0 * u_cam));
return elementary + (beta * dphi - beta * beta * sch_i2(u_cam, u_R, beta));
}
/* ------------------------------------------------------------------------- *
* Rotation and canonical <-> backend bridging.
* ------------------------------------------------------------------------- */
static void rotate_axis(const double v[3], const double axis[3], double angle,
double out[3]) {
const double c = cos(angle), s = sin(angle);
double cross[3];
cross3(axis, v, cross);
const double adotv = dot3(axis, v);
for (int i = 0; i < 3; ++i)
out[i] = v[i] * c + cross[i] * s + axis[i] * adotv * (1.0 - c);
}
/* The algebraic monopole formulas below assume the asymptotic frame axes are
* the backend Cartesian axes; a rotated frame would require rotating the
* momentum and the worldtube. */
static int sch_frame_is_aligned(const SpacetimeAsymptoticEnd *end) {
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j) {
const double expected = i == j ? 1.0 : 0.0;
if (fabs(end->frame_axes[i][j] - expected) > 1e-12)
return 0;
}
return 1;
}
int asymptotic_schwarzschild_canonical_from_state(
const SpacetimeAsymptoticEnd *end, const MetricData *metric, double t,
const double x[3], const double Pi[3], double log_alpha_p0,
SchwarzschildCanonical *out) {
if (end == NULL || metric == NULL || out == NULL || end->mass <= 0.0 ||
!sch_frame_is_aligned(end))
return -1;
double beta_dot_pi = 0.0;
for (int i = 0; i < 3; ++i)
beta_dot_pi += metric->beta[i] * Pi[i];
const double energy = exp(log_alpha_p0) * (metric->alpha - beta_dot_pi);
if (!isfinite(energy) || energy <= 0.0)
return -1;
double rel[3];
for (int i = 0; i < 3; ++i)
rel[i] = x[i] - end->frame_origin[i];
const double radius = sqrt(dot3(rel, rel));
if (!(radius > 0.0))
return -1;
double Lvec[3];
cross3(rel, Pi, Lvec);
const double Lmag = sqrt(dot3(Lvec, Lvec));
const double denom = metric->alpha - beta_dot_pi;
out->end_id = end->end_id;
out->t = t;
out->rho = radius / end->mass;
out->energy = energy;
for (int i = 0; i < 3; ++i)
out->rhat[i] = rel[i] / radius;
if (Lmag > 0.0) {
out->beta = (Lmag / denom) / end->mass;
for (int i = 0; i < 3; ++i)
out->Lhat[i] = Lvec[i] / Lmag;
} else {
out->beta = 0.0;
out->Lhat[0] = out->Lhat[1] = out->Lhat[2] = 0.0;
}
double inv[3][3];
const double det =
metric->gamma[0][0] * (metric->gamma[1][1] * metric->gamma[2][2] -
metric->gamma[1][2] * metric->gamma[2][1]) -
metric->gamma[0][1] * (metric->gamma[1][0] * metric->gamma[2][2] -
metric->gamma[1][2] * metric->gamma[2][0]) +
metric->gamma[0][2] * (metric->gamma[1][0] * metric->gamma[2][1] -
metric->gamma[1][1] * metric->gamma[2][0]);
inv[0][0] = (metric->gamma[1][1] * metric->gamma[2][2] -
metric->gamma[1][2] * metric->gamma[2][1]) / det;
inv[0][1] = (metric->gamma[0][2] * metric->gamma[2][1] -
metric->gamma[0][1] * metric->gamma[2][2]) / det;
inv[0][2] = (metric->gamma[0][1] * metric->gamma[1][2] -
metric->gamma[0][2] * metric->gamma[1][1]) / det;
inv[1][0] = (metric->gamma[1][2] * metric->gamma[2][0] -
metric->gamma[1][0] * metric->gamma[2][2]) / det;
inv[1][1] = (metric->gamma[0][0] * metric->gamma[2][2] -
metric->gamma[0][2] * metric->gamma[2][0]) / det;
inv[1][2] = (metric->gamma[0][2] * metric->gamma[1][0] -
metric->gamma[0][0] * metric->gamma[1][2]) / det;
inv[2][0] = (metric->gamma[1][0] * metric->gamma[2][1] -
metric->gamma[1][1] * metric->gamma[2][0]) / det;
inv[2][1] = (metric->gamma[0][1] * metric->gamma[2][0] -
metric->gamma[0][0] * metric->gamma[2][1]) / det;
inv[2][2] = (metric->gamma[0][0] * metric->gamma[1][1] -
metric->gamma[0][1] * metric->gamma[1][0]) / det;
double dxdt[3];
for (int i = 0; i < 3; ++i) {
dxdt[i] = -metric->beta[i];
for (int j = 0; j < 3; ++j)
dxdt[i] += metric->alpha * inv[i][j] * Pi[j];
}
double radial = 0.0;
for (int i = 0; i < 3; ++i)
radial += -dxdt[i] * out->rhat[i];
out->radial_sign = radial > 0.0 ? 1 : (radial < 0.0 ? -1 : 0);
return 0;
}
int asymptotic_schwarzschild_state_from_canonical(
const SpacetimeAsymptoticEnd *end, const SchwarzschildCanonical *c,
double x[3], double Pi[3], double *log_alpha_p0) {
if (end == NULL || c == NULL || x == NULL || Pi == NULL ||
!sch_frame_is_aligned(end))
return -1;
const double rho = c->rho;
if (!(rho > 2.0))
return -1;
const double Q = 1.0 - c->beta * c->beta * (1.0 - 2.0 / rho) / (rho * rho);
if (!(Q >= 0.0))
return -1;
const double sqrtQ = sqrt(Q);
double e_phi[3] = {0.0, 0.0, 0.0};
if (c->beta > 0.0)
cross3(c->Lhat, c->rhat, e_phi);
const double s_aff = -(double)c->radial_sign; /* physical (future) radial */
const double kr = s_aff * c->energy * sqrtQ;
const double ktang = c->beta * c->energy / rho;
double kvec[3];
for (int i = 0; i < 3; ++i)
kvec[i] = kr * c->rhat[i] + ktang * e_phi[i];
const double kt_s = c->energy / (1.0 - 2.0 / rho);
const double kt_ks = kt_s + (2.0 / (rho - 2.0)) * kr;
const double alpha = 1.0 / sqrt(1.0 + 2.0 / rho);
const double ak0 = alpha * kt_ks;
if (!isfinite(ak0) || ak0 <= 0.0)
return -1;
for (int i = 0; i < 3; ++i) {
x[i] = end->frame_origin[i] + end->mass * rho * c->rhat[i];
const double kcov = kvec[i] + (2.0 / rho) * c->rhat[i] * (kr + kt_ks);
Pi[i] = kcov / ak0;
}
if (log_alpha_p0 != NULL)
*log_alpha_p0 = log(ak0);
return 0;
}
int asymptotic_schwarzschild_finish(const SpacetimeAsymptoticEnd *end,
const SchwarzschildCanonical *canonical,
double n_infinity[3],
double *frequency_ratio) {
if (end == NULL || canonical == NULL || n_infinity == NULL)
return -1;
const double phi = asymptotic_schwarzschild_phi(canonical->rho,
canonical->beta);
if (!isfinite(phi))
return -1;
if (canonical->beta > 0.0) {
double e_phi[3];
cross3(canonical->Lhat, canonical->rhat, e_phi);
for (int i = 0; i < 3; ++i)
n_infinity[i] = cos(phi) * canonical->rhat[i] -
sin(phi) * e_phi[i];
} else {
const double s = canonical->radial_sign >= 0 ? 1.0 : -1.0;
for (int i = 0; i < 3; ++i)
n_infinity[i] = s * canonical->rhat[i];
}
normalize3(n_infinity);
if (frequency_ratio != NULL)
*frequency_ratio = 1.0 / canonical->energy;
return 0;
}
int asymptotic_schwarzschild_preroute(
const SpacetimeAsymptoticEnd *end, double worldtube_radius,
const SchwarzschildCanonical *camera, SchwarzschildRouteKind *kind,
double *activate_t, double x[3], double Pi[3], double *log_alpha_p0,
double n_infinity[3], double *frequency_ratio) {
if (end == NULL || camera == NULL || kind == NULL)
return -1;
const double R = worldtube_radius;
if (!(R > 2.0) || !(camera->rho >= R))
return -1;
if (R / 1.0 < 64.0) {
*kind = SCH_ROUTE_UNSUPPORTED;
return 0;
}
const double beta_R = R / sqrt(1.0 - 2.0 / R);
if (camera->radial_sign >= 0) {
/* Past propagation is outward or tangent: no entry, immediate infinity
* endpoint. (radial_sign == 0 means the camera is on the boundary with a
* tangent ray.) */
const double phi = asymptotic_schwarzschild_phi(camera->rho, camera->beta);
if (!isfinite(phi)) {
*kind = SCH_ROUTE_UNSUPPORTED;
return 0;
}
if (camera->beta > 0.0) {
double e_phi[3];
cross3(camera->Lhat, camera->rhat, e_phi);
for (int i = 0; i < 3; ++i)
n_infinity[i] = cos(phi) * camera->rhat[i] - sin(phi) * e_phi[i];
} else {
for (int i = 0; i < 3; ++i)
n_infinity[i] = camera->rhat[i];
}
normalize3(n_infinity);
*frequency_ratio = 1.0 / camera->energy;
*kind = SCH_ROUTE_ESCAPED;
return 0;
}
if (camera->beta < beta_R) {
const double dphi = asymptotic_schwarzschild_phi(R, camera->beta) -
asymptotic_schwarzschild_phi(camera->rho, camera->beta);
double rhat_entry[3];
if (camera->beta > 0.0)
rotate_axis(camera->rhat, camera->Lhat, -dphi, rhat_entry);
else
for (int i = 0; i < 3; ++i)
rhat_entry[i] = camera->rhat[i];
SchwarzschildCanonical entry = *camera;
entry.rho = R;
for (int i = 0; i < 3; ++i)
entry.rhat[i] = rhat_entry[i];
entry.radial_sign = -1;
if (asymptotic_schwarzschild_state_from_canonical(end, &entry, x, Pi,
log_alpha_p0))
return -1;
const double T =
sch_time_transfer(camera->rho, R, camera->beta);
*activate_t = camera->t - end->mass * T;
*kind = SCH_ROUTE_ENTRY;
return 0;
}
/* Inward but misses: turn before R and escape. */
const double rho_turn = asymptotic_schwarzschild_turning_rho(camera->beta);
if (!isfinite(rho_turn) || rho_turn > camera->rho) {
*kind = SCH_ROUTE_UNSUPPORTED;
return 0;
}
const double total =
2.0 * asymptotic_schwarzschild_phi(rho_turn, camera->beta) -
asymptotic_schwarzschild_phi(camera->rho, camera->beta);
if (!isfinite(total)) {
*kind = SCH_ROUTE_UNSUPPORTED;
return 0;
}
if (camera->beta > 0.0) {
double e_phi[3];
cross3(camera->Lhat, camera->rhat, e_phi);
for (int i = 0; i < 3; ++i)
n_infinity[i] = cos(total) * camera->rhat[i] - sin(total) * e_phi[i];
} else {
for (int i = 0; i < 3; ++i)
n_infinity[i] = camera->rhat[i];
}
normalize3(n_infinity);
*frequency_ratio = 1.0 / camera->energy;
*kind = SCH_ROUTE_ESCAPED;
return 0;
}
src/asymptotic_schwarzschild.h
+66
View File
@@ -0,0 +1,66 @@
#ifndef ASYMPTOTIC_SCHWARZSCHILD_H
#define ASYMPTOTIC_SCHWARZSCHILD_H
#include "geodesic.h"
#include "spacetime.h"
/* Canonical photon state for a fixed, concentric Schwarzschild monopole
* exterior. All radial quantities are in units of the mass: rho = r / M.
* `Lhat` is the (unit) conserved angular-momentum direction = normalize(x x
* Pi); `beta` is the impact parameter b/M > 0. `radial_sign` is the sign of
* dr/ds along the renderer's past propagation (s = t_camera - t): +1 outward
* into the past, -1 inward into the past. `energy` is E = -p_t with the
* camera normalization E_camera = 1. */
typedef struct {
SpacetimeEndId end_id;
double t;
double rho;
double rhat[3];
double Lhat[3];
double beta;
double energy;
int radial_sign;
} SchwarzschildCanonical;
/* Angular primitive Phi(rho, beta): the azimuth swept on the outward branch
* from radius rho to infinity. Returns NAN outside the supported domain. */
double asymptotic_schwarzschild_phi(double rho, double beta);
/* Larger positive turning radius for the given impact parameter, or INFINITY
* when no turning point exists (beta <= 3 sqrt(3)). */
double asymptotic_schwarzschild_turning_rho(double beta);
/* Convert a backend state into the canonical form. `metric` must be the
* Schwarzschild Kerr-Schild metric at (t, x). */
int asymptotic_schwarzschild_canonical_from_state(
const SpacetimeAsymptoticEnd *end, const MetricData *metric, double t,
const double x[3], const double Pi[3], double log_alpha_p0,
SchwarzschildCanonical *out);
/* Rebuild the backend state at the stored radius / radial directions. */
int asymptotic_schwarzschild_state_from_canonical(
const SpacetimeAsymptoticEnd *end, const SchwarzschildCanonical *canonical,
double x[3], double Pi[3], double *log_alpha_p0);
/* Infinity endpoint for an outward crossing at the canonical radius. */
int asymptotic_schwarzschild_finish(const SpacetimeAsymptoticEnd *end,
const SchwarzschildCanonical *canonical,
double n_infinity[3],
double *frequency_ratio);
/* Pre-route a camera ray outside the worldtube. Fills one of the route
* kinds. `worldtube_radius` is R/M. */
typedef enum {
SCH_ROUTE_ENTRY,
SCH_ROUTE_ESCAPED,
SCH_ROUTE_TIME_RANGE_EXHAUSTED,
SCH_ROUTE_UNSUPPORTED
} SchwarzschildRouteKind;
int asymptotic_schwarzschild_preroute(
const SpacetimeAsymptoticEnd *end, double worldtube_radius,
const SchwarzschildCanonical *camera, SchwarzschildRouteKind *kind,
double *activate_t, double x[3], double Pi[3], double *log_alpha_p0,
double n_infinity[3], double *frequency_ratio);
#endif
src/frame.c
+14 -1
View File
@@ -116,6 +116,7 @@ int frame_lens_mesh_trace(FrameLensMesh *mesh, const SpacetimeSource *spacetime,
RayEndpoint endpoint = geodesic_trace_past(spacetime, observer,
vertex->camera_direction, trace);
vertex->status = endpoint.status;
vertex->end_id = endpoint.end_id;
vertex->traced = 1;
if (endpoint.status == RAY_ENDPOINT_ESCAPED) {
for (int axis = 0; axis < 3; ++axis)
@@ -437,11 +438,20 @@ static int all_vertices_traced(const FrameLensMesh *mesh) {
static int terminal_mismatch(const LensVertex *a, const LensVertex *b,
const LensVertex *c) {
int escaped = 0, captured = 0;
int escaped = 0, captured = 0, have_end = 0;
SpacetimeEndId end = SPACETIME_END_NONE;
const LensVertex *vertices[] = {a, b, c};
for (size_t i = 0; i < 3; ++i) {
escaped |= vertices[i]->status == RAY_ENDPOINT_ESCAPED;
captured |= vertices[i]->status == RAY_ENDPOINT_CAPTURED;
if (vertices[i]->status == RAY_ENDPOINT_ESCAPED) {
if (!have_end) {
end = vertices[i]->end_id;
have_end = 1;
} else if (vertices[i]->end_id != end) {
return 1; /* two different infinity ends must not be interpolated */
}
}
}
return escaped && captured;
}
@@ -576,6 +586,7 @@ int frame_lens_mesh_install_sample(FrameLensMesh *mesh, size_t sample_id,
? &mesh->vertices[sample->vertex_id]
: &sample->vertex;
vertex->status = endpoint->status;
vertex->end_id = endpoint->end_id;
vertex->traced = 1;
if (endpoint->status == RAY_ENDPOINT_ESCAPED) {
for (int axis = 0; axis < 3; ++axis)
@@ -610,6 +621,8 @@ static int discrete_jacobian(const FrameLensMesh *mesh,
if (a->status != RAY_ENDPOINT_ESCAPED || b->status != RAY_ENDPOINT_ESCAPED ||
c->status != RAY_ENDPOINT_ESCAPED)
return 0;
if (a->end_id != b->end_id || a->end_id != c->end_id)
return 0;
const double image_area = spherical_signed_area(
a->camera_direction, b->camera_direction, c->camera_direction);
if (!isfinite(image_area) || fabs(image_area) <= 1e-15)
src/frame.h
+3
View File
@@ -15,6 +15,9 @@ typedef struct {
double n_infinity[3];
double log_frequency_ratio;
RayEndpointStatus status;
/* Asymptotic end this escaped vertex belongs to; a triangle must not
* interpolate across two different ends. */
SpacetimeEndId end_id;
int traced;
} LensVertex;
src/geodesic.c
+220 -31
View File
@@ -1,4 +1,5 @@
#include "geodesic.h"
#include "asymptotic.h"
#include <math.h>
#include <stddef.h>
@@ -105,15 +106,14 @@ static int rk4(const MetricSlab *slab, double t, double h, State *s) {
return 0;
}
int geodesic_initialize_past_ray(const MetricSlab *slab,
const ObserverState *o, const double n[3],
State *s) {
MetricData m;
int geodesic_initialize_past_ray_metric(const MetricData *metric,
const ObserverState *o,
const double n[3], State *s) {
const MetricData *m = metric;
if (m->alpha <= 0)
return -1;
double k[4] = {o->tetrad[0][0], o->tetrad[0][1], o->tetrad[0][2],
o->tetrad[0][3]};
if (spacetime_slab_eval(slab, o->coordinate_time, o->coordinate_position, &m) ||
m.alpha <= 0)
return -1;
for (int a = 0; a < 3; a++)
for (int mu = 0; mu < 4; mu++)
k[mu] -= n[a] * o->tetrad[a + 1][mu];
@@ -123,15 +123,24 @@ int geodesic_initialize_past_ray(const MetricSlab *slab,
s->x[i] = o->coordinate_position[i];
s->Pi[i] = 0;
for (int j = 0; j < 3; j++)
s->Pi[i] += m.gamma[i][j] * (k[j + 1] + m.beta[j] * k[0]);
s->Pi[i] /= m.alpha * k[0];
s->Pi[i] += m->gamma[i][j] * (k[j + 1] + m->beta[j] * k[0]);
s->Pi[i] /= m->alpha * k[0];
}
s->log_alpha_p0 = log(m.alpha * k[0]);
s->log_alpha_p0 = log(m->alpha * k[0]);
s->coordinate_time = o->coordinate_time;
s->steps = 0;
return isfinite(s->log_alpha_p0) ? 0 : -1;
}
int geodesic_initialize_past_ray(const MetricSlab *slab,
const ObserverState *o, const double n[3],
State *s) {
MetricData m;
if (spacetime_slab_eval(slab, o->coordinate_time, o->coordinate_position, &m))
return -1;
return geodesic_initialize_past_ray_metric(&m, o, n, s);
}
static int escaped_direction(const MetricSlab *slab, double t,
const State *s, double n[3]) {
MetricData m;
@@ -151,6 +160,82 @@ static int escaped_direction(const MetricSlab *slab, double t,
return 0;
}
/* Three-state lifecycle selection. Only a backend that declares no ends at
* all may use the legacy region test; a declared but inconsistent or
* unsupported end is an explicit protocol error, never a silent fallback. */
typedef enum {
ASYM_LIFECYCLE_NONE,
ASYM_LIFECYCLE_READY,
ASYM_LIFECYCLE_PROTOCOL_ERROR
} AsymLifecycleMode;
static AsymLifecycleMode asym_lifecycle_mode(const SpacetimeSource *source) {
const size_t count = spacetime_asymptotic_end_count(source);
if (count == 0)
return ASYM_LIFECYCLE_NONE;
for (size_t i = 0; i < count; ++i) {
SpacetimeAsymptoticEnd end;
if (spacetime_asymptotic_end(source, i, &end))
return ASYM_LIFECYCLE_PROTOCOL_ERROR;
if (end.exterior_kind != ASYMPTOTIC_EXTERIOR_MINKOWSKI &&
end.exterior_kind != ASYMPTOTIC_EXTERIOR_SCHWARZSCHILD_MONOPOLE)
return ASYM_LIFECYCLE_PROTOCOL_ERROR;
}
return ASYM_LIFECYCLE_READY;
}
/* Bounded root localization of the first inside->outside worldtube crossing
* within one accepted step. Re-integrates from `before` with fractional step
* sizes; `after` lands on the outside end of the bracket. Returns
* ASYMPTOTIC_OK, ASYMPTOTIC_TIME_RANGE_EXHAUSTED (a midpoint fell into a
* history hole), or ASYMPTOTIC_INVALID. */
static AsymptoticStatus localize_worldtube_crossing(const MetricSlab *slab,
SpacetimeEndId end_id,
const State *before,
double h, State *after) {
double f_lo = 0.0, f_hi = 1.0;
for (int iteration = 0; iteration < 64; ++iteration) {
const double f = 0.5 * (f_lo + f_hi);
State mid = *before;
if (rk4(slab, before->coordinate_time, h * f, &mid))
return ASYMPTOTIC_INVALID;
mid.coordinate_time = before->coordinate_time + h * f;
double value;
const int status = asymptotic_worldtube_value(
slab->source, end_id, mid.coordinate_time, mid.x, &value);
if (status != ASYMPTOTIC_OK)
return (AsymptoticStatus)status;
if (value >= 0.0)
f_hi = f;
else
f_lo = f;
}
*after = *before;
if (rk4(slab, before->coordinate_time, h * f_hi, after))
return ASYMPTOTIC_INVALID;
after->coordinate_time = before->coordinate_time + h * f_hi;
after->steps = before->steps + 1;
return ASYMPTOTIC_OK;
}
static GeodesicAdvanceResult legacy_escape_or_capture(const MetricSlab *slab,
const State *s,
SpacetimeRayStatus status,
RayEndpoint *out) {
out->status =
status == SPACETIME_RAY_ESCAPED ? RAY_ENDPOINT_ESCAPED
: RAY_ENDPOINT_CAPTURED;
if (out->status == RAY_ENDPOINT_ESCAPED) {
if (escaped_direction(slab, s->coordinate_time, s, out->n_infinity) == 0)
out->frequency_ratio = exp(-s->log_alpha_p0);
else
out->status = RAY_ENDPOINT_INTEGRATION_FAILURE;
}
return out->status == RAY_ENDPOINT_INTEGRATION_FAILURE
? GEODESIC_ADVANCE_FAILED
: GEODESIC_ADVANCE_TERMINATED;
}
GeodesicAdvanceResult geodesic_advance_past_ray(
const MetricSlab *slab, State *s, double slab_left_time,
const GeodesicTraceConfig *config, RayEndpoint *out) {
@@ -158,36 +243,105 @@ GeodesicAdvanceResult geodesic_advance_past_ray(
!config->max_steps || !isfinite(slab_left_time) ||
slab_left_time > s->coordinate_time)
return GEODESIC_ADVANCE_FAILED;
const AsymLifecycleMode mode = asym_lifecycle_mode(slab->source);
if (mode == ASYM_LIFECYCLE_PROTOCOL_ERROR) {
out->status = RAY_ENDPOINT_INVALID;
return GEODESIC_ADVANCE_FAILED;
}
const int directed = mode == ASYM_LIFECYCLE_READY;
const size_t end_count =
directed ? spacetime_asymptotic_end_count(slab->source) : 0;
while (s->coordinate_time > slab_left_time) {
if (config->capture_log_alpha_p0 > 0.0 &&
s->log_alpha_p0 >= config->capture_log_alpha_p0) {
out->status = RAY_ENDPOINT_CAPTURED;
return GEODESIC_ADVANCE_TERMINATED;
}
SpacetimeRayStatus status =
const SpacetimeRayStatus status =
spacetime_slab_classify(slab, s->coordinate_time, s->x);
if (status != SPACETIME_RAY_ACTIVE) {
out->status = status == SPACETIME_RAY_ESCAPED ? RAY_ENDPOINT_ESCAPED
: RAY_ENDPOINT_CAPTURED;
if (out->status == RAY_ENDPOINT_ESCAPED &&
escaped_direction(slab, s->coordinate_time, s, out->n_infinity) == 0)
out->frequency_ratio = exp(-s->log_alpha_p0);
else if (out->status == RAY_ENDPOINT_ESCAPED)
out->status = RAY_ENDPOINT_INTEGRATION_FAILURE;
return out->status == RAY_ENDPOINT_INTEGRATION_FAILURE
? GEODESIC_ADVANCE_FAILED
: GEODESIC_ADVANCE_TERMINATED;
if (status == SPACETIME_RAY_CAPTURED) {
out->status = RAY_ENDPOINT_CAPTURED;
return GEODESIC_ADVANCE_TERMINATED;
}
if (!directed && status != SPACETIME_RAY_ACTIVE)
return legacy_escape_or_capture(slab, s, status, out);
if (s->steps >= config->max_steps) {
out->status = RAY_ENDPOINT_MAX_STEPS;
return GEODESIC_ADVANCE_TERMINATED;
}
const double h = -fmin(config->coordinate_time_step,
s->coordinate_time - slab_left_time);
const State before = *s;
if (rk4(slab, s->coordinate_time, h, s))
return GEODESIC_ADVANCE_FAILED;
s->coordinate_time += h;
++s->steps;
if (!directed)
continue;
if (spacetime_slab_classify(slab, s->coordinate_time, s->x) ==
SPACETIME_RAY_CAPTURED) {
out->status = RAY_ENDPOINT_CAPTURED;
return GEODESIC_ADVANCE_TERMINATED;
}
for (size_t i = 0; i < end_count; ++i) {
SpacetimeAsymptoticEnd end;
if (spacetime_asymptotic_end(slab->source, i, &end)) {
out->status = RAY_ENDPOINT_INVALID;
return GEODESIC_ADVANCE_FAILED;
}
double f_before, f_after;
const int before_status = asymptotic_worldtube_value(
slab->source, end.end_id, before.coordinate_time, before.x,
&f_before);
const int after_status = asymptotic_worldtube_value(
slab->source, end.end_id, s->coordinate_time, s->x, &f_after);
if (before_status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED ||
after_status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) {
/* A first-class terminal reason, matching pre-route exhaustion:
* preserve the end id and install it as terminated provenance. */
out->status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED;
out->end_id = end.end_id;
return GEODESIC_ADVANCE_TERMINATED;
}
if (before_status != ASYMPTOTIC_OK || after_status != ASYMPTOTIC_OK) {
/* The backend cannot describe its own worldtube; this is an explicit
* failure, not a physical escape. */
out->status = RAY_ENDPOINT_INVALID;
return GEODESIC_ADVANCE_FAILED;
}
/* Strict inside->outside: the step must end strictly outside, so a
* single touch at F == 0 (a tangent) is not accepted as a crossing.
* A crossing whose root lands exactly on a step boundary is picked up
* on the following step as f_before == 0, f_after > 0. */
if (f_before > 0.0 || f_after <= 0.0)
continue;
State crossing;
const AsymptoticStatus localized = localize_worldtube_crossing(
slab, end.end_id, &before, h, &crossing);
if (localized == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) {
out->status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED;
out->end_id = end.end_id;
return GEODESIC_ADVANCE_TERMINATED;
}
if (localized != ASYMPTOTIC_OK) {
out->status = RAY_ENDPOINT_INVALID;
return GEODESIC_ADVANCE_FAILED;
}
const AsymptoticStatus transfer = asymptotic_finish_escape(
slab->source, end.end_id, crossing.coordinate_time, crossing.x,
crossing.Pi, crossing.log_alpha_p0, out);
if (transfer == ASYMPTOTIC_OK)
return GEODESIC_ADVANCE_TERMINATED;
out->status = transfer == ASYMPTOTIC_TIME_RANGE_EXHAUSTED
? RAY_ENDPOINT_TIME_RANGE_EXHAUSTED
: RAY_ENDPOINT_INVALID;
if (transfer == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) {
out->end_id = end.end_id;
return GEODESIC_ADVANCE_TERMINATED;
}
return GEODESIC_ADVANCE_FAILED;
}
}
return GEODESIC_ADVANCE_ACTIVE;
}
@@ -198,22 +352,57 @@ RayEndpoint geodesic_trace_past(const SpacetimeSource *source,
const GeodesicTraceConfig *config) {
RayEndpoint out = {.frequency_ratio = 0,
.magnification = 1,
.end_id = SPACETIME_END_NONE,
.status = RAY_ENDPOINT_INTEGRATION_FAILURE};
State state;
MetricSlab *slab = NULL;
if (!source || !observer || !config || config->coordinate_time_step <= 0 ||
!config->max_steps || fabs(dot(n, n) - 1) > 1e-10)
return out;
if (spacetime_load_slab(source, observer->coordinate_time,
observer->coordinate_time -
config->coordinate_time_step * config->max_steps - 1.0,
&slab) ||
geodesic_initialize_past_ray(slab, observer, n, &state)) {
spacetime_free_slab(slab);
AsymptoticRoute route;
const AsymptoticStatus route_status =
asymptotic_route_camera(source, observer, n, &route);
if (route_status == ASYMPTOTIC_UNSUPPORTED) {
out.status = RAY_ENDPOINT_INVALID;
return out;
}
if (route_status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) {
out.status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED;
out.end_id = route.end_id;
return out;
}
if (route_status != ASYMPTOTIC_OK) {
out.status = RAY_ENDPOINT_INVALID;
return out;
}
if (route.kind == ASYMPTOTIC_ROUTE_ESCAPED) {
for (int i = 0; i < 3; ++i)
out.n_infinity[i] = route.n_infinity[i];
out.frequency_ratio = route.frequency_ratio;
out.end_id = route.end_id;
out.status = RAY_ENDPOINT_ESCAPED;
return out;
}
if (route.kind == ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED) {
out.status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED;
out.end_id = route.end_id;
return out;
}
if (route.kind != ASYMPTOTIC_ROUTE_INSIDE &&
route.kind != ASYMPTOTIC_ROUTE_ENTRY) {
out.status = RAY_ENDPOINT_INVALID;
return out;
}
State state = {.coordinate_time = route.activate_t,
.x = {route.x[0], route.x[1], route.x[2]},
.Pi = {route.Pi[0], route.Pi[1], route.Pi[2]},
.log_alpha_p0 = route.log_alpha_p0,
.steps = 0};
const double last_time =
observer->coordinate_time - config->coordinate_time_step * config->max_steps;
route.activate_t - config->coordinate_time_step * config->max_steps;
MetricSlab *slab = NULL;
if (spacetime_load_slab(source, route.activate_t, last_time - 1.0, &slab)) {
out.status = RAY_ENDPOINT_INVALID;
return out;
}
if (geodesic_advance_past_ray(slab, &state, last_time, config, &out) ==
GEODESIC_ADVANCE_ACTIVE)
out.status = RAY_ENDPOINT_MAX_STEPS;
src/geodesic.h
+12 -1
View File
@@ -8,13 +8,18 @@ typedef enum {
RAY_ENDPOINT_ESCAPED,
RAY_ENDPOINT_CAPTURED,
RAY_ENDPOINT_MAX_STEPS,
RAY_ENDPOINT_INTEGRATION_FAILURE
RAY_ENDPOINT_INTEGRATION_FAILURE,
RAY_ENDPOINT_TIME_RANGE_EXHAUSTED,
RAY_ENDPOINT_INVALID
} RayEndpointStatus;
typedef struct {
double n_infinity[3];
double frequency_ratio; /* E_camera / E_infinity */
double magnification; /* Filled by the future local inverse lens map. */
/* Meaningful for RAY_ENDPOINT_ESCAPED and for
* RAY_ENDPOINT_TIME_RANGE_EXHAUSTED; SPACETIME_END_NONE otherwise. */
SpacetimeEndId end_id;
RayEndpointStatus status;
} RayEndpoint;
@@ -52,6 +57,12 @@ int geodesic_initialize_past_ray(const MetricSlab *slab,
const ObserverState *observer,
const double camera_direction[3],
GeodesicRayState *state);
/* Metric-based core of the initialization above; used by the asymptotic
* pre-route, which evaluates the metric at the camera event directly. */
int geodesic_initialize_past_ray_metric(const MetricData *metric,
const ObserverState *observer,
const double camera_direction[3],
GeodesicRayState *state);
GeodesicAdvanceResult geodesic_advance_past_ray(
const MetricSlab *slab, GeodesicRayState *state,
double slab_left_time, const GeodesicTraceConfig *config,
src/main.c
+9 -6
View File
@@ -870,8 +870,12 @@ static double alcubierre_step_budget(const Settings *s) {
static GeodesicTraceConfig trace_config(const Settings *s) {
#ifdef SPACETIME_SCHWARZSCHILD
(void)s;
/* The directed worldtube crossing makes an escaping ray traverse the
* interior as a round trip from the entry sphere (in, turn, back out),
* rather than the old one-way stop at the first outside sample. The step
* budget must cover roughly twice the escape sphere plus margin. */
return (GeodesicTraceConfig){.coordinate_time_step = 0.1,
.max_steps = 4096,
.max_steps = 65536,
.capture_log_alpha_p0 = 8.0};
#elif defined(SPACETIME_ALCUBIERRE)
const double step = alcubierre_time_step(s);
@@ -949,11 +953,9 @@ static int build_observer(const Settings *s, const SpacetimeSource *spacetime,
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;
}
/* A camera outside the escape sphere is supported by the asymptotic
* exterior module for every declared end kind; unsupported exteriors are
* reported through the ray endpoints instead. */
MetricData metric;
if (spacetime_eval(spacetime, camera.coordinate_time, camera.position, &metric)) {
fputs("Could not evaluate metric at the camera event.\n", stderr);
@@ -1133,6 +1135,7 @@ static int trace_movie_generation(Movie *movie, const Settings *s,
return -1;
}
}
ray_pool_preroute(&rays, spacetime);
double slab_hi = movie->frames[movie->frame_count - 1].coordinate_time;
size_t slab_id = 0;
while (ray_pool_has_live(&rays)) {
src/ray.c
+68 -15
View File
@@ -1,5 +1,8 @@
#include "ray.h"
#include "asymptotic.h"
#include <math.h>
#include <omp.h>
#include <stdlib.h>
@@ -11,8 +14,9 @@ int ray_pool_init(RayPool *p, size_t capacity) {
if (!(RAY_ALLOC(t) && RAY_ALLOC(x0) && RAY_ALLOC(x1) && RAY_ALLOC(x2) &&
RAY_ALLOC(p0) && RAY_ALLOC(p1) && RAY_ALLOC(p2) && RAY_ALLOC(observer) &&
RAY_ALLOC(direction0) && RAY_ALLOC(direction1) && RAY_ALLOC(direction2) &&
RAY_ALLOC(log_alpha_p0) && RAY_ALLOC(steps) && RAY_ALLOC(frame_id) &&
RAY_ALLOC(vertex_id) && RAY_ALLOC(status) && RAY_ALLOC(endpoint))) {
RAY_ALLOC(log_alpha_p0) && RAY_ALLOC(activate_t) && RAY_ALLOC(steps) &&
RAY_ALLOC(frame_id) && RAY_ALLOC(vertex_id) && RAY_ALLOC(status) &&
RAY_ALLOC(endpoint))) {
ray_pool_destroy(p);
return -1;
}
@@ -29,6 +33,7 @@ int ray_pool_append(RayPool *p, const ObserverState *observer,
if (observer == NULL || direction == NULL)
return -1;
p->t[i] = observer->coordinate_time;
p->activate_t[i] = observer->coordinate_time;
p->observer[i] = observer;
p->direction0[i] = direction[0];
p->direction1[i] = direction[1];
@@ -37,28 +42,75 @@ int ray_pool_append(RayPool *p, const ObserverState *observer,
p->vertex_id[i] = vertex_id;
p->status[i] = RAY_POOL_PENDING;
p->endpoint[i] = (RayEndpoint){.magnification = 1.0,
.end_id = SPACETIME_END_NONE,
.status = RAY_ENDPOINT_INTEGRATION_FAILURE};
++p->count;
return 0;
}
void ray_pool_activate_in_time_range(RayPool *p, const MetricSlab *slab) {
void ray_pool_preroute(RayPool *p, const SpacetimeSource *source) {
if (p == NULL || source == NULL)
return;
#pragma omp parallel for schedule(static)
for (size_t i = 0; i < p->count; ++i) {
if (p->status[i] != RAY_POOL_PENDING || p->t[i] > slab->t_hi ||
p->t[i] <= slab->t_lo)
if (p->status[i] != RAY_POOL_PENDING)
continue;
GeodesicRayState state;
if (geodesic_initialize_past_ray(
slab, p->observer[i],
(double[]){p->direction0[i], p->direction1[i], p->direction2[i]},
&state)) {
AsymptoticRoute route;
const AsymptoticStatus status = asymptotic_route_camera(
source, p->observer[i],
(double[]){p->direction0[i], p->direction1[i], p->direction2[i]},
&route);
if (status == ASYMPTOTIC_UNSUPPORTED || status == ASYMPTOTIC_INVALID) {
p->endpoint[i].status = RAY_ENDPOINT_INVALID;
p->endpoint[i].end_id = route.end_id;
p->status[i] = RAY_POOL_FAILED;
continue;
}
p->x0[i] = state.x[0]; p->x1[i] = state.x[1]; p->x2[i] = state.x[2];
p->p0[i] = state.Pi[0]; p->p1[i] = state.Pi[1]; p->p2[i] = state.Pi[2];
p->log_alpha_p0[i] = state.log_alpha_p0;
p->steps[i] = state.steps;
if (status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) {
p->endpoint[i].status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED;
p->endpoint[i].end_id = route.end_id;
p->status[i] = RAY_POOL_TERMINATED;
continue;
}
if (route.kind == ASYMPTOTIC_ROUTE_ESCAPED) {
for (int axis = 0; axis < 3; ++axis)
p->endpoint[i].n_infinity[axis] = route.n_infinity[axis];
p->endpoint[i].frequency_ratio = route.frequency_ratio;
p->endpoint[i].end_id = route.end_id;
p->endpoint[i].status = RAY_ENDPOINT_ESCAPED;
p->status[i] = RAY_POOL_TERMINATED;
continue;
}
if (route.kind == ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED) {
p->endpoint[i].status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED;
p->endpoint[i].end_id = route.end_id;
p->status[i] = RAY_POOL_TERMINATED;
continue;
}
if (route.kind != ASYMPTOTIC_ROUTE_INSIDE &&
route.kind != ASYMPTOTIC_ROUTE_ENTRY) {
p->endpoint[i].status = RAY_ENDPOINT_INVALID;
p->status[i] = RAY_POOL_FAILED;
continue;
}
p->activate_t[i] = route.activate_t;
p->x0[i] = route.x[0];
p->x1[i] = route.x[1];
p->x2[i] = route.x[2];
p->p0[i] = route.Pi[0];
p->p1[i] = route.Pi[1];
p->p2[i] = route.Pi[2];
p->log_alpha_p0[i] = route.log_alpha_p0;
}
}
void ray_pool_activate_in_time_range(RayPool *p, const MetricSlab *slab) {
for (size_t i = 0; i < p->count; ++i) {
if (p->status[i] != RAY_POOL_PENDING || p->activate_t[i] > slab->t_hi ||
p->activate_t[i] <= slab->t_lo)
continue;
p->t[i] = p->activate_t[i];
p->steps[i] = 0;
p->status[i] = RAY_POOL_ACTIVE;
}
}
@@ -107,7 +159,8 @@ void ray_pool_destroy(RayPool *p) {
free(p->t); free(p->x0); free(p->x1); free(p->x2); free(p->observer);
free(p->direction0); free(p->direction1); free(p->direction2);
free(p->p0); free(p->p1); free(p->p2); free(p->log_alpha_p0);
free(p->steps); free(p->frame_id); free(p->vertex_id); free(p->status);
free(p->activate_t); free(p->steps); free(p->frame_id); free(p->vertex_id);
free(p->status);
free(p->endpoint);
*p = (RayPool){0};
}
src/ray.h
+6
View File
@@ -15,6 +15,10 @@ typedef enum {
typedef struct {
double *t, *x0, *x1, *x2, *p0, *p1, *p2, *log_alpha_p0;
/* Coordinate time at which the pre-routed interior state becomes valid.
* For a camera inside a worldtube this equals the camera time; for an
* exterior hit it is the earlier entry time. */
double *activate_t;
const ObserverState **observer;
double *direction0, *direction1, *direction2;
unsigned int *steps;
@@ -28,6 +32,8 @@ int ray_pool_init(RayPool *pool, size_t capacity);
int ray_pool_append(RayPool *pool, const ObserverState *observer,
const double direction[3],
size_t frame_id, size_t vertex_id);
/* Pre-route every still-PENDING ray once, before the slab sweep. */
void ray_pool_preroute(RayPool *pool, const SpacetimeSource *source);
void ray_pool_activate_in_time_range(RayPool *pool, const MetricSlab *slab);
void ray_pool_advance_active(RayPool *pool, const MetricSlab *slab,
const GeodesicTraceConfig *config);
src/spacetime.h
+66
View File
@@ -1,6 +1,9 @@
#ifndef SPACETIME_H
#define SPACETIME_H
#include <stddef.h>
#include <stdint.h>
typedef struct {
double alpha;
double beta[3];
@@ -17,6 +20,42 @@ typedef enum {
SPACETIME_RAY_CAPTURED
} SpacetimeRayStatus;
/* Stable identifier for one asymptotic end (infinity) of a backend. Backends
* may describe more than one; the current analytic backends expose one. */
typedef uint32_t SpacetimeEndId;
#define SPACETIME_END_NONE ((SpacetimeEndId)0xffffffffu)
typedef enum {
ASYMPTOTIC_EXTERIOR_MINKOWSKI,
ASYMPTOTIC_EXTERIOR_SCHWARZSCHILD_MONOPOLE
} AsymptoticExteriorKind;
/* Declared asymptotic end. `frame_origin` and the columns of `frame_axes`
* express the asymptotic reference frame in backend coordinates; spatial
* `n_infinity` values use the same coordinate axes as the observer tetrad. */
typedef struct {
SpacetimeEndId end_id;
AsymptoticExteriorKind exterior_kind;
double mass;
double frame_origin[3];
double frame_axes[3][3];
} SpacetimeAsymptoticEnd;
/* Escape worldtube sample at one coordinate time. A zero `radius_rate` and a
* time-independent `velocity` describe the fixed/constant-velocity cases used
* in this phase. `valid == 0` means the backend cannot describe the worldtube
* at this time (history exhausted); callers must not treat that as a miss. */
typedef struct {
double center[3];
double velocity[3];
double radius;
double radius_rate;
/* Nonzero when `velocity` and `radius_rate` are exact throughout the
* current motion segment, so the first entry has a closed quadratic form. */
int velocity_constant;
int valid;
} SpacetimeEscapeWorldtubeSample;
typedef struct SpacetimeSource SpacetimeSource;
typedef struct MetricSlab MetricSlab;
@@ -38,6 +77,20 @@ typedef struct {
MetricData *metric);
SpacetimeRayStatus (*classify_slab)(const MetricSlab *slab, double t,
const double x[3]);
/* Declared asymptotic ends and their moving escape worldtubes. Backends
* without an escape sphere may leave these NULL. */
size_t (*asymptotic_end_count)(const SpacetimeSource *source);
int (*asymptotic_end)(const SpacetimeSource *source, size_t index,
SpacetimeAsymptoticEnd *out);
int (*escape_worldtube_sample)(const SpacetimeSource *source,
SpacetimeEndId end_id, double t,
SpacetimeEscapeWorldtubeSample *out);
/* Coordinate time of the next motion-segment boundary reached while
* integrating backward in time, i.e. the largest boundary strictly less
* than `t`. Return NAN when the worldtube description has a single open
* segment. */
double (*escape_worldtube_next_segment)(const SpacetimeSource *source,
SpacetimeEndId end_id, double t);
/* Analytic backends have negligible per-ray metric state. A numerical
* backend must opt in once its metric slabs and evaluator workspaces need
* to reserve memory alongside the private HDR render buffers. */
@@ -76,6 +129,19 @@ int spacetime_slab_eval(const MetricSlab *slab, double t, const double x[3],
MetricData *metric);
SpacetimeRayStatus spacetime_slab_classify(const MetricSlab *slab, double t,
const double x[3]);
size_t spacetime_asymptotic_end_count(const SpacetimeSource *source);
int spacetime_asymptotic_end(const SpacetimeSource *source, size_t index,
SpacetimeAsymptoticEnd *out);
int spacetime_escape_worldtube_sample(const SpacetimeSource *source,
SpacetimeEndId end_id, double t,
SpacetimeEscapeWorldtubeSample *out);
double spacetime_escape_worldtube_next_segment(const SpacetimeSource *source,
SpacetimeEndId end_id, double t);
/* Common structural validation that every successful constructor must pass
* before returning. A source that passes is a promise that it can safely
* enter ray tracing; backend-specific history/segment validation stays in the
* backend constructor. On failure the constructor must destroy its context. */
int spacetime_source_finalize(SpacetimeSource *source);
int spacetime_limits_render_workers_by_memory(const SpacetimeSource *source);
#endif
src/spacetime_alcubierre.c
+43
View File
@@ -127,9 +127,48 @@ static void alcubierre_destroy(SpacetimeSource *source) {
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,
};
@@ -156,6 +195,10 @@ int spacetime_create_alcubierre(SpacetimeSource *source, double vs,
context->escape_radius = escape_radius;
source->ops = &alcubierre_ops;
source->context = context;
if (spacetime_source_finalize(source)) {
alcubierre_destroy(source);
return -1;
}
return 0;
}
src/spacetime_common.c
+74
View File
@@ -1,5 +1,6 @@
#include "spacetime.h"
#include <math.h>
#include <stddef.h>
#include <stdlib.h>
@@ -64,6 +65,79 @@ SpacetimeRayStatus spacetime_slab_classify(const MetricSlab *slab, double t,
return spacetime_classify(slab->source, t, x);
}
size_t spacetime_asymptotic_end_count(const SpacetimeSource *source) {
return source == NULL || source->ops == NULL ||
source->ops->asymptotic_end_count == NULL
? 0
: source->ops->asymptotic_end_count(source);
}
int spacetime_asymptotic_end(const SpacetimeSource *source, size_t index,
SpacetimeAsymptoticEnd *out) {
return source == NULL || source->ops == NULL || out == NULL ||
source->ops->asymptotic_end == NULL
? -1
: source->ops->asymptotic_end(source, index, out);
}
int spacetime_escape_worldtube_sample(const SpacetimeSource *source,
SpacetimeEndId end_id, double t,
SpacetimeEscapeWorldtubeSample *out) {
if (source == NULL || source->ops == NULL || out == NULL ||
source->ops->escape_worldtube_sample == NULL)
return -1;
return source->ops->escape_worldtube_sample(source, end_id, t, out);
}
double spacetime_escape_worldtube_next_segment(const SpacetimeSource *source,
SpacetimeEndId end_id,
double t) {
if (source == NULL || source->ops == NULL ||
source->ops->escape_worldtube_next_segment == NULL)
return NAN;
return source->ops->escape_worldtube_next_segment(source, end_id, t);
}
int spacetime_source_finalize(SpacetimeSource *source) {
if (source == NULL || source->ops == NULL || source->context == NULL)
return -1;
const SpacetimeOps *ops = source->ops;
if (ops->eval == NULL || ops->classify == NULL || ops->destroy == NULL)
return -1;
const size_t count = spacetime_asymptotic_end_count(source);
if (count == 0)
return 0; /* legacy backend without asymptotic ends */
if (ops->asymptotic_end == NULL || ops->escape_worldtube_sample == NULL)
return -1;
if (count > 64)
return -1;
SpacetimeEndId ids[64];
for (size_t i = 0; i < count; ++i) {
SpacetimeAsymptoticEnd end;
if (spacetime_asymptotic_end(source, i, &end))
return -1;
if (end.end_id == SPACETIME_END_NONE)
return -1;
for (size_t j = 0; j < i; ++j)
if (ids[j] == end.end_id)
return -1;
ids[i] = end.end_id;
if (end.exterior_kind != ASYMPTOTIC_EXTERIOR_MINKOWSKI &&
end.exterior_kind != ASYMPTOTIC_EXTERIOR_SCHWARZSCHILD_MONOPOLE)
return -1;
if (!isfinite(end.mass) || end.mass < 0.0)
return -1;
for (int k = 0; k < 3; ++k) {
if (!isfinite(end.frame_origin[k]))
return -1;
for (int l = 0; l < 3; ++l)
if (!isfinite(end.frame_axes[k][l]))
return -1;
}
}
return 0;
}
int spacetime_limits_render_workers_by_memory(const SpacetimeSource *source) {
return source != NULL && source->ops != NULL &&
source->ops->limit_render_workers_by_memory;
src/spacetime_minkowski.c
+42
View File
@@ -33,9 +33,47 @@ static void minkowski_destroy(SpacetimeSource *source) {
source->ops = NULL;
}
static size_t minkowski_asymptotic_end_count(const SpacetimeSource *source) {
(void)source;
return 1;
}
static int minkowski_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 minkowski_escape_worldtube_sample(
const SpacetimeSource *source, SpacetimeEndId end_id, double t,
SpacetimeEscapeWorldtubeSample *out) {
const MinkowskiContext *context = source->context;
if (end_id != 0)
return -1;
*out = (SpacetimeEscapeWorldtubeSample){.center = {0.0, 0.0, 0.0},
.velocity = {0.0, 0.0, 0.0},
.radius = context->escape_radius,
.radius_rate = 0.0,
.velocity_constant = 1,
.valid = 1};
(void)t;
return 0;
}
static const SpacetimeOps minkowski_ops = {
.eval = minkowski_eval,
.classify = minkowski_classify,
.asymptotic_end_count = minkowski_asymptotic_end_count,
.asymptotic_end = minkowski_asymptotic_end,
.escape_worldtube_sample = minkowski_escape_worldtube_sample,
.destroy = minkowski_destroy,
};
@@ -48,6 +86,10 @@ int spacetime_create_minkowski(SpacetimeSource *source, double escape_radius) {
context->escape_radius = escape_radius;
source->ops = &minkowski_ops;
source->context = context;
if (spacetime_source_finalize(source)) {
minkowski_destroy(source);
return -1;
}
return 0;
}
src/spacetime_schwarzschild.c
+43
View File
@@ -104,9 +104,48 @@ static void schwarzschild_ks_destroy(SpacetimeSource *source) {
source->ops = NULL;
}
static size_t schwarzschild_ks_asymptotic_end_count(
const SpacetimeSource *source) {
(void)source;
return 1;
}
static int schwarzschild_ks_asymptotic_end(
const SpacetimeSource *source, size_t index, SpacetimeAsymptoticEnd *out) {
const SchwarzschildKsContext *context = source->context;
if (index != 0)
return -1;
*out = (SpacetimeAsymptoticEnd){
.end_id = 0,
.exterior_kind = ASYMPTOTIC_EXTERIOR_SCHWARZSCHILD_MONOPOLE,
.mass = context->mass,
.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 schwarzschild_ks_escape_worldtube_sample(
const SpacetimeSource *source, SpacetimeEndId end_id, double t,
SpacetimeEscapeWorldtubeSample *out) {
const SchwarzschildKsContext *context = source->context;
if (end_id != 0)
return -1;
*out = (SpacetimeEscapeWorldtubeSample){.center = {0.0, 0.0, 0.0},
.velocity = {0.0, 0.0, 0.0},
.radius = context->escape_radius,
.radius_rate = 0.0,
.velocity_constant = 1,
.valid = 1};
(void)t;
return 0;
}
static const SpacetimeOps schwarzschild_ks_ops = {
.eval = schwarzschild_ks_eval,
.classify = schwarzschild_ks_classify,
.asymptotic_end_count = schwarzschild_ks_asymptotic_end_count,
.asymptotic_end = schwarzschild_ks_asymptotic_end,
.escape_worldtube_sample = schwarzschild_ks_escape_worldtube_sample,
.destroy = schwarzschild_ks_destroy,
};
@@ -123,6 +162,10 @@ int spacetime_create_schwarzschild_ks(SpacetimeSource *source, double mass,
*context = (SchwarzschildKsContext){mass, escape_radius, capture_radius};
source->ops = &schwarzschild_ks_ops;
source->context = context;
if (spacetime_source_finalize(source)) {
schwarzschild_ks_destroy(source);
return -1;
}
return 0;
}