# HG changeset patch
# User Christof Kaufmann <christofkaufmann.dev@gmail.com>
# Date 1643025149 -3600
# Node ID 79d6280fb00a3a49c1149472f9191f16b25521af
# Parent  c3494f4088959fdb193d8a6baebdb4e3118db9c2
interp2: Fix issues with complex input (bug #61903).

* scripts/polynomial/pchip.m: Fix interpolation of imaginary part of complex
numbers.
* scripts/general/interp1.m: Fix value for extrapolation of complex values.
Add BISTs.

diff -r c3494f408895 -r 79d6280fb00a scripts/general/interp1.m
--- a/scripts/general/interp1.m	Thu Jan 27 18:06:01 2022 +0100
+++ b/scripts/general/interp1.m	Mon Jan 24 12:52:29 2022 +0100
@@ -123,7 +123,7 @@
   endif
 
   method = "linear";
-  extrap = NA;
+  extrap = [];
   xi = [];
   ispp = false;
   have_xi = false;
@@ -167,6 +167,14 @@
     endif
   endif
 
+  if (isempty (extrap))
+    if (iscomplex (y))
+      extrap = NA + 1i*NA;
+    else
+      extrap = NA;
+    endif
+  endif
+
   ## reshape matrices for convenience
   x = x(:);
   nx = rows (x);
@@ -503,6 +511,9 @@
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> ...
+%!  (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!   interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -520,6 +531,9 @@
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> ...
+%!  (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!   interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -540,6 +554,9 @@
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> ...
+%!  (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!   interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 ## This test is expected to fail, so commented out.
 ## "previous" and "next" options are not symmetric w.r.t to flipping xp,yp
 #%!assert (interp1 (xp,yp,xi,style),...
@@ -559,6 +576,9 @@
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> ...
+%!  (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!   interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 #%!assert (interp1 (xp,yp,xi,style),...
 #%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -579,6 +599,9 @@
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> ...
+%!  (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!   interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 #%!assert (interp1 (xp,yp,xi,style),...
 #%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -596,6 +619,9 @@
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> ...
+%!  (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!   interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 #%!assert (interp1 (xp,yp,xi,style),...
 #%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -616,6 +642,9 @@
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> ...
+%!  (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!   interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -633,6 +662,9 @@
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> ...
+%!  (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!   interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -655,6 +687,9 @@
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> ...
+%!  (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!   interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -672,6 +707,9 @@
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> ...
+%!  (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!   interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -692,6 +730,9 @@
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> ...
+%!  (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!   interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -709,6 +750,9 @@
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> ...
+%!  (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!   interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -729,6 +773,9 @@
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> ...
+%!  (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!   interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -746,6 +793,9 @@
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> ...
+%!  (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!   interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
diff -r c3494f408895 -r 79d6280fb00a scripts/polynomial/pchip.m
--- a/scripts/polynomial/pchip.m	Thu Jan 27 18:06:01 2022 +0100
+++ b/scripts/polynomial/pchip.m	Mon Jan 24 12:52:29 2022 +0100
@@ -111,6 +111,11 @@
 
   ## Compute derivatives.
   d = __pchip_deriv__ (x, y, 2);
+  if (iscomplex (y))
+    ## __pchip_deriv__ ignores imaginary part.  Do it again for imag part.
+    ## FIXME: Adapt __pchip_deriv__ to correctly handle complex input.
+    d += 1i * __pchip_deriv__ (x, imag (y), 2);
+  endif
   d1 = d(:, 1:n-1);
   d2 = d(:, 2:n);