fixed some math, added neutron decay

src/reactor/control/control_rod.cpp

@@ -39,7 +39,7 @@ void control_rod::set_reactivity(double a)
 
 void control_rod::update(double secs)
 {
-	update_rod();
+	update_rod(secs);
 	
 	double m = 1 - std::pow(0.5, (1 - absorbed / limit) * inserted * max);
 	double r_fast = vals[val_t::N_FAST] * m;

src/reactor/coolant/pipe.cpp

@@ -23,7 +23,7 @@ void pipe::update(double secs)
 {
 	double v;
 	
-	update_rod();
+	update_rod(secs);
 	
 	v = vessel->add_steam(vals[val_t::HEAT]);
 	steam = vals[val_t::HEAT] - v;

src/reactor/fuel/fuel_rod.cpp

@@ -32,7 +32,7 @@ double fuel_rod::get_k(val_t type) const
 
 void fuel_rod::update(double secs)
 {
-	update_rod();
+	update_rod(secs);
 	
 	s.add_slow_neutrons(vals[val_t::N_SLOW]);
 

src/reactor/fuel/sample.cpp

@@ -66,7 +66,7 @@ void sample::update(double secs)
 	// sim::reactor::fuelulate fuel use
 	fuel -= neutrons_fuel;
 	energy += neutrons_fuel;
-	waste.add_fissile(neutrons_fuel * 2);
+	waste.add_fissile(neutrons_fuel * 5);
 
 	// do the poison
 	te_135 += neutrons_fuel * (1.0 / 8.0);

src/reactor/fuel/waste.cpp

@@ -4,13 +4,10 @@
 
 using namespace sim::reactor::fuel;
 
-static const double SLOWDOWN_M = 1.0 / 4.0;
-
 void waste::update(double secs)
 {
 	double hl = 1;
-	double next_h[waste::N - 1] = {0};
+	double next[waste::N - 1] = {0};
-	double next_l[waste::N - 1] = {0};
 
 	for(int i = 0; i < waste::N - 1; i++)
 	{
@@ -18,8 +15,7 @@ void waste::update(double secs)
 		double h = high[i] * m;
 		double l = low[i] * m;
 
-		next_h[i] += h * (1 - SLOWDOWN_M);
+		next[i] += l + h;
-		next_l[i] += l + h * SLOWDOWN_M;
 		high[i] -= h;
 		low[i] -= l;
 
@@ -30,8 +26,7 @@ void waste::update(double secs)
 
 	for(int i = 0; i < waste::N - 1; i++)
 	{
-		high[i + 1] += next_h[i];
+		low[i + 1] += next[i];
-		low[i + 1] += next_l[i];
 	}
 }
 

src/reactor/rod.cpp

@@ -15,9 +15,9 @@ void rod::add(val_t type, double v)
 	vals_in[type] += v;
 }
 
-double rod::extract(val_t type, double k, double o)
+double rod::extract(val_t type, double s, double k, double o)
 {
-	double v = (1 - std::pow(0.5, k * get_k(type))) * 0.5 * (get(type) - o);
+	double v = (1 - std::pow(0.5, s * -std::log2(1 - k * get_k(type)))) * 0.5 * (get(type) - o);
 	vals_in[type] -= v;
 	return v;
 }
@@ -27,11 +27,11 @@ void rod::interact(rod* o, double secs)
 	for(int i = 0; i < rod::VAL_N; i++)
 	{
 		val_t v = (val_t)i;
-		add(v, o->extract(v, secs * get_k(v), get(v)));
+		add(v, o->extract(v, secs, get_k(v), get(v)));
 	}
 }
 
-void rod::update_rod()
+void rod::update_rod(double secs)
 {
 	for(int i = 0; i < rod::VAL_N; i++)
 	{
@@ -39,5 +39,10 @@ void rod::update_rod()
 		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;
 }
 

src/reactor/rod.hpp

@@ -22,7 +22,7 @@ public:
 	virtual void interact(rod* o, double secs);
 	virtual void update(double secs) { };
 	virtual void add(val_t type, double v);
-	virtual double extract(val_t type, double k, double o);
+	virtual double extract(val_t type, double s, double k, double o);
 	virtual double get(val_t type) const;
 
 	virtual bool should_display() const { return false; }
@@ -55,7 +55,7 @@ protected:
 	virtual double get_k(val_t type) const { return 0; }
 	virtual const char* get_name() const { return "Empty"; }
 
-	void update_rod();
+	void update_rod(double secs);
 };
 
 }