demo-tube.rs

1use mrlycore::errors::Result;2use mrlyfig::board::{Board, Frame};3use mrlyfig::{ink, plot, save, Grid};4use mrlymath::two;56const CODE: u128 = 495;7const SIDE: usize = 3;8const LEVEL: usize = 5;9const SPAN: usize = 243;10const EPS: f64 = 7.0;11const CELL: f64 = 2.0;12const STEPS: usize = 1200;1314fn lower(f: &[f64], d: &mut [f64], hull: &mut [usize], edge: &mut [f64]) {15    let n = f.len();16    let mut k = 0;17    hull[0] = 0;18    edge[0] = f64::NEG_INFINITY;19    edge[1] = f64::INFINITY;20    for q in 1..n {21        let mut cut;22        loop {23            let p = hull[k];24            cut = ((f[q] + (q * q) as f64) - (f[p] + (p * p) as f64)) / (2 * (q - p)) as f64;25            if k > 0 && cut <= edge[k] {26                k -= 1;27            } else {28                break;29            }30        }31        k += 1;32        hull[k] = q;33        edge[k] = cut;34        edge[k + 1] = f64::INFINITY;35    }36    k = 0;37    for (q, slot) in d.iter_mut().enumerate() {38        while edge[k + 1] < q as f64 {39            k += 1;40        }41        let gap = q as f64 - hull[k] as f64;42        *slot = gap * gap + f[hull[k]];43    }44}4546fn transform(types: &[u8], side: usize) -> Vec<f64> {47    let far = (4 * side * side) as f64;48    let mut square: Vec<f64> = types49        .iter()50        .map(|&b| if b != 0 { 0.0 } else { far })51        .collect();52    let mut lane = vec![0.0; side];53    let mut out = vec![0.0; side];54    let mut hull = vec![0usize; side];55    let mut edge = vec![0.0; side + 1];56    for col in 0..side {57        for row in 0..side {58            lane[row] = square[row * side + col];59        }60        lower(&lane, &mut out, &mut hull, &mut edge);61        for row in 0..side {62            square[row * side + col] = out[row];63        }64    }65    for row in 0..side {66        lane.copy_from_slice(&square[row * side..(row + 1) * side]);67        lower(&lane, &mut out, &mut hull, &mut edge);68        square[row * side..(row + 1) * side].copy_from_slice(&out);69    }70    square.iter().map(|v| v.sqrt()).collect()71}7273fn limit(t: f64) -> f64 {74    let dimension = 8f64.ln() / 3f64.ln();75    let (flat, slope, bend) = if t < 0.5 {76        (1.0, 4.0 / 5.0, -4.0 / 7.0)77    } else {78        (9.0 / 8.0, 3.0 / 10.0, -1.0 / 14.0)79    };80    t.powf(dimension - 2.0) * (flat + slope * t + bend * t * t)81}8283fn main() -> Result<()> {84    let cells = two::create(CODE, SIDE, LEVEL, 0, SIDE)?;85    let types = cells.types().bytes().to_vec();86    assert_eq!(types.len(), SPAN * SPAN);87    assert_eq!(types.iter().filter(|&&b| b != 0).count(), 32768);88    let dist = transform(&types, SPAN);89    assert_eq!(dist.iter().filter(|&&v| v > 0.0 && v <= EPS).count(), 20440);9091    let seam = limit(1.0 / 3.0);92    assert!((seam - 379.0 / 280.0).abs() < 1e-12);93    assert!((limit(1.0) - seam).abs() < 1e-12);94    let phases: Vec<f64> = (0..=STEPS)95        .map(|i| 1.0 / 3.0 + (2.0 / 3.0) * i as f64 / STEPS as f64)96        .collect();97    let profile: Vec<f64> = phases.iter().map(|&t| limit(t)).collect();98    let high = profile.iter().copied().fold(f64::MIN, f64::max);99    let low = profile.iter().copied().fold(f64::MAX, f64::min);100    assert!(((high - low) / low - 0.003662).abs() < 1e-5);101102    let mut board = Board::square();103    let margin = (board.width as f64 * 0.08).round();104    let plate = CELL * SPAN as f64;105    let sheet = Frame::new(margin, margin, plate, plate);106    let lattice = Grid::new(sheet, SPAN, SPAN, 0.0);107    for row in 0..SPAN {108        for col in 0..SPAN {109            let at = row * SPAN + col;110            if types[at] != 0 {111                lattice.fill(&mut board, col, row, ink::fg());112            } else if dist[at] <= EPS {113                lattice.fill(&mut board, col, row, ink::blue());114            }115        }116    }117    let tall = 320.0;118    let panel = Frame::new(119        board.width as f64 - margin - plate,120        board.height as f64 - margin - tall,121        plate,122        tall,123    );124    plot::axis(&mut board, panel, ink::line());125    plot::curve(126        &mut board,127        panel.inset(24.0),128        &phases,129        &profile,130        4.0,131        ink::yellow(),132    );133    save("demo-tube", &board)?;134    Ok(())135}