diff --git a/src/geodesic.c b/src/geodesic.c index 1a554ea..3adeb89 100644 --- a/src/geodesic.c +++ b/src/geodesic.c @@ -1,5 +1,6 @@ #include "geodesic.h" #include +#include typedef GeodesicRayState State; typedef struct { @@ -29,11 +30,11 @@ static int invert(double a[3][3], double b[3][3]) { } /* Equation (4) and (5) of Bohn et al., arXiv:1410.7775. */ -static int rhs(const SpacetimeSource *source, double t, const State *s, +static int rhs(const MetricSlab *slab, double t, const State *s, Derivative *out) { MetricData m; double inv[3][3], up[3] = {0}, da_pi = 0, k_pi_pi = 0; - if (spacetime_eval(source, t, s->x, &m) || m.alpha <= 0 || + if (spacetime_slab_eval(slab, t, s->x, &m) || m.alpha <= 0 || invert(m.gamma, inv)) return -1; for (int i = 0; i < 3; i++) @@ -79,19 +80,19 @@ static State add(const State *s, const Derivative *d, double h) { return r; } -static int rk4(const SpacetimeSource *source, double t, double h, State *s) { +static int rk4(const MetricSlab *slab, double t, double h, State *s) { Derivative a, b, c, d; State q; - if (rhs(source, t, s, &a)) + if (rhs(slab, t, s, &a)) return -1; q = add(s, &a, h / 2); - if (rhs(source, t + h / 2, &q, &b)) + if (rhs(slab, t + h / 2, &q, &b)) return -1; q = add(s, &b, h / 2); - if (rhs(source, t + h / 2, &q, &c)) + if (rhs(slab, t + h / 2, &q, &c)) return -1; q = add(s, &c, h); - if (rhs(source, t + h, &q, &d)) + if (rhs(slab, t + h, &q, &d)) return -1; for (int i = 0; i < 3; i++) { s->x[i] += h * (a.x[i] + 2 * b.x[i] + 2 * c.x[i] + d.x[i]) / 6; @@ -104,13 +105,13 @@ static int rk4(const SpacetimeSource *source, double t, double h, State *s) { return 0; } -int geodesic_initialize_past_ray(const SpacetimeSource *source, +int geodesic_initialize_past_ray(const MetricSlab *slab, const ObserverState *o, const double n[3], State *s) { MetricData m; double k[4] = {o->tetrad[0][0], o->tetrad[0][1], o->tetrad[0][2], o->tetrad[0][3]}; - if (spacetime_eval(source, o->coordinate_time, o->coordinate_position, &m) || + if (spacetime_slab_eval(slab, o->coordinate_time, o->coordinate_position, &m) || m.alpha <= 0) return -1; for (int a = 0; a < 3; a++) @@ -131,11 +132,11 @@ int geodesic_initialize_past_ray(const SpacetimeSource *source, return isfinite(s->log_alpha_p0) ? 0 : -1; } -static int escaped_direction(const SpacetimeSource *source, double t, +static int escaped_direction(const MetricSlab *slab, double t, const State *s, double n[3]) { MetricData m; double inv[3][3], norm = 0; - if (spacetime_eval(source, t, s->x, &m) || invert(m.gamma, inv)) + if (spacetime_slab_eval(slab, t, s->x, &m) || invert(m.gamma, inv)) return -1; for (int i = 0; i < 3; i++) { n[i] = 0; @@ -151,9 +152,9 @@ static int escaped_direction(const SpacetimeSource *source, double t, } GeodesicAdvanceResult geodesic_advance_past_ray( - const SpacetimeSource *source, State *s, double slab_left_time, + const MetricSlab *slab, State *s, double slab_left_time, const GeodesicTraceConfig *config, RayEndpoint *out) { - if (!source || !s || !config || !out || config->coordinate_time_step <= 0 || + if (!slab || !s || !config || !out || config->coordinate_time_step <= 0 || !config->max_steps || !isfinite(slab_left_time) || slab_left_time > s->coordinate_time) return GEODESIC_ADVANCE_FAILED; @@ -164,12 +165,12 @@ GeodesicAdvanceResult geodesic_advance_past_ray( return GEODESIC_ADVANCE_TERMINATED; } SpacetimeRayStatus status = - spacetime_classify(source, s->coordinate_time, s->x); + 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(source, s->coordinate_time, s, out->n_infinity) == 0) + 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; @@ -183,7 +184,7 @@ GeodesicAdvanceResult geodesic_advance_past_ray( } const double h = -fmin(config->coordinate_time_step, s->coordinate_time - slab_left_time); - if (rk4(source, s->coordinate_time, h, s)) + if (rk4(slab, s->coordinate_time, h, s)) return GEODESIC_ADVANCE_FAILED; s->coordinate_time += h; ++s->steps; @@ -199,14 +200,23 @@ RayEndpoint geodesic_trace_past(const SpacetimeSource *source, .magnification = 1, .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 || - geodesic_initialize_past_ray(source, observer, n, &state)) + !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); + return out; + } const double last_time = observer->coordinate_time - config->coordinate_time_step * config->max_steps; - if (geodesic_advance_past_ray(source, &state, last_time, config, &out) == + if (geodesic_advance_past_ray(slab, &state, last_time, config, &out) == GEODESIC_ADVANCE_ACTIVE) out.status = RAY_ENDPOINT_MAX_STEPS; + spacetime_free_slab(slab); return out; } diff --git a/src/geodesic.h b/src/geodesic.h index 5631315..c8da3bd 100644 --- a/src/geodesic.h +++ b/src/geodesic.h @@ -48,12 +48,12 @@ RayEndpoint geodesic_trace_past(const SpacetimeSource *source, const ObserverState *observer, const double camera_direction[3], const GeodesicTraceConfig *config); -int geodesic_initialize_past_ray(const SpacetimeSource *source, +int geodesic_initialize_past_ray(const MetricSlab *slab, const ObserverState *observer, const double camera_direction[3], GeodesicRayState *state); GeodesicAdvanceResult geodesic_advance_past_ray( - const SpacetimeSource *source, GeodesicRayState *state, + const MetricSlab *slab, GeodesicRayState *state, double slab_left_time, const GeodesicTraceConfig *config, RayEndpoint *endpoint); #endif diff --git a/src/main.c b/src/main.c index e59f500..2a836a3 100644 --- a/src/main.c +++ b/src/main.c @@ -255,15 +255,19 @@ static int render_movie(const Settings *s, const StarCatalog *catalog, goto done; for (size_t f = 0; f < movie.frame_count; ++f) for (size_t v = 0; v < movie.frames[f].mesh.vertex_count; ++v) - if (ray_pool_append(&rays, spacetime, &movie.frames[f].observer, + if (ray_pool_append(&rays, &movie.frames[f].observer, movie.frames[f].mesh.vertices[v].camera_direction, f, v)) goto done; double slab_hi = movie.frames[movie.frame_count - 1].coordinate_time; while (ray_pool_has_live(&rays)) { const double slab_lo = slab_hi - s->slab_duration; - ray_pool_activate_in_time_range(&rays, slab_hi, slab_lo); - ray_pool_advance_active(&rays, spacetime, slab_lo, &trace); + MetricSlab *slab = NULL; + if (spacetime_load_slab(spacetime, slab_hi, slab_lo, &slab)) + goto done; + ray_pool_activate_in_time_range(&rays, slab); + ray_pool_advance_active(&rays, slab, &trace); + spacetime_free_slab(slab); slab_hi = slab_lo; } for (size_t i = 0; i < rays.count; ++i) { diff --git a/src/ray.c b/src/ray.c index 7a4aeb4..7b09715 100644 --- a/src/ray.c +++ b/src/ray.c @@ -9,7 +9,8 @@ int ray_pool_init(RayPool *p, size_t capacity) { *p = (RayPool){.capacity = capacity}; #define RAY_ALLOC(field) (p->field = calloc(capacity, sizeof *p->field)) if (!(RAY_ALLOC(t) && RAY_ALLOC(x0) && RAY_ALLOC(x1) && RAY_ALLOC(x2) && - RAY_ALLOC(p0) && RAY_ALLOC(p1) && RAY_ALLOC(p2) && + 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_pool_destroy(p); @@ -19,20 +20,19 @@ int ray_pool_init(RayPool *p, size_t capacity) { return 0; } -int ray_pool_append(RayPool *p, const SpacetimeSource *source, - const ObserverState *observer, const double direction[3], +int ray_pool_append(RayPool *p, const ObserverState *observer, + const double direction[3], size_t frame_id, size_t vertex_id) { if (p == NULL || p->count == p->capacity) return -1; - GeodesicRayState s; const size_t i = p->count; - if (geodesic_initialize_past_ray(source, observer, direction, &s)) + if (observer == NULL || direction == NULL) return -1; - p->t[i] = s.coordinate_time; - p->x0[i] = s.x[0]; p->x1[i] = s.x[1]; p->x2[i] = s.x[2]; - p->p0[i] = s.Pi[0]; p->p1[i] = s.Pi[1]; p->p2[i] = s.Pi[2]; - p->log_alpha_p0[i] = s.log_alpha_p0; - p->steps[i] = s.steps; + p->t[i] = observer->coordinate_time; + p->observer[i] = observer; + p->direction0[i] = direction[0]; + p->direction1[i] = direction[1]; + p->direction2[i] = direction[2]; p->frame_id[i] = frame_id; p->vertex_id[i] = vertex_id; p->status[i] = RAY_POOL_PENDING; @@ -42,14 +42,29 @@ int ray_pool_append(RayPool *p, const SpacetimeSource *source, return 0; } -void ray_pool_activate_in_time_range(RayPool *p, double t_hi, double t_lo) { - for (size_t i = 0; i < p->count; ++i) - if (p->status[i] == RAY_POOL_PENDING && p->t[i] <= t_hi && p->t[i] > t_lo) - p->status[i] = RAY_POOL_ACTIVE; +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->t[i] > slab->t_hi || + p->t[i] <= slab->t_lo) + continue; + GeodesicRayState state; + if (geodesic_initialize_past_ray( + slab, p->observer[i], + (double[]){p->direction0[i], p->direction1[i], p->direction2[i]}, + &state)) { + 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; + p->status[i] = RAY_POOL_ACTIVE; + } } -void ray_pool_advance_active(RayPool *p, const SpacetimeSource *source, - double t_lo, const GeodesicTraceConfig *config) { +void ray_pool_advance_active(RayPool *p, const MetricSlab *slab, + const GeodesicTraceConfig *config) { #pragma omp parallel for schedule(static) for (size_t i = 0; i < p->count; ++i) { if (p->status[i] != RAY_POOL_ACTIVE) @@ -60,7 +75,7 @@ void ray_pool_advance_active(RayPool *p, const SpacetimeSource *source, .log_alpha_p0 = p->log_alpha_p0[i], .steps = p->steps[i]}; const GeodesicAdvanceResult result = - geodesic_advance_past_ray(source, &s, t_lo, config, &p->endpoint[i]); + geodesic_advance_past_ray(slab, &s, slab->t_lo, config, &p->endpoint[i]); p->t[i] = s.coordinate_time; p->x0[i] = s.x[0]; p->x1[i] = s.x[1]; p->x2[i] = s.x[2]; p->p0[i] = s.Pi[0]; p->p1[i] = s.Pi[1]; p->p2[i] = s.Pi[2]; @@ -85,7 +100,8 @@ int ray_pool_has_live(const RayPool *p) { void ray_pool_destroy(RayPool *p) { if (p == NULL) return; - free(p->t); free(p->x0); free(p->x1); free(p->x2); + 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->endpoint); diff --git a/src/ray.h b/src/ray.h index 0248ef1..baf780a 100644 --- a/src/ray.h +++ b/src/ray.h @@ -15,6 +15,8 @@ typedef enum { typedef struct { double *t, *x0, *x1, *x2, *p0, *p1, *p2, *log_alpha_p0; + const ObserverState **observer; + double *direction0, *direction1, *direction2; unsigned int *steps; size_t *frame_id, *vertex_id; uint8_t *status; @@ -23,12 +25,12 @@ typedef struct { } RayPool; int ray_pool_init(RayPool *pool, size_t capacity); -int ray_pool_append(RayPool *pool, const SpacetimeSource *source, - const ObserverState *observer, const double direction[3], +int ray_pool_append(RayPool *pool, const ObserverState *observer, + const double direction[3], size_t frame_id, size_t vertex_id); -void ray_pool_activate_in_time_range(RayPool *pool, double t_hi, double t_lo); -void ray_pool_advance_active(RayPool *pool, const SpacetimeSource *source, - double t_lo, const GeodesicTraceConfig *config); +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); int ray_pool_has_live(const RayPool *pool); void ray_pool_destroy(RayPool *pool); diff --git a/src/spacetime.h b/src/spacetime.h index 01a4051..30e77a0 100644 --- a/src/spacetime.h +++ b/src/spacetime.h @@ -18,12 +18,26 @@ typedef enum { } SpacetimeRayStatus; typedef struct SpacetimeSource SpacetimeSource; +typedef struct MetricSlab MetricSlab; + +struct MetricSlab { + const SpacetimeSource *source; + double t_hi, t_lo; + void *context; +}; typedef struct { int (*eval)(const SpacetimeSource *source, double t, const double x[3], MetricData *metric); SpacetimeRayStatus (*classify)(const SpacetimeSource *source, double t, const double x[3]); + int (*load_slab)(const SpacetimeSource *source, double t_hi, double t_lo, + MetricSlab **out); + void (*free_slab)(MetricSlab *slab); + int (*eval_slab)(const MetricSlab *slab, double t, const double x[3], + MetricData *metric); + SpacetimeRayStatus (*classify_slab)(const MetricSlab *slab, double t, + const double x[3]); void (*destroy)(SpacetimeSource *source); } SpacetimeOps; @@ -44,5 +58,12 @@ int spacetime_eval(const SpacetimeSource *source, double t, const double x[3], MetricData *metric); SpacetimeRayStatus spacetime_classify(const SpacetimeSource *source, double t, const double x[3]); +int spacetime_load_slab(const SpacetimeSource *source, double t_hi, double t_lo, + MetricSlab **out); +void spacetime_free_slab(MetricSlab *slab); +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]); #endif diff --git a/src/spacetime_common.c b/src/spacetime_common.c index be121b5..ba625c9 100644 --- a/src/spacetime_common.c +++ b/src/spacetime_common.c @@ -1,6 +1,7 @@ #include "spacetime.h" #include +#include void spacetime_destroy(SpacetimeSource *source) { if (source != NULL && source->ops != NULL) @@ -20,3 +21,45 @@ SpacetimeRayStatus spacetime_classify(const SpacetimeSource *source, double t, ? SPACETIME_RAY_CAPTURED : source->ops->classify(source, t, x); } + +int spacetime_load_slab(const SpacetimeSource *source, double t_hi, double t_lo, + MetricSlab **out) { + if (source == NULL || source->ops == NULL || out == NULL || t_lo >= t_hi) + return -1; + if (source->ops->load_slab != NULL) + return source->ops->load_slab(source, t_hi, t_lo, out); + MetricSlab *slab = malloc(sizeof *slab); + if (slab == NULL) + return -1; + *slab = (MetricSlab){.source = source, .t_hi = t_hi, .t_lo = t_lo}; + *out = slab; + return 0; +} + +void spacetime_free_slab(MetricSlab *slab) { + if (slab == NULL) + return; + if (slab->source != NULL && slab->source->ops != NULL && + slab->source->ops->free_slab != NULL) + slab->source->ops->free_slab(slab); + else + free(slab); +} + +int spacetime_slab_eval(const MetricSlab *slab, double t, const double x[3], + MetricData *metric) { + if (slab == NULL || t < slab->t_lo || t > slab->t_hi) + return -1; + if (slab->source->ops->eval_slab != NULL) + return slab->source->ops->eval_slab(slab, t, x, metric); + return spacetime_eval(slab->source, t, x, metric); +} + +SpacetimeRayStatus spacetime_slab_classify(const MetricSlab *slab, double t, + const double x[3]) { + if (slab == NULL || t < slab->t_lo || t > slab->t_hi) + return SPACETIME_RAY_CAPTURED; + if (slab->source->ops->classify_slab != NULL) + return slab->source->ops->classify_slab(slab, t, x); + return spacetime_classify(slab->source, t, x); +} diff --git a/tests/test_geodesic.c b/tests/test_geodesic.c index 8822a5a..164be8e 100644 --- a/tests/test_geodesic.c +++ b/tests/test_geodesic.c @@ -27,8 +27,14 @@ static int check_ray(const SpacetimeSource *source, int main(void) { SpacetimeSource source = {0}; + MetricSlab *slab = NULL; + MetricData metric; const ObserverState observer = observer_fixed_at_origin(); - if (spacetime_create_minkowski(&source, 10.0)) + if (spacetime_create_minkowski(&source, 10.0) || + spacetime_load_slab(&source, 0.0, -1.0, &slab) || + spacetime_slab_eval(slab, -0.5, (double[]){0.0, 0.0, 0.0}, &metric) || + metric.alpha != 1.0 || + !spacetime_slab_eval(slab, 0.25, (double[]){0.0, 0.0, 0.0}, &metric)) return 1; int result = check_ray(&source, &observer, (double[]){1.0, 0.0, 0.0}, (double[]){0.0, 0.0, -1.0}) || @@ -57,6 +63,7 @@ int main(void) { } } observer_track_destroy(&accelerated); + spacetime_free_slab(slab); spacetime_destroy(&source); return result; }