diff --git a/.gitignore b/.gitignore
index bf89451..c228cbc 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,6 @@
 /target
+/tmp
+/test.mp4
 
 # CLion
 **/.idea
\ No newline at end of file
diff --git a/Cargo.lock b/Cargo.lock
index f1a5cda..67456c1 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -320,9 +320,9 @@ dependencies = [
 
 [[package]]
 name = "image"
-version = "0.23.13"
+version = "0.23.14"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "293f07a1875fa7e9c5897b51aa68b2d8ed8271b87e1a44cb64b9c3d98aabbc0d"
+checksum = "24ffcb7e7244a9bf19d35bf2883b9c080c4ced3c07a9895572178cdb8f13f6a1"
 dependencies = [
  "bytemuck",
  "byteorder",
@@ -346,6 +346,7 @@ dependencies = [
  "console",
  "lazy_static",
  "number_prefix",
+ "rayon",
  "regex",
 ]
 
@@ -399,9 +400,9 @@ checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646"
 
 [[package]]
 name = "libc"
-version = "0.2.86"
+version = "0.2.91"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "b7282d924be3275cec7f6756ff4121987bc6481325397dde6ba3e7802b1a8b1c"
+checksum = "8916b1f6ca17130ec6568feccee27c156ad12037880833a3b842a823236502e7"
 
 [[package]]
 name = "libm"
@@ -769,9 +770,9 @@ checksum = "388a1df253eca08550bef6c72392cfe7c30914bf41df5269b68cbd6ff8f570a3"
 
 [[package]]
 name = "serde"
-version = "1.0.123"
+version = "1.0.125"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "92d5161132722baa40d802cc70b15262b98258453e85e5d1d365c757c73869ae"
+checksum = "558dc50e1a5a5fa7112ca2ce4effcb321b0300c0d4ccf0776a9f60cd89031171"
 
 [[package]]
 name = "serde_cbor"
diff --git a/Cargo.toml b/Cargo.toml
index c83a630..6f02c8c 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -5,11 +5,11 @@ authors = ["mindv0rtex <mindv0rtex@users.noreply.github.com>"]
 edition = "2018"
 
 [dependencies]
-chrono = "0.4"
-image = "0.23"
-indicatif = "0.15"
+chrono = "0.4.19"
+image = "0.23.14"
+indicatif = {version = "0.15.0", features = ["rayon"]}
 itertools = "0.10"
-rand = "0.8"
+rand = "0.8.3"
 rand_distr = "0.4"
 rayon = "1.5"
 
@@ -19,3 +19,15 @@ criterion = "0.3"
 [[bench]]
 name = "trig"
 harness = false
+
+[profile.dev]
+codegen-units = 1
+opt-level = 3
+target-cpu = "native"
+lto = "thin"
+
+[profile.release]
+codegen-units = 1
+opt-level = 3
+target-cpu = "native"
+lto = "fat"
\ No newline at end of file
diff --git a/TODO.md b/TODO.md
new file mode 100644
index 0000000..b98e02a
--- /dev/null
+++ b/TODO.md
@@ -0,0 +1,3 @@
+### Todo:
+- Auto create a mp4 from generate images
+    - Instead of using the command `ffmpeg -r 20 -i tmp/out_%d.png -vcodec  libx264 -crf 25 test.mp4` maybe use a rust library to do the same (more research needed)
\ No newline at end of file
diff --git a/src/blur.rs b/src/blur.rs
index fd56d8e..b3bed17 100644
--- a/src/blur.rs
+++ b/src/blur.rs
@@ -6,6 +6,14 @@ pub struct Blur {
     row_buffer: Vec<f32>,
 }
 
+impl Clone for Blur {
+    fn clone(&self) -> Blur {
+        return Blur {
+            row_buffer: self.row_buffer.clone(),
+        }
+    }
+}
+
 impl Blur {
     pub fn new(width: usize) -> Self {
         Blur {
diff --git a/src/grid.rs b/src/grid.rs
index 8b94d05..2979f8a 100644
--- a/src/grid.rs
+++ b/src/grid.rs
@@ -16,6 +16,19 @@ pub struct PopulationConfig {
     deposition_amount: f32,
 }
 
+impl Clone for PopulationConfig {
+    fn clone(&self) -> PopulationConfig {
+        return PopulationConfig {
+            sensor_distance: self.sensor_distance,
+            step_distance: self.step_distance,
+            sensor_angle: self.sensor_angle,
+            rotation_angle: self.rotation_angle,
+            decay_factor: self.decay_factor,
+            deposition_amount: self.deposition_amount,
+        }
+    }
+}
+
 impl Display for PopulationConfig {
     fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
         write!(
@@ -78,6 +91,20 @@ pub struct Grid {
     blur: Blur,
 }
 
+impl Clone for Grid {
+    fn clone(&self) -> Grid {
+        return Grid {
+            config: self.config.clone(), 
+            width: self.width.clone(), 
+            height: self.height.clone(), 
+            data:  self.data.clone(), 
+            buf: self.buf.clone(), 
+            blur: self.blur.clone(),
+        }
+        // return Grid::new();
+    }
+}
+
 impl Grid {
     /// Create a new grid filled with random floats in the [0.0..1.0) range.
     pub fn new<R: Rng + ?Sized>(width: usize, height: usize, rng: &mut R) -> Self {
diff --git a/src/imgdata.rs b/src/imgdata.rs
new file mode 100644
index 0000000..dc7eec0
--- /dev/null
+++ b/src/imgdata.rs
@@ -0,0 +1,53 @@
+use crate::{
+    grid::{combine, Grid, PopulationConfig},
+    palette::{random_palette, Palette},
+};
+
+
+use rayon::iter::{ParallelIterator, IntoParallelIterator};
+
+// for file stuff
+use std::fs;
+use std::io::{BufRead, Write, BufReader};
+use std::fs::File;
+use std::path::Path;
+use std::fs::OpenOptions;
+
+/*
+fn get_resumed_primes(file_path: &str) -> Vec<i32> {
+    let path = Path::new(file_path);
+    let lines = lines_from_file(path);
+
+    let resumed_primes = lines.par_iter().map(|x| {
+        return str::replace(&str::replace(x, "Invalid: ", ""), "Prime: ", "").parse::<i32>().unwrap();
+    }).collect();
+    return resumed_primes;
+}
+*/
+
+// Class for storing data that will be used to create images
+pub struct ImgData {
+    pub grids: Vec<Grid>,
+    pub palette: Palette,
+    pub iteration: i32,
+}
+
+impl Clone for ImgData {
+    fn clone(&self) -> ImgData {
+        return ImgData {
+            grids: self.grids.clone(),
+            palette: self.palette.clone(),
+            iteration: self.iteration.clone(),
+        }
+    }
+}
+
+impl ImgData {
+    pub fn new(in_grids: Vec<Grid>, in_palette: Palette, in_iteration: i32) -> Self {
+        ImgData {
+            grids: in_grids, 
+            palette: in_palette,
+            iteration: in_iteration,
+        }
+    }
+}
\ No newline at end of file
diff --git a/src/lib.rs b/src/lib.rs
index a2f9cd2..c743b7b 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -4,3 +4,4 @@ pub mod model;
 mod palette;
 pub mod trig; // for benchmarking
 mod util;
+mod imgdata; // for storing image data
\ No newline at end of file
diff --git a/src/main.rs b/src/main.rs
index e62bbce..7694354 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -4,33 +4,35 @@ use physarum::model;
 use rand::Rng;
 
 fn main() {
-    let n_iterations = 400;
-    let (width, height) = (1024, 1024);
+    let n_iterations = 16384;
+    // let (width, height) = (512, 512);
+    // let (width, height) = (1024, 1024);
+    let (width, height) = (2048, 2048);
+
     let n_particles = 1 << 22;
+    println!("n_particles: {}", n_particles);
     let diffusivity = 1;
     let mut rng = rand::thread_rng();
 
-    loop {
-        let pb = ProgressBar::new(n_iterations);
-        pb.set_style(
-            ProgressStyle::default_bar()
-                .template(
-                    "{spinner:.green} [{elapsed_precise}] [{bar:40.cyan/blue}] {pos}/{len} ({eta})",
-                )
-                .progress_chars("#>-"),
-        );
+    let pb = ProgressBar::new(n_iterations);
+    pb.set_style(
+        ProgressStyle::default_bar()
+            .template(
+                "{spinner:.green} [{elapsed_precise}] [{bar:40.cyan/blue}] {pos}/{len} ({eta} {percent}%, {per_sec})",
+            )
+            .progress_chars("#>-"),
+    );
 
-        let n_populations = 1 + rng.gen_range(1..4);
-        let mut model = model::Model::new(width, height, n_particles, n_populations, diffusivity);
-        model.print_configurations();
+    let n_populations = 1 + rng.gen_range(1..4);
+    let mut model = model::Model::new(width, height, n_particles, n_populations, diffusivity);
+    model.print_configurations();
 
-        for i in 0..n_iterations {
-            model.step();
-            pb.set_position(i);
-        }
-        pb.finish();
-
-        let now: DateTime<Utc> = Utc::now();
-        model.save_to_image(format!("out_{}.png", now.timestamp()).as_str());
+    for i in 0..n_iterations {
+        model.step();
+        pb.set_position(i);
     }
+    pb.finish();
+
+    model.render_all_imgdata();
+    model.flush_image_data();
 }
diff --git a/src/model.rs b/src/model.rs
index d6b80f5..968514a 100644
--- a/src/model.rs
+++ b/src/model.rs
@@ -1,8 +1,11 @@
 use crate::{
     grid::{combine, Grid, PopulationConfig},
     palette::{random_palette, Palette},
+    imgdata::ImgData,
 };
 
+
+
 use rand::{seq::SliceRandom, Rng};
 use rand_distr::{Distribution, Normal};
 use rayon::prelude::*;
@@ -10,6 +13,11 @@ use rayon::prelude::*;
 use itertools::multizip;
 use std::f32::consts::TAU;
 
+use std::time::{Duration, Instant};
+use rayon::iter::{ParallelIterator, IntoParallelIterator};
+
+use indicatif::{ParallelProgressIterator, ProgressBar, ProgressStyle};
+
 /// A single Physarum agent. The x and y positions are continuous, hence we use floating point
 /// numbers instead of integers.
 #[derive(Debug)]
@@ -67,6 +75,9 @@ pub struct Model {
     iteration: i32,
 
     palette: Palette,
+
+    // List of ImgData to be processed post-simulation into images
+    img_data_vec: Vec<ImgData>,
 }
 
 impl Model {
@@ -123,29 +134,33 @@ impl Model {
             diffusivity,
             iteration: 0,
             palette: random_palette(),
+            img_data_vec: Vec::new(),
         }
     }
 
     fn pick_direction<R: Rng + ?Sized>(center: f32, left: f32, right: f32, rng: &mut R) -> f32 {
         if (center > left) && (center > right) {
-            0.0
+            return 0.0;
         } else if (center < left) && (center < right) {
-            *[-1.0, 1.0].choose(rng).unwrap()
+            return *[-1.0, 1.0].choose(rng).unwrap();
         } else if left < right {
-            1.0
+            return 1.0;
         } else if right < left {
-            -1.0
-        } else {
-            0.0
+            return -1.0;
         }
+        return 0.0;
     }
 
     /// Perform a single simulation step.
     pub fn step(&mut self) {
+        let save_image: bool = true;
+
         // Combine grids
         let grids = &mut self.grids;
         combine(grids, &self.attraction_table);
 
+        println!("Starting tick for all agents...");
+        let agents_tick_time = Instant::now();
         self.agents.par_iter_mut().for_each(|agent| {
             let grid = &grids[agent.population_id];
             let PopulationConfig {
@@ -159,10 +174,14 @@ impl Model {
 
             let xc = agent.x + agent.angle.cos() * sensor_distance;
             let yc = agent.y + agent.angle.sin() * sensor_distance;
-            let xl = agent.x + (agent.angle - sensor_angle).cos() * sensor_distance;
-            let yl = agent.y + (agent.angle - sensor_angle).sin() * sensor_distance;
-            let xr = agent.x + (agent.angle + sensor_angle).cos() * sensor_distance;
-            let yr = agent.y + (agent.angle + sensor_angle).sin() * sensor_distance;
+            
+            let agent_add_sens = agent.angle + sensor_angle;
+            let agent_sub_sens = agent.angle - sensor_angle;
+
+            let xl = agent.x + agent_sub_sens.cos() * sensor_distance;
+            let yl = agent.y + agent_sub_sens.sin() * sensor_distance;
+            let xr = agent.x + agent_add_sens.cos() * sensor_distance;
+            let yr = agent.y + agent_add_sens.sin() * sensor_distance;
 
             // Sense. We sense from the buffer because this is where we previously combined data
             // from all the grid.
@@ -176,6 +195,10 @@ impl Model {
             agent.rotate_and_move(direction, rotation_angle, step_distance, width, height);
         });
 
+        let agents_tick_elapsed = agents_tick_time.elapsed().as_millis();
+        let ms_per_agent: f64 = (agents_tick_elapsed as f64) / (self.agents.len() as f64);
+        println!("Finished tick for all agents. took {}ms\nTime peragent: {}ms", agents_tick_time.elapsed().as_millis(), ms_per_agent);
+
         // Deposit
         for agent in self.agents.iter() {
             self.grids[agent.population_id].deposit(agent.x, agent.y);
@@ -186,15 +209,56 @@ impl Model {
         self.grids.par_iter_mut().for_each(|grid| {
             grid.diffuse(diffusivity);
         });
+
+        /*
+        println!("Saving image...");
+        let image_save_time = Instant::now();
+        self.save_to_image(format!("./tmp/out_{}.png", self.iteration).as_str());
+        println!("Saved image took {}", image_save_time.elapsed().as_millis());
+        */
+        println!("Saving imgdata...");
+        let image_save_time = Instant::now();
+        self.save_image_data();
+        println!("Saved imgdata, took {}", image_save_time.elapsed().as_millis());
+        
         self.iteration += 1;
     }
 
-    /// Output the current trail layer as a grayscale image.
-    pub fn save_to_image(&self, name: &str) {
-        let (width, height) = (self.grids[0].width, self.grids[0].height);
+
+    fn save_image_data(&mut self) {
+        let grids = self.grids.clone();
+        self.img_data_vec.push(ImgData::new(grids, self.palette, self.iteration));
+    }
+
+    pub fn flush_image_data(&mut self) {
+        self.img_data_vec.clear();
+    }
+
+    pub fn render_all_imgdata(&self) {
+        let pb = ProgressBar::new(self.img_data_vec.len() as u64);
+        pb.set_style(ProgressStyle::default_bar().template(
+            "{spinner:.green} [{elapsed_precise}] [{bar:40.cyan/blue}] ({pos}/{len}, {percent}%, {per_sec})",
+        ));
+
+        for img in &self.img_data_vec {
+            Self::save_to_image(img.to_owned());
+            pb.inc(1);
+        }
+        pb.finish();
+
+        /*
+        img_data_list.par_iter().progress_with(pb)
+            .foreach(|&img| {
+                save_to_image(img);
+            });
+        */
+    }
+
+    pub fn save_to_image(imgdata: ImgData) {
+        let (width, height) = (imgdata.grids[0].width, imgdata.grids[0].height);
         let mut img = image::RgbImage::new(width as u32, height as u32);
 
-        let max_values: Vec<_> = self
+        let max_values: Vec<_> = imgdata
             .grids
             .iter()
             .map(|grid| grid.quantile(0.999) * 1.5)
@@ -205,8 +269,7 @@ impl Model {
                 let i = y * width + x;
                 let (mut r, mut g, mut b) = (0.0_f32, 0.0_f32, 0.0_f32);
                 for (grid, max_value, color) in
-                    multizip((&self.grids, &max_values, &self.palette.colors))
-                {
+                    multizip((&imgdata.grids, &max_values, &imgdata.palette.colors)) {
                     let mut t = (grid.data()[i] / max_value).clamp(0.0, 1.0);
                     t = t.powf(1.0 / 2.2); // gamma correction
                     r += color.0[0] as f32 * t;
@@ -220,6 +283,7 @@ impl Model {
             }
         }
 
-        img.save(name).unwrap();
+    
+        img.save(format!("./tmp/out_{}.png", imgdata.iteration).as_str()).unwrap();
     }
 }
diff --git a/src/palette.rs b/src/palette.rs
index 8c6006c..92ca1ae 100644
--- a/src/palette.rs
+++ b/src/palette.rs
@@ -9,14 +9,14 @@ pub fn random_palette() -> Palette {
     let mut rng = thread_rng();
     let mut palette = PALETTES[rng.gen_range(0..PALETTES.len())];
     palette.colors.shuffle(&mut rng);
-    palette
+    return palette;
 }
 
 const fn hex_to_color(c: usize) -> image::Rgb<u8> {
     let r = (c >> 16) & 0xff;
     let g = (c >> 8) & 0xff;
     let b = (c >> 0) & 0xff;
-    image::Rgb::<u8>([r as u8, g as u8, b as u8])
+    return image::Rgb::<u8>([r as u8, g as u8, b as u8]);
 }
 
 const PALETTES: [Palette; 8] = [