fix for floating point issues

src/coolant/condenser.cpp

@@ -19,11 +19,6 @@ condenser::condenser(fluid_t type, double height, double diameter, double level)
 	this->level = level;
 }
 
-double condenser::get_bubble_hl()
-{
-	return (level / volume) * diameter * 0.5 / fluid.bubble_speed;
-}
-
 void condenser::update(double secs)
 {
 	((sim::coolant::fluid_holder*)this)->update(secs);

src/coolant/condenser.hpp

@@ -10,8 +10,6 @@ class condenser : public fluid_holder
 {
 	const double height;
 	const double diameter;
-	
-	virtual double get_bubble_hl();
 
 public:
 

src/coolant/fluid_holder.cpp

@@ -5,6 +5,7 @@
 #include "../reactor/fuel/half_life.hpp"
 
 #include <cmath>
+#include <iostream>
 
 using namespace sim::coolant;
 
@@ -19,8 +20,12 @@ double fluid_holder::add_heat(double m1, double t1)
 	double t = t1 - t2;
 	double m2 = (fluid.l_to_g(level) + steam) * fluid.jPgk;
 	double m = m1 + m2;
+
+	if(m1 == 0 || m2 == 0)
+		return t1;
 	
    	heat = t1 - t * m2 / m;
+
 	return heat;
 }
 
@@ -124,7 +129,5 @@ void fluid_holder::update(double secs)
 	level = fluid.g_to_l(l);
 	heat -= diff * fluid.jPg / (fluid.l_to_g(level) + steam) / fluid.jPgk;
 
-	if(diff > 0) steam_suspended += diff;
-	steam_suspended *= reactor::fuel::half_life::get(secs, get_bubble_hl());
 }
 

src/coolant/fluid_holder.hpp

@@ -15,7 +15,6 @@ protected:
 
 	double level = 0; // litres
 	double steam = 0; // grams
-	double steam_suspended = 0; // grams
 	double heat = 0; // celsius
 	
 	fluid_holder(fluid_t fluid, double volume);
@@ -33,10 +32,6 @@ public:
 	constexpr double get_steam() const { return steam; } // grams
 	constexpr double get_mass() const { return fluid.l_to_g(level) + steam; } // grams
 	constexpr double get_steam_density() const { return steam / (volume - level); } // g/L
-	constexpr double get_steam_suspended() const { return steam_suspended; } // grams
-	constexpr double get_void_ratio() const { double s = steam_suspended / get_steam_density(); return s / (level + s); }
-
-	virtual double get_bubble_hl() = 0;
 
 	double get_pressure() const; // pascals
 	void update(double dt);

src/electric/turbine.cpp

@@ -19,11 +19,6 @@ turbine::turbine(coolant::fluid_t type, double height, double diameter, double l
 	this->level = level;
 }
 
-double turbine::get_bubble_hl()
-{
-	return (level / volume) * diameter * 0.5 / fluid.bubble_speed;
-}
-
 void turbine::update(double secs)
 {
 	((sim::coolant::fluid_holder*)this)->update(secs);

src/electric/turbine.hpp

@@ -11,8 +11,6 @@ class turbine : public sim::coolant::fluid_holder
 	const double height;
 	const double diameter;
 	
-	virtual double get_bubble_hl();
-
 public:
 
 	turbine(coolant::fluid_t type, double height, double diameter, double level);

src/graphics/mesh/mesh.cpp

@@ -178,7 +178,11 @@ static bool calc_intercept_vert(vec3 v[3], vec3 pos, vec3& path, vec3& path_n, v
 	path -= normal * d;
 	normal_last = normal;
 	l = glm::length(path);
-	path_n = path / l;
+
+	if(l > 0)
+	{
+		path_n = path / l;
+	}
 
 	return true;
 }
@@ -207,6 +211,11 @@ vec3 mesh::calc_intersect(vec3 pos, vec3 path, vec3& normal_last) const
 		{
 			i_found = i;
 		}
+
+		if(l == 0)
+		{
+			return path;
+		}
 	}
 
 	for(unsigned int i = 0; i < i_found; i += 3)
@@ -218,6 +227,11 @@ vec3 mesh::calc_intersect(vec3 pos, vec3 path, vec3& normal_last) const
 		};
 		
 		calc_intercept_vert(v, pos, path, path_n, normal_last, l);
