From 3ec85791874d5f874cd531c5eca928ac0d85b9b3 Mon Sep 17 00:00:00 2001 From: graham Date: Wed, 11 Nov 2009 13:23:26 +0000 Subject: [PATCH] Stabilising division by p. --- .../derived/pitchForkRing/pitchForkRing.C | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.C index 8a3a081f6f..b94a4f3a10 100644 --- a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.C +++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.C @@ -71,9 +71,10 @@ scalar pitchForkRing::energy(const vector r) const { scalar p = sqrt(r.x()*r.x() + r.y()*r.y()); - scalar pMinusRSqr = (p - rOrbit_)*(p - rOrbit_); + scalar pMinusRSqr = sqr(p - rOrbit_); - return -0.5 * mu_ * pMinusRSqr + return + -0.5 * mu_ * pMinusRSqr + 0.25 * pMinusRSqr * pMinusRSqr + 0.5 * alpha_ * r.z() * r.z(); } @@ -87,9 +88,9 @@ vector pitchForkRing::force(const vector r) const return vector ( - (mu_ - pMinusR * pMinusR) * pMinusR * r.x()/p, - (mu_ - pMinusR * pMinusR) * pMinusR * r.y()/p, - -alpha_ * r.z() + (mu_ - sqr(pMinusR)) * pMinusR * r.x()/(p + VSMALL), + (mu_ - sqr(pMinusR)) * pMinusR * r.y()/(p + VSMALL), + - alpha_ * r.z() ); }