use crate::config::parameters::C;
use sciforge::hub::prelude::astronomy::stellar;
#[derive(Debug, Clone, Copy)]
pub struct SpacetimePoint {
pub t: f64,
pub r: f64,
pub theta: f64,
pub phi: f64,
}
impl SpacetimePoint {
pub fn new(t: f64, r: f64, theta: f64, phi: f64) -> Self {
Self { t, r, theta, phi }
}
}
#[derive(Debug, Clone)]
pub struct KerrSpacetime {
pub mass: f64,
pub spin: f64,
rg: f64,
a: f64,
}
impl KerrSpacetime {
pub fn new(mass: f64, spin: f64) -> Self {
let rs = stellar::schwarzschild_radius(mass);
let rg = rs / 2.0;
let a = spin * rg;
Self { mass, spin, rg, a }
}
fn sigma(&self, r: f64, theta: f64) -> f64 {
r * r + self.a * self.a * theta.cos().powi(2)
}
fn delta(&self, r: f64) -> f64 {
r * r - 2.0 * self.rg * r + self.a * self.a
}
fn big_a(&self, r: f64, theta: f64) -> f64 {
let a2 = self.a * self.a;
let del = self.delta(r);
(r * r + a2).powi(2) - a2 * del * theta.sin().powi(2)
}
pub fn event_horizon(&self) -> f64 {
self.rg + (self.rg * self.rg - self.a * self.a).sqrt()
}
pub fn cauchy_horizon(&self) -> f64 {
self.rg - (self.rg * self.rg - self.a * self.a).sqrt()
}
pub fn ergosphere(&self, theta: f64) -> f64 {
self.rg + (self.rg * self.rg - self.a * self.a * theta.cos().powi(2)).sqrt()
}
pub fn metric_tensor(&self, p: &SpacetimePoint) -> [[f64; 4]; 4] {
let sig = self.sigma(p.r, p.theta);
let del = self.delta(p.r);
let big_a = self.big_a(p.r, p.theta);
let sin2 = p.theta.sin().powi(2);
let mut g = [[0.0f64; 4]; 4];
g[0][0] = -(1.0 - 2.0 * self.rg * p.r / sig) * C * C;
g[0][3] = -2.0 * self.rg * p.r * self.a * sin2 * C / sig;
g[3][0] = g[0][3];
g[1][1] = sig / del;
g[2][2] = sig;
g[3][3] = big_a * sin2 / sig;
g
}
pub fn christoffel_radial_tt(&self, r: f64, theta: f64) -> f64 {
let sig = self.sigma(r, theta);
let del = self.delta(r);
self.rg * del * (r * r - self.a * self.a * theta.cos().powi(2)) / (sig * sig * sig)
}
pub fn geodesic_step_rk4(&self, state: &mut [f64; 8], d_lambda: f64) {
let a = self.spin;
let compute_derivs = |s: &[f64; 8]| -> [f64; 8] {
let r = s[1];
let theta = s[2];
let a2 = a * a;
let sig = r * r + a2 * theta.cos().powi(2);
let del = r * r - 2.0 * r + a2;
let dr_sig = 2.0 * r;
let dtheta_sig = -2.0 * a2 * theta.sin() * theta.cos();
let dr_del = 2.0 * r - 2.0;
let mut d = [0.0f64; 8];
d[0] = s[4];
d[1] = s[5];
d[2] = s[6];
d[3] = s[7];
let sig2 = sig * sig;
let sin2 = theta.sin().powi(2);
let cos_th = theta.cos();
let sin_th = theta.sin();
d[4] = -(dr_del * sig - del * dr_sig) / sig2 * s[4] * s[5] / del.abs().max(1e-30)
+ a * sin2 * dr_sig / sig2 * s[5] * s[7]
+ r * a * 2.0 * sin_th * cos_th / sig * s[6] * s[7];
d[5] = -del / sig2 * s[4] * s[4]
+ (dr_del / (2.0 * del.abs().max(1e-30)) - dr_sig / (2.0 * sig))
* del
* s[5]
* s[5]
- dtheta_sig / (2.0 * sig) * del * s[6] * s[6]
+ del / sig * r * sin2 * s[7] * s[7];
d[6] = -dr_sig / (2.0 * sig) * s[5] * s[6] - dtheta_sig / (2.0 * sig) * s[6] * s[6]
+ sin_th * cos_th * s[7] * s[7];
d[7] = -2.0 * s[5] * s[7] / r.abs().max(1e-30)
- 2.0 * cos_th / sin_th.abs().max(1e-30) * s[6] * s[7];
d
};
let k1 = compute_derivs(state);
let mut s2 = [0.0f64; 8];
for i in 0..8 {
s2[i] = state[i] + 0.5 * d_lambda * k1[i];
}
let k2 = compute_derivs(&s2);
for i in 0..8 {
s2[i] = state[i] + 0.5 * d_lambda * k2[i];
}
let k3 = compute_derivs(&s2);
for i in 0..8 {
s2[i] = state[i] + d_lambda * k3[i];
}
let k4 = compute_derivs(&s2);
for i in 0..8 {
state[i] += d_lambda / 6.0 * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
}
}
pub fn trace_photon(
&self,
initial: SpacetimePoint,
direction: [f64; 3],
n_steps: usize,
d_lambda: f64,
) -> Vec<SpacetimePoint> {
let rg = self.rg;
let a = self.spin;
let r_n = initial.r / rg;
let vr_n = direction[0] / C;
let vth_n = direction[1] * rg / C;
let vphi_n = direction[2] * rg / C;
let sig_n = r_n * r_n + a * a * initial.theta.cos().powi(2);
let del_n = r_n * r_n - 2.0 * r_n + a * a;
let g_tt_n = -(1.0 - 2.0 * r_n / sig_n);
let g_rr_n = sig_n / del_n.abs().max(1e-30);
let g_thth_n = sig_n;
let big_a_n = (r_n * r_n + a * a).powi(2) - a * a * del_n * initial.theta.sin().powi(2);
let g_phph_n = big_a_n * initial.theta.sin().powi(2) / sig_n;
let g_tph_n = -2.0 * r_n * a * initial.theta.sin().powi(2) / sig_n;
let spatial_n =
g_rr_n * vr_n * vr_n + g_thth_n * vth_n * vth_n + g_phph_n * vphi_n * vphi_n;
let qa = g_tt_n;
let qb = 2.0 * g_tph_n * vphi_n;
let qc = spatial_n;
let disc = qb * qb - 4.0 * qa * qc;
let dt_n = ((-qb + disc.max(0.0).sqrt()) / (2.0 * qa)).abs();
let metric_scale = (sig_n / del_n.abs().max(1e-30)).sqrt();
let adaptive_dl = d_lambda * metric_scale.clamp(1e-6, 1e6);
let mut state = [
initial.t * C / rg,
r_n,
initial.theta,
initial.phi,
dt_n,
vr_n,
vth_n,
vphi_n,
];
let rh_n = self.event_horizon() / rg;
let mut trajectory = Vec::with_capacity(n_steps);
trajectory.push(initial);
let mut dl = adaptive_dl;
for _ in 0..n_steps {
self.geodesic_step_rk4(&mut state, dl);
if state[1] < rh_n * 1.01 {
trajectory.push(SpacetimePoint::new(
state[0] * rg / C,
state[1] * rg,
state[2],
state[3],
));
break;
}
if state[1] > 1e6 {
break;
}
let local_del = state[1] * state[1] - 2.0 * state[1] + a * a;
let local_sig = state[1] * state[1] + a * a * state[2].cos().powi(2);
dl = d_lambda
* (local_sig / local_del.abs().max(1e-30))
.sqrt()
.clamp(1e-6, 1e6);
trajectory.push(SpacetimePoint::new(
state[0] * rg / C,
state[1] * rg,
state[2],
state[3],
));
}
trajectory
}
}