From 2b438504e0fc68ea8224e88675247e555ee6a6e6 Mon Sep 17 00:00:00 2001
From: Tavian Barnes <tavianator@gmail.com>
Date: Wed, 6 Oct 2010 23:56:11 -0400
Subject: Handle edge case in AdaptiveIntegrator::step().

---
 src/vZ/Adaptive.hpp | 16 ++++++++++------
 1 file changed, 10 insertions(+), 6 deletions(-)

diff --git a/src/vZ/Adaptive.hpp b/src/vZ/Adaptive.hpp
index 256bda6..05fc439 100644
--- a/src/vZ/Adaptive.hpp
+++ b/src/vZ/Adaptive.hpp
@@ -76,30 +76,34 @@ namespace vZ
   GenericAdaptiveIntegrator<Y>::step()
   {
     static const Scalar S = Scalar(19)/Scalar(20); // Arbitrary saftey factor
-    Scalar newH;
+    Scalar newH = this->h();
     Y y;
 
     // Attempt the integration step in a loop
-    bool rejected = true;
-    while (rejected) {
+    while (true) {
       KVector k = calculateK(m_a);
       y         = calculateY(k, m_b);
       Y yStar   = calculateY(k, m_bStar);
 
       // Get an error estimate
+
       using std::abs;
       using std::pow;
       Scalar delta = abs(y - yStar);
-      Scalar scale = m_atol + std::max(abs(y), abs(this->y()))*m_rtol;
 
-      newH = S*this->h()*pow(scale/delta, Scalar(1)/m_order);
+      if (delta == Scalar(0)) {
+        break;
+      }
+
+      Scalar scale = m_atol + std::max(abs(y), abs(this->y()))*m_rtol;
+      newH         = S*this->h()*pow(scale/delta, Scalar(1)/m_order);
 
       if (delta > scale) {
         // Reject the step
         this->h(newH);
         ++m_rejections;
       } else {
-        rejected = false;
+        break;
       }
     }
 
-- 
cgit v1.2.3