planetsfactory 0.0.4

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

pub fn rk4_step<F: Fn(f64, f64) -> f64>(f: &F, t: f64, y: f64, dt: f64) -> f64 {
    let k1 = f(t, y);
    let k2 = f(t + 0.5 * dt, y + 0.5 * dt * k1);
    let k3 = f(t + 0.5 * dt, y + 0.5 * dt * k2);
    let k4 = f(t + dt, y + dt * k3);
    y + dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
}

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;
        let w = if i % 2 == 0 { 2.0 } else { 4.0 };
        sum += w * 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 format_si(value: f64) -> String {
    let prefixes = [
        (1e24, "Y"),
        (1e21, "Z"),
        (1e18, "E"),
        (1e15, "P"),
        (1e12, "T"),
        (1e9, "G"),
        (1e6, "M"),
        (1e3, "k"),
        (1.0, ""),
        (1e-3, "m"),
        (1e-6, "µ"),
        (1e-9, "n"),
        (1e-12, "p"),
    ];
    let abs = value.abs();
    for &(thresh, prefix) in &prefixes {
        if abs >= thresh {
            return format!("{:.3} {}", value / thresh, prefix);
        }
    }
    format!("{:.3e}", value)
}