From 67773f14bd99a848479d68e434dbd69b185a33a2 Mon Sep 17 00:00:00 2001
From: Jay Robson <jsrobson10@gmail.com>
Date: Thu, 18 Jan 2024 18:00:39 +1100
Subject: [PATCH] improved simulation stuff

---
 src/coolant/valve.hpp                |  56 +++++++++++
 src/main.cpp                         | 114 +++++++++++++---------
 src/reactor/control/control_rod.cpp  |  17 ++--
 src/reactor/control/control_rod.hpp  |  12 +--
 src/reactor/control/graphite_rod.cpp |   5 +-
 src/reactor/control/graphite_rod.hpp |   2 +-
 src/reactor/coolant/heater.cpp       |   7 +-
 src/reactor/coolant/heater.hpp       |   2 +-
 src/reactor/coolant/pipe.cpp         |  18 ++--
 src/reactor/coolant/pipe.hpp         |   6 +-
 src/reactor/coolant/vessel.cpp       |  32 +++----
 src/reactor/fuel/fuel_rod.cpp        |  18 ++--
 src/reactor/fuel/fuel_rod.hpp        |   2 +-
 src/reactor/fuel/sample.cpp          |  77 +++++----------
 src/reactor/fuel/sample.hpp          |  19 ++--
 src/reactor/fuel/waste.cpp           |   5 +-
 src/reactor/reactor.hpp              | 137 +++++++++++++++------------
 src/reactor/rod.cpp                  |  45 ++++++---
 src/reactor/rod.hpp                  |  15 ++-
 19 files changed, 335 insertions(+), 254 deletions(-)
 create mode 100644 src/coolant/valve.hpp

diff --git a/src/coolant/valve.hpp b/src/coolant/valve.hpp
new file mode 100644
index 0000000..8498832
--- /dev/null
+++ b/src/coolant/valve.hpp
@@ -0,0 +1,56 @@
+
+#pragma once
+
+namespace sim::coolant
+{
+
+template <class A>
+class valve
+{
+	A* a;
+	
+	const double max;
+	const double pressure;
+
+	double state = 0;
+	double rate = 0;
+
+public:
+
+	constexpr valve(A& a, double max, double pressure) : a(&a), max(max), pressure(pressure)
+	{
+
+	}
+
+	constexpr double get_state() const
+	{
+		return state;
+	}
+
+	constexpr void set_state(double v)
+	{
+		if(v > 1) v = 1;
+		if(v < 0) v = 0;
+		state = v;
+	}
+
+	constexpr void open(double v)
+	{
+		set_state(state + v);
+	}
+
+	constexpr void update(double secs)
+	{
+		rate = a->extract_steam(secs, state * max, pressure) / secs;
+	}
+
+	friend std::ostream& operator<<(std::ostream& o, const valve& v)
+	{
+		o << "Opened: " << (v.state * 100) << " %\n";
+		o << "Rate: " << v.rate << " g/s\n";
+		return o;
+	}
+};
+
+};
+
diff --git a/src/main.cpp b/src/main.cpp
index 5a64e57..691bdbf 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -6,6 +6,7 @@
 #include "reactor/coolant/heater.hpp"
 #include "reactor/coolant/vessel.hpp"
 #include "coolant/fluid_t.hpp"
+#include "coolant/valve.hpp"
 #include "display.hpp"
 
 #include <cmath>
@@ -14,6 +15,9 @@
 #include <curses.h>
 #include <sys/time.h>
 
+const bool do_graphics = true;
+const bool do_clock = true;
+
 unsigned long get_now()
 {
 	struct timeval tv;
@@ -23,29 +27,38 @@ unsigned long get_now()
 
 int main()
 {
-	initscr();
-	cbreak();
-	noecho();
-	keypad(stdscr, TRUE);
-	nodelay(stdscr, TRUE);
-	curs_set(0);
+	std::random_device rd;
+	std::mt19937 rand(rd());
+	
+	if(do_graphics)
+	{
+		initscr();
+		cbreak();
+		noecho();
+		keypad(stdscr, TRUE);
+		nodelay(stdscr, TRUE);
+		curs_set(0);
+	}
 	
 	sim::reactor::coolant::vessel vessel(200, 400, sim::coolant::WATER);
 	sim::reactor::reactor<5, 5> reactor = sim::reactor::builder<5, 5>(
 		sim::reactor::fuel::fuel_rod(2000, 4000),
-		sim::reactor::control::control_rod(10000, 1),
+		sim::reactor::control::control_rod(vessel, 10000, 1),
 		sim::reactor::coolant::pipe(vessel), {
-			"#  ##",
-			"#FCF ",
-			" CHC ",
-			" FCF#",
-			"##  #"
+			"#C#C#",
+			"CFCFC",
+			"#C#C#",
+			"CFCFC",
+			"#C#C#"
 		});
 
+	sim::coolant::valve<sim::reactor::coolant::vessel> valve(vessel, 1, 101000);
+
 	double secs = 0;
 	long clock = get_now();
 	double speed = 1;
 	int framerate = 100;
+	int steps_extra = 1;
 
 	for(;;)
 	{
@@ -80,29 +93,23 @@ secs:		ss << s << "s\n";
 			ss << "Speed: " << speed << "x\n\n";
 		}
 
-		int skip = 1;
-
-		while(speed / framerate / skip > 1)
+		for(int i = 0; i < steps_extra; i++)
 		{
-			skip *= 2;
-		}
-
-		for(int i = 0; i < skip; i++)
-		{
-			double dt = speed / framerate / skip;
-
-			reactor.update(dt);
+			double dt = speed / framerate / steps_extra;
+			reactor.update(rand, dt);
 			vessel.update();
-
-			vessel.extract_steam(dt, 0.01, 101000);
+			valve.update(dt);
+			secs += dt;
 		}
 		
-		secs += speed / framerate;
-		
 		ss << "Vessel\n" << vessel << "\n";
+		ss << "Steam Valve\n" << valve << "\n";
 
-		erase();
-		display::draw_text(1, 0, ss.str().c_str());
+		if(do_graphics)
+		{
+			erase();
+			display::draw_text(1, 0, ss.str().c_str());
+		}
 
 		const int X = 1, Y = 30;
 		const int W = 36, H = 11;
@@ -110,7 +117,8 @@ secs:		ss << s << "s\n";
 		for(int x = 0; x < reactor.width; x++)
 		for(int y = 0; y < reactor.height; y++)
 		{
-			sim::reactor::rod* r = reactor.rods[x][y];
+			int id = y * reactor.width + x;
+			sim::reactor::rod* r = reactor.rods[id];
 
 			if(!r->should_display())
 			{
@@ -123,23 +131,30 @@ secs:		ss << s << "s\n";
 			int px = X + (H - 1) * y;
 			int py = Y + (W - 1) * x;
 
-			display::draw_text(px + 1, py + 2, ss.str().c_str());
-			display::draw_box(px, py, H, W);
+			if(do_graphics)
+			{
+				display::draw_text(px + 1, py + 2, ss.str().c_str());
+				display::draw_box(px, py, H, W);
+			}
 
-			if(r->should_select() && x == reactor.cursor_x && y == reactor.cursor_y)
+			if(do_graphics && r->should_select() && id == reactor.cursor)
 			{
 				display::draw_text(px + 1, py + W - 5, "[ ]");
 			}
 
-			if(r->is_selected())
+			if(do_graphics && r->is_selected())
 			{
 				display::draw_text(px + 1, py + W - 4, "#");
 			}
 		}
 
-		refresh();
+		int c = 0;
 
-		int c = getch();
+		if(do_graphics)
+		{
+			refresh();
+			c = getch();
+		}
 
 		switch(c)
 		{
@@ -150,10 +165,10 @@ secs:		ss << s << "s\n";
 			reactor.move_cursor(1);
 			break;
 		case KEY_UP:
-			reactor.update_selected(0.01);
+			reactor.update_selected(1);
 			break;
 		case KEY_DOWN:
-			reactor.update_selected(-0.01);
+			reactor.update_selected(-1);
 			break;
 		case ' ':
 			reactor.toggle_selected();
@@ -164,17 +179,26 @@ secs:		ss << s << "s\n";
 		case 'g':
 			speed /= 10;
 			break;
+		case 'r':
+			valve.open(0.001);
+			break;
+		case 'f':
+			valve.open(-0.001);
+			break;
 		}
 
-		long now = get_now();
-
-		while(clock + 1e6 / framerate > now)
+		if(do_clock)
 		{
-			usleep(clock + 1e6 / framerate - now);
-			now = get_now();
-		}
+			long now = get_now();
 
-		clock += 1e6 / framerate;
+			while(clock + 1e6 / framerate > now)
+			{
+				usleep(clock + 1e6 / framerate - now);
+				now = get_now();
+			}
+
+			clock += 1e6 / framerate;
+		}
 	}
 
 	return 0;
diff --git a/src/reactor/control/control_rod.cpp b/src/reactor/control/control_rod.cpp
index 2d3142c..e955b98 100644
--- a/src/reactor/control/control_rod.cpp
+++ b/src/reactor/control/control_rod.cpp
@@ -5,10 +5,9 @@
 
 using namespace sim::reactor::control;
 
-control_rod::control_rod(double limit, double max)
+control_rod::control_rod(coolant::vessel& v, double limit, double max) : coolant::pipe(v), limit(limit), max(max)
 {
-	this->limit = limit;
-	this->max = max;
+	
 }
 
 void control_rod::display(std::ostream& o) const
@@ -17,11 +16,6 @@ void control_rod::display(std::ostream& o) const
 	o << "Use: " << absorbed << " / " << limit << " mol\n";
 };
 
-double control_rod::get_k(val_t type) const
-{
-	return 0.5;
-}
-
 void control_rod::set_reactivity(double a)
 {
 	inserted = 1 - a;
@@ -30,14 +24,17 @@ void control_rod::set_reactivity(double a)
 void control_rod::update(double secs)
 {
 	update_rod(secs);
-	
-	double m = 1 - std::pow(0.5, (1 - absorbed / limit) * inserted * max);
+
+	double k = (1 - absorbed / limit) * inserted * max;
+	double m = 1 - std::pow(0.5, secs * -std::log2(1 - k));
 	double r_fast = vals[val_t::N_FAST] * m;
 	double r_slow = vals[val_t::N_SLOW] * m;
 
 	vals[val_t::N_FAST] -= r_fast;
 	vals[val_t::N_SLOW] -= r_slow;
 	absorbed += r_fast + r_slow;
+	
+	update_pipe(secs);
 }
 
 void control_rod::update_selected(double a)
diff --git a/src/reactor/control/control_rod.hpp b/src/reactor/control/control_rod.hpp
index d83c877..806a47d 100644
--- a/src/reactor/control/control_rod.hpp
+++ b/src/reactor/control/control_rod.hpp
@@ -1,26 +1,26 @@
 
 #pragma once
 
-#include "../rod.hpp"
+#include "../coolant/pipe.hpp"
 
 namespace sim::reactor::control
 {
 
-class control_rod : public sim::reactor::rod
+class control_rod : public sim::reactor::coolant::pipe
 {
+	const double limit;
+	const double max;
+
 	double inserted = 1;
 	double absorbed = 0;
-	double limit;
-	double max;
 	
-	virtual double get_k(sim::reactor::rod::val_t type) const;
 	virtual void display(std::ostream& o) const;
 
 	virtual const char* get_name() const { return "Control Rod"; }
 
 public:
 
-	control_rod(double limit, double max);
+	control_rod(coolant::vessel& v, double limit, double max);
 
 	virtual void update(double secs);
 	void set_reactivity(double a);
diff --git a/src/reactor/control/graphite_rod.cpp b/src/reactor/control/graphite_rod.cpp
index 129889a..19f8907 100644
--- a/src/reactor/control/graphite_rod.cpp
+++ b/src/reactor/control/graphite_rod.cpp
@@ -26,9 +26,8 @@ void graphite_rod::update(double secs)
 {
 	update_rod(secs);
 	
-	double v = vals[val_t::N_FAST];
-	vals_in[val_t::N_FAST] -= v;
-	vals_in[val_t::N_SLOW] += v;
+	vals[val_t::N_SLOW] = vals[val_t::N_FAST];
+	vals[val_t::N_FAST] = 0;
 }
 
 void graphite_rod::update_selected(double a)
diff --git a/src/reactor/control/graphite_rod.hpp b/src/reactor/control/graphite_rod.hpp
index a0583ef..202381c 100644
--- a/src/reactor/control/graphite_rod.hpp
+++ b/src/reactor/control/graphite_rod.hpp
@@ -10,10 +10,10 @@ class graphite_rod : public sim::reactor::rod
 {
 	double inserted = 0;
 	
-	virtual double get_k(sim::reactor::rod::val_t type) const;
 	virtual void display(std::ostream& o) const;
 
 	virtual const char* get_name() const { return "Graphite Rod"; }
+	virtual double get_k(val_t type) const;
 
 public:
 
diff --git a/src/reactor/coolant/heater.cpp b/src/reactor/coolant/heater.cpp
index 42acf00..1d3eecb 100644
--- a/src/reactor/coolant/heater.cpp
+++ b/src/reactor/coolant/heater.cpp
@@ -7,7 +7,7 @@ void heater::update(double secs)
 {
 	update_rod(secs);
 
-	vals_in[val_t::HEAT] += rate * secs;
+	vals[val_t::HEAT] += rate * secs;
 }
 
 void heater::update_selected(double a)
@@ -15,11 +15,6 @@ void heater::update_selected(double a)
 	rate += a;
 }
 
-double heater::get_k(sim::reactor::rod::val_t type) const
-{
-	return 0.5;
-}
-
 void heater::display(std::ostream& o) const
 {
 	o << "Rate: " << rate << " C/s\n";
diff --git a/src/reactor/coolant/heater.hpp b/src/reactor/coolant/heater.hpp
index c6d90f4..05cab44 100644
--- a/src/reactor/coolant/heater.hpp
+++ b/src/reactor/coolant/heater.hpp
@@ -10,10 +10,10 @@ class heater : public sim::reactor::rod
 {
 	double rate = 0;
 
-	virtual double get_k(sim::reactor::rod::val_t type) const;
 	virtual void display(std::ostream& o) const;
 
 	virtual const char* get_name() const { return "Heater"; }
+	virtual double get_k(val_t type) const { return 0.5; }
 
 public:
 
diff --git a/src/reactor/coolant/pipe.cpp b/src/reactor/coolant/pipe.cpp
index b0edd86..3e7824e 100644
--- a/src/reactor/coolant/pipe.cpp
+++ b/src/reactor/coolant/pipe.cpp
@@ -16,16 +16,14 @@ double pipe::get_k(val_t type) const
 
 void pipe::update(double secs)
 {
-	double v;
-	
 	update_rod(secs);
-	
-	v = vessel->add_heat(vals[val_t::HEAT]);
-	steam = vals[val_t::HEAT] - v;
-	vals[val_t::HEAT] = v;
-
-	v = vals[val_t::N_FAST];
-	vals_in[val_t::N_FAST] -= v;
-	vals_in[val_t::N_SLOW] += v;
+	update_pipe(secs);
+}
+
+void pipe::update_pipe(double secs)
+{
+	vals[val_t::HEAT] = vessel->add_heat(vals[val_t::HEAT]);
+	vals[val_t::N_SLOW] += vals[val_t::N_FAST];
+	vals[val_t::N_FAST] = 0;
 }
 
diff --git a/src/reactor/coolant/pipe.hpp b/src/reactor/coolant/pipe.hpp
index e0035c2..1e877bd 100644
--- a/src/reactor/coolant/pipe.hpp
+++ b/src/reactor/coolant/pipe.hpp
@@ -9,20 +9,22 @@ namespace sim::reactor::coolant
 
 class pipe : public sim::reactor::rod
 {
+protected:
+	
 	coolant::vessel* vessel;
 	double steam;
 	
 	virtual double get_k(sim::reactor::rod::val_t type) const;
-	
 	virtual const char* get_name() const { return "Coolant"; }
+	void update_pipe(double secs);
 
 public:
 
 	pipe(coolant::vessel& v);
 
 	virtual bool should_display() const { return true; }
-
 	virtual void update(double secs);
+
 };
 
 }
diff --git a/src/reactor/coolant/vessel.cpp b/src/reactor/coolant/vessel.cpp
index 84d8a8b..0c154a7 100644
--- a/src/reactor/coolant/vessel.cpp
+++ b/src/reactor/coolant/vessel.cpp
@@ -21,6 +21,7 @@ void vessel::update()
 
 	double s = steam + n;
 	double l = fluid.l_to_g(level) - n;
+	double v = fluid.l_to_g(volume);
 
 	if(l < 0)
 	{
@@ -28,12 +29,18 @@ void vessel::update()
 		l = 0;
 	}
 
-	if(fluid.g_to_l(s) > volume)
+	if(s > v)
 	{
-		s = fluid.l_to_g(volume);
+		s = v;
 		l = 0;
 	}
 
+	if(l > v)
+	{
+		l = v;
+		s = 0;
+	}
+
 	double diff = s - steam;
 
 	steam = s;
@@ -56,33 +63,22 @@ double vessel::extract_steam(double dt, double a, double p2)
 {
 	// calculate the mass moved
 	double p1 = get_pressure();
-	double p = p1 - p2;
-	double m = 1;
+	double p = (p1 - p2) * 0.001; // mPa or g/m/s^2
 
 	if(p == 0)
 	{
 		return 0;
 	}
 
-	if(p < 0)
-	{
-		m = -1;
-		p = -p;
+	double V = (volume - level) * 0.001; // m^3
+	double mass = std::min(dt * a * p / std::sqrt( V * std::abs(p) / steam ), steam);
 
+	if(std::isnan(mass))
+	{
 		return 0;
 	}
 
-	double V = (volume - level) * 0.001;
-	double v = std::sqrt( V * p / steam );
-	double mass = m * dt * a * p / v;
-
-	if(mass > steam)
-	{
-		mass = steam;
-	}
-
 	steam -= mass;
-
 	return mass;
 }
 
diff --git a/src/reactor/fuel/fuel_rod.cpp b/src/reactor/fuel/fuel_rod.cpp
index 1169061..fc96c2d 100644
--- a/src/reactor/fuel/fuel_rod.cpp
+++ b/src/reactor/fuel/fuel_rod.cpp
@@ -12,26 +12,22 @@ void fuel_rod::display(std::ostream& o) const
 {
 	o << "Fuel: " << s.get_fuel() << " / " << s.get_mass() << " mol\n";
 	o << "Efficiency: " << (s.get_efficiency() * 100) << " %\n";
-	o << "Energy: +" << s.get_energy() << " C\n";
+	o << "Output: " << s.get_energy() << " W/s\n";
 	o << "Iodine: " << s.get_i_135() << " mol\n";
 	o << "Xenon: " << s.get_xe_135() << " mol\n";
 }
 
-double fuel_rod::get_k(val_t type) const
-{
-	return 0.5;
-}
-
 void fuel_rod::update(double secs)
 {
 	update_rod(secs);
 	
+	s.clear_energy();
+	s.clear_fast_neutrons();
 	s.add_slow_neutrons(vals[val_t::N_SLOW]);
-
-	vals_in[val_t::HEAT] += s.extract_energy();
-	vals_in[val_t::N_FAST] += s.extract_fast_neutrons();
-	vals_in[val_t::N_SLOW] -= vals[val_t::N_SLOW];
-	
 	s.update(secs);
+
+	vals[val_t::HEAT] += s.get_energy() * secs;
+	vals[val_t::N_FAST] += s.get_fast_neutrons();
+	vals[val_t::N_SLOW] = 0;
 }
 
diff --git a/src/reactor/fuel/fuel_rod.hpp b/src/reactor/fuel/fuel_rod.hpp
index 3c08b49..892cd55 100644
--- a/src/reactor/fuel/fuel_rod.hpp
+++ b/src/reactor/fuel/fuel_rod.hpp
@@ -11,7 +11,7 @@ class fuel_rod : public sim::reactor::rod
 {
 	sample s;
 	
-	virtual double get_k(sim::reactor::rod::val_t type) const;
+	virtual double get_k(val_t type) const { return 0.5; }
 	virtual void display(std::ostream& o) const;
 	
 	virtual const char* get_name() const { return "Fuel"; }
diff --git a/src/reactor/fuel/sample.cpp b/src/reactor/fuel/sample.cpp
index 43270a8..8ded7d5 100644
--- a/src/reactor/fuel/sample.cpp
+++ b/src/reactor/fuel/sample.cpp
@@ -2,6 +2,8 @@
 #include "sample.hpp"
 #include "half_life.hpp"
 
+#include <iostream>
+
 using namespace sim::reactor::fuel;
 
 static const double NEUTRON_BG = 1e-10;
@@ -12,44 +14,11 @@ sample::sample(double fuel, double mass)
 	this->mass = mass;
 }
 
-void sample::absorb_slow_neutrons()
-{
-	// absorb neutrons
-	double volume = get_volume();
-	double neutrons = slow_neutrons;
-	double neutrons_fuel = neutrons * (fuel / volume);
-	double neutrons_iodine = neutrons * (i_135 / volume);
-	double neutrons_xenon = neutrons * ((xe_135 * Xe_135_M) / volume);
-	double neutrons_total = neutrons;
-	slow_neutrons = 0;
-	
-	xe_135 -= neutrons_xenon;
-	i_135 -= neutrons_iodine;
-	fuel -= neutrons_fuel;
-
-	// do the poison
-	te_135 += neutrons_fuel * (1.0 / 50.0);
-	xe_135 -= neutrons_xenon;
-	i_135 -= neutrons_iodine;
-	
-	// deal with these edge cases
-	if(xe_135 < 0) { slow_neutrons -= xe_135; xe_135 = 0; }
-	if(i_135 < 0) { slow_neutrons -= i_135; i_135 = 0; }
-	if(fuel < 0) { slow_neutrons -= fuel; fuel = 0; }
-
-	efficiency = neutrons_fuel / neutrons_total;
-	
-	// simulate fuel use
-	energy += neutrons_fuel;
-	waste.add_fissile(neutrons_fuel * 6);
-
-}
-
 void sample::update(double secs)
 {
 	double m;
 
-	// sim::reactor::fuelulate waste and extract products
+	// decay waste and extract products
 	waste.update(secs);
 	fast_neutrons += waste.extract_neutrons();
 	energy += waste.extract_energy();
@@ -69,27 +38,29 @@ void sample::update(double secs)
 	te_135 *= m;
 
 	// absorb neutrons
-	slow_neutrons += NEUTRON_BG * secs;
-	absorb_slow_neutrons();
-}
+	double volume = get_volume();
+	double neutrons = slow_neutrons + NEUTRON_BG * secs;
+	double neutrons_fuel = neutrons * (fuel / volume);
+	double neutrons_iodine = neutrons * (i_135 / volume);
+	double neutrons_xenon = neutrons * ((xe_135 * Xe_135_M) / volume);
+	double neutrons_total = neutrons;
+	slow_neutrons = 0;
 
-double sample::extract_energy()
-{
-	double v = energy;
-	energy = 0;
-	return v;
-}
+	te_135 += neutrons_fuel * (1.0 / 50.0);
+	xe_135 -= neutrons_xenon;
+	i_135 -= neutrons_iodine;
+	fuel -= neutrons_fuel;
+	
+	// deal with these edge cases
+	if(xe_135 < 0) { slow_neutrons += xe_135; xe_135 = 0; }
+	if(i_135 < 0) { slow_neutrons += i_135; i_135 = 0; }
+	if(fuel < 0) { slow_neutrons += fuel; fuel = 0; }
 
-double sample::extract_fast_neutrons()
-{
-	double v = fast_neutrons;
-	fast_neutrons = 0;
-	return v;
-}
-
-void sample::add_slow_neutrons(double a)
-{
-	slow_neutrons += a;
+	efficiency = neutrons_fuel / neutrons_total;
+	
+	// simulate fuel use
+	energy += neutrons_fuel / secs;
+	waste.add_fissile(neutrons_fuel * 6);
 }
 
 void sample::display(std::ostream& o) const
diff --git a/src/reactor/fuel/sample.hpp b/src/reactor/fuel/sample.hpp
index 937d3af..c69cacc 100644
--- a/src/reactor/fuel/sample.hpp
+++ b/src/reactor/fuel/sample.hpp
@@ -14,37 +14,40 @@ class sample
 
 	sim::reactor::fuel::waste waste;
 	
+	// mol
 	double fuel = 0;
 	double i_135 = 0;
 	double xe_135 = 0;
 	double te_135 = 0;
 	double mass = 0;
 
-	double energy = 0;
-	double fast_neutrons = 0;
-	double slow_neutrons = 0;
+	double energy = 0; // W/s
+	double fast_neutrons = 0; // mol/s
+	double slow_neutrons = 0; // mol/s
 	double efficiency = 0;
 	
 	void display(std::ostream& o) const;
-	void absorb_slow_neutrons();
+	void absorb_slow_neutrons(double secs);
 
 public:
 
 	sample(double fuel, double mass);
 
 	void update(double secs);
-	double extract_energy();
-	double extract_fast_neutrons();
-	void add_slow_neutrons(double a);
-	
+
 	constexpr double get_fuel() const { return fuel; }
 	constexpr double get_mass() const { return mass; }
 	constexpr double get_energy() const { return energy; }
+	constexpr double get_fast_neutrons() const { return fast_neutrons; }
 	constexpr double get_volume() const { return mass + xe_135 * Xe_135_M; }
 	constexpr double get_efficiency() const { return efficiency; }
 	constexpr double get_te_135() const { return te_135; }
 	constexpr double get_i_135() const { return i_135; }
 	constexpr double get_xe_135() const { return xe_135; }
+	
+	constexpr void clear_energy() { energy = 0; }
+	constexpr void clear_fast_neutrons() { fast_neutrons = 0; }
+	constexpr void add_slow_neutrons(double a) { slow_neutrons += a; }
 
 	friend std::ostream& operator<<(std::ostream& o, const sample& s)
 	{
diff --git a/src/reactor/fuel/waste.cpp b/src/reactor/fuel/waste.cpp
index cc0f919..291d516 100644
--- a/src/reactor/fuel/waste.cpp
+++ b/src/reactor/fuel/waste.cpp
@@ -6,10 +6,10 @@ using namespace sim::reactor::fuel;
 
 void waste::update(double secs)
 {
-	double hl = 1;
 	double next[waste::N - 1] = {0};
+	double hl = 1;
 
-	for(int i = 0; i < waste::N - 1; i++)
+	for(int i = 0; i < waste::N - 1; hl *= 2, i++)
 	{
 		double m = 1 - half_life::get(secs, hl);
 		double h = high[i] * m;
@@ -21,7 +21,6 @@ void waste::update(double secs)
 
 		neutrons += h;
 		energy += h + l;
-		hl *= 2;
 	}
 
 	for(int i = 0; i < waste::N - 1; i++)
diff --git a/src/reactor/reactor.hpp b/src/reactor/reactor.hpp
index d1854b0..856fba4 100644
--- a/src/reactor/reactor.hpp
+++ b/src/reactor/reactor.hpp
@@ -4,6 +4,9 @@
 #include "rod.hpp"
 
 #include <array>
+#include <random>
+#include <algorithm>
+#include <iostream>
 
 namespace sim::reactor
 {
@@ -11,99 +14,113 @@ namespace sim::reactor
 template <int W, int H>
 struct reactor
 {
-	const static int width = W;
-	const static int height = H;
+	constexpr const static int width = W;
+	constexpr const static int height = H;
+	constexpr const static int size = W*H;
 	
-	std::array<std::array<rod*, H>, W> rods;
+	rod* rods[size];
 
-	int cursor_x = 0;
-	int cursor_y = 0;
+	int cursor = 0;
 
-	reactor(std::array<rod*, W * H> rods)
+	reactor(std::array<rod*, size> rods)
 	{
-		for(int y = 0; y < H; y++)
-		for(int x = 0; x < W; x++)
+		for(int i = 0; i < size; i++)
 		{
-			this->rods[x][y] = rods[y * W + x];
+			this->rods[i] = rods[i];
 		}
 	}
 
-	void update(double secs)
+	void update(std::mt19937& rand, double secs)
 	{
-		// do interactions
-		for(int x = 1; x < W; x++)
+		int rods_lookup[size];
+
+		for(int i = 0; i < size; i++)
 		{
-			rods[x][0]->interact(rods[x - 1][0], secs);
+			rods_lookup[i] = i;
 		}
 
-		for(int y = 1; y < H; y++)
+		for(int i = 0; i < size; i++)
 		{
-			rods[0][y]->interact(rods[0][y - 1], secs);
+			rods[i]->update(secs);
 		}
 
-		for(int y = 1; y < H; y++)
-		for(int x = 1; x < W; x++)
-		{
-			rod* r = rods[x][y];
-			r->interact(rods[x - 1][y], secs);
-			r->interact(rods[x][y - 1], secs);
-		}
+		update_interactions(rand, rods_lookup, secs / 2);
+	}
 
-		// do updates
-		for(int y = 0; y < H; y++)
-		for(int x = 0; x < W; x++)
+	void update_selected(int v)
+	{
+		for(int i = 0; i < size; i++)
 		{
-			rods[x][y]->update(secs);
+			rod* r = rods[i];
+
+			if(r->is_selected())
+			{
+				r->update_rod_selected(v);
+			}
 		}
 	}
 
-	void update_selected(double v)
+	int move_cursor(int d)
 	{
-		for(int y = 0; y < H; y++)
-		for(int x = 0; x < W; x++)
-		if(rods[x][y]->is_selected())
+		for(int i = 0; i < size; i++)
 		{
-			rods[x][y]->update_selected(v);
-		}
-	}
-
-	void move_cursor(int d)
-	{
-		for(;;)
-		{
-			cursor_x += d;
-
-			while(cursor_x >= W)
+			cursor = (cursor + d) % size;
+			
+			if(cursor < 0)
 			{
-				cursor_x -= W;
-				cursor_y += 1;
+				cursor += size;
 			}
 
-			while(cursor_x < 0)
+			if(rods[cursor]->should_select())
 			{
-				cursor_x += W;
-				cursor_y -= 1;
-			}
-
-			cursor_y %= H;
-
-			if(cursor_y < 0)
-			{
-				cursor_y += H;
-			}
-
-			if(rods[cursor_x][cursor_y]->should_select())
-			{
-				return;
+				return cursor;
 			}
 		}
+
+		return 0;
 	}
 
 	void toggle_selected()
 	{
-		if(rods[cursor_x][cursor_y]->should_select())
+		if(rods[cursor]->should_select())
 		{
-			rods[cursor_x][cursor_y]->toggle_selected();
+			rods[cursor]->toggle_selected();
+		}
+	}
+
+private:
+
+	void update_tile(std::mt19937& rand, double secs, int i, int x, int y)
+	{
+		std::array<int, 2> nb_lookup[4] = {{-1, 0}, {1, 0}, {0, -1}, {0, 1}};
+		std::shuffle(nb_lookup, &nb_lookup[3], rand);
+
+		for(int j = 0; j < 4; j++)
+		{
+			int xp = x + nb_lookup[j][0];
+			int yp = y + nb_lookup[j][1];
+			
+			if(xp >= 0 && yp >= 0 && xp < width && yp < height)
+			{
+				rods[i]->interact(rods[yp * width + xp], secs / 2);
+			}
+		}
+	}
+
+	void update_interactions(std::mt19937& rand, int* rods_lookup, double secs)
+	{
+		std::shuffle(rods_lookup, &rods_lookup[size - 1], rand);
+
+		for(int id = 0; id < size; id++)
+		{
+			int i = rods_lookup[id];
+			int x = i % width;
+			int y = i / width;
+
+			for(int j = 0; j < 4; j++)
+			{
+				update_tile(rand, secs, i, x, y);
+			}
 		}
 	}
 };
diff --git a/src/reactor/rod.cpp b/src/reactor/rod.cpp
index 7842167..8b5ed22 100644
--- a/src/reactor/rod.cpp
+++ b/src/reactor/rod.cpp
@@ -12,7 +12,7 @@ double rod::get(val_t type) const
 
 void rod::add(val_t type, double v)
 {
-	vals_in[type] += v;
+	vals[type] += v;
 }
 
 double rod::extract(val_t type, double s, double k, double o)
@@ -20,14 +20,16 @@ double rod::extract(val_t type, double s, double k, double o)
 	k *= get_k(type);
 
 	double m = 1;
-
-	if(k < 1)
+	k = 1 - k * get_k(type);
+   
+	if(k > 0)
 	{
-		m = 1 - std::pow(0.5, s * -std::log2(1 - k * get_k(type)));
+		m = 1 - std::pow(0.5, s * -std::log2(k));
 	}
-	
+
 	double v = m * 0.5 * (get(type) - o);
-	vals_in[type] -= v;
+
+	vals[type] -= v;
 	return v;
 }
 
@@ -40,18 +42,35 @@ void rod::interact(rod* o, double secs)
 	}
 }
 
+double rod::get_speed() const
+{
+	int m = motion < 0 ? -1 : 1;
+	return motion == 0 ? 0 : (std::pow(10, std::abs(motion)) * 1e-10 * m);
+}
+
 void rod::update_rod(double secs)
 {
-	for(int i = 0; i < rod::VAL_N; i++)
-	{
-		val_t v = (val_t)i;
-		vals[v] += vals_in[v];
-		vals_in[v] = 0;
-	}
-
 	// decay the free neutrons
 	double m = std::pow(0.5, secs / 879.4);
 	vals[val_t::N_FAST] *= m;
 	vals[val_t::N_SLOW] *= m;
+
+	if(motion != 0 && !is_selected())
+	{
+		motion = 0;
+	}
+
+	if(motion != 0)
+	{
+		update_selected(get_speed() * secs);
+	}
+}
+
+void rod::update_rod_selected(int m)
+{
+	motion += m;
+
+	if(motion > 10) motion = 10;
+	if(motion < -10) motion = -10;
 }
 
diff --git a/src/reactor/rod.hpp b/src/reactor/rod.hpp
index 2c6da60..eac3b50 100644
--- a/src/reactor/rod.hpp
+++ b/src/reactor/rod.hpp
@@ -27,17 +27,24 @@ public:
 
 	virtual bool should_display() const { return false; }
 	virtual bool should_select() const { return false; }
-	virtual void update_selected(double a) { }
+	void update_rod_selected(int m);
 	
 	constexpr void toggle_selected() { selected = !selected; }
-	constexpr bool is_selected() { return selected; }
+	constexpr bool is_selected() const { return selected; }
 
 	friend std::ostream& operator<<(std::ostream& o, const rod& r)
 	{
 		if(!r.should_display()) return o;
 
 		o << r.get_name() << "\n";
+
+		if(r.is_selected())
+		{
+			o << "Speed: " << r.get_speed() << "\n";
+		}
+
 		r.display(o);
+
 		o << "Heat: " << r.get(val_t::HEAT) << " C\n";
 		o << "Fast: " << r.get(val_t::N_FAST) << " mol\n";
 		o << "Slow: " << r.get(val_t::N_SLOW) << " mol\n";
@@ -47,15 +54,17 @@ public:
 
 protected:
 
-	double vals_in[VAL_N] = {0};
 	double vals[VAL_N] = {0};
 	bool selected = false;
+	int motion = 0;
 
 	virtual void display(std::ostream& o) const { };
 	virtual double get_k(val_t type) const { return 0; }
 	virtual const char* get_name() const { return "Empty"; }
+	virtual void update_selected(double a) { }
 
 	void update_rod(double secs);
+	double get_speed() const;
 };
 
 }