# HG changeset patch
# User Christof Kaufmann <christofkaufmann.dev@gmail.com>
# Date 1642847104 -3600
# Node ID 27566803b3b1a9ef6e0b15684241fa2c49e20591
# Parent  ffed0ae35edaca6e022576e8d381878e05b29746
interp2: Fix issues with complex input (bug #61863).

* scripts/general/interp2.m: Fix interpolation of imaginary part of complex
numbers for "pchip". Fix value for extrapolation of complex values. Add BIST.
* doc/interpreter/contributors.in: Add name to list of contributors.

diff -r ffed0ae35eda -r 27566803b3b1 doc/interpreter/contributors.in
--- a/doc/interpreter/contributors.in	Fri Jan 21 15:23:02 2022 +0100
+++ b/doc/interpreter/contributors.in	Sat Jan 22 11:25:04 2022 +0100
@@ -183,6 +183,7 @@
 Lute Kamstra
 Fotios Kasolis
 Thomas Kasper
+Christof Kaufmann
 Joel Keay
 Mumit Khan
 Paul Kienzle
diff -r ffed0ae35eda -r 27566803b3b1 scripts/general/interp2.m
--- a/scripts/general/interp2.m	Fri Jan 21 15:23:02 2022 +0100
+++ b/scripts/general/interp2.m	Sat Jan 22 11:25:04 2022 +0100
@@ -254,9 +254,21 @@
       ## first order derivatives
       DX = __pchip_deriv__ (X, Z, 2);
       DY = __pchip_deriv__ (Y, Z, 1);
-      ## Compute mixed derivatives row-wise and column-wise, use the average.
+      ## Compute mixed derivatives row-wise and column-wise.  Use the average.
       DXY = (__pchip_deriv__ (X, DY, 2) + __pchip_deriv__ (Y, DX, 1))/2;
 
+      if (iscomplex (Z))
+        ## __pchip_deriv__ works only on real part.  Do it again for imag part.
+        ## FIXME: Adapt __pchip_deriv__ to correctly handle complex input.
+
+        ## first order derivatives
+        DX += 1i * __pchip_deriv__ (X, imag (Z), 2);
+        DY += 1i * __pchip_deriv__ (Y, imag (Z), 1);
+        ## Compute mixed derivatives row-wise and column-wise.  Use the average.
+        DXY += 1i * (__pchip_deriv__ (X, imag (DY), 2)
+                     + __pchip_deriv__ (Y, imag (DX), 1))/2;
+      endif
+
       ## do the bicubic interpolation
       hx = diff (X); hx = hx(xidx);
       hy = diff (Y); hy = hy(yidx);
@@ -327,7 +339,11 @@
 
   ## extrapolation 'extrap'
   if (isempty (extrap))
-    extrap = NA;
+    if (iscomplex (Z))
+      extrap = complex (NA, NA);
+    else
+      extrap = NA;
+    endif
   endif
 
   if (X(1) < X(end))
@@ -480,6 +496,20 @@
 %!
 %! assert (result, expected, 1000*eps);
 
+## Test that interpolating a complex matrix is equivalent to interpolating its
+## real and imaginary parts separately.
+%!test <61863>
+%! xi = [2.5, 3.5];
+%! yi = [0.5, 1.5]';
+%! orig = rand (2, 3) + 1i * rand (2, 3);
+%! for method = {"nearest", "linear", "pchip", "spline"}
+%!   interp_complex = interp2 (orig, xi, yi, method{1});
+%!   interp_real = interp2 (real (orig), xi, yi, method{1});
+%!   interp_imag = interp2 (imag (orig), xi, yi, method{1});
+%!   assert (real (interp_complex), interp_real)
+%!   assert (imag (interp_complex), interp_imag)
+%! endfor
+
 %!test  # 2^n refinement form
 %! x = [1,2,3];
 %! y = [4,5,6,7];