still a buggy mess

nasty_function.c++

@@ -1,51 +0,0 @@
-
-#include <cmath>
-
-double calc(double& val1, double& val2)
-{
-	double J_m2 = J_m*J_m;
-	double L_m2 = L_m*L_m;
-	double V_t2 = V_t*V_t;
-	double n_g2 = n_g*n_g;
-	double n_l2 = n_l*n_l;
-	double K2 = K*K;
-	double R2 = R*R;
-
-	double part1 = (
-			-J_m*R*n_g +
-			2*K*R*n_g +
-			K*R*n_l -
-			L_m*n_g -
-			V_t
-	);
-
-	double part2 = std::sqrt(
-			J_m2*R2*n_g2 +
-			2*J_m*K*R2*n_g*n_l +
-			2*J_m*L_m*R*n_g2 +
-			4*J_m*L_m*R*n_g*n_l +
-			4*J_m*L_m*R*n_l2 -
-			2*J_m*R*V_t*n_g -
-			4*J_m*R*V_t*n_l +
-			K2*R2*n_l2 -
-			2*K*L_m*R*n_g*n_l -
-			4*K*L_m*R*n_l2 +
-			2*K*R*V_t*n_l +
-			L_m2*n_g2 +
-			4*L_m2*n_g*n_l +
-			4*L_m2*n_l2 -
-			2*L_m*V_t*n_g -
-			4*L_m*V_t*n_l +
-			V_t2
-	);
-
-	double part3 = (2*(
-			J_m*R -
-			K*R +
-			L_m
-	));
-	
-	val1 = (part1 - part2) / part3;
-	val2 = (part1 + part2) / part3;
-}
-

src/coolant/fluid_holder.cpp

@@ -6,9 +6,12 @@
 
 #include <cmath>
 #include <iostream>
+#include <complex>
 
 using namespace sim::coolant;
 
+typedef std::complex<double> complex;
+
 fluid_holder::fluid_holder(fluid_t fluid, double volume, double extra_mass) : fluid(fluid), volume(volume), extra_mass(extra_mass)
 {
 
@@ -118,28 +121,37 @@ void fluid_holder::update(double secs)
 
 	if(mass > 0)
 	{
-		// use ideal gas law to get target steam density in mol/L
+		// use ideal gas law to get target steam pressure
 		double heat_k = conversions::temperature::c_to_k(heat);
 		double target_pressure = fluid.vapor_pressure.calc_p(heat);
-		double density = target_pressure / (constants::R * heat_k) / 1000;
-		
-		double m_c = fluid.l_to_mol(1);
-		double n_t = fluid.l_to_mol(level) + fluid.g_to_mol(steam);
-		double v_l = (n_t - density * volume) / (m_c - density);
-		double n_l = fluid.l_to_mol(v_l);
 
-		if(n_l < 0)
+		double K = heat_k;
+		double R = 1000 * constants::R;
+		double P = target_pressure;
+		double J_m = fluid.jPg * fluid.gPmol;
+		double L_m = fluid.gPmol / fluid.gPl;
+		double n_g = fluid.g_to_mol(steam);
+		double n_l = fluid.l_to_mol(level);
+		double V_t = volume;
+
+		double n = (-K*R*n_g - L_m*P*n_l + P*V_t)/(K*R - L_m*P) * 0.5;
+
+		/*if(std::abs(n.imag()) > std::numeric_limits<double>::epsilon())
 		{
-			v_l = 0;
+			throw std::runtime_error("Nonzero imaginary component");
-			n_l = 0;
+		}*/
+
+		double l = level - fluid.mol_to_l(n);
+
+		if(l < 0)
+		{
+			n -= fluid.l_to_mol(l);
+			l = 0;
 		}
 
-		double n_diff = n_l - fluid.l_to_mol(level);
+		level = l;
-		double steam_add = -fluid.mol_to_g(n_diff);
+		steam += fluid.mol_to_g(n);
-
+		heat -= n * J_m / mass;
-		level += fluid.mol_to_l(n_diff);
-		steam += steam_add;
-		heat -= steam_add * fluid.jPg / mass;
 	}
 }
 

src/tests.cpp

@@ -24,11 +24,9 @@ void tests::run()
 {
 	fluid_holder fhs[] = {
 		fluid_holder(WATER, 75398, 0),
-		fluid_holder(WATER, 75398, 0),
 	};
 
 	valve vs[] = {
-		valve(&fhs[0], &fhs[1], 1, 0.1),
 	};
 
 	fhs[0].level = 100;