research-dimensions.rs

1use mrlycore::errors::Result;2use mrlyfig::{ink, plot, save, Board};34const RE_LO: f64 = -1.0;5const RE_HI: f64 = 1.0;6const IM_REACH: f64 = 40.0;78fn cexp(z: (f64, f64)) -> (f64, f64) {9    let e = z.0.exp();10    (e * z.1.cos(), e * z.1.sin())11}1213fn cdiv(a: (f64, f64), b: (f64, f64)) -> (f64, f64) {14    let d = b.0 * b.0 + b.1 * b.1;15    ((a.0 * b.0 + a.1 * b.1) / d, (a.1 * b.0 - a.0 * b.1) / d)16}1718fn power(base: f64, s: (f64, f64)) -> (f64, f64) {19    let l = base.ln();20    cexp((-s.0 * l, -s.1 * l))21}2223fn residual(s: (f64, f64)) -> (f64, f64) {24    let a = power(3.0, s);25    let b = power(5.0, s);26    (a.0 + b.0 - 1.0, a.1 + b.1)27}2829fn slope(s: (f64, f64)) -> (f64, f64) {30    let a = power(3.0, s);31    let b = power(5.0, s);32    let (l3, l5) = (3f64.ln(), 5f64.ln());33    (-l3 * a.0 - l5 * b.0, -l3 * a.1 - l5 * b.1)34}3536fn newton(seed: (f64, f64)) -> Option<(f64, f64)> {37    let mut s = seed;38    for _ in 0..80 {39        if s.0.abs() > 6.0 || s.1.abs() > 200.0 {40            return None;41        }42        let d = slope(s);43        if d.0.hypot(d.1) < 1e-14 {44            return None;45        }46        let step = cdiv(residual(s), d);47        s = (s.0 - step.0, s.1 - step.1);48    }49    let r = residual(s);50    if r.0.hypot(r.1) > 1e-10 {51        return None;52    }53    if s.0 < RE_LO || s.0 > RE_HI || s.1.abs() > IM_REACH {54        return None;55    }56    Some(s)57}5859fn control_poles() -> Vec<(f64, f64)> {60    let mut out: Vec<(f64, f64)> = Vec::new();61    for i in 0..=80 {62        for j in 0..=320 {63            let seed = (64                RE_LO + (RE_HI - RE_LO) * i as f64 / 80.0,65                -IM_REACH + 2.0 * IM_REACH * j as f64 / 320.0,66            );67            if let Some(root) = newton(seed) {68                if !out69                    .iter()70                    .any(|p| (p.0 - root.0).hypot(p.1 - root.1) < 1e-6)71                {72                    out.push(root);73                }74            }75        }76    }77    out78}7980fn main() -> Result<()> {81    let real = 2f64.ln() / 3f64.ln();82    let omega = 2.0 * std::f64::consts::PI / 3f64.ln();83    let lattice: Vec<(f64, f64)> = (-6..=6).map(|m| (real, m as f64 * omega)).collect();84    let control = control_poles();85    assert_eq!(lattice.len(), 13);86    assert_eq!(control.len(), 21);8788    let mut board = Board::square();89    let frame = board.frame(0.08);90    plot::axis(&mut board, frame, ink::line());91    let at = |re: f64, im: f64| {92        (93            frame.x + frame.w * (re - RE_LO) / (RE_HI - RE_LO),94            frame.y + frame.h * (1.0 - (im + IM_REACH) / (2.0 * IM_REACH)),95        )96    };97    board.segment(at(0.0, -IM_REACH), at(0.0, IM_REACH), 1.6, ink::line());98    board.segment(at(RE_LO, 0.0), at(RE_HI, 0.0), 1.6, ink::line());99    board.segment(100        at(real, -IM_REACH),101        at(real, IM_REACH),102        1.8,103        ink::fade(ink::blue(), 0.35),104    );105    for pole in &control {106        let (x, y) = at(pole.0, pole.1);107        board.disc(x, y, 7.0, ink::orange());108    }109    for pole in &lattice {110        let (x, y) = at(pole.0, pole.1);111        board.disc(x, y, 9.0, ink::blue());112    }113    save("research-dimensions", &board)?;114    Ok(())115}