paper-half-interval-mobius.rs

1use figures::{ink, save, Board};2use mrlyrs::core::error::Result;3use std::f64::consts::PI;45const NAME: &str = "paper-half-interval-mobius";6const GAMMA: f64 = 0.577_215_664_901_532_9;7const LOW: f64 = 2.0;8const HIGH: f64 = 7.0;9const SAMPLES: usize = 1600;10const FIFTH: u64 = 94939;11const QUARTER: u64 = 3789;1213// CERTIFICATE1415fn rate(base: u64) -> f64 {16    let b = base as f64;17    let fill = (b + 1.0) / 2.0;18    let cp = 2f64.sqrt() - 4.0 / PI;19    let tan = (3.0 * PI / 8.0 + PI / (4.0 * b)).tan();20    let x0 = (b / PI) * ((b + 3.0).ln() + GAMMA + tan.ln()) + cp * (b + 1.0).powi(2) / (8.0 * b);21    (fill + x0 + 0.5 / (PI / (2.0 * b)).sin()) / fill22}2324fn margin(base: u64, bar: f64) -> f64 {25    bar - rate(base).ln() / (base as f64).ln()26}2728fn odd(x: f64) -> u64 {29    let n = x.round() as u64;30    if n % 2 == 1 {31        n32    } else {33        n + 134    }35}3637// FIGURE3839fn main() -> Result<()> {40    assert!(margin(FIFTH - 2, 0.2) < 0.0 && margin(FIFTH, 0.2) > 0.0);41    assert!(margin(QUARTER - 2, 0.25) < 0.0 && margin(QUARTER, 0.25) > 0.0);42    assert!((1.3678e-8..1.3679e-8).contains(&margin(FIFTH, 0.2)));43    assert!((7.9625e-6..7.9626e-6).contains(&margin(QUARTER, 0.25)));44    assert!((0.199_387_1..=0.199_387_2).contains(&(0.2 - margin(100_003, 0.2))));4546    let mut bases: Vec<u64> = (0..=SAMPLES)47        .map(|i| odd(10f64.powf(LOW + (HIGH - LOW) * i as f64 / SAMPLES as f64)))48        .collect();49    bases.extend([QUARTER - 2, QUARTER, FIFTH - 2, FIFTH]);50    bases.sort_unstable();51    bases.dedup();52    for &b in &bases {53        assert_eq!(margin(b, 0.2) > 0.0, b >= FIFTH);54        assert_eq!(margin(b, 0.25) > 0.0, b >= QUARTER);55    }5657    let mut board = Board::square();58    let frame = board.frame(0.08);59    let top = margin(*bases.last().unwrap(), 0.25);60    let bottom = margin(bases[0], 0.2);61    let pad = 0.04 * (top - bottom);62    let (hi, lo) = (top + pad, bottom - pad);63    let at = |b: u64, m: f64| {64        let u = ((b as f64).log10() - LOW) / (HIGH - LOW);65        (66            frame.x + frame.w * u,67            frame.y + frame.h * (hi - m) / (hi - lo),68        )69    };7071    for decade in 3..7 {72        let x = frame.x + frame.w * (decade as f64 - LOW) / (HIGH - LOW);73        board.segment((x, frame.y), (x, frame.y + frame.h), 1.4, ink::line());74    }75    let zero = at(bases[0], 0.0).1;76    board.segment((frame.x, zero), (frame.x + frame.w, zero), 2.2, ink::dim());7778    for (bar, wall, color) in [(0.25, QUARTER, ink::orange()), (0.2, FIFTH, ink::blue())] {79        let pts: Vec<(f64, f64)> = bases.iter().map(|&b| at(b, margin(b, bar))).collect();80        board.polyline(&pts, 5.0, color);81        let (x, y) = at(wall, 0.0);82        board.disc(x, y, 15.0, ink::ground());83        board.disc(x, y, 9.5, color);84    }85    save(NAME, &board)?;86    Ok(())87}