Use field lines to display the electric field

Cargo.toml

@@ -11,6 +11,7 @@ ansi_term = "0.12"
 decorum = "0.3"
 enum_dispatch = "0.3"
 image = "0.23"
+imageproc = "0.21"
 ndarray = "0.13"
 piecewise-linear = "0.1"
 y4m = "0.7"

src/geom.rs

@@ -1,5 +1,5 @@
-use decorum::R64;
-use std::ops::{Add, Neg};
+use decorum::{R64, Real as _};
+use std::ops::{Add, AddAssign, Mul, Neg, Sub};
 
 
 fn real(f: f64) -> R64 {
@@ -15,10 +15,16 @@ pub struct Point {
 impl Add for Point {
     type Output = Point;
     fn add(self, other: Point) -> Point {
-        Point {
-            x: self.x + other.x,
-            y: self.y + other.y,
-        }
+        let mut ret = self.clone();
+        ret += other;
+        ret
+    }
+}
+
+impl AddAssign for Point {
+    fn add_assign(&mut self, other: Point) {
+        self.x += other.x;
+        self.y += other.y;
     }
 }
 
@@ -32,6 +38,23 @@ impl Neg for Point {
     }
 }
 
+impl Sub for Point {
+    type Output = Point;
+    fn sub(self, other: Point) -> Point {
+        self + (-other)
+    }
+}
+
+impl Mul<f64> for Point {
+    type Output = Point;
+    fn mul(self, s: f64) -> Point {
+        Point {
+            x: self.x * s,
+            y: self.y * s,
+        }
+    }
+}
+
 impl Into<(f64, f64)> for Point {
     fn into(self) -> (f64, f64) {
         (self.x(), self.y())
@@ -53,6 +76,21 @@ impl Point {
     pub fn y(&self) -> f64 {
         self.y.into()
     }
+
+    pub fn round(&self) -> Point {
+        Point {
+            x: self.x.round(),
+            y: self.y.round(),
+        }
+    }
+
+    pub fn mag_sq(&self) -> f64 {
+        (self.x*self.x + self.y*self.y).into()
+    }
+
+    pub fn mag(&self) -> f64 {
+        self.mag_sq().sqrt()
+    }
 }
 
 #[derive(Copy, Clone, Debug, Default)]

src/render.rs

@@ -1,6 +1,8 @@
 use ansi_term::Color::RGB;
+use crate::geom::Point;
 use crate::{Material as _, SimState};
 use image::{RgbImage, Rgb};
+use imageproc::{pixelops, drawing};
 use std::convert::TryInto as _;
 use std::fs::File;
 use std::path::PathBuf;
@@ -18,35 +20,94 @@ use y4m;
 //         -(-c).sqrt()
 //     }
 // }
+//
+
+trait SimStateRenderExt {
+    fn to_image(&self) -> RgbImage;
+    fn e_vector(&self, xidx: u32, yidx: u32, size: u32) -> Point;
+}
+
+impl SimStateRenderExt for SimState {
+    fn to_image(&self) -> RgbImage {
+        let w = self.width().try_into().unwrap();
+        let h = self.height().try_into().unwrap();
+        let evec_spacing = 10;
+        let mut image = RgbImage::new(w, h);
+        for y in 0..h {
+            for x in 0..w {
+                let cell = self.get(x as usize, y as usize);
+                //let r = norm_color(cell.bz() * consts::C);
+                //let r = 0;
+                let r = norm_color(cell.mat().mz()*1.0e-2);
+                let b = (55.0*cell.mat().conductivity()).min(255.0) as u8;
+                //let b = 0;
+                //let b = norm_color(cell.ey());
+                //let g = 0;
+                //let g = norm_color(cell.ex());
+                //let g = norm_color(curl(cell.ex(), cell.ey()));
+                let g = norm_color((cell.bz() * 1.0e4).into());
+                image.put_pixel(x, y, Rgb([r, g, b]));
+            }
+        }
+
+        for y in 0..h {
+            for x in 0..w {
+                if x % evec_spacing == 0 && y % evec_spacing == 0 {
+                    let mut vec = self.e_vector(x, y, evec_spacing) * 0.01;
+                    if vec.mag() > evec_spacing as f64 {
+                        vec = vec * (evec_spacing as f64 / vec.mag());
+                    }
+                    let center = Point::new(x as _, y as _) + Point::new(evec_spacing as _, evec_spacing as _)*0.5;
+                    image.draw_field_arrow(center, vec, Rgb([0xff, 0xff, 0xff]));
+                }
+            }
+        }
+        image
+    }
+
+    fn e_vector(&self, xidx: u32, yidx: u32, size: u32) -> Point {
+        let mut e = Point::default();
+        let w = self.width().try_into().unwrap();
+        let h = self.height().try_into().unwrap();
+        let xstart = xidx.min(w);
+        let ystart = yidx.min(h);
+        let xend = (xstart + size).min(w);
+        let yend = (ystart + size).min(h);
+        for y in ystart..yend {
+            for x in xstart..xend {
+                e += self.get(x as usize, y as usize).e();
+            }
+        }
+        let xw = xend - xstart;
+        let yw = yend - ystart;
+        if xw == 0 || yw == 0 {
+            // avoid division by zero
+            Point::new(0.0, 0.0)
+        } else {
+            e * (1.0 / ((xw*yw) as f64))
+        }
+    }
+}
+
+trait ImageRenderExt {
+    fn draw_field_arrow(&mut self, center: Point, rel: Point, color: Rgb<u8>);
+}
+
+impl ImageRenderExt for RgbImage {
+    fn draw_field_arrow(&mut self, center: Point, rel: Point, color: Rgb<u8>) {
+        let start = (center - rel * 0.5).round();
+        let end = (center + rel * 0.5).round();
+        let i_start = (start.x() as _, start.y() as _);
+        let i_end = (end.x() as _, end.y() as _);
+        drawing::draw_antialiased_line_segment_mut(self, i_start, i_end, color, pixelops::interpolate);
+        //drawing::draw_line_segment_mut(self, i_start, i_end, color);
+    }
+}
 
 fn norm_color(v: f64) -> u8 {
     (v * 64.0 + 128.0).max(0.0).min(255.0) as u8
 }
 
-fn simstate_to_image(state: &SimState) -> RgbImage {
-    let w = state.width().try_into().unwrap();
-    let h = state.height().try_into().unwrap();
-    let mut image = RgbImage::new(w, h);
-    for y in 0..h {
-        for x in 0..w {
-            let cell = state.get(x as usize, y as usize);
-            //let r = norm_color(cell.bz() * consts::C);
-            //let r = 0;
-            let r = norm_color(cell.mat().mz()*1.0e-2);
-            let b = (55.0*cell.mat().conductivity()).min(255.0) as u8;
-            //let b = 0;
-            //let b = norm_color(cell.ey());
-            //let g = 0;
-            //let g = norm_color(cell.ex());
-            //let g = norm_color(curl(cell.ex(), cell.ey()));
-            let g = norm_color((cell.bz() * 1.0e4).into());
-            image.put_pixel(x, y, Rgb([r, g, b]));
-        }
-    }
-    image
-}
-
-
 pub trait Renderer {
     fn render(&mut self, state: &SimState);
 }
@@ -76,7 +137,7 @@ pub struct ColorTermRenderer;
 impl Renderer for ColorTermRenderer {
     fn render(&mut self, state: &SimState) {
         let square = "█";
-        let im = simstate_to_image(state);
+        let im = state.to_image();
         let buf: String = im
             .enumerate_rows()
             .map(|(_y, row)| {
@@ -106,7 +167,7 @@ impl Y4MRenderer {
 
 impl Renderer for Y4MRenderer {
     fn render(&mut self, state: &SimState) {
-        let im = simstate_to_image(state);
+        let im = state.to_image();
         if self.encoder.is_none() {
             let writer = File::create(&self.out_path).unwrap();
             self.encoder = Some(y4m::encode(im.width() as usize, im.height() as usize, y4m::Ratio::new(30, 1))

src/sim.rs

@@ -1,4 +1,5 @@
 use crate::{consts};
+use crate::geom::Point;
 use crate::mat::{GenericMaterial, Material};
 
 use decorum::R64;
@@ -115,6 +116,9 @@ impl<M> Cell<M> {
     pub fn ey(&self) -> f64 {
         self.state.ey()
     }
+    pub fn e(&self) -> Point {
+        Point::new(self.ex(), self.ey())
+    }
     pub fn hz(&self) -> f64 {
         self.state.hz()
     }