From 89e66616be3108437756fc7160c151e38d3e7c9a Mon Sep 17 00:00:00 2001
From: ericway1024 <ericway1024@gmail.com>
Date: Thu, 7 Jul 2022 20:57:11 +0100
Subject: [PATCH 1/3] Add gradient to linear interpolant

---
 nevis/_interpolation.py | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)

diff --git a/nevis/_interpolation.py b/nevis/_interpolation.py
index 7311341..5782a9b 100644
--- a/nevis/_interpolation.py
+++ b/nevis/_interpolation.py
@@ -67,6 +67,27 @@ def __call__(self, x, y):
         # Final result
         return np.where(f1 == f2, f1, (y2 - y) * f1 + (y - y1) * f2)
 
+    def grad(self, x, y):
+        """
+        Returns the gradient of the linear interpolant at the given point.
+        """
+        ny, nx = self._heights.shape
+        x, y = x / self._resolution - 0.5, y / self._resolution - 0.5
+
+        # Find nearest grid points
+        # x1 (left), x2 (right), y1 (bottom), y2 (top).
+        x1 = np.minimum(nx - 2, np.maximum(0, int(x)))
+        y1 = np.minimum(ny - 2, np.maximum(0, int(y)))
+        x2 = x1 + 1
+        y2 = y1 + 1
+
+        # Heights at nearest grid points (subscripts are x_y)
+        h11 = self._heights[y1, x1]
+        h12 = self._heights[y2, x1]
+        h21 = self._heights[y1, x2]
+
+        return (h21 - h11) / self._resolution, (h12 - h11) / self._resolution
+
 
 def spline(verbose=False):
     """

From 53638a0085ba76c97407248d4f976763de2b068c Mon Sep 17 00:00:00 2001
From: ericway1024 <ericway1024@gmail.com>
Date: Fri, 8 Jul 2022 10:51:04 +0100
Subject: [PATCH 2/3] Make gradient part of linear interpolation

---
 nevis/_interpolation.py | 49 +++++++++++++++++++----------------------
 1 file changed, 23 insertions(+), 26 deletions(-)

diff --git a/nevis/_interpolation.py b/nevis/_interpolation.py
index 5782a9b..66a060f 100644
--- a/nevis/_interpolation.py
+++ b/nevis/_interpolation.py
@@ -13,10 +13,17 @@
 
 class linear_interpolant(object):
     """
-    Returns a linear interpolation over the full GB data set.
+    Returns a linear interpolation and optionally its gradient over the full
+    GB data set.
 
-    The returned function takes two arguments ``x`` and ``y`` (both in metres)
-    and returns an interpolated height ``z`` (in meters).
+    When ``grad`` is set to ``False`` (by default), the returned function
+    takes two arguments ``x`` and ``y`` (both in metres) and returns an
+    interpolated height ``z`` (in meters).
+
+    When ``grad`` is set to ``True``, the returned function takes two
+    arguments ``x`` and ``y`` (both in metres) and returns a tuple ``(z, g)``,
+    where ``z`` is an interpolated height (in meters) and ``g`` is a tuple
+    ``(dz/dx, dz/dy)``.
 
     The height for each grid point ``(i, j)`` is assumed to be in the center of
     the square from ``(i, j)`` to ``(i + 1, j + 1)``.
@@ -26,13 +33,18 @@ class linear_interpolant(object):
         f = linear_interpolation()
         print(f(1000, 500))
 
+        f_grad = linear_interpolation(grad=True)
+        z, (gx, gy) = f_grad(1000, 500)
+        print(f"Height: {z:.2f} m, gradient: ({gx:.2f} m/m, {gy:.2f} m/m)")
+
     """
     # Note: This is technically a class, but used as a function here so
     # following the underscore naming convention.
 
-    def __init__(self):
+    def __init__(self, grad=False):
         self._heights = nevis.gb()
         self._resolution = nevis.spacing()
+        self.grad = grad
 
     def __call__(self, x, y):
         ny, nx = self._heights.shape
@@ -65,28 +77,13 @@ def __call__(self, x, y):
         f2 = np.where(h12 == h22, h12, (x2 - x) * h12 + (x - x1) * h22)
 
         # Final result
-        return np.where(f1 == f2, f1, (y2 - y) * f1 + (y - y1) * f2)
-
-    def grad(self, x, y):
-        """
-        Returns the gradient of the linear interpolant at the given point.
-        """
-        ny, nx = self._heights.shape
-        x, y = x / self._resolution - 0.5, y / self._resolution - 0.5
-
-        # Find nearest grid points
-        # x1 (left), x2 (right), y1 (bottom), y2 (top).
-        x1 = np.minimum(nx - 2, np.maximum(0, int(x)))
-        y1 = np.minimum(ny - 2, np.maximum(0, int(y)))
-        x2 = x1 + 1
-        y2 = y1 + 1
-
-        # Heights at nearest grid points (subscripts are x_y)
-        h11 = self._heights[y1, x1]
-        h12 = self._heights[y2, x1]
-        h21 = self._heights[y1, x2]
-
-        return (h21 - h11) / self._resolution, (h12 - h11) / self._resolution
+        f = np.where(f1 == f2, f1, (y2 - y) * f1 + (y - y1) * f2)
+        if self.grad:
+            # Gradient of the interpolant
+            g = (h21 - h11) / self._resolution, (h12 - h11) / self._resolution
+            return f, g
+        else:
+            return f
 
 
 def spline(verbose=False):

From 26ca5ff47c7ed8246969e530c24774055445a74e Mon Sep 17 00:00:00 2001
From: ericway1024 <ericway1024@gmail.com>
Date: Fri, 8 Jul 2022 11:54:38 +0100
Subject: [PATCH 3/3] Make interpolation a float

---
 nevis/_interpolation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/nevis/_interpolation.py b/nevis/_interpolation.py
index 66a060f..6d6cfcd 100644
--- a/nevis/_interpolation.py
+++ b/nevis/_interpolation.py
@@ -77,7 +77,7 @@ def __call__(self, x, y):
         f2 = np.where(h12 == h22, h12, (x2 - x) * h12 + (x - x1) * h22)
 
         # Final result
-        f = np.where(f1 == f2, f1, (y2 - y) * f1 + (y - y1) * f2)
+        f = float(np.where(f1 == f2, f1, (y2 - y) * f1 + (y - y1) * f2))
         if self.grad:
             # Gradient of the interpolant
             g = (h21 - h11) / self._resolution, (h12 - h11) / self._resolution