From 01f26d607143fa1b81cefadeb932d8457f7bab8a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 2 May 2022 22:09:06 -0500 Subject: [PATCH] Fix bug with torus distance when particle is coincident --- src/surface.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/surface.cpp b/src/surface.cpp index dbfa0bf3d..c0b3f0a59 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1026,9 +1026,11 @@ double torus_distance(double x1, double x2, double x3, double u1, double u2, double c1p = 2 * four_A2 * (u1 * x1 + u2 * x2); double c0p = four_A2 * (x1 * x1 + x2 * x2); - // Coefficient for equation: a t^4 + b t^3 + c t^2 + d t + e = 0 + // Coefficient for equation: a t^4 + b t^3 + c t^2 + d t + e = 0. If the point + // is coincident, the 'e' coefficient should be zero. Explicitly setting it to + // zero helps avoid numerical issues below with root finding. double coeff[5]; - coeff[0] = c0 * c0 - c0p; + coeff[0] = coincident ? 0.0 : c0 * c0 - c0p; coeff[1] = 2 * c0 * c1 - c1p; coeff[2] = c1 * c1 + 2 * c0 * c2 - c2p; coeff[3] = 2 * c1 * c2;