From 8909dbbfd38c3bbc89e96ae6b904111946864f51 Mon Sep 17 00:00:00 2001 From: "Robert C. Helling" Date: Sun, 23 Aug 2015 10:48:38 +0200 Subject: 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 Signed-off-by: Dirk Hohndel --- deco.c | 31 +++++++++++++++++++++++++------ 1 file changed, 25 insertions(+), 6 deletions(-) (limited to 'deco.c') diff --git a/deco.c b/deco.c index d802ca539..198ffe85b 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) -- cgit v1.2.3-70-g09d2