From efda5dd343cff2f8935409d96e1be1ccd5189445 Mon Sep 17 00:00:00 2001 From: "Robert C. Helling" Date: Wed, 1 Jul 2015 12:27:42 +0200 Subject: [PATCH] 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 --- deco.c | 45 ++++++++++++++++++++------------------------- 1 file changed, 20 insertions(+), 25 deletions(-) diff --git a/deco.c b/deco.c index 6dce93d..107ad05 100644 --- a/deco.c +++ b/deco.c @@ -281,37 +281,32 @@ void nuclear_regeneration(double time) // 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; - // 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); + 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); - // 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); + // A*r^3 + B*r^2 + C == 0 + // Solved with the help of mathematica - // According to the algorithm's authors... - low_bound = B / A; - high_bound = onset_radius; + 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); - valH = high_bound * high_bound * (A * high_bound - B) - C; - valL = low_bound * low_bound * (A * low_bound - B) - C; + double BA = B/A; + double CA = C/A; - 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; + double discriminant = CA * (4 * BA * BA * BA + 27 * CA); + + // 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 -- 2.3.2 (Apple Git-55)