VPM-B: Add simple Boyle's law compensation.

This is a very basic implementation that uses bin search for solving the cubic.
It's not called anywhere at this stage, needs to be changed to analytic solution.

Signed-off-by: Jan Darowski <jan.darowski@gmail.com>

1 files changed, 48 insertions, 9 deletions

diff --git a/deco.c b/deco.c
index a8012213c..3baefc286 100644
--- a/deco.c
+++ b/deco.c
@@ -114,6 +114,11 @@ double allowable_n2_gradient[16];
 double allowable_he_gradient[16];
 double total_gradient[16];
 
+double bottom_n2_gradient[16];
+double bottom_he_gradient[16];
+
+double initial_n2_gradient[16];
+double initial_he_gradient[16];
 
 static double tissue_tolerance_calc(const struct dive *dive)
 {
@@ -231,8 +236,8 @@ void vpmb_start_gradient()
 	double gradient_n2, gradient_he;
 
 	for (ci = 0; ci < 16; ++ci) {
-		allowable_n2_gradient[ci] = 2.0 * (vpmb_config.surface_tension_gamma / vpmb_config.skin_compression_gammaC) * ((vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma) / n2_regen_radius[ci]);
-		allowable_he_gradient[ci] = 2.0 * (vpmb_config.surface_tension_gamma / vpmb_config.skin_compression_gammaC) * ((vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma) / he_regen_radius[ci]);
+		initial_n2_gradient[ci] = bottom_n2_gradient[ci] = allowable_n2_gradient[ci] = 2.0 * (vpmb_config.surface_tension_gamma / vpmb_config.skin_compression_gammaC) * ((vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma) / n2_regen_radius[ci]);
+		initial_he_gradient[ci] = bottom_he_gradient[ci] = allowable_he_gradient[ci] = 2.0 * (vpmb_config.surface_tension_gamma / vpmb_config.skin_compression_gammaC) * ((vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma) / he_regen_radius[ci]);
 
 		total_gradient[ci] = ((allowable_n2_gradient[ci] * tissue_n2_sat[ci]) + (allowable_he_gradient[ci] * tissue_he_sat[ci])) / (tissue_n2_sat[ci] + tissue_he_sat[ci]);
 	}
@@ -244,19 +249,53 @@ void vpmb_next_gradient(double deco_time)
 	double gradient_n2, gradient_he;
 	double n2_b, n2_c;
 	double he_b, he_c;
-	deco_time /= 60.0 ;
+	double desat_time = deco_time / 60.0;
 
 	for (ci = 0; ci < 16; ++ci) {
-		n2_b = allowable_n2_gradient[ci] + ((vpmb_config.crit_volume_lambda * vpmb_config.surface_tension_gamma) / (vpmb_config.skin_compression_gammaC * (deco_time + buehlmann_N2_t_halflife[ci] * 60.0 / log(2.0))));
-		he_b = allowable_he_gradient[ci] + ((vpmb_config.crit_volume_lambda * vpmb_config.surface_tension_gamma) / (vpmb_config.skin_compression_gammaC * (deco_time + buehlmann_He_t_halflife[ci] * 60.0 / log(2.0))));
+		n2_b = initial_n2_gradient[ci] + (vpmb_config.crit_volume_lambda * vpmb_config.surface_tension_gamma) / (vpmb_config.skin_compression_gammaC * desat_time);
+		he_b = initial_he_gradient[ci] + (vpmb_config.crit_volume_lambda * vpmb_config.surface_tension_gamma) / (vpmb_config.skin_compression_gammaC * desat_time);
 
 		n2_c = vpmb_config.surface_tension_gamma * vpmb_config.surface_tension_gamma * vpmb_config.crit_volume_lambda * max_n2_crushing_pressure[ci];
-		n2_c = n2_c / (vpmb_config.skin_compression_gammaC * vpmb_config.skin_compression_gammaC * (deco_time + buehlmann_N2_t_halflife[ci] * 60.0 / log(2.0)));
+		n2_c = n2_c / (vpmb_config.skin_compression_gammaC * vpmb_config.skin_compression_gammaC * desat_time);
 		he_c = vpmb_config.surface_tension_gamma * vpmb_config.surface_tension_gamma * vpmb_config.crit_volume_lambda * max_he_crushing_pressure[ci];
-		he_c = he_c / (vpmb_config.skin_compression_gammaC * vpmb_config.skin_compression_gammaC * (deco_time + buehlmann_He_t_halflife[ci] * 60.0 / log(2.0)));
+		he_c = he_c / (vpmb_config.skin_compression_gammaC * vpmb_config.skin_compression_gammaC * desat_time);
+
+		bottom_n2_gradient[ci] = allowable_n2_gradient[ci] = 0.5 * ( n2_b + sqrt(n2_b * n2_b - 4.0 * n2_c));
+		bottom_he_gradient[ci] = allowable_he_gradient[ci] = 0.5 * ( he_b + sqrt(he_b * he_b - 4.0 * he_c));
+
+		total_gradient[ci] = ((allowable_n2_gradient[ci] * tissue_n2_sat[ci]) + (allowable_he_gradient[ci] * tissue_he_sat[ci])) / (tissue_n2_sat[ci] + tissue_he_sat[ci]);
+	}
+}
 
-		allowable_n2_gradient[ci] = 0.5 * ( n2_b + sqrt(n2_b * n2_b - 4.0 * n2_c));
-		allowable_he_gradient[ci] = 0.5 * ( he_b + sqrt(he_b * he_b - 4.0 * he_c));
+double update_gradient(double first_stop_pressure, double next_stop_pressure, double first_gradient)
+{
+	double first_radius = 2.0 * vpmb_config.surface_tension_gamma / first_gradient;
+	double A = next_stop_pressure;
+	double B = -2.0 * vpmb_config.surface_tension_gamma;
+	double C = (first_stop_pressure + 2.0 * vpmb_config.surface_tension_gamma / first_radius) * pow(first_radius, 3.0);
+
+	double low = first_radius;
+	double high = first_radius * pow(first_stop_pressure / next_stop_pressure, (1.0/3.0));
+	double next_radius;
+	double value;
+	int ci;
+	for (ci = 0; ci < 100; ++ci){
+		next_radius = (high + low) /2.0;
+		value = A * pow(next_radius, 3.0) - B * next_radius * next_radius - C;
+		if (value < 0)
+			low = next_radius;
+		else
+			high = next_radius;
+	}
+	return 2.0 * vpmb_config.surface_tension_gamma / next_radius;
+}
+
+void boyles_law(double first_stop_pressure, double next_stop_pressure)
+{
+	int ci;
+	for (ci = 0; ci < 16; ++ci) {
+		allowable_n2_gradient[ci] = update_gradient(first_stop_pressure, next_stop_pressure, bottom_n2_gradient[ci]);
+		allowable_he_gradient[ci] = update_gradient(first_stop_pressure, next_stop_pressure, bottom_he_gradient[ci]);
 
 		total_gradient[ci] = ((allowable_n2_gradient[ci] * tissue_n2_sat[ci]) + (allowable_he_gradient[ci] * tissue_he_sat[ci])) / (tissue_n2_sat[ci] + tissue_he_sat[ci]);
 	}