VPM-B: use an analytic solution for nucleon inner pressure instead of binary root search

According to mathematica

In[4]:= f[x_] := x^3 - b x^2 - c

In[18]:= Solve[f[x] == 0, x]

Out[18]= {{x ->
   1/3 (b + (
      2^(1/3) b^2)/(2 b^3 + 27 c + 3 Sqrt[3] Sqrt[4 b^3 c + 27 c^2])^(
      1/3) + (2 b^3 + 27 c + 3 Sqrt[3] Sqrt[4 b^3 c + 27 c^2])^(1/3)/
      2^(1/3))}, {x ->
   b/3 - ((1 + I Sqrt[3]) b^2)/(
    3 2^(2/3) (2 b^3 + 27 c + 3 Sqrt[3] Sqrt[4 b^3 c + 27 c^2])^(
     1/3)) - ((1 - I Sqrt[3]) (2 b^3 + 27 c +
       3 Sqrt[3] Sqrt[4 b^3 c + 27 c^2])^(1/3))/(6 2^(1/3))}, {x ->
   b/3 - ((1 - I Sqrt[3]) b^2)/(
    3 2^(2/3) (2 b^3 + 27 c + 3 Sqrt[3] Sqrt[4 b^3 c + 27 c^2])^(
     1/3)) - ((1 + I Sqrt[3]) (2 b^3 + 27 c +
       3 Sqrt[3] Sqrt[4 b^3 c + 27 c^2])^(1/3))/(6 2^(1/3))}}

For the values of b and c encounterd in the algorithm, the first solution is in fact the
only real one that we are after. So we can use this solution instead of doing a binary
search for the root of the cubic.

Signed-off-by: Robert C. Helling <helling@atdotde.de>
Signed-off-by: Jan Darowski <jan.darowski@gmail.com>

1 files changed, 18 insertions, 26 deletions

diff --git a/deco.c b/deco.c
index b8aae27a0..70430bb66 100644
--- a/deco.c
+++ b/deco.c
@@ -209,38 +209,30 @@ double he_factor(int period_in_seconds, int ci)
 // Calculates the nucleons inner pressure during the impermeable period
 double calc_inner_pressure(double crit_radius, double onset_tension, double current_ambient_pressure)
 {
-	double onset_radius;
-	double current_radius;
-	double A, B, C, low_bound, high_bound, root;
-	double valH, valL;
-	int ci;
-	const int max_iters = 10;
+	double onset_radius = 1.0 / (vpmb_config.gradient_of_imperm / (2.0 * (vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma)) + 1.0 / crit_radius);
 
-	// const, depends only on config.
-	onset_radius = 1.0 / (vpmb_config.gradient_of_imperm / (2.0 * (vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma)) + 1.0 / crit_radius);
+	// A*r^3 + B*r^2 + C == 0
+	// Solved with the help of mathematica
 
-	// A*r^3 + B*r^2 + C = 0
-	A = current_ambient_pressure - vpmb_config.gradient_of_imperm + (2.0 * (vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma)) / onset_radius;
-	B = 2.0 * (vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma);
-	C = onset_tension * pow(onset_radius, 3);
+	double A = current_ambient_pressure - vpmb_config.gradient_of_imperm + (2.0 * (vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma)) / onset_radius;
+	double B = 2.0 * (vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma);
+	double C = onset_tension * pow(onset_radius, 3);
 
-	// According to the algorithm's authors...
-	low_bound = B / A;
-	high_bound = onset_radius;
+	double BA = B/A;
+	double CA = C/A;
 
-	valH = high_bound * high_bound * (A * high_bound - B) - C;
-	valL = low_bound * low_bound * (A * low_bound - B) - C;
+	double discriminant = CA * (4 * BA * BA * BA + 27 * CA);
 
-	// Binary search for equations root.
-	for (ci = 0; ci < max_iters; ++ci) {
-		current_radius = (high_bound + low_bound) *0.5;
-		root = (current_radius * current_radius * (A * current_radius - B)) - C;
-		if (root >= 0.0)
-			high_bound = current_radius;
-		else
-			low_bound = current_radius;
+	// Let's make sure we have a real solution:
+	if (discriminant < 0.0) {
+		// This should better not happen
+		report_error("Complex solution for inner pressure encountered!\n A=%f\tB=%f\tC=%f\n", A, B, C);
+		return 0.0;
 	}
-	return onset_tension * (pow(onset_radius, 3) / pow(current_radius, 3));
+	double denominator = pow(BA * BA * BA + 1.5 * (9 * CA + sqrt(3.0) * sqrt(discriminant)), 1/3.0);
+	double current_radius = (BA + BA * BA / denominator + denominator) / 3.0;
+
+	return onset_tension * onset_radius * onset_radius * onset_radius / (current_radius * current_radius * current_radius);
 }
 
 // Calculates the crushing pressure in the given moment. Updates crushing_onset_tension and critical radius if needed