double root = sqrtpos( r1*r1 - 4. * r2 * r0 );
if (abs(r2) < TINY || root < TINY) return Vec4(0., 0., 0., -1.);
xInvHad = 0.5 * (root / abs(r2) - r1 / r2);
- xDirHad = (cM0 - cM3 * xInvHad) / (cM2 + cM4 * xInvHad);
-
+ if (abs(cM2 + cM4 * xInvHad) > 1.e-10) {
+ xDirHad = (cM0 - cM3 * xInvHad) / (cM2 + cM4 * xInvHad);
+ } else {
+ xDirHad = 1.e10 * (cM0 - cM3 * xInvHad);
+ }
// Define position of new trial vertex.
xDirNew = (iDirNew == iDirOld) ? xDirOld - xDirHad : 1. - xDirHad;
xInvNew = (iInvNew == iInvOld) ? xInvOld + xInvHad : xInvHad;