wiki-nodal-domains.rs

1use mrlycore::errors::Result;2use mrlyfig::{field, ink, save, Board, Frame, Ramp};3use std::f64::consts::PI;45const RES: usize = 256;6const REACH: usize = 2;78fn main() -> Result<()> {9    let mut board = Board::square();10    let frame = board.frame(0.08);11    let gap = frame.w * 0.06;12    let side = (frame.w - gap) / 2.0;13    let panels = [14        Frame::new(frame.x, frame.y, side, side),15        Frame::new(frame.x + side + gap, frame.y, side, side),16        Frame::new(frame.x, frame.y + side + gap, side, side),17        Frame::new(frame.x + side + gap, frame.y + side + gap, side, side),18    ];19    let modes: [(Box<dyn Fn(f64, f64) -> f64>, u32, usize); 4] = [20        (Box::new(|x, y| mode(1, 1, x, y)), 2, 1),21        (Box::new(|x, y| mode(2, 1, x, y)), 5, 2),22        (Box::new(|x, y| mode(2, 2, x, y)), 8, 4),23        (Box::new(|x, y| mode(1, 3, x, y) + mode(3, 1, x, y)), 10, 2),24    ];25    assert_eq!([index(2), index(5), index(8), index(10)], [1, 2, 4, 5]);26    let ramp = Ramp::new(vec![ink::orange(), ink::ground(), ink::blue()]);27    for (panel, (f, lambda, want)) in panels.iter().zip(modes.iter()) {28        let signs = signs(f);29        let domains = domains(&signs);30        assert_eq!(domains, *want);31        assert!(domains <= index(*lambda));32        let values = nodal(&signs);33        field::draw_range(&mut board, *panel, RES, RES, &values, (-1.0, 1.0), &ramp);34    }35    save("wiki-nodal-domains", &board)?;36    Ok(())37}3839fn mode(m: u32, n: u32, x: f64, y: f64) -> f64 {40    (m as f64 * PI * x).sin() * (n as f64 * PI * y).sin()41}4243fn index(lambda: u32) -> usize {44    let mut below = 0;45    for m in 1..=10 {46        for n in 1..=10 {47            if m * m + n * n < lambda {48                below += 1;49            }50        }51    }52    below + 153}5455fn signs(f: &dyn Fn(f64, f64) -> f64) -> Vec<i8> {56    let mut out = Vec::with_capacity(RES * RES);57    for row in 0..RES {58        for col in 0..RES {59            let x = (col as f64 + 0.5) / RES as f64;60            let y = (row as f64 + 0.5) / RES as f64;61            out.push(if f(x, y) >= 0.0 { 1 } else { -1 });62        }63    }64    out65}6667fn neighbours(i: usize) -> impl Iterator<Item = usize> {68    let (row, col) = (i / RES, i % RES);69    let mut out = Vec::with_capacity(4);70    if row > 0 {71        out.push(i - RES);72    }73    if row + 1 < RES {74        out.push(i + RES);75    }76    if col > 0 {77        out.push(i - 1);78    }79    if col + 1 < RES {80        out.push(i + 1);81    }82    out.into_iter()83}8485fn domains(signs: &[i8]) -> usize {86    let mut seen = vec![false; signs.len()];87    let mut count = 0;88    for start in 0..signs.len() {89        if seen[start] {90            continue;91        }92        count += 1;93        let mut stack = vec![start];94        seen[start] = true;95        while let Some(i) = stack.pop() {96            for j in neighbours(i) {97                if !seen[j] && signs[j] == signs[i] {98                    seen[j] = true;99                    stack.push(j);100                }101            }102        }103    }104    count105}106107fn nodal(signs: &[i8]) -> Vec<f64> {108    (0..signs.len())109        .map(|i| {110            let (row, col) = (i / RES, i % RES);111            let rows = row.saturating_sub(REACH)..(row + REACH + 1).min(RES);112            let cols = col.saturating_sub(REACH)..(col + REACH + 1).min(RES);113            let crossed = rows114                .flat_map(|r| cols.clone().map(move |c| r * RES + c))115                .any(|j| signs[j] != signs[i]);116            if crossed {117                0.0118            } else {119                signs[i] as f64120            }121        })122        .collect()123}