blackholesfactory 0.0.4

Black hole factory — create and simulate stellar, intermediate-mass, and supermassive black holes with full Kerr physics
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

pub fn lerp(a: f64, b: f64, t: f64) -> f64 {
    a + (b - a) * t
}

pub fn clamp(val: f64, lo: f64, hi: f64) -> f64 {
    val.clamp(lo, hi)
}

pub fn rk4_step<F: Fn(f64, &[f64]) -> Vec<f64>>(f: &F, t: f64, y: &[f64], dt: f64) -> Vec<f64> {
    let k1 = f(t, y);
    let y2: Vec<f64> = y
        .iter()
        .zip(&k1)
        .map(|(yi, ki)| yi + 0.5 * dt * ki)
        .collect();
    let k2 = f(t + 0.5 * dt, &y2);
    let y3: Vec<f64> = y
        .iter()
        .zip(&k2)
        .map(|(yi, ki)| yi + 0.5 * dt * ki)
        .collect();
    let k3 = f(t + 0.5 * dt, &y3);
    let y4: Vec<f64> = y.iter().zip(&k3).map(|(yi, ki)| yi + dt * ki).collect();
    let k4 = f(t + dt, &y4);
    y.iter()
        .enumerate()
        .map(|(i, yi)| yi + dt / 6.0 * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]))
        .collect()
}

pub fn simpson_integrate<F: Fn(f64) -> f64>(f: &F, a: f64, b: f64, n: usize) -> f64 {
    let n = if n.is_multiple_of(2) { n } else { n + 1 };
    let h = (b - a) / n as f64;
    let mut sum = f(a) + f(b);
    for i in 1..n {
        let x = a + i as f64 * h;
        sum += if i % 2 == 0 { 2.0 } else { 4.0 } * f(x);
    }
    sum * h / 3.0
}

pub fn log_range(start: f64, end: f64, n: usize) -> Vec<f64> {
    let ls = start.ln();
    let le = end.ln();
    (0..n)
        .map(|i| (ls + (le - ls) * i as f64 / (n - 1).max(1) as f64).exp())
        .collect()
}

pub fn spherical_to_cartesian(r: f64, theta: f64, phi: f64) -> (f64, f64, f64) {
    (
        r * theta.sin() * phi.cos(),
        r * theta.sin() * phi.sin(),
        r * theta.cos(),
    )
}

pub fn cartesian_to_spherical(x: f64, y: f64, z: f64) -> (f64, f64, f64) {
    let r = (x * x + y * y + z * z).sqrt();
    let theta = if r > 0.0 { (z / r).acos() } else { 0.0 };
    let phi = y.atan2(x);
    (r, theta, phi)
}