+
+		if(l == 0)
+		{
+			return path;
+		}
 	}
 
 	return path;

src/main.cpp

@@ -4,6 +4,7 @@
 #include <random>
 #include <sstream>
 #include <cmath>
+#include <cfenv>
 
 #include "reactor/coolant/vessel.hpp"
 #include "coolant/fluid_t.hpp"
@@ -29,6 +30,8 @@ unsigned long get_now()
 
 int main()
 {
+	feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
+	
 	graphics::window::create();
 
 	long clock = get_now();

src/reactor/coolant/vessel.cpp

@@ -4,6 +4,7 @@
 #include "../../conversions/temperature.hpp"
 #include "../fuel/half_life.hpp"
 
+#include <iostream>
 #include <cmath>
 
 using namespace sim::reactor::coolant;
@@ -22,13 +23,50 @@ vessel::vessel(sim::coolant::fluid_t fluid, double height, double diameter, doub
 	this->level = level;
 }
 
-double vessel::get_bubble_hl()
+double vessel::get_steam_suspended() const
+{
+	return steam_suspended;
+}
+
+double vessel::get_void_ratio() const
+{
+	double s = steam_suspended / get_steam_density();
+	
+	if(s == 0)
+	{
+		return 0;
+	}
+
+	return s / (level + s);
+}
+
+double vessel::get_bubble_hl() const
 {
 	return (level / volume) * height * 0.5 / fluid.bubble_speed;
 }
 
 void vessel::update(double secs)
 {
+	double s = steam;
+
 	((sim::coolant::fluid_holder*)this)->update(secs);
+
+	double diff = steam - s;
+	double hl = get_bubble_hl();
+
+	if(diff > 0)
+	{
+		steam_suspended += diff;
+	}
+
+	if(hl > 0)
+	{
+		steam_suspended *= reactor::fuel::half_life::get(secs, get_bubble_hl());
+	}
+
+	else
+	{
+		steam_suspended = 0;
+	}
 }
 

src/reactor/coolant/vessel.hpp

@@ -14,10 +14,14 @@ public:
 	
 	const double height; // meters
 	const double diameter; // meters
+	double steam_suspended = 0; // grams
 
 	vessel(sim::coolant::fluid_t fluid, double height, double diameter, double level);
 
-	virtual double get_bubble_hl();
+	double get_steam_suspended() const; // grams
+	double get_void_ratio() const;
+	double get_bubble_hl() const;
+
 	void update(double secs);
 
 	friend std::ostream& operator<<(std::ostream& o, const vessel& v)

src/reactor/reactor.hpp

@@ -14,10 +14,11 @@ struct reactor
 {
 	const double cell_width;
 	const double cell_height;
+	
 	const int width;
 	const int height;
 	const int size;
-
+	
 	std::vector<std::unique_ptr<rod>> rods;
 	double rod_speed = 0;
 	int cursor;
@@ -38,7 +39,7 @@ struct reactor
 	double get_total(rod::val_t type);
 	int move_cursor(int d);
 	void toggle_selected();
-
+	
 private:
 
 	void update_tile(double secs, int i, int x, int y);

src/reactor/rod.cpp

@@ -31,6 +31,7 @@ double rod::extract(val_t type, double s, double k, double o)
 		m = 1 - std::pow(0.5, s * -std::log2(k));
 	}
 
+
 	double v = m * 0.5 * (get(type) - o);
 	vals[type] -= v;
 

src/system.cpp

@@ -37,7 +37,7 @@ system::system()
 	
 	vessel = std::make_unique<reactor::coolant::vessel>(sim::coolant::WATER, 8, 10, 300);
 	reactor = std::make_unique<reactor::reactor>(sim::reactor::builder(19, 19, 1.0 / 4.0, 4, reactor::fuel::fuel_rod(0.5), vessel.get(), layout));
-	condenser = std::make_unique<coolant::condenser>(sim::coolant::WATER, 8, 6, 200);
+	condenser = std::make_unique<coolant::condenser>(sim::coolant::WATER, 8, 6, 20);
 	turbine = std::make_unique<electric::turbine>(sim::coolant::WATER, 6, 3, 20);
 	
 	turbine_inlet_valve = std::make_unique<coolant::valve>(vessel.get(), turbine.get(), 1);