Hi,
here is another attempt. This at least computes the correct roots and avoids the division by zero thanks to some algebraic simplifications. |
From dda7655dc6c184990bcad172309695dea0ce3acb Mon Sep 17 00:00:00 2001 From: "Robert C. Helling" <[email protected]> Date: Sun, 23 Aug 2015 10:48:38 +0200 Subject: [PATCH] Simplified Boyle's law correction and new cubic root
This uses a bit of algebra to simplify the formula for Boyle's law correction and introduces another algebraic cubic root solver that can also deal with negative discriminants thanks to a trigonometric formula from Wikipedia. Signed-off-by: Robert C. Helling <[email protected]> --- deco.c | 31 +++++++++++++++++++++++++------ 1 file changed, 25 insertions(+), 6 deletions(-) diff --git a/deco.c b/deco.c index d802ca5..198ffe8 100644 --- a/deco.c +++ b/deco.c @@ -346,16 +346,35 @@ double solve_cubic(double A, double B, double C) } +// Solve another cubic equation, this time +// x^3 - B x - C == 0 +// Use trigonometric formula for negative discriminants (see Wikipedia for details) + +double solve_cubic2(double B, double C) +{ + double discriminant = 27 * C * C - 4 * cube(B); + if (discriminant < 0.0) { + return 2.0 * sqrt(B / 3.0) * cos(acos(3.0 * C * sqrt(3.0 / B) / (2.0 * B)) / 3.0); + } + + double denominator = pow(9 * C + sqrt(3 * discriminant), 1 / 3.0); + + return pow(2.0 / 3.0, 1.0 / 3.0) * B / denominator + denominator / pow(18.0, 1.0 / 3.0); +} + +// This is a simplified formula avoiding radii. It uses the fact that Boyle's law says +// pV = (G + P_amb) / G^3 is constant to solve for the new gradient G. + double update_gradient(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.mbar / 1000.0 + 2.0 * vpmb_config.surface_tension_gamma / first_radius) * cube(first_radius); + double B = cube(first_gradient) / (first_stop_pressure.mbar / 1000.0 + first_gradient); + double C = next_stop_pressure * B; - double next_radius = solve_cubic(A, B, C); + double new_gradient = solve_cubic2(B, C); - return 2.0 * vpmb_config.surface_tension_gamma / next_radius; + if (new_gradient < 0.0) + report_error("Negative gradient encountered!"); + return new_gradient; } void boyles_law(double next_stop_pressure) -- 1.9.5 (Apple Git-50.3)
Patch 2 I sent remains valid. This definitely deals with the symptoms. But there might still be a problem with the profile ceiling calculations as when calculating the 60min@100m dive there is a ceiling violation shown. This needs some further care. Will do this later. Best Robert |
signature.asc
Description: Message signed with OpenPGP using GPGMail
_______________________________________________ subsurface mailing list [email protected] http://lists.subsurface-divelog.org/cgi-bin/mailman/listinfo/subsurface
