From f3efb298461abe2287491465fe22572d4d131322 Mon Sep 17 00:00:00 2001
From: Don Cross <cosinekitty@gmail.com>
Date: Wed, 26 Jun 2019 17:06:14 -0400
Subject: [PATCH] Python: Implemented quadratic interpolation, to be used by
 Search.

---
 generate/template/astronomy.py | 40 ++++++++++++++++++++++++++++++++++
 source/python/astronomy.py     | 40 ++++++++++++++++++++++++++++++++++
 2 files changed, 80 insertions(+)

diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py
index 1e45306d..de59bdf7 100644
--- a/generate/template/astronomy.py
+++ b/generate/template/astronomy.py
@@ -969,6 +969,46 @@ def _CalcChebyshev(model, time):
 
 # END CHEBYSHEV
 #----------------------------------------------------------------------------
+# BEGIN Search
+
+def _QuadInterp(tm, dt, fa, fm, fb):
+    Q = (fb + fa)/2 - fm
+    R = (fb - fa)/2
+    S = fm
+
+    if Q == 0:
+        # This is a line, not a parabola.
+        if R == 0:
+            # This is a HORIZONTAL line... can't make progress!
+            return None
+        x = -S / R
+        if not (-1 <= x <= +1):
+            return None  # out of bounds
+    else:
+        # It really is a parabola. Find roots x1, x2.
+        u = R*R - 4*Q*S
+        if u <= 0:
+            return None
+        ru = math.sqrt(u)
+        x1 = (-R + ru) / (2 * Q)
+        x2 = (-R - ru) / (2 * Q)
+
+        if -1 <= x1 <= +1:
+            if -1 <= x2 <= +1:
+                # Two solutions... so parabola intersects twice.
+                return None
+            x = x1
+        elif -1 <= x2 <= +1:
+            x = x2
+        else:
+            return None
+
+    t = tm + x*dt
+    df_dt = (2*Q*x + R) / dt
+    return (x, t, df_dt)
+
+# END Search
+#----------------------------------------------------------------------------
 
 def HelioVector(body, time):
     if body == BODY_PLUTO:
diff --git a/source/python/astronomy.py b/source/python/astronomy.py
index 352f6c69..c13e830c 100644
--- a/source/python/astronomy.py
+++ b/source/python/astronomy.py
@@ -1807,6 +1807,46 @@ def _CalcChebyshev(model, time):
 
 # END CHEBYSHEV
 #----------------------------------------------------------------------------
+# BEGIN Search
+
+def _QuadInterp(tm, dt, fa, fm, fb):
+    Q = (fb + fa)/2 - fm
+    R = (fb - fa)/2
+    S = fm
+
+    if Q == 0:
+        # This is a line, not a parabola.
+        if R == 0:
+            # This is a HORIZONTAL line... can't make progress!
+            return None
+        x = -S / R
+        if not (-1 <= x <= +1):
+            return None  # out of bounds
+    else:
+        # It really is a parabola. Find roots x1, x2.
+        u = R*R - 4*Q*S
+        if u <= 0:
+            return None
+        ru = math.sqrt(u)
+        x1 = (-R + ru) / (2 * Q)
+        x2 = (-R - ru) / (2 * Q)
+
+        if -1 <= x1 <= +1:
+            if -1 <= x2 <= +1:
+                # Two solutions... so parabola intersects twice.
+                return None
+            x = x1
+        elif -1 <= x2 <= +1:
+            x = x2
+        else:
+            return None
+
+    t = tm + x*dt
+    df_dt = (2*Q*x + R) / dt
+    return (x, t, df_dt)
+
+# END Search
+#----------------------------------------------------------------------------
 
 def HelioVector(body, time):
     if body == BODY_PLUTO: