* Bug 16629 fixed: interp1() fixed + complex + extended extrapolations 77/21677/6
Samuel GOUGEON [Mon, 18 Jan 2021 06:25:45 +0000 (07:25 +0100)]
  http://bugzilla.scilab.org/16629
  interp1 page fixed & updated (PDF): http://bugzilla.scilab.org/attachment.cgi?id=5217

Change-Id: Iec59c6c2ffa312d6d282c52828b40966bc8871d9

13 files changed:
scilab/CHANGES.md
scilab/modules/helptools/etc/images_md5.txt
scilab/modules/helptools/images/cell_view.png [new file with mode: 0644]
scilab/modules/helptools/images/interp1_1.png
scilab/modules/helptools/images/library_view.png [new file with mode: 0644]
scilab/modules/helptools/images/list_view.png [new file with mode: 0644]
scilab/modules/helptools/images/structure_view.png [new file with mode: 0644]
scilab/modules/interpolation/help/en_US/interp1.xml
scilab/modules/interpolation/help/ja_JP/interp1.xml [deleted file]
scilab/modules/interpolation/help/pt_BR/interp1.xml [deleted file]
scilab/modules/interpolation/help/ru_RU/interp1.xml [new file with mode: 0644]
scilab/modules/interpolation/macros/interp1.sci
scilab/modules/interpolation/tests/unit_tests/interp1.tst [new file with mode: 0644]

index 3c27968..136da4a 100644 (file)
@@ -200,6 +200,9 @@ Feature changes and additions
 * `unique` is enabled for any 2D sparse arrays, in simple, 'c' and 'r' modes.
 * `%chars` constant added, to easily access to some selected sets of unicode symbols.
 * Lists are displayed in a more compact and comprehensive way.
+* `interp1` is upgraded:
+  - complex numbers `y` now supported: the real and imaginary parts are interpolated separately.
+  - extrapolation option extended: `edgevalue` mode added for all interpolations; `periodic` mode added for linear and spline interpolations. `linear` mode added for the spline interpolations.
 
 Help pages:
 -----------
@@ -354,10 +357,14 @@ Bug Fixes
 * [#16553](https://bugzilla.scilab.org/16553): `unique(["" ""])` returned `["" ""]`.
 * [#16557](https://bugzilla.scilab.org/16557): `macr2tree` + `tree2code` translated `e={2}` into `"e=1"` and `e={2,"ab"}` into `"e=[2,"ab"]"`.
 * [#16559](https://bugzilla.scilab.org/16553): `isempty(A)` was true for sparse matrix of dimension 2^16 or larger.
+<<<<<<< HEAD
 * [#16622](https://bugzilla.scilab.org/16622): `inv` could no longer be overloaded for hypermatrices of decimal or complex numbers.
 * [#16623](https://bugzilla.scilab.org/16623): `rand(2,2,2)^2` yielded a wrong result instead of trying to call the `%s_p_s` overload for input hypermatrices.
 * [#16644](https://bugzilla.scilab.org/16644): `input("message:")` yielded a wrong error message about `mprintf` in case of non-interpretable input.
 * [#16654](https://bugzilla.scilab.org/16654): `interp` was leaking memory.
+=======
+* [#16629](https://bugzilla.scilab.org/16629): `interp1`'s documentation did not tell the spline edges conditions ; extrapolation modes were poorly explained ; the description of the result's size was completely wrong ; x as an option was not documented. A wrong extrapolation value could silently return a wrong result. There was some dead code like `if varargin(5)==%nan`. A bugged error message yielded its own error. When x is implicit, the argument index in error messages could be wrong. `periodic` and `edgevalue` extrapolation modes were not available. `linear` extrapolation was not available for splines. When `xp` is an hypermatrix with `size(xp,1)==1`, the size of the result was irregular/wrong.
+>>>>>>> ae66d0ee461... * Bug 16629 fixed: interp1() fixed + complex + extended extrapolations
 
 
 ### Bugs fixed in 6.1.0:
index 459a3ae..674e87c 100644 (file)
@@ -802,7 +802,7 @@ hsv2rgb_1.png=7d20c259e94301d9763fbddb7bff4784
 hsvcolormap_1.png=11918d88bcc793200af0b9f1b58b0554
 iir_1.png=e675c2755f68ddc4611c849895b63012
 intdec_1.png=da3896e7d2e8468dd33edf57ccbe4480
-interp1_1.png=0e9a3f4319b2818ce4921b9bc3008d80
+interp1_1.png=e76a8e0a6070fbe078060844707c49fe
 interp2d_1.png=f4af61bc3faf493d778169ec7decc7ae
 interp2d_2.png=a62363849a3a9a3571dc1942fcebfbd3
 interp2d_3.png=a9f2a0b0e59ac8cedae3ad22e0b2cc46
diff --git a/scilab/modules/helptools/images/cell_view.png b/scilab/modules/helptools/images/cell_view.png
new file mode 100644 (file)
index 0000000..2905e6c
Binary files /dev/null and b/scilab/modules/helptools/images/cell_view.png differ
index 7002a47..2b0659b 100644 (file)
Binary files a/scilab/modules/helptools/images/interp1_1.png and b/scilab/modules/helptools/images/interp1_1.png differ
diff --git a/scilab/modules/helptools/images/library_view.png b/scilab/modules/helptools/images/library_view.png
new file mode 100644 (file)
index 0000000..d3df3ad
Binary files /dev/null and b/scilab/modules/helptools/images/library_view.png differ
diff --git a/scilab/modules/helptools/images/list_view.png b/scilab/modules/helptools/images/list_view.png
new file mode 100644 (file)
index 0000000..0a9d969
Binary files /dev/null and b/scilab/modules/helptools/images/list_view.png differ
diff --git a/scilab/modules/helptools/images/structure_view.png b/scilab/modules/helptools/images/structure_view.png
new file mode 100644 (file)
index 0000000..106a0a9
Binary files /dev/null and b/scilab/modules/helptools/images/structure_view.png differ
index 1408549..2b6b370 100644 (file)
@@ -2,8 +2,8 @@
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2008 - INRIA - Farid BELAHCENE
- *
  * Copyright (C) 2012 - 2016 - Scilab Enterprises
+ * Copyright (C) 2021 - Samuel GOUGEON - Le Mans Université
  *
  * This file is hereby licensed under the terms of the GNU GPL v2.0,
  * pursuant to article 5.3.4 of the CeCILL v.2.1.
           xmlns:scilab="http://www.scilab.org" xml:id="interp1" xml:lang="en">
     <refnamediv>
         <refname>interp1</refname>
-        <refpurpose>one_dimension interpolation function</refpurpose>
+        <refpurpose>1D interpolation in nearest, linear or spline mode</refpurpose>
     </refnamediv>
     <refsynopsisdiv>
         <title>Syntax</title>
         <synopsis>
+          yp = interp1(y, xp)
           yp = interp1(x, y, xp)
-          yp = interp1(x, y, xp, method)
-          yp = interp1(x, y, xp, method, extrapolation)
+          yp = interp1(.., xp, method)
+          yp = interp1(.., xp, method, extrapolation)
         </synopsis>
     </refsynopsisdiv>
     <refsection>
         <title>Arguments</title>
         <variablelist>
             <varlistentry>
-                <term>xp</term>
+                <term>x</term>
                 <listitem>
-                    <para>reals scalar, vector or matrix (or hypermatrix)</para>
+                    vector of at least 2 real numbers: Abscissas of known interpolation nodes,
+                    without duplicates. By default,
+                    <itemizedlist>
+                        <listitem>
+                            if <varname>y</varname> is a vector: <literal>x=1:length(y)</literal>.
+                        </listitem>
+                        <listitem>
+                            if <varname>y</varname> is a matrix or an hypermatrix:
+                            <literal>x=1:size(y,1)</literal>.
+                        </listitem>
+                    </itemizedlist>
+                    <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>x</term>
+                <term>y</term>
                 <listitem>
-                    <para>reals vector</para>
+                    vector, matrix or hypermatrix of real or complex numbers: values at known
+                    interpolation nodes, at the corresponding <varname>x</varname> abscissas.
+                    <itemizedlist>
+                        <listitem>
+                            if <varname>y</varname> is a vector, <varname>x</varname> and
+                            <varname>y</varname> must have the same length.
+                        </listitem>
+                        <listitem>
+                            if <varname>y</varname> is a matrix or an hypermatrix, we must have
+                            <literal>length(x)==size(y,1)</literal>. Each column of
+                            <varname>y</varname> is then interpolated versus the same
+                            <varname>x</varname> abscissas, for the given <varname>xp</varname>.
+                        </listitem>
+                    </itemizedlist>
+                    <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>method</term>
+                <term>xp</term>
                 <listitem>
-                    <para>(optional) string defining the interpolation method</para>
+                    scalar, vector, matrix or hypermatrix or decimal numbers: abscissas of
+                    points whose values <varname>yp</varname> must be computed according
+                    to data of interpolation nodes.
+                    <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>extrapolation</term>
+                <term>yp</term>
                 <listitem>
-                    <para>(optional) string, or real value defining the yp(j) components
-                        for xp(j) values outside [x1,xn] interval.
-                    </para>
+                    vector, matrix, or hypermatrix of numbers: interpolated <varname>y</varname>
+                    values at the given <varname>xp</varname>.
+                    <itemizedlist>
+                        <listitem>
+                            if <varname>y</varname> is a vector: <varname>yp</varname> has the
+                            size of <varname>xp</varname>.
+                            <para/>
+                        </listitem>
+                        <listitem>
+                            if <varname>y</varname> is a matrix or an hypermatrix:
+                            <itemizedlist>
+                                <listitem>
+                                    if <varname>xp</varname> is a scalar or a vector:
+                                    <literal>size(yp)</literal> is <literal>[length(xp) size(y)(2:$)]</literal>
+                                </listitem>
+                                <listitem>
+                                    if <varname>xp</varname> is a matrix or an hypermatrix:
+                                    <literal>size(yp)</literal> is <literal>[size(xp) size(y)(2:$)]</literal>
+                                </listitem>
+                            </itemizedlist>
+                        </listitem>
+                    </itemizedlist>
+                    <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>yp</term>
+                <term>method</term>
                 <listitem>
-                    <para>vector, or matrix (or hypermatrix)</para>
+                    string defining the interpolation method. Possible values and processing are:
+                    <para/>
+                    <table>
+                        <tr>
+                            <td><emphasis role="bold">"linear"</emphasis>:</td>
+                            <td>
+                                linear interpolation between consecutive nodes, used by default.
+                            </td>
+                        </tr>
+                        <tr>
+                            <td><emphasis role="bold">"spline"</emphasis>:</td>
+                            <td>
+                                interpolation by cubic splines
+                            </td>
+                        </tr>
+                        <tr>
+                            <td><emphasis role="bold">"nearest"</emphasis>:</td>
+                            <td>
+                                <para>for each value xp(j), yp(j) takes the value or y(i)
+                                    corresponding to x(i) the nearest neighbor of xp(j)
+                                </para>
+                            </td>
+                        </tr>
+                    </table>
+                    <para/>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>extrapolation</term>
+                <listitem>
+                    string or number defining the yp(j) components for xp(j) values outside the
+                    [x(1)=min(x),x($)=max(x)] interval. We suppose here-below that <varname>x</varname>
+                    and <varname>y</varname> have already been sorted accordingly.
+                    <para/>
+                    <table>
+                        <tr>
+                            <td><emphasis role="bold">"extrap"</emphasis>:</td>
+                            <td>
+                                <literal>interp1(x,y,xp, method, "extrap")</literal> is equivalent
+                                to <literal>interp1(x,y,xp, method, method)</literal>.
+                            </td>
+                        </tr>
+                        <tr>
+                            <td><emphasis role="bold">"linear"</emphasis>:</td>
+                            <td>
+                                Can be used with the <literal>"spline"</literal> (and obviously
+                                <literal>"linear"</literal>) interpolation methods.
+                            </td>
+                        </tr>
+                        <tr>
+                            <td><emphasis role="bold">"periodic"</emphasis>:</td>
+                            <td>
+                                This extrapolation type can be used with the
+                                <literal>"linear"</literal> or <literal>"spline"</literal>
+                                interpolation methods. Then: if <varname>y</varname> is a
+                                vector, <literal>y(1)==y($)</literal> is required ; otherwise
+                                <literal>y(1,:)==y($,:)</literal> is required.
+                            </td>
+                        </tr>
+                        <tr>
+                            <td><emphasis role="bold">"edgevalue"</emphasis>:</td>
+                            <td>
+                                Then <literal>yp(i)=y(1)</literal> for every
+                                <literal>xp(i)&lt;x(1)</literal>, and
+                                <literal>yp(i)=y($)</literal> for every <literal>xp(i)>x($)</literal>.
+                            </td>
+                        </tr>
+                        <tr>
+                            <td><emphasis role="bold">padding</emphasis>:</td>
+                            <td>
+                                <varname>padding</varname> is a decimal or complex number
+                                used to set <literal>yp(i)=padding</literal> for every
+                                <literal>xp(i) ∉ [min(x),max(x)]</literal>. Example:
+                                <literal>yi=interp1(x,y,xp,method, 0)</literal>.
+                            </td>
+                        </tr>
+                        <tr>
+                            <td><emphasis role="bold">(none)</emphasis>:</td>
+                            <td>
+                                By default, the extrapolation is performed by splines when
+                                splines are used for the interpolation, and by padding
+                                with %nan when the interpolation is linear or by "nearest" node.
+                            </td>
+                        </tr>
+                    </table>
                 </listitem>
             </varlistentry>
         </variablelist>
     <refsection>
         <title>Description</title>
         <para>
-            Given <literal>(x,y,xp)</literal>, this function performs the yp
-            components corresponding to xp by the interpolation(linear by default)
-            defined by x and y.
+            Given <literal>(x,y,xp)</literal>, this function computes the yp
+            components corresponding to xp by the interpolation between known
+            data provided by (x,y) nodes.
         </para>
         <para>
-            If <literal>yp</literal> is a vector then the length of xp is equal
-            to the length of <literal>yp,</literal> if <literal>yp</literal> is a
-            matrix then <literal>xp</literal> have the same length than the length of
-            each columns of yp, if <literal>yp</literal> is a hypermatrix then the
-            length of <literal>xp</literal> have the same length than the first
-            dimension of <literal>yp. </literal>
-        </para>
-        <para>If size(y)=[C,M1,M2,M3,...,Mj] and size(xp)=[N1,N2,N3,...,Nk] then
-            size(yp)=[N1,N2,..,Nk,M1,M2,...Mj] and length of x must be equal to
-            size(y,1)
+            <varname>x</varname> is priorly sorted in ascending order, and
+            <varname>y</varname> values or per column are then sorted accordingly.
         </para>
         <para>
-            The <literal>method</literal> parameter sets the evaluation rule for
-            interpolation
+            <emphasis role="bold">Interpolation of complex values</emphasis>: When
+            <varname>y</varname> is complex, its real and imaginary parts are interpolated
+            separately, and then added to build the complex <varname>yp</varname>.
         </para>
-        <variablelist>
-            <varlistentry>
-                <term>"linear"</term>
-                <listitem>
-                    <para>
-                        the interpolation is defined by linear method (see <link linkend="interpln">interpln</link>)
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>"spline"</term>
-                <listitem>
-                    <para>the interpolation is defined by cubic spline interpolation (
-                        see <link linkend="splin">splin</link> , <link linkend="interp">interp</link>)
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>"nearest"</term>
-                <listitem>
-                    <para>for each value xp(j), yp(j) takes the value or y(i)
-                        corresponding to x(i) the nearest neighbor of xp(j)
-                    </para>
-                </listitem>
-            </varlistentry>
-        </variablelist>
         <para>
-            The <literal>extrapolation</literal> parameter sets the evaluation
-            rule for extrapolation, i.e for <literal>xp(i) </literal>not in [x1,xn]
-            interval
+            <emphasis role="bold">interp1(x,y,xp,"nearest")</emphasis>: For any xp
+            at the middle of an <literal>[x(i),x(i+1)]</literal> interval, the upper
+            bound <literal>x(i+1)</literal> is considered as the nearest
+            <varname>x</varname> value, and <literal>yp=y(i+1)</literal> is assigned.
         </para>
-        <variablelist>
-            <varlistentry>
-                <term>"extrap"</term>
-                <listitem>
-                    <para>the extrapolation is performed by the defined method.
-                        yp=interp1(x,y,xp,method,"extrap")
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>real value</term>
-                <listitem>
-                    <para>you can choose a real value for extrapolation, in this way
-                        yp(i) takes this value for xp(i) not in [x1,xn] interval, for
-                        example 0 (but also nan or inf). yi=interp1(x,y,xp,method, 0)
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>by default</term>
-                <listitem>
-                    <para>the extrapolation is performed by the defined method (for
-                        spline method), and by nan for linear and nearest methods.
-                        yp=interp1(x,y,xp,method)
-                    </para>
-                </listitem>
-            </varlistentry>
-        </variablelist>
+        <refsect3>
+            <title>linear interpolations</title>
+            They are performed through the <literal>linear_interpn(..)</literal> function,
+            with the corresponding <literal>"edgevalue"→"C0"</literal>,
+            <literal>"linear"→"natural"</literal>, <literal>"periodic"→"periodic"</literal>
+            extrapolation option.
+        </refsect3>
+        <refsect3>
+            <title>spline interpolations</title>
+            <para>
+                <emphasis role="bold">interp1(..,xp,"spline")</emphasis> or
+                <emphasis role="bold">interp1(..,xp,"spline","spline")</emphasis> or
+                <emphasis role="bold">interp1(..,xp,"spline","extrap")</emphasis>
+                use <literal>not_a_knot</literal> edges conditions. Extrapolation is performed
+                by using both spline polynomials computed at the <literal>(x,y)</literal> edges.
+            </para>
+            <para>
+                <emphasis role="bold">interp1(..,xp,"spline","edgevalue")</emphasis>
+                uses <literal>not_a_knot</literal> edges conditions and then calls
+                <literal>interp(..,"C0")</literal> to perform the actual interpolation
+                and extrapolation.
+            </para>
+            <para>
+                <emphasis role="bold">interp1(..,xp,"spline","periodic")</emphasis>
+                calls both <literal>splin(..)</literal> and then <literal>interp(..)</literal>
+                with their <literal>"periodic"</literal> option.
+            </para>
+            <para>
+                <emphasis role="bold">interp1(..,xp,"spline","linear")</emphasis>
+                calls <literal>splin(..,"natural")</literal> for linear edges conditions,
+                and then feeds <literal>interp(..,"linear")</literal>.
+            </para>
+        </refsect3>
     </refsection>
     <refsection>
         <title>Examples</title>
         <programlisting role="example"><![CDATA[
-x = linspace(0, 3, 20);
-y = x.^2;
-xx = linspace(0, 3, 100);
-yy1 = interp1(x, y, xx, 'linear');
-yy2 = interp1(x, y, xx, 'spline');
-yy3 = interp1(x, y, xx, 'nearest');
-plot(xx', [yy1; yy2; yy3]', x, y, '*')
-xtitle('interpolation of square function')
-legend(['linear','spline','nearest'], "in_upper_left");
+x = linspace(0, 10, 11)';
+y = sin(x);
+xx = linspace(0,10,1000)';
+yy2 = interp1(x, y, xx, 'linear');
+yy1 = interp1(x, y, xx, 'nearest');
+yy3 = interp1(x, y, xx, 'spline');
+clf
+h = plot(xx, [yy1 yy2 yy3], x, y, '.')
+h(1).mark_size = 8;
+title "Interpolation of a poorly sampled sin() function" fontsize 3
+legend(['nearest','linear','spline','nodes'], "in_lower_left");
  ]]></programlisting>
         <scilab:image>
-            x=linspace(0,3,20);
-            y=x.^2;
-            xx=linspace(0,3,100);
-            yy1=interp1(x,y,xx,'linear');
-            yy2=interp1(x,y,xx,'spline');
-            yy3=interp1(x,y,xx,'nearest');
-            plot(xx',[yy1;yy2;yy3]',x,y,'*')
-            xtitle('interpolation of square function')
-            legend(['linear','spline','nearest'],a=2)
+            x = linspace(0, 10, 11)';
+            y = sin(x);
+            xx = linspace(0,10,1000)';
+            yy2 = interp1(x, y, xx, 'linear');
+            yy1 = interp1(x, y, xx, 'nearest');
+            yy3 = interp1(x, y, xx, 'spline');
+            clf
+            h = plot(xx, [yy1 yy2 yy3], x, y, '.')
+            h(1).mark_size = 8;
+            title "Interpolation of a poorly sampled sin() function" fontsize 3
+            legend(['nearest','linear','spline','nodes'], "in_lower_left");
+            gcf().axes_size = [600 400];
         </scilab:image>
     </refsection>
     <refsection role="see also">
@@ -181,10 +292,10 @@ legend(['linear','spline','nearest'], "in_upper_left");
                 <link linkend="interp">interp</link>
             </member>
             <member>
-                <link linkend="interpln">interpln</link>
+                <link linkend="splin">splin</link>
             </member>
             <member>
-                <link linkend="splin">splin</link>
+                <link linkend="linear_interpn">linear_interpn</link>
             </member>
         </simplelist>
     </refsection>
@@ -192,8 +303,34 @@ legend(['linear','spline','nearest'], "in_upper_left");
         <title>History</title>
         <revhistory>
             <revision>
-                <revnumber>5.4.0</revnumber>
-                <revremark>previously, imaginary part of input arguments were implicitly ignored.</revremark>
+                <revnumber>6.1.1</revnumber>
+                <revremark>
+                    <itemizedlist>
+                        <listitem>
+                            For complex <varname>y</varname> values, <literal>imag(y)</literal> is
+                            no longer ignored: <literal>real(y)</literal> and <literal>imag(y)</literal>
+                            parts are separately interpolated.
+                        </listitem>
+                        <listitem>
+                            <literal>"periodic"</literal> extrapolation added for the linear and
+                            spline interpolations.
+                        </listitem>
+                        <listitem>
+                            <literal>"edgevalue"</literal> extrapolation added for all nearest,
+                            linear and spline interpolations.
+                        </listitem>
+                        <listitem>
+                            <literal>"linear"</literal> extrapolation added for the spline
+                            interpolation.
+                        </listitem>
+                        <listitem>
+                            When <varname>xp</varname> is an hypermatrix and
+                            <literal>size(xp,1)==1</literal>, <literal>size(yp)</literal>
+                            is now always <literal>[size(xp) size(y)(2,$)</literal> instead of
+                            <literal>[size(xp,2:$), size(y)(2,$)</literal>.
+                        </listitem>
+                    </itemizedlist>
+                </revremark>
             </revision>
         </revhistory>
     </refsection>
diff --git a/scilab/modules/interpolation/help/ja_JP/interp1.xml b/scilab/modules/interpolation/help/ja_JP/interp1.xml
deleted file mode 100644 (file)
index 547ad0c..0000000
+++ /dev/null
@@ -1,200 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<!--
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2008 - INRIA - Farid BELAHCENE
- *
- * Copyright (C) 2012 - 2016 - Scilab Enterprises
- *
- * This file is hereby licensed under the terms of the GNU GPL v2.0,
- * pursuant to article 5.3.4 of the CeCILL v.2.1.
- * This file was originally licensed under the terms of the CeCILL v2.1,
- * and continues to be available under such terms.
- * For more information, see the COPYING file which you should have received
- * along with this program.
- *
- -->
-<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink"
-          xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns5="http://www.w3.org/1999/xhtml"
-          xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook"
-          xmlns:scilab="http://www.scilab.org" xml:id="interp1" xml:lang="ja">
-    <refnamediv>
-        <refname>interp1</refname>
-        <refpurpose>一次元補間関数</refpurpose>
-    </refnamediv>
-    <refsynopsisdiv>
-        <title>呼び出し手順</title>
-        <synopsis>
-          yp = interp1(x, y, xp)
-          yp = interp1(x, y, xp, method)
-          yp = interp1(x, y, xp, method, extrapolation)
-        </synopsis>
-    </refsynopsisdiv>
-    <refsection>
-        <title>パラメータ</title>
-        <variablelist>
-            <varlistentry>
-                <term>xp</term>
-                <listitem>
-                    <para>実数スカラー, ベクトルまたは行列 (またはハイパー行列)</para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>x</term>
-                <listitem>
-                    <para>実数ベクトル</para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>method</term>
-                <listitem>
-                    <para>(オプションの) 補間方法を定義する文字列</para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>extrapolation</term>
-                <listitem>
-                    <para>(オプションの) 文字列, または領域 [x1,xn] の外側の xp(j)の値について
-                        yp(j) の要素を定義する実数値.
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>yp</term>
-                <listitem>
-                    <para>ベクトルまたは行列 (またはハイパー行列)</para>
-                </listitem>
-            </varlistentry>
-        </variablelist>
-    </refsection>
-    <refsection>
-        <title>説明</title>
-        <para>
-            <literal>(x,y,xp)</literal>を指定すると, この関数は
-            xp に対応する yp の要素を x および y により定義された
-            (デフォルトでは線形)補間により計算します.
-        </para>
-        <para>
-            <literal>yp</literal> がベクトルの場合,
-            xp の長さは <literal>yp</literal> の長さに等しくなります.
-            <literal>yp</literal> が行列の場合,
-            <literal>xp</literal> は yp の各列の長さと同じ長さとなります.
-            <literal>yp</literal> がハイパー行列の場合,
-            <literal>xp</literal> の長さは<literal>yp</literal> の最初の次元と同じ
-            になります.
-        </para>
-        <para>size(y)=[C,M1,M2,M3,...,Mj] かつ size(xp)=[N1,N2,N3,...,Nk] の場合,
-            size(yp)=[N1,N2,..,Nk,M1,M2,...Mj] となり, x の長さは size(y,1) に等しくなります.
-        </para>
-        <para>
-            <literal>method</literal> パラメータは補間の計算手法を設定します.
-        </para>
-        <variablelist>
-            <varlistentry>
-                <term>"linear"</term>
-                <listitem>
-                    <para>
-                        補間は線形手法により定義されます ( <link linkend="interpln">interpln</link>参照)
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>"spline"</term>
-                <listitem>
-                    <para>補間は3次スプライン補間により定義されます (
-                        <link linkend="splin">splin</link> , <link linkend="interp">interp</link>参照)
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>"nearest"</term>
-                <listitem>
-                    <para>各 xp(j), yp(j) の値は,
-                        xp(j)の最近傍にある x(i) に対応する y(i) の値を選択します.
-                    </para>
-                </listitem>
-            </varlistentry>
-        </variablelist>
-        <para>
-            <literal>extrapolation</literal> パラメータは捕外用,すなわち
-            <literal>xp(i) </literal> が領域[x1,xn]内にない時
-            の計算手法を設定します.
-        </para>
-        <variablelist>
-            <varlistentry>
-                <term>"extrap"</term>
-                <listitem>
-                    <para>捕外が定義された手法により実行されます.
-                        yp=interp1(x,y,xp,method,"extrap")
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>実数値</term>
-                <listitem>
-                    <para>捕外用の実数値を選択できます.
-                        この場合, 区間 [x1,xn] の中にない xp(i) について yp(i) は
-                        この値とします.
-                        yi=interp1(x,y,xp,method, 0)
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>デフォルトの動作</term>
-                <listitem>
-                    <para>捕外は(スプライン手法の場合)定義された手法により行われます
-                        線形または最近傍手法の場合は nan となります.
-                        yp=interp1(x,y,xp,method)
-                    </para>
-                </listitem>
-            </varlistentry>
-        </variablelist>
-    </refsection>
-    <refsection>
-        <title>例</title>
-        <programlisting role="example"><![CDATA[
-x = linspace(0, 3, 20);
-y = x^2;
-xx = linspace(0, 3, 100);
-yy1 = interp1(x, y, xx, 'linear');
-yy2 = interp1(x, y, xx, 'spline');
-yy3 = interp1(x, y, xx, 'nearest');
-plot(xx, [yy1; yy2; yy3], x, y, '*')
-xtitle('interpolation of square function')
-legend(['linear','spline','nearest'], "in_upper_left");
- ]]></programlisting>
-        <scilab:image>
-            x=linspace(0,3,20);
-            y=x.^2;
-            xx=linspace(0,3,100);
-            yy1=interp1(x,y,xx,'linear');
-            yy2=interp1(x,y,xx,'spline');
-            yy3=interp1(x,y,xx,'nearest');
-            plot(xx',[yy1;yy2;yy3]',x,y,'*')
-            xtitle('interpolation of square function')
-            legend(['linear','spline','nearest'],a=2)
-        </scilab:image>
-    </refsection>
-    <refsection role="see also">
-        <title>参照</title>
-        <simplelist type="inline">
-            <member>
-                <link linkend="interp">interp</link>
-            </member>
-            <member>
-                <link linkend="interpln">interpln</link>
-            </member>
-            <member>
-                <link linkend="splin">splin</link>
-            </member>
-        </simplelist>
-    </refsection>
-    <refsection>
-        <title>履歴</title>
-        <revhistory>
-            <revision>
-                <revnumber>5.4.0</revnumber>
-                <revremark>以前では, 入力引数の虚部は暗黙のうちに無視されていました.</revremark>
-            </revision>
-        </revhistory>
-    </refsection>
-</refentry>
diff --git a/scilab/modules/interpolation/help/pt_BR/interp1.xml b/scilab/modules/interpolation/help/pt_BR/interp1.xml
deleted file mode 100644 (file)
index 2bae0d1..0000000
+++ /dev/null
@@ -1,182 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<!--
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2008 - INRIA - Farid BELAHCENE
- *
- * Copyright (C) 2012 - 2016 - Scilab Enterprises
- *
- * This file is hereby licensed under the terms of the GNU GPL v2.0,
- * pursuant to article 5.3.4 of the CeCILL v.2.1.
- * This file was originally licensed under the terms of the CeCILL v2.1,
- * and continues to be available under such terms.
- * For more information, see the COPYING file which you should have received
- * along with this program.
- *
- -->
-<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink"
-          xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns5="http://www.w3.org/1999/xhtml"
-          xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook"
-          xmlns:scilab="http://www.scilab.org" xml:id="interp1" xml:lang="pt">
-    <refnamediv>
-        <refname>interp1</refname>
-        <refpurpose>função de interpolação 1d</refpurpose>
-    </refnamediv>
-    <refsynopsisdiv>
-        <title>Seqüência de Chamamento</title>
-        <synopsis>
-          yp = interp1(x, y, xp)
-          yp = interp1(x, y, xp, method)
-          yp = interp1(x, y, xp, method, extrapolation)
-        </synopsis>
-    </refsynopsisdiv>
-    <refsection>
-        <title>Parâmetros</title>
-        <variablelist>
-            <varlistentry>
-                <term>xp</term>
-                <listitem>
-                    <para>escalar real, vetor ou matriz (ou hipermatriz) de reais
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>x</term>
-                <listitem>
-                    <para>vetor de reais</para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>method</term>
-                <listitem>
-                    <para>(opcional) string definindo o método de interpolação </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>extrapolation</term>
-                <listitem>
-                    <para>(opcional) string ou valor real defindo os componentes yp(j)
-                        para os valores xp(j) fora do intervalo [x1,xn].
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>yp</term>
-                <listitem>
-                    <para>vetor ou matriz (ou hipermatriz) </para>
-                </listitem>
-            </varlistentry>
-        </variablelist>
-    </refsection>
-    <refsection>
-        <title>Descrição</title>
-        <para>
-            Dados <literal>(x,y,xp)</literal>, esta função faz corresponder os
-            yp componentes aos xp por interpolação (linear por padrão) definida por x
-            e y.
-        </para>
-        <para>
-            Se <literal>yp</literal> é um vetor, então o comprimento de xp é
-            igual ao comprimento de <literal>yp</literal>. Se <literal>yp</literal> é
-            uma matriz, então <literal>xp</literal> tem o mesmo comprimento que cada
-            uma das colunas de yp. Se <literal>yp</literal> é uma hipermatriz, então o
-            comprimento de <literal>xp</literal> é o mesmo da primeira dimensão de
-            <literal>yp. </literal>
-        </para>
-        <para>Se size(y)=[C,M1,M2,M3,...,Mj] e size(xp)= [N1,N2,N3,...,Nk] então
-            size(yp)= [N1,N2,..,Nk,M1,M2,...Mj] e o comprimento de x deve ser igual a
-            size(y,1)
-        </para>
-        <para>
-            O parâmetro <literal>method</literal> ajusta a regra avaliação para
-            interpolação.
-        </para>
-        <variablelist>
-            <varlistentry>
-                <term>"linear"</term>
-                <listitem>
-                    <para>
-                        a interpolação é definida pelo método linear (ver <link linkend="interpln">interpln</link>)
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>"spline"</term>
-                <listitem>
-                    <para>
-                        definição de interpolação por spline cúbico (ver <link linkend="splin">splin</link> , <link linkend="interp">interp</link>)
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>"nearest"</term>
-                <listitem>
-                    <para>para cada valor xp(j), yp(j) toma o valor ou y(i)
-                        correspondente ao x(i), o vízinho mais próximo de xp(j)
-                    </para>
-                    <para/>
-                </listitem>
-            </varlistentry>
-        </variablelist>
-        <para>
-            O parâmetro <literal>extrapolation</literal> ajusta a regra de
-            avaliação para extrapolação, i.e para <literal>xp(i) </literal>fora do
-            intervalo [x1,xn]
-        </para>
-        <variablelist>
-            <varlistentry>
-                <term>"extrap"</term>
-                <listitem>
-                    <para>a extrapolação é realizada pelo método definido.
-                        yp=interp1(x,y,xp,method,"extrap")
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>valor real</term>
-                <listitem>
-                    <para>você pode escolher um valor real para extrapolação. Deste
-                        modo, yp(i) toma este valor para xp(i) fora do intervalo [x1,xn],
-                        por exemplo 0 (mas também nan ou inf). yi=interp1(x,y,xp,method, 0)
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>por padrão</term>
-                <listitem>
-                    <para>a extrapolação é realizada pelo método definido (para método
-                        spline), e por nan para os métodos "linear" e "nearest".
-                        yp=interp1(x,y,xp,method)
-                    </para>
-                </listitem>
-            </varlistentry>
-        </variablelist>
-    </refsection>
-    <refsection>
-        <title>Exemplos</title>
-        <programlisting role="example"><![CDATA[
-x = linspace(0, 3, 20);
-y = x^2;
-xx = linspace(0, 3, 100);
-yy1 = interp1(x, y, xx, 'linear');
-yy2 = interp1(x, y, xx, 'spline');
-yy3 = interp1(x, y, xx, 'nearest');
-plot(xx, [yy1; yy2; yy3], x, y, '*')
-xtitle('interpolation of square function')
-legend(['linear','spline','nearest'], "in_upper_left");
- ]]></programlisting>
-    </refsection>
-    <refsection role="see also">
-        <title>Ver Também</title>
-        <simplelist type="inline">
-            <member>
-                <link linkend="interp">interp</link>
-            </member>
-            <member>
-                <link linkend="interpln">interpln</link>
-            </member>
-            <member>
-                <link linkend="splin">splin</link>
-            </member>
-        </simplelist>
-    </refsection>
-</refentry>
diff --git a/scilab/modules/interpolation/help/ru_RU/interp1.xml b/scilab/modules/interpolation/help/ru_RU/interp1.xml
new file mode 100644 (file)
index 0000000..d0c573f
--- /dev/null
@@ -0,0 +1,356 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!--
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) 2008 - INRIA - Farid BELAHCENE
+ * Copyright (C) 2012 - 2016 - Scilab Enterprises
+ * Copyright (C) 2021 - Samuel GOUGEON - Le Mans Université
+ *
+ * This file is hereby licensed under the terms of the GNU GPL v2.0,
+ * pursuant to article 5.3.4 of the CeCILL v.2.1.
+ * This file was originally licensed under the terms of the CeCILL v2.1,
+ * and continues to be available under such terms.
+ * For more information, see the COPYING file which you should have received
+ * along with this program.
+ *
+ -->
+<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink"
+          xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns5="http://www.w3.org/1999/xhtml"
+          xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook"
+          xmlns:scilab="http://www.scilab.org" xml:id="interp1" xml:lang="ru">
+    <refnamediv>
+        <refname>interp1</refname>
+        <refpurpose>одномерная интерполяция в режиме nearest, linear или spline</refpurpose>
+    </refnamediv>
+    <refsynopsisdiv>
+        <title>Синтаксис</title>
+        <synopsis>
+          yp = interp1(y, xp)
+          yp = interp1(x, y, xp)
+          yp = interp1(.., xp, method)
+          yp = interp1(.., xp, method, extrapolation)
+        </synopsis>
+    </refsynopsisdiv>
+    <refsection>
+        <title>Аргументы</title>
+        <variablelist>
+            <varlistentry>
+                <term>x</term>
+                <listitem>
+                    вектор по меньшей мере из двух вещественных чисел: ось абсцисс известных узлов
+                    интерполяции, без повторов. По умолчанию
+                    <itemizedlist>
+                        <listitem>
+                            если <varname>y</varname> является вектором: <literal>x=1:length(y)</literal>.
+                        </listitem>
+                        <listitem>
+                            если <varname>y</varname> является матрицей или гиперматрицей:
+                            <literal>x=1:size(y,1)</literal>.
+                        </listitem>
+                    </itemizedlist>
+                    <para/>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>y</term>
+                <listitem>
+                    вектор, матрица или гиперматрица вещественных или комплексных
+                    чисел: значения в узлах интерполяции в соответствующих значениях
+                    оси абсцисс <varname>x</varname>.
+                    <itemizedlist>
+                        <listitem>
+                            если <varname>y</varname> является вектором, то <varname>x</varname>
+                            и <varname>y</varname> должны быть одной длины.
+                        </listitem>
+                        <listitem>
+                            если <varname>y</varname> является матрицей или гиперматрицей, то мы
+                            должны иметь <literal>length(x)==size(y,1)</literal>. Каждый столбец
+                            <varname>y</varname> тогда интерполирован в зависимости от тех же
+                            значений оси абсцисс <varname>x</varname>, для указанного
+                            <varname>xp</varname>.
+                        </listitem>
+                    </itemizedlist>
+                    <para/>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>xp</term>
+                <listitem>
+                    скаляр, вектор, матрица или гиперматрица десятичных чисел:
+                    ось абсцисс точек, чьи значения <varname>yp</varname> должны быть
+                    вычислены в соответствии с данными узлов интерполяции.
+                    <para/>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>yp</term>
+                <listitem>
+                    вектор, матрица или гиперматрица чисел: интерполированные значения
+                    <varname>y</varname> в указанном <varname>xp</varname>.
+                    <itemizedlist>
+                        <listitem>
+                            если <varname>y</varname> является вектором:
+                            <varname>yp</varname> имеет размер <varname>xp</varname>.
+                            <para/>
+                        </listitem>
+                        <listitem>
+                            если <varname>y</varname> является матрицей или гиперматрицей:
+                            <itemizedlist>
+                                <listitem>
+                                    если <varname>xp</varname> является скаляром или вектором:
+                                    <literal>size(yp)</literal> равен
+                                    <literal>[length(xp) size(y)(2:$)]</literal>
+                                </listitem>
+                                <listitem>
+                                    если <varname>xp</varname> является матрицей или
+                                    гиперматрицей: <literal>size(yp)</literal> равен
+                                    <literal>[size(xp) size(y)(2:$)]</literal>
+                                </listitem>
+                            </itemizedlist>
+                        </listitem>
+                    </itemizedlist>
+                    <para/>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>method</term>
+                <listitem>
+                    строковое значение, определяющее метод интерполяции. Возможными
+                    значениями и обработкой являются:
+                    <para/>
+                    <table>
+                        <tr>
+                            <td><emphasis role="bold">"linear"</emphasis>:</td>
+                            <td>
+                                линейная интерполяция между последовательными узлами, используется
+                                по умолчанию.
+                            </td>
+                        </tr>
+                        <tr>
+                            <td><emphasis role="bold">"spline"</emphasis>:</td>
+                            <td>
+                                интерполяция кубическими сплайнами
+                            </td>
+                        </tr>
+                        <tr>
+                            <td><emphasis role="bold">"nearest"</emphasis>:</td>
+                            <td>
+                                <para>
+                                    для каждого значения <literal>xp(j)</literal> <literal>yp(j)</literal>
+                                    принимает значение <literal>y(i)</literal>, соответствующее
+                                    <literal>x(i)</literal>, ближайшего соседа
+                                    <literal>xp(j)</literal>
+                                </para>
+                            </td>
+                        </tr>
+                    </table>
+                    <para/>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>extrapolation</term>
+                <listitem>
+                    строковое значение или число, определяющее элементы <literal>yp(j)</literal>
+                    для значений <literal>xp(j)</literal> за пределами интервала
+                    <literal>[x(1)=min(x),x($)=max(x)]</literal>. Мы полагаем здесь и далее,
+                    что <varname>x</varname> и <varname>y</varname> уже соответственно
+                    отсортированы.
+                    <para/>
+                    <table>
+                        <tr>
+                            <td><emphasis role="bold">"extrap"</emphasis>:</td>
+                            <td>
+                                <literal>interp1(x,y,xp, method, "extrap")</literal> эквивалентен
+                                <literal>interp1(x,y,xp, method, method)</literal>.
+                            </td>
+                        </tr>
+                        <tr>
+                            <td><emphasis role="bold">"linear"</emphasis>:</td>
+                            <td>
+                                Может использоваться с методами интерполяции
+                                <literal>"spline"</literal> (и, очевидно,
+                                <literal>"linear"</literal>).
+                            </td>
+                        </tr>
+                        <tr>
+                            <td><emphasis role="bold">"periodic"</emphasis>:</td>
+                            <td>
+                                Этот тип экстраполяции может использоваться с методами
+                                интерполяции <literal>"linear"</literal> или
+                                <literal>"spline"</literal>. Тогда: если <varname>y</varname>
+                                является вектором, то требуется, чтобы
+                                <literal>y(1)==y($)</literal>; в противном случае требуется,
+                                чтобы <literal>y(1,:)==y($,:)</literal>.
+                            </td>
+                        </tr>
+                        <tr>
+                            <td><emphasis role="bold">"edgevalue"</emphasis>:</td>
+                            <td>
+                                Тогда <literal>yp(i)=y(1)</literal> для каждого
+                                <literal>xp(i)&lt;x(1)</literal>, и
+                                <literal>yp(i)=y($)</literal> для каждого
+                                <literal>xp(i)>x($)</literal>.
+                            </td>
+                        </tr>
+                        <tr>
+                            <td><emphasis role="bold">padding</emphasis>:</td>
+                            <td>
+                                <varname>padding</varname> является десятичным или
+                                комплексным числом, используемым для установки
+                                <literal>yp(i)=padding</literal> для каждого
+                                <literal>xp(i) ∉ [min(x),max(x)]</literal>. Например:
+                                <literal>yi=interp1(x,y,xp,method, 0)</literal>.
+                            </td>
+                        </tr>
+                        <tr>
+                            <td><emphasis role="bold">(none)</emphasis>:</td>
+                            <td>
+                                По умолчанию экстраполяция выполняется сплайнами, когда
+                                сплайны используются для интерполяции и дополнением
+                                значениями <literal>%nan</literal>, если интерполяция
+                                линейна или по "ближайшему" узлу.
+                            </td>
+                        </tr>
+                    </table>
+                </listitem>
+            </varlistentry>
+        </variablelist>
+    </refsection>
+    <refsection>
+        <title>Описание</title>
+        <para>
+            Указывая <literal>(x,y,xp)</literal>, данная функция вычисляет элементы
+            <varname>yp</varname>, соответствующие <varname>xp</varname> с помощью интерполяции
+            между известными данными, указанными в узлах <literal>(x,y)</literal>.
+        </para>
+        <para>
+            <varname>x</varname> предварительно отсортирован в порядке возрастания, а значения
+            <varname>y</varname> либо по столбцам тогда сортируются соответственно.
+        </para>
+        <para>
+            <emphasis role="bold">Интерполяция комплексных значений</emphasis>:
+            Если <varname>y</varname> является комплексным, то его вещественная и мнимая
+            части интерполируются отдельно, и потом суммируются для формирования
+            комплексного <varname>yp</varname>.
+        </para>
+        <para>
+            <emphasis role="bold">interp1(x,y,xp,"nearest")</emphasis>:
+            Для любого <varname>xp</varname> в середине интервала
+            <literal>[x(i),x(i+1)]</literal> верхняя граница <literal>x(i+1)</literal>
+            рассматривается как ближайшее значение <varname>x</varname>, и присваивается
+            <literal>yp=y(i+1)</literal>.
+        </para>
+        <refsect3>
+            <title>линейные интерполяции</title>
+            Они выполняются через функцию <literal>linear_interpn(..)</literal>,
+            с соответствующей опцией интерполяции <literal>"edgevalue"→"C0"</literal>,
+            <literal>"linear"→"natural"</literal>, <literal>"periodic"→"periodic"</literal>.
+        </refsect3>
+        <refsect3>
+            <title>интерполяции сплайнами</title>
+            <para>
+                <emphasis role="bold">interp1(..,xp,"spline")</emphasis> или
+                <emphasis role="bold">interp1(..,xp,"spline","spline")</emphasis> или
+                <emphasis role="bold">interp1(..,xp,"spline","extrap")</emphasis>
+                используют условия граней <literal>не_узел</literal>.
+                Экстраполяция выполняется с помощью обоих сплайновых полиномов,
+                вычисленных на гранях <literal>(x,y)</literal>.
+            </para>
+            <para>
+                <emphasis role="bold">interp1(..,xp,"spline","edgevalue")</emphasis>
+                использует условия граней <literal>не_узел</literal>, а затем вызывает
+                <literal>interp(..,"C0")</literal>, чтобы выполнить фактическую
+                интерполяцию и экстраполяцию.
+            </para>
+            <para>
+                <emphasis role="bold">interp1(..,xp,"spline","periodic")</emphasis>
+                вызывает оба <literal>splin(..)</literal>, а затем <literal>interp(..)</literal>
+                с их опцией <literal>"periodic"</literal>.
+            </para>
+            <para>
+                <emphasis role="bold">interp1(..,xp,"spline","linear")</emphasis>
+                вызывает <literal>splin(..,"natural")</literal> для условий линейных граней, а затем
+                передаётся в <literal>interp(..,"linear")</literal>.
+            </para>
+        </refsect3>
+    </refsection>
+    <refsection>
+        <title>Примеры</title>
+        <programlisting role="example"><![CDATA[
+x = linspace(0, 10, 11)';
+y = sin(x);
+xx = linspace(0,10,1000)';
+yy2 = interp1(x, y, xx, 'linear');
+yy1 = interp1(x, y, xx, 'nearest');
+yy3 = interp1(x, y, xx, 'spline');
+clf
+h = plot(xx, [yy1 yy2 yy3], x, y, '.')
+h(1).mark_size = 8;
+title "Interpolation of a poorly sampled sin() function" fontsize 3
+legend(['nearest','linear','spline','nodes'], "in_lower_left");
+ ]]></programlisting>
+        <scilab:image>
+            x = linspace(0, 10, 11)';
+            y = sin(x);
+            xx = linspace(0,10,1000)';
+            yy2 = interp1(x, y, xx, 'linear');
+            yy1 = interp1(x, y, xx, 'nearest');
+            yy3 = interp1(x, y, xx, 'spline');
+            clf
+            h=plot(xx, [yy1 yy2 yy3],x, y, '.')
+            h(1).mark_size = 8;
+            title "Interpolation of a poorly sampled sin() function" fontsize 3
+            legend(['nearest','linear','spline','nodes'], "in_lower_left");
+            gcf().axes_size = [600 400];
+        </scilab:image>
+    </refsection>
+    <refsection role="see also">
+        <title>Смотрите также</title>
+        <simplelist type="inline">
+            <member>
+                <link linkend="interp">interp</link>
+            </member>
+            <member>
+                <link linkend="splin">splin</link>
+            </member>
+            <member>
+                <link linkend="linear_interpn">linear_interpn</link>
+            </member>
+        </simplelist>
+    </refsection>
+    <refsection>
+        <title>История</title>
+        <revhistory>
+            <revision>
+                <revnumber>6.1.1</revnumber>
+                <revremark>
+                    <itemizedlist>
+                        <listitem>
+                            Для комплексных значений <varname>y</varname>,
+                            <literal>imag(y)</literal> более не игнорируется:
+                            <literal>real(y)</literal> и <literal>imag(y)</literal>
+                            части интерполируются раздельно.
+                        </listitem>
+                        <listitem>
+                            Экстраполяция <literal>"periodic"</literal> добавлена для
+                            линейной и сплайновой интерполяций.
+                        </listitem>
+                        <listitem>
+                            Экстраполяция <literal>"edgevalue"</literal> добавлена для
+                            всех интерполяций ближайшего соседа, линейной и сплайновой.
+                        </listitem>
+                        <listitem>
+                            Экстраполяция <literal>"linear"</literal> добавлена для
+                            сплайновой интерполяции.
+                        </listitem>
+                        <listitem>
+                            Когда <varname>xp</varname> является гиперматрицей и
+                            <literal>size(xp,1)==1</literal>, <literal>size(yp)</literal>
+                            теперь всегда <literal>[size(xp) size(y)(2,$)</literal> вместо
+                            <literal>[size(xp,2:$), size(y)(2,$)</literal>.
+                        </listitem>
+                    </itemizedlist>
+                </revremark>
+            </revision>
+        </revhistory>
+    </refsection>
+</refentry>
index 56dac19..1709911 100644 (file)
@@ -1,7 +1,7 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) INRIA - Farid BELAHCENE
-//
 // Copyright (C) 2012 - 2016 - Scilab Enterprises
+// Copyright (C) 2021 - Samuel GOUGEON - Le Mans Université
 //
 // This file is hereby licensed under the terms of the GNU GPL v2.0,
 // pursuant to article 5.3.4 of the CeCILL v.2.1.
 // For more information, see the COPYING file which you should have received
 // along with this program.
 
-function yi=interp1(varargin)
-
-    // yi=interp1(x,y,xi[,method[,interpolation)
-    // This function performs the yi values corresponding to xi by interpolation defined by x and y.
-    // Inputs :
-    // x , y : given data, x is a reals vector, y is a vector, matrix, or hypermatrix of reals
-    // if y is a vector, the length of x must be equal to the length of y,
-    // else the size of the first dimension of y must be equal to length of x.
-    // xi : a vector, matrix, or hypermatrix of reals
-    // Output
-    // yi : reals vector, matrix or hypermatrix, the values corresponding to xi by interpolation defined by x and y
-    // if size(y)=[C,N1,N2,N3,....] and size(xi)=[M1,M2,M3,M4] then size(xi)=[M1,M2,M3,M4,N1,N2,N3,N4,..], and length of x must be equal to C.
-    // Several kind of intepolations may be computed by selecting the appropriate method  parameter:
-    // The methods are:
-    // linear : this is the default method (using the interp Scilab function)
-    // spline : this is the cubic spline interpolation  (using interpln and splin Scilab functions)
-    // nearest : yi take the values corresponding to the nearest neighbor of xi
+function yi = interp1(varargin)
+    // yi = interp1(x, y, xi)
+    // yi = interp1(x, y, xi, method)
+    // yi = interp1(x, y, xi, method, extrapolation)
+    // yi = interp1(y, xi)
+    // yi = interp1(y, xi, method)
+    // yi = interp1(y, xi, method, extrapolation)
     //
-    // Several kind of extrapolations may be computed :
-    // 'extrap' : the extrapolation points is performed by the defined method
-    //  real value : you can choose a real value for extrapolation, in this way yp(i) takes this value for xp(i) not in [x1,xn] interval, for example 0 (but also nan or inf).
-    // by default the extrapolation is performed by the defined method (for spline method), and by nan for linear and nearest method.
-    // F.B
+    // x: real vector
+    // y: real or complex vector of length(x) or array of size(y,1)==length(x)
+    // xi: real scalar, vector, matrix or hypermatrix of any size.
+    // method = "linear" (default) | "nearest" | "spline"
+    // extrapolation = "extrap" (= method) | real or complex padding scalar |
+    //  * if method = "linear": "linear" | "periodic" | "edgevalue"
+    //  * if method = "spline": "linear" | "periodic" | "edgevalue" | "spline"
+    //  * if method = "nearest": "nearest" = "edgevalue"
 
-    rhs=size(varargin)
+    // ========================
+    // CHECKING INPUT ARGUMENTS
+    // ========================
+    rhs = size(varargin)
     // 2 <= Number of input arguments <= 5
     if rhs<2 | rhs>5 then
-        error(msprintf(gettext("%s: Wrong number of input arguments: Must be between %d and %d.\n"),"interp1",2,5));
+        msg = _("%s: Wrong number of input arguments: Must be between %d and %d.\n")
+        error(msprintf(msg, "interp1", 2, 5));
     end
 
-    //if yi=interp1(x,y,xi,..) not change
-    //if yi=interp1(y,xi,...) replace input argument by yi=interp1(x,y,xi,..),whith x=1:size(y,1) by default
-    if rhs==2 | (rhs>2 & type(varargin(3))==10) then
-        if isvector(varargin(1)) then
-            X=1:size(varargin(1),"*")
-        elseif size(size(varargin(1)),"*")==2 then
-            if (size(varargin(1),1)>1 & size(varargin(1),2)>1) then
-                X=1:size(varargin(1),1)
-            end
+    // Checking x
+    // ----------
+    iy = 2  // index of y argin, for error messages
+    if rhs >= 3 & type(varargin(1))<>0 & type(varargin(3))<>10 then
+        // x is explicit: interp1(x, y, xi,..)
+        x = varargin(1)
+        if type(x)<>1 | ~isreal(x) then
+            msg = _("%s: Argument #%d: Real numbers expected.\n")
+            error(msprintf(msg,"interp1",1))
+        elseif ~isvector(x)
+            msg = _("%s: Argument #%d: Vector expected.\n")
+            error(msprintf(msg, "interp1", 1))
+        elseif or(isnan(x)) then
+            msg = _("%s: Argument #%d: Nan value forbidden.\n")
+            error(msprintf(msg, "interp1", 1))
+        end
+        x = matrix(x,1,-1)
+        y = varargin(2)
+    else // Default x: 
+        if type(varargin(1))<>0 // interp1(y, xi,..)
+            iy = 1
+            y = varargin(1)
+        else                    // interp1(, y, xi,..)
+            y = varargin(2)
+        end
+        if isvector(y) then
+            x = 1:size(y,"*")
+            y = y(:)
         else
-            error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector or matrix expected.\n"),"interp1",2));
+            x = 1:size(y,1)
         end
-        for i=rhs:-1:1
-            varargin(i+1)=varargin(i)
+    end
+    if length(x) <= 1
+        if iy==2
+            msg = _("%s: Argument #%d: At least 2 elements expected.\n")
+            error(msprintf(msg, "interp1", 1))
+        else
+            msg = _("%s: Argument #%d: vector or at least 2 rows expected.\n")
+            error(msprintf(msg, "interp1", 2))
         end
-        varargin(1)=X
     end
+    // Sorting x in increasing order:
+    [?, k] = gsort(x,"g","i")
+    x = x(k)    // should anyway be a row
 
-    //............................
-    //ininialisation of xi
-    //............................
-    //xi components must be reals
-    xi=varargin(3)
-    if type(xi)<>1 then
-        error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of reals expected.\n"),"interp1",3));
-    end
-    //delete the dimension of xi equal to one after the second dimension
-    //or the first dimension
-    xisize=size(xi);
-    while size(xisize,"*")>=2 & xisize($)==1
-        xisize=xisize(1:$-1);
+    // Checking y
+    // ----------
+    if type(y) <> 1 then
+        msg = _("%s: Argument #%d: Decimal or complex numbers expected.\n")
+        error(msprintf(msg, "interp1", iy));
     end
-    xisizetemp=xisize
-    if size(xisize,"*")>=2 then
-        if xisize(1)==1 then
-            xisize=xisize(2:$);
+    if isvector(y) then
+        y = y(:)
+        if length(y) <> length(x)
+            msg = _("%s: Arguments #%d and #%d: Same numbers of elements expected.\n")
+            error(msprintf(msg, "interp1", 1, 2))
+        end
+    else
+        if size(y,1) <> length(x) then
+            msg = _("%s: Arguments #%d and #%d: Same numbers of rows expected.\n")
+            error(msprintf(msg, "interp1", 1, 2))
         end
     end
+    // Sorting y according to x:
+    cmd = "y = y(k," + strcat(repmat(":",[1 ndims(y)-1]),",")+")"
+    execstr(cmd)
 
-    //-------------------------
-    //Initialisation of x, y
-    //-------------------------
-    x=varargin(1);
-    y=varargin(2);
-    //x must be real vector
-    if type(x)<>1 then
-        error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of reals expected.\n"),"interp1",1));
+    // Checking xi
+    // -----------
+    if type(varargin(iy+1))<>1 | ~isreal(varargin(iy+1)) then
+        msg = _("%s: Argument #%d: Decimal numbers expected.\n")
+        error(msprintf(msg, "interp1", iy+1))
     end
-    //y components must be reals
-    if type(y)<>1 then
-        error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of reals expected.\n"),"interp1",2));
+    xi = varargin(iy+1)
+    if xi == [] then
+        yi = []
+        return
     end
-    //verification of x,y line/column
-    if isvector(x) then
-        if find(isnan(x))<>[] then
-            error(msprintf(gettext("%s: Wrong value for input argument #%d: Reals expected but some NaN found.\n"),"interp1",1));
-        end
-        if isvector(y) then
-            if size(x,"*")<>size(y,"*") then
-                error(msprintf(gettext("%s: Wrong size for input arguments #%d and #%d: Same size expected.\n"),"interp1",1,2));
-            end
-        elseif size(size(y),"*")>=2 then
-            if size(x,"*")<>size(y,1) then
-                error(msprintf(gettext("%s: Wrong size for input arguments #%d and #%d: Same size expected.\n"),"interp1",1,2));
-            end
-        else
-            error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector or matrix expected.\n"),"interp1",2));
+    xisize = size(xi)
+    if length(size(y))==2 & or(size(y)==1) then  // y = 2D vector
+        if iscolumn(xisize) then
+            xisize = length(xi)
         end
     else
-        error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector expected.\n"),"interp1",1));
-    end
-
-    // xi : increase order sorting (for xi)
-    [xtemp,p]=gsort(matrix(x,1,-1),"c","i")
-    x=matrix(xtemp,size(x))
-    x=matrix(x,1,-1)
-    if isvector(y) then
-        y=y(p)
-    elseif size(size(y),"*") then
-        for l=1:size(y,"*")/size(y,1)
-            y(:,l)=y(p,l)
+        if length(size(xi))==2 & or(size(xi)==1) then // xi = scalar or 2D vector
+            xisize = length(xi)
         end
-    else
-        error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector or matrix expected.\n"),"interp1",2));
-    end
-
-    //-------------------------------------------------
-    // CASE : 3 inputs arguments : yi=interp1(x,y,xi)
-    //-------------------------------------------------
-
-    //default method : linear method is used
-    if size(varargin)==3 then
-        yi=interp1(x,y,xi,"linear",%nan)
     end
 
-    //--------------------------------------------------
-    // CASE : 4 inputs arguments : yi=interp1(x,y,xi,method)
-    //--------------------------------------------------
-
-    if size(varargin)==4 then
-        select part(varargin(4),1)
-            //-------------------------------------------
-            // Linear method : yi=linear(x,y,xi,'linear')
-            //-------------------------------------------
-            // the values of extrapolation points are nan for linear method
-        case "l"
-
-            yi=interp1(x,y,xi,"linear",%nan)
-
-            //-------------------------------------------
-            // Spline method  yi=interp1(x,y,xi,'spline')
-            //-------------------------------------------
-            // the extrapolation used the spline method
-        case "s"
-            if xi==[] then
-                yi=[]
-                return
-            end
-
-            yi=interp1(x,y,xi,"spline","extrap")
-            //----------------------------------------------
-            // Nearest method  yi=interp1(x,y,xi,'nearest')
-            //----------------------------------------------
-            // the values of extrapolation points are nan for nearest method
-        case "n"
-            if xi==[] then
-                yi=[]
-                return
-            end
-            yi=interp1(x,y,xi,"nearest",%nan)
-        else
-            error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or ''%s'' or ''%s'' expected.\n"),"interp1",4,"linear","nearest"));
+    // Checking method
+    // ---------------
+    if rhs > iy+1 & type(varargin(iy+2))<>0 then
+        method = varargin(iy+2)
+        if size(method,"*") <> 1
+            msg = _("%s: Argument #%d: Scalar (1 element) expected.\n")
+            error(msprintf(msg, "interp1", iy+2))
+        end
+        if type(method) <> 10 | ~or(convstr(method)==["linear","nearest","spline"])
+            msg = _("%s: Argument #%d: Must be in the set {%s}.\n")
+            error(msprintf(msg, "interp1", iy+2, """linear"",""spline"",""nearest"""))
         end
+        method = convstr(method)
+    else
+        method = "linear"
     end
 
-    //-------------------------------------------------------------------------------------------------
-    // CASE : 5 inputs arguments :  yi=interp1(x,y,xi,method,'extrap')  or  yi=interp1(x,y,xi,method,extrapval)
-    //-------------------------------------------------------------------------------------------------
-
-    if size(varargin)==5 then
-        select part(varargin(4),1)
-            //-----------------------------------------------------------------------------------
-            // Linear method : yi=linear(x,y,xi,'linear','extrap') or  yi=interp1(x,y,xi,method,extrapval)
-            //------------------------------------------------------------------------------------
-        case "l"
-            xitemp=matrix(xi,-1,1)
-            // y is a vector
-            if isvector(y) then
-                yi = resize_matrix(0, size(xitemp))
-                [x,ind]=gsort(matrix(x,1,-1),"c","i")
-                if varargin(5)==%nan then
-                    yi=linear_interpn(xitemp,x,y(ind),"by_nan");
-                end
-                if type(varargin(5))==10 then
-                    if varargin(5)<>"extrap" then
-                        error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
-                    else
-                        yi=linear_interpn(xitemp,x,y(ind),"natural");
+    // Checking extrapolation method
+    // -----------------------------
+    extraVal = []
+    if rhs > iy+2 then
+        extra = varargin(iy+3)
+        err = %f
+        if type(extra)==10 & size(extra,"*")==1
+            extra = convstr(extra)
+            if extra=="extrap"
+                extra = "natural"
+            elseif or(method==["linear" "spline"])
+                if extra=="linear" then
+                    if method=="linear"
+                        extra = "natural"
+                    end
+                elseif extra=="edgevalue"
+                    extra = "C0"
+                elseif extra=="periodic"
+                    msg = []
+                    if isvector(y) & y($)<>y(1)
+                        msg = _("%s: Argument #%d: periodicity y($)==y(1) expected.\n")
+                    elseif ~isvector(y) & or(y($,:)<>y(1,:))
+                        msg = _("%s: Argument #%d: periodicity y($,:)==y(1,:) expected.\n")
                     end
-                elseif type(varargin(5))==1  then
-                    yi=linear_interpn(xitemp,x,y(ind),"by_nan");
-                    if ~isnan(varargin(5)) then
-                        k=find(xitemp>max(x)|xitemp<min(x))
-                        yi(k)=varargin(5)
+                    if msg <> []
+                        error(msprintf(msg, "interp1",iy))
                     end
+                elseif extra<>"natural"
+                    err = %t
                 end
-                if size(xisize,"*")>=2
-                    yi=matrix(yi,xisize)
+            else // nearest
+                if or(extra==["edgevalue" "nearest" "natural"])
+                    extra = "natural"
                 else
-                    yi=matrix(yi,xisizetemp)
-                end
-
-                // y is matrix or hypermatrix
-            elseif size(size(y),"*")>=2 then
-                ysize=size(y)
-                ky=ysize(2:$)
-                yi = resize_matrix(0, [size(xitemp) ky])
-                [x,ind]=gsort(matrix(x,1,-1),"c","i")
-                //extrapolation
-                if type(varargin(5))==10 then
-                    if varargin(5)<>"extrap" then
-                        error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
-                    else
-                        if xitemp==[] then
-                            yi=[]
-                            return
-                        end
-                        for l=1:size(y,"*")/size(y,1)
-                            ytemp=y(:,l)
-                            yi(:,l)=matrix(linear_interpn(xitemp,x,ytemp(ind),"natural"),size(xitemp))
-                        end
-                    end
-                elseif type(varargin(5))==1 then
-                    if xitemp==[] then
-                        yi=[]
-                        return
-                    end
-                    for l=1:size(y,"*")/size(y,1)
-                        ytemp=y(:,l)
-                        yi(:,l)=matrix(linear_interpn(xitemp,x,ytemp(ind),"by_nan"),size(xitemp))
-                    end
-                    if ~isnan(varargin(5)) then
-                        k=find(xitemp>max(x)|xitemp<min(x))
-                        yi(k,:)=varargin(5)
-                    end
+                    err = %t
                 end
-                yi=matrix(yi,[xisize,ky])
+            end
+        else
+            if type(extra)<>1 | length(extra)<>1
+                err = %t
             else
-                error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector or matrix expected.\n"),"interp1",2));
+                extraVal = extra
+                extra = "natural"
             end
+        end
+        if err
+            if method=="nearest"
+                tmp = """extrap"",""edgevalue"",""nearest"""
+            else
+                tmp = """extrap"",""linear"",""periodic"",""edgevalue"""
+            end
+            msg = _("%s: Argument #%d: scalar number or one of {%s} expected.\n")
+            error(msprintf(msg, "interp1", iy+3, tmp))
+        end
+    else
+        extra = "natural"
+        if method <> "spline"
+            extraVal = %nan
+        end
+    end
 
-            //-------------------------------------------------------------------------------------
-            // Spline method  yi=interp1(x,y,xi,'spline','extrap') or  yi=interp1(x,y,xi,'spline',extrapval)
-            //-------------------------------------------------------------------------------------
-        case "s"
-            if xi==[] then
-                if varargin(5)=="extrap"|type(varargin(5))==1 then
-                    yi=[]
-                    return
-                else
-                    error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
+    // ==========
+    // PROCESSING
+    // ==========
+    if ~isreal(y) then
+        if extraVal==[]
+            [extraR, extraI] = (extra, extra);
+        else
+            [extraR, extraI] = (real(extraVal), imag(extraVal));
+        end
+        yir = interp1(x, real(y), xi, method, extraR)
+        yii = interp1(x, imag(y), xi, method, extraI)
+        yi = complex(yir, yii)
+        return
+    end
+    // REAL CASE
+    // =========
+    select part(method, 1)
+    case "l"
+        // LINEAR
+        // ------
+        xitemp = matrix(xi,-1,1)
+        // y is a vector
+        if isvector(y) then
+            yi = resize_matrix(0, size(xitemp))
+            if extraVal <> [] then
+                yi = linear_interpn(xitemp, x, y, "by_nan");
+                if ~isnan(extraVal) then
+                    k = find(xitemp<min(x) | xitemp>max(x))
+                    yi(k) = extraVal
                 end
+            else
+                yi = linear_interpn(xitemp, x, y, extra);
             end
-            xitemp=matrix(xi,-1,1)
-            //y is a vector
-            if isvector(y) then
-                yi = resize_matrix(0, size(xitemp))
-                yi=interp(xitemp,matrix(x,1,-1),matrix(y,1,-1),splin(matrix(x,1,-1),matrix(y,1,-1)),"natural");
-                if type(varargin(5))==10 then
-                    if varargin(5)<>"extrap" then
-                        error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
-                    end
-                elseif type(varargin(5))==1  then
-                    k=find(xitemp>max(x)|xitemp<min(x))
-                    yi(k)=varargin(5)
+            yi = matrix(yi, xisize)
+
+        else    // y is matrix or hypermatrix
+            ky = size(y)(2:$)
+            ncolumns = prod(ky)
+            yi = resize_matrix(0, [size(xitemp) ky])
+            if extraVal <> [] then
+                for ic = 1:ncolumns
+                    tmp = linear_interpn(xitemp, x, y(:,ic), "by_nan")
+                    yi(:,ic) = matrix(tmp, size(xitemp))
                 end
-                if size(xisize,"*")>=2
-                    yi=matrix(yi,xisize)
-                else
-                    yi=matrix(yi,xisizetemp)
+                if ~isnan(extraVal) then
+                    k = find(xitemp < min(x) | xitemp > max(x))
+                    yi(k,:) = extraVal
                 end
-                //y is a matrix or a hypermatrix
-            elseif size(size(y),"*")>=2 then
-                ky=size(y)
-                ky=ky(2:$)
-                yi = resize_matrix(0, [size(xitemp) ky])
-                for l=1:size(y,"*")/size(y,1)
-                    yi(:,l)=matrix(interp(matrix(xi,-1,1),matrix(x,-1,1),y(:,l),splin(matrix(x,-1,1),y(:,l)),"natural"),size(xitemp))//les composante de yi
+            else // extra=="linear"=="natural" | "periodic" | "edgevalue"=="C0"
+                for ic = 1:ncolumns
+                    tmp = linear_interpn(xitemp, x, y(:,ic), extra)
+                    yi(:,ic) = matrix(tmp, size(xitemp))
                 end
-                //extrapolation
-                if type(varargin(5))==10 then
-                    if varargin(5)<>"extrap" then
-                        error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
-                    end
-                elseif type(varargin(5))==1  then
-                    k=find(xitemp>max(x)|xitemp<min(x))
-                    yi(k,:)=varargin(5)
+            end
+            yi = matrix(yi, [xisize, ky])
+        end
+
+    case "s"
+        // SPLINE
+        // ------
+        xitemp = matrix(xi,-1,1)
+        // y is a vector
+        if isvector(y) then
+            yi = resize_matrix(0, size(xitemp))
+            yrow = matrix(y,1,-1)
+            if extra=="linear"
+                yi = interp(xitemp, x, yrow, splin(x, yrow, "natural"), extra)
+            elseif extra=="periodic"
+                yi = interp(xitemp, x, yrow, splin(x, yrow, extra), extra)
+            else    // default not_a_knot
+                yi = interp(xitemp, x, yrow, splin(x, yrow), extra)
+            end
+            if extraVal <> []  then
+                k = find(xitemp<min(x) | xitemp>max(x))
+                yi(k) = extraVal
+            end
+            yi = matrix(yi, xisize)
+
+        else    // y is a matrix or a hypermatrix
+            ky = size(y)(2:$)
+            ncolumns = prod(ky)
+            yi = resize_matrix(0, [size(xitemp) ky])
+            x  = matrix(x,-1,1)
+            xi = matrix(xi,-1,1)
+            for ic = 1:ncolumns
+                tmp = y(:,ic)
+                tmp = interp(xi, x, tmp, splin(x, tmp), extra)
+                yi(:,ic) = matrix(tmp, size(xitemp))
+            end
+            // extrapolation
+            if extraVal <> []  then
+                k = find(xitemp<min(x) | xitemp>max(x))
+                yi(k,:) = extraVal
+            end
+            yi = matrix(yi,[xisize,ky])
+        end
+
+    case "n"
+        // NEAREST
+        // -------
+        // if all xi values are nan, retuns nan values for yi
+        if and(isnan(xi)) then
+            yi = xi
+            return
+        end
+
+        xitemp = matrix(xi,1,-1)
+        knan = isnan(xitemp)
+        [xitemp,p] = gsort(xitemp(~knan),"g","i")
+        // For each xi, we search for the k index of the x nearest to xi:
+        k = zeros(xitemp);
+        nx = length(x);
+        j = length(xitemp);
+        i = nx;
+        while j>=1 & i>=1
+            if xitemp(j) >= x(i)  then
+                if i <> nx then
+                    k(j) = i;
                 end
-                yi=matrix(yi,[xisize,ky])
+                j = j - 1;
             else
-                error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector or matrix expected.\n"),"interp1",2));
+                i = i - 1;
             end
+        end
+        k(xitemp < x(1)) = 1;
+        k(xitemp >= x(nx)) = nx - 1; 
+        i = find(xitemp >= matrix((x(k)+x(k+1))/2,size(k)));
+        if i~=[]
+            k(i) = k(i)+1;
+        end
 
-            //---------------------------------------------------------------------------------------
-            // Nearest method  yi=interp1(x,y,xi,'nearest','extrap') or yi=interp1(x,y,xi,'nearest',extrapval)
-            //---------------------------------------------------------------------------------------
-        case "n"
-            //if all xi values are nan, retuns nan values for yi
-            if size(find(isnan(xi)),"*")==size(xi,"*") then
-                if varargin(5)=="extrap"|type(varargin(5))==1 then
-                    yi=xi
-                    return
-                else
-                    error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
-                end
+        // y is vector
+        if isvector(y) then
+            yi = y(k)
+            yi = matrix(yi,1,-1)
+            // extrapolation
+            if extraVal <> [] then
+                n = find(xitemp<min(x) | xitemp>max(x) | isnan(xitemp))
+                yi(n) = extraVal
             end
-            if xi==[] then
-                if varargin(5)=="extrap"|type(varargin(5))==1 then
-                    yi=[]
-                    return
-                else
-                    error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
-                end
+            yitemp = yi
+            yi(p)  = yitemp
+            ytemp  = yi
+            yi = matrix(xi,1,-1)
+            yi(knan) = %nan
+            yi(~knan) = ytemp
+            yi = matrix(yi, xisize)
+
+        else  // y is a matrix or a hypermatrix
+            ky = size(y)(2:$)
+            yi = resize_matrix(0, [size(xitemp,"*") ky])
+            ncolumns = prod(ky)
+            for ic = 1:ncolumns
+                yi(:,ic) = y(k,ic)
             end
-            //y is vector
-            if isvector(y) then
-                xitemp=matrix(xi,1,-1)
-                knan=find(isnan(xitemp))
-                knotnan=find(~isnan(xitemp))
-                [xitemp,p]=gsort(matrix(xitemp(knotnan),1,-1),"c","i")
-                yi=matrix(xi,1,-1)
-                k=zeros(xitemp)
-                x_size=size(x,"*")
-                j=size(xitemp,"*")
-                i=x_size
-                while j>=1 & i>=1
-                    if xitemp(j)>=x(i)  then
-                        if i<>x_size then
-                            k(j)=i
-                        end
-                        j=j-1
-                    else
-                        i=i-1
-                    end
-                end
-                k(xitemp<x(1)) = 1;
-                k(xitemp>=x(x_size)) = x_size-1;
-                i = find(xitemp >= matrix((x(k)+x(k+1))/2,size(k)));
-                if i~=[]
-                    k(i) = k(i)+1;
-                end
-                yi=y(k)
-                yi=matrix(yi,1,-1)
-                //extrapolation
-                if type(varargin(5))==10 then
-                    if varargin(5)<>"extrap" then
-                        error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
-                    end
-                elseif type(varargin(5))==1 then
-                    n=find(xitemp>max(x)|xitemp<min(x)|isnan(xitemp)|isnan(xitemp))
-                    yi(n)=varargin(5)
-                else
-                    error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
-                end
-                yitemp=yi
-                yi(p)=yitemp
-                ytemp=yi
-                yi=matrix(xi,1,-1)
-                yi(knan)=%nan
-                yi(knotnan)=ytemp
-                if size(xisize,"*")>=2
-                    yi=matrix(yi,xisize)
-                else
-                    yi=matrix(yi,xisizetemp)
-                end
-                //y is a matrix or a hypermatrix
-            elseif size(size(y),"*")>=2 then
-                xitemp=matrix(xi,1,-1)
-                knan=find(isnan(xitemp))
-                knotnan=find(~isnan(xitemp))
-                [xitemp,p]=gsort(xitemp(knotnan),"c","i")
-                ind=size(y)
-                ind=ind(2:$)
-                yi = resize_matrix(0, [size(xitemp,"*") ind])
-                k=zeros(xitemp)
-                x_size=size(x,"*")
-                j=size(xitemp,"*")
-                i=x_size
-                while j>=1 & i>=1
-                    if xitemp(j)>=x(i)  then
-                        if i<>x_size then
-                            k(j)=i
-                        end
-                        j=j-1
-                    else
-                        i=i-1
-                    end
+            // extrapolation
+            if extraVal <> [] then
+                n = find(xitemp < min(x) | xitemp > max(x))
+                for ic = 1:ncolumns
+                    yi(n,ic) = extraVal
                 end
-                k(xitemp<x(1)) = 1;
-                k(xitemp>=x(x_size)) = x_size-1;
-                i = find(xitemp >= matrix((x(k)+x(k+1))/2,size(k)));
-                if i~=[]
-                    k(i) = k(i)+1;
-                end
-                for l=1:size(y,"*")/size(y,1)
-                    ytemp=matrix(y(:,l),1,-1)
-                    yi(:,l) =ytemp(k)
-                end
-                //extrapolation
-                if type(varargin(5))==10 then
-                    if varargin(5)<>"extrap" then
-                        error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
-                    end
-                elseif type(varargin(5))==1 then
-                    n=find(xitemp>max(x)|xitemp<min(x))
-                    for l=1:size(y,"*")/size(y,1)
-                        yi(n,l)=varargin(5)
-                    end
-                else
-                    error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
-                end
-                yitemp=yi
-                for l=1:size(y,"*")/size(y,1)
-                    yi(p,l)=yitemp(:,l)
-                end
-                yitemp=yi
-                yi = resize_matrix(0, [size(xi,"*") ind])
-                for l=1:size(y,"*")/size(y,1)
-                    yi(knan,l)=%nan
-                    yi(knotnan,l)=yitemp(:,l)
-                end
-                yi=matrix(yi,[xisize,ind])
-            else
-                error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector or matrix expected.\n"),"interp1",2));
             end
-        else
-            error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or real expected.\n"),"interp1",5,"extrap"));
+            yitemp = yi
+            for ic = 1:ncolumns
+                yi(p,ic) = yitemp(:,ic)
+            end
+            yitemp = yi
+            yi = resize_matrix(0, [size(xi,"*") ky])
+            for ic = 1:ncolumns
+                yi(knan,ic) = %nan
+                yi(~knan,ic) = yitemp(:,ic)
+            end
+            yi = matrix(yi,[xisize,ky])
         end
     end
 endfunction
diff --git a/scilab/modules/interpolation/tests/unit_tests/interp1.tst b/scilab/modules/interpolation/tests/unit_tests/interp1.tst
new file mode 100644 (file)
index 0000000..3a62bc7
--- /dev/null
@@ -0,0 +1,654 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2021 - Samuel GOUGEON - Le Mans Université
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- CLI SHELL MODE -->
+// <-- NO CHECK REF -->
+
+// Checking the size of the result
+// ===============================
+for m = ["linear", "spline", "nearest"]
+    for x = list(1:0.5:4, (1:0.5:4)')
+        // Y vector => size(result) == size(xi)
+        Xi = list(%pi, 1.2:0.2:2, (1.2:0.2:2)', 1+rand(2,5)*4, ..
+             1+rand(1,2,5)*4, 1+rand(2,1,5)*4, 1+rand(2,5,6)*4);
+        Y = list(5:11, (5:11)');
+        for y = Y
+            for xi = Xi
+                r = interp1(x, y, xi, m);
+                assert_checkequal(size(r), size(xi));
+            end
+        end
+
+        // y = matrix or hypermatrix, xi = scalar or vector
+        Y = list(rand(7,3), rand(7,1,3), rand(7,3,4));
+        Xi = list(%pi, 1.2:0.2:2, (1.2:0.2:2)');
+        for y = Y
+            for xi = Xi
+                r = interp1(x, y, xi, m);
+                assert_checkequal(size(r), [length(xi) size(y)(2:$)]);
+            end
+        end
+
+        // Y and xi are matrices or hypermatrices
+        Y = list(rand(7,3), rand(7,1,3), rand(7,3,4));
+        Xi = list(1+rand(2,5)*4, 1+rand(1,2,5)*4, 1+rand(2,1,5)*4, 1+rand(2,5,6)*4);
+        for y = Y
+            for xi = Xi
+                r = interp1(x, y, xi, m);
+                assert_checkequal(size(r), [size(xi) size(y)(2:$)]);
+            end
+        end
+    end
+end
+
+
+// =============
+// y is a vector
+// =============
+x = 1:4;
+// method = "linear" (default)
+// ---------------------------
+// extrapolation "by_nan" (default)
+for y = list(2:5, (2:5)')
+    assert_checkequal(interp1(x,y, 0), %nan);
+    assert_checkequal(interp1(x,y,1.5), 2.5);
+    assert_checkequal(interp1(x,y,[0.5:5]), [%nan 2.5 3.5 4.5 %nan]);
+    assert_checkequal(interp1(x,y,[0.5:5]'), [%nan 2.5 3.5 4.5 %nan]');
+    assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5]), [%nan 2.5 3.5;4.5 %nan %nan]);
+    Ref = cat(3, [%nan 2.5 3.5], [4.5 %nan %nan]);
+    assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5])), Ref);
+end
+// extrapolation "extrap" = "linear"
+for y = list(2:5, (2:5)')
+    assert_checkequal(interp1(x,y, 0, ,"extrap"), 1);
+    assert_checkequal(interp1(x,y,1.5, ,"extrap"), 2.5);
+    assert_checkequal(interp1(x,y,[0.5:5], ,"extrap"), [1.5 2.5 3.5 4.5 5.5]);
+    assert_checkequal(interp1(x,y,[0.5:5]', ,"extrap"), [1.5 2.5 3.5 4.5 5.5]');
+    assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], ,"extrap"), ..
+                      [1.5 2.5 3.5;4.5 5.5 6.5]);
+    Ref = cat(3, [1.5 2.5 3.5], [4.5 5.5 6.5]);
+    assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), ,"extrap"), Ref);
+end
+// extrapolation by padding
+for y = list(2:5, (2:5)')
+    assert_checkequal(interp1(x,y, 0, ,%pi), %pi);
+    assert_checkequal(interp1(x,y,1.5, ,%pi), 2.5);
+    assert_checkequal(interp1(x,y,[0.5:5], ,%pi), [%pi 2.5 3.5 4.5 %pi]);
+    assert_checkequal(interp1(x,y,[0.5:5]', ,%pi), [%pi 2.5 3.5 4.5 %pi]');
+    assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], ,%pi), ..
+                      [%pi 2.5 3.5;4.5 %pi %pi]);
+    Ref = cat(3, [%pi 2.5 3.5], [4.5 %pi %pi]);
+    assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), ,%pi), Ref);
+end
+// periodic extrapolation
+for y = list([2 3 4 2], [2 3 4 2]')
+    assert_checkequal(interp1(x,y, 0, ,"periodic"), 4);
+    assert_checkequal(interp1(x,y,1.5, ,"periodic"), 2.5);
+    assert_checkequal(interp1(x,y,[0.5:5], ,"periodic"), [3 2.5 3.5 3 2.5]);
+    assert_checkequal(interp1(x,y,[0.5:5]', ,"periodic"), [3 2.5 3.5 3 2.5]');
+    assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], ,"periodic"), ..
+                      [3 2.5 3.5 ; 3 2.5 3.5]);
+    Ref = cat(3, [3 2.5 3.5], [3 2.5 3.5]);
+    assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), ,"periodic"), Ref);
+end
+// extrapolation by edgevalue
+e = "edgevalue";
+for y = list(2:5, (2:5)')
+    assert_checkequal(interp1(x,y, 0, , e), 2);
+    assert_checkequal(interp1(x,y,1.5, , e), 2.5);
+    assert_checkequal(interp1(x,y,[0.5:5], , e), [2 2.5 3.5 4.5 5]);
+    assert_checkequal(interp1(x,y,[0.5:5]', , e), [2 2.5 3.5 4.5 5]');
+    assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], , e), [2 2.5 3.5 ; 4.5 5 5]);
+    Ref = cat(3, [2 2.5 3.5], [4.5 5 5]);
+    assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), , e), Ref);
+end
+
+// method = "nearest"
+// ------------------
+m = "nearest";
+// extrapolation "by_nan" (default)
+for y = list(2:5, (2:5)')
+    assert_checkequal(interp1(x,y, [0 5.4], m), [%nan %nan]);
+    assert_checkequal(interp1(x,y,1.5, m), 3);
+    assert_checkequal(interp1(x,y,[0.5:5], m), [%nan 3 4 5 %nan]);
+    assert_checkequal(interp1(x,y,[0.5:5]', m), [%nan 3 4 5 %nan]');
+    assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m), [%nan 3 4;5 %nan %nan]);
+    Ref = cat(3, [%nan 3 4], [5 %nan %nan]);
+    assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m), Ref);
+end
+// extrapolation "extrap" = "nearest" = "edgevalue"
+for e = ["extrap" "edgevalue" "nearest"]
+    for y = list(2:5, (2:5)')
+        assert_checkequal(interp1(x,y, [0 5.4], m, e), [2 5]);
+        assert_checkequal(interp1(x,y,1.5, m, e), 3);
+        assert_checkequal(interp1(x,y,[0.5:5], m, e), [2 3 4 5 5]);
+        assert_checkequal(interp1(x,y,[0.5:5]', m, e), [2 3 4 5 5]');
+        assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, e), ..
+                          [2 3 4 ; 5 5 5]);
+        Ref = cat(3, [2 3 4], [5 5 5]);
+        assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, e), Ref);
+    end
+end
+// extrapolation by padding
+for y = list(2:5, (2:5)')
+    assert_checkequal(interp1(x,y, [0 5.4], m, %pi), [%pi %pi]);
+    assert_checkequal(interp1(x,y, 1.5, m, %pi), 3);
+    assert_checkequal(interp1(x,y, 0.5:5, m, %pi), [%pi 3 4 5 %pi]);
+    assert_checkequal(interp1(x,y,[0.5:5]', m, %pi), [%pi 3 4 5 %pi]');
+    assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, %pi), ..
+                      [%pi 3 4 ; 5 %pi %pi]);
+    Ref = cat(3, [%pi 3 4], [5 %pi %pi]);
+    assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, %pi), Ref);
+end
+
+// method = "spline"
+// -----------------
+// extrapolation default == "spline" (=="natural")
+m = "spline";
+for y = list(2:5, (2:5)')
+    assert_checkalmostequal(interp1(x,y, [0 5.4], m), [1 6.4], 5*%eps);
+    assert_checkequal(interp1(x,y,1.5, m), 2.5);
+    assert_checkalmostequal(interp1(x,y,[0.5:5], m), 1.5:5.5, %eps);
+    assert_checkalmostequal(interp1(x,y,[0.5:5]', m), (1.5:5.5)', %eps);
+    assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m), [1.5 2.5 3.5;4.5 5.5 6.5], 5*%eps);
+    Ref = cat(3, [1.5 2.5 3.5], [4.5 5.5 6.5]);
+    assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m), Ref, 5*%eps);
+end
+// extrapolation "extrap" == "spline"
+for y = list(2:5, (2:5)')
+    assert_checkalmostequal(interp1(x,y, [0 5.4], m, "extrap"), [1 6.4], 5*%eps);
+    assert_checkequal(interp1(x,y,1.5, m, "extrap"), 2.5);
+    assert_checkalmostequal(interp1(x,y,[0.5:5], m, "extrap"), 1.5:5.5, %eps);
+    assert_checkalmostequal(interp1(x,y,[0.5:5]', m, "extrap"), (1.5:5.5)', %eps);
+    Ref = [1.5 2.5 3.5 ; 4.5 5.5 6.5];
+    assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, "extrap"), Ref, 5*%eps);
+    Ref = cat(3, [1.5 2.5 3.5], [4.5 5.5 6.5]);
+    assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, "extrap"), Ref, 5*%eps);
+end
+// extrapolation by padding
+for y = list(2:5, (2:5)')
+    assert_checkequal(interp1(x,y, [0 5.4], m, %pi), [%pi %pi]);
+    assert_checkequal(interp1(x,y, 1.5, m, %pi), 2.5);
+    assert_checkequal(interp1(x,y, 0.5:5, m, %pi), [%pi 2.5 3.5 4.5 %pi]);
+    assert_checkequal(interp1(x,y,[0.5:5]', m, %pi), [%pi 2.5 3.5 4.5 %pi]');
+    assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, %pi), ..
+                      [%pi 2.5 3.5 ; 4.5 %pi %pi]);
+    Ref = cat(3, [%pi 2.5 3.5], [4.5 %pi %pi]);
+    assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, %pi), Ref);
+end
+// periodic extrapolation
+for y = list([2 3 4 2], [2 3 4 2]')
+    assert_checkalmostequal(interp1(x,y, [0 5.4], m, "periodic"), [4 3.736], %eps);
+    assert_checkequal(interp1(x,y, 1.5, m, "periodic"), 2.125);
+    assert_checkalmostequal(interp1(x,y, 0.5:5, m, "periodic"), [3 2.125 3.875 3 2.125], %eps);
+    assert_checkalmostequal(interp1(x,y,[0.5:5]', m, "periodic"), [3 2.125 3.875 3 2.125]', %eps);
+    assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, "periodic"), ..
+                      [3 2.125 3.875 ; 3 2.125 3.875], %eps);
+    Ref = cat(3, [3 2.125 3.875], [3 2.125 3.875]);
+    assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, "periodic"), Ref, %eps);
+end
+// extrapolation by edgevalue
+e = "edgevalue";
+for y = list(2:5, (2:5)')
+    assert_checkequal(interp1(x,y, [0 5.4], m, e), [2 5]);
+    assert_checkequal(interp1(x,y, 1.5, m, e), 2.5);
+    assert_checkequal(interp1(x,y, 0.5:5, m, e), [2 2.5 3.5 4.5 5]);
+    assert_checkequal(interp1(x,y,[0.5:5]', m, e), [2 2.5 3.5 4.5 5]');
+    assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, e), ..
+                                  [2 2.5 3.5 ; 4.5 5 5]);
+    Ref = cat(3, [2 2.5 3.5], [4.5 5 5]);
+    assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, e), Ref);
+end
+
+// =============
+// y is a matrix
+// =============
+x = 1:4;
+y = [0 2 4 6 ; 0.5 1.5 2.5 3.5]';
+// method = "linear" (default)
+// ---------------------------
+// extrapolation "by_nan" (default)
+assert_checkequal(interp1(x,y, 0), [%nan %nan]);
+assert_checkequal(interp1(x,y,1.5), [1 1]);
+assert_checkequal(interp1(x,y, 0.5:5  ), [%nan 1 3 5 %nan; %nan 1 2 3 %nan]');
+assert_checkequal(interp1(x,y,(0.5:5)'), [%nan 1 3 5 %nan; %nan 1 2 3 %nan]');
+Ref = cat(3,[%nan 1 3 ; 5 %nan %nan], [%nan 1 2 ; 3 %nan %nan]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5]), Ref);
+Ref = matrix([%nan 1 3 5 %nan %nan %nan 1 2 3 %nan %nan],[1 3 2 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5])), Ref);
+
+// extrapolation "extrap" = "linear"
+assert_checkequal(interp1(x,y, 0, ,"extrap"), [-2 -0.5]);
+assert_checkequal(interp1(x,y,1.5, ,"extrap"), [1 1]);
+assert_checkequal(interp1(x,y, 0.5:4,   ,"extrap"), [-1 1 3 5 ; 0 1 2 3]');
+assert_checkequal(interp1(x,y,(0.5:4)', ,"extrap"), [-1 1 3 5 ; 0 1 2 3]');
+Ref = cat(3, [-1 1 3 ; 5 7 9], [0 1 2 ; 3 4 5]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], ,"extrap"), Ref);
+Ref = matrix([-1 1 3 5 7 9 0 1 2 3 4 5],[1 3 2 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), ,"extrap"), Ref);
+
+// extrapolation by padding
+assert_checkequal(interp1(x,y, 0, ,%pi), [%pi %pi]);
+assert_checkequal(interp1(x,y,1.5, ,%pi), [1 1]);
+assert_checkequal(interp1(x,y, 0.5:5, ,%pi  ), [%pi 1 3 5 %pi; %pi 1 2 3 %pi]');
+assert_checkequal(interp1(x,y,(0.5:5)', ,%pi), [%pi 1 3 5 %pi; %pi 1 2 3 %pi]');
+Ref = cat(3,[%pi 1 3 ; 5 %pi %pi], [%pi 1 2 ; 3 %pi %pi]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], ,%pi), Ref);
+Ref = matrix([%pi 1 3 5 %pi %pi %pi 1 2 3 %pi %pi],[1 3 2 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), ,%pi), Ref);
+
+// periodic extrapolation
+y = [0 2 4 0 ; 0.5 1.5 2.5 0.5]';
+assert_checkalmostequal(interp1(x,y, 5, ,"periodic"), [2 1.5], 2*%eps);
+assert_checkalmostequal(interp1(x,y,1.5, ,"periodic"), [1 1], 2*%eps);
+assert_checkalmostequal(interp1(x,y, 0.5:5, ,"periodic"  ), [2 1 3 2 1; 1.5 1 2 1.5 1]', 2*%eps);
+assert_checkalmostequal(interp1(x,y,(0.5:5)', ,"periodic"), [2 1 3 2 1; 1.5 1 2 1.5 1]', 2*%eps);
+Ref = cat(3,[2 1 3 ; 2 1 3], [1.5 1 2 ; 1.5 1 2]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], ,"periodic"), Ref, 2*%eps);
+Ref = matrix([2 1 3 2 1 3 1.5 1 2 1.5 1 2],[1 3 2 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), ,"periodic"), Ref, 2*%eps);
+
+// extrapolation by edgevalue
+e = "edgevalue";
+y = [0 2 4 6 ; 0.5 1.5 2.5 3.5]';
+assert_checkequal(interp1(x,y, 0, , e), [0 0.5]);
+assert_checkequal(interp1(x,y,1.5, , e), [1 1]);
+assert_checkequal(interp1(x,y, 0.5:5, , e  ), [0 1 3 5 6 ; 0.5 1 2 3 3.5]');
+assert_checkequal(interp1(x,y,(0.5:5)', , e), [0 1 3 5 6 ; 0.5 1 2 3 3.5]');
+Ref = cat(3,[0 1 3 ; 5 6 6], [0.5 1 2 ; 3 3.5 3.5]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], , e), Ref);
+Ref = matrix([0 1 3 5 6 6 0.5 1 2 3 3.5 3.5],[1 3 2 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), ,e), Ref);
+
+// method = "nearest"
+// ------------------
+m = "nearest";
+// extrapolation "by_nan" (default)
+assert_checkequal(interp1(x,y, [0 5.4], m), [%nan %nan ; %nan %nan]);
+assert_checkequal(interp1(x,y,1.5, m), [2 1.5]);
+assert_checkequal(interp1(x,y,[0.5:4], m), [%nan %nan ; 2 1.5 ; 4 2.5 ; 6 3.5]);
+assert_checkequal(interp1(x,y,[0.5:4]', m),[%nan %nan ; 2 1.5 ; 4 2.5 ; 6 3.5]);
+Ref = cat(3, [%nan 2 4 ; 6 %nan %nan], [%nan 1.5 2.5 ; 3.5 %nan %nan]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m), Ref);  // OErr
+Ref = matrix([%nan 2 4 6 %nan %nan %nan 1.5 2.5 3.5 %nan %nan],[1 3 2 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m), Ref);
+
+// extrapolation "extrap" = "nearest" = "edgevalue"
+for e = ["extrap" "nearest" "edgevalue"]
+    assert_checkequal(interp1(x,y, [0 5.4], m, e), [0 0.5 ; 6 3.5]);
+    assert_checkequal(interp1(x,y,1.5, m, e), [2 1.5]);
+    assert_checkequal(interp1(x,y,[0.5:4], m, e), y);
+    assert_checkequal(interp1(x,y,[0.5:4]', m, e), y);
+    Ref = cat(3, [0 2 4 ; 6 6 6], [0.5 1.5 2.5 ; 3.5 3.5 3.5]);
+    assert_checkequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m, e), Ref);
+    Ref = matrix([0 2 4 6 6 6 0.5 1.5 2.5 3.5 3.5 3.5],[1 3 2 2]);
+    assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, e), Ref);
+end
+
+// extrapolation by padding
+assert_checkequal(interp1(x,y, [0 5.4], m, %pi), [%pi %pi ; %pi %pi]);
+assert_checkequal(interp1(x,y,1.5, m, %pi), [2 1.5]);
+assert_checkequal(interp1(x,y,[0.5:4], m, %pi), [%pi %pi ; 2 1.5 ; 4 2.5 ; 6 3.5]);
+assert_checkequal(interp1(x,y,[0.5:4]', m, %pi),[%pi %pi ; 2 1.5 ; 4 2.5 ; 6 3.5]);
+Ref = cat(3, [%pi 2 4 ; 6 %pi %pi], [%pi 1.5 2.5 ; 3.5 %pi %pi]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m, %pi), Ref); // OErr
+Ref = matrix([%pi 2 4 6 %pi %pi %pi 1.5 2.5 3.5 %pi %pi],[1 3 2 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, %pi), Ref); // OErr
+
+// method = "spline"
+// -----------------
+m = "spline";
+// extrapolation default == "spline" (=="natural")
+assert_checkalmostequal(interp1(x,y, [0 5.4], m), [-2 -0.5 ; 8.8 4.9], 5*%eps);
+assert_checkequal(interp1(x,y,1.5, m), [1 1]);
+assert_checkequal(interp1(x,y,[0.5:4], m), [-1 1 3 5 ; 0 1 2 3]');
+assert_checkequal(interp1(x,y,[0.5:4]', m), [-1 1 3 5 ; 0 1 2 3]');
+Ref = cat(3,[-1 1 3 ; 5 7 9], [0 1 2 ; 3 4 5]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m), Ref, 5*%eps);
+Ref = matrix([-1  1  3  5  7  9  0  1  2  3  4  5 ],[1 3 2 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m), Ref, 5*%eps);
+
+// extrapolation "extrap" == "spline"
+assert_checkalmostequal(interp1(x,y, [0 5.4], m, "extrap"), [-2 -0.5 ; 8.8 4.9], 5*%eps);
+assert_checkequal(interp1(x,y,1.5, m, "extrap"), [1 1]);
+assert_checkequal(interp1(x,y,[0.5:4], m, "extrap"), [-1 1 3 5 ; 0 1 2 3]');
+assert_checkequal(interp1(x,y,[0.5:4]', m, "extrap"), [-1 1 3 5 ; 0 1 2 3]');
+Ref = cat(3,[-1 1 3 ; 5 7 9], [0 1 2 ; 3 4 5]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, "extrap"), Ref, 5*%eps);
+Ref = matrix([-1  1  3  5  7  9  0  1  2  3  4  5 ],[1 3 2 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, "extrap"), Ref, 5*%eps);
+
+// extrapolation by padding
+assert_checkalmostequal(interp1(x,y, [0 5.4], m, %pi), [%pi %pi ; %pi %pi], 5*%eps);
+assert_checkequal(interp1(x,y,1.5, m, %pi), [1 1]);
+assert_checkequal(interp1(x,y,[0.5:4], m, %pi),  [%pi 1 3 5 ; %pi 1 2 3]');
+assert_checkequal(interp1(x,y,[0.5:4]', m, %pi), [%pi 1 3 5 ; %pi 1 2 3]');
+Ref = cat(3,[%pi 1 3 ; 5 %pi %pi], [%pi 1 2 ; 3 %pi %pi]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, %pi), Ref, 5*%eps);
+Ref = matrix([%pi 1  3  5  %pi %pi %pi  1  2  3  %pi %pi],[1 3 2 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, %pi), Ref);
+
+// extrapolation by edgevalue
+e = "edgevalue";
+assert_checkequal(interp1(x,y, [0 5.4], m, e), [0 0.5 ; 6 3.5]);
+assert_checkequal(interp1(x,y,1.5, m, e), [1 1]);
+assert_checkequal(interp1(x,y,[0.5:4], m, e),  [0 1 3 5 ; 0.5 1 2 3]');
+assert_checkequal(interp1(x,y,[0.5:4]', m, e), [0 1 3 5 ; 0.5 1 2 3]');
+Ref = cat(3, [0 1 3 ; 5 6 6], [0.5 1 2 ; 3 3.5 3.5]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, e), Ref);
+Ref = matrix([0 1 3 5 6 6 0.5 1 2 3 3.5 3.5], [1 3 2 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, e), Ref);
+
+// periodic extrapolation
+y = [0 2 4 0 ; 0.5 1.5 2.5 0.5]';
+assert_checkalmostequal(interp1(x,y, [0 5.4], m, "periodic"), [4 2.5 ; 3.136 2.068], 2*%eps);
+assert_checkalmostequal(interp1(x,y,1.5, m, "periodic"), [0.625 0.8125], 5*%eps);
+Ref = [3.125  2.0625 ; 0.625  0.8125 ; 3.375  2.1875 ; 3.125  2.0625];
+assert_checkequal(interp1(x,y,[0.5:4], m, "periodic"),  Ref);
+assert_checkequal(interp1(x,y,[0.5:4]', m, "periodic"), Ref);
+Ref = cat(3, [3.125 0.625 3.375 ; 3.125 0.625 3.375], ..
+             [2.0625 0.8125 2.1875 ; 2.0625 0.8125 2.1875]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], m, "periodic"), Ref, 5*%eps);
+Ref = matrix([3.125 0.625 3.375 3.125 0.625 3.375 2.0625 0.8125 2.1875 2.0625 0.8125 2.1875],[1 3 2 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, "periodic"), Ref, 5*%eps);
+
+
+// ===================
+// y is an hypermatrix
+// ===================
+y = cat(4,[0 2 4 6 ; 0.5 1.5 2.5 3.5]', [1:4 ; 2:5]');
+// method = "linear" (default)
+// ---------------------------
+// extrapolation "by_nan" (default)
+assert_checkequal(interp1(x,y, 0), cat(4,[%nan %nan],[%nan %nan]));    // OWrong size and values
+assert_checkequal(interp1(x,y,1.5), cat(4,[1 1],[1.5 2.5]));           // OWrong size
+Ref = cat(4,[%nan 1 3 5 ; %nan 1 2 3]', [%nan 1.5:3.5 ; %nan 2.5:4.5]');
+assert_checkequal(interp1(x,y, 0.5:4), Ref);
+assert_checkequal(interp1(x,y,(0.5:4)'), Ref);
+Ref = matrix([%nan 5 1 %nan 3 %nan %nan 3 1 %nan 2 %nan %nan 3.5 1.5 %nan ..
+   2.5 %nan %nan 4.5 2.5 %nan 3.5 %nan],[2 3 2 1 2]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5]), Ref);
+Ref = matrix([%nan 1 3 5 %nan %nan %nan 1 2 3 %nan %nan %nan 1.5 2.5 3.5 ..
+              %nan %nan %nan 2.5 3.5 4.5 %nan %nan],[1 3 2 2 1 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5])), Ref);
+
+// extrapolation "extrap" = "linear"
+assert_checkequal(interp1(x,y, 0, ,"extrap"), cat(4, [-2 -0.5],[0 1]));   // OWrong size
+assert_checkequal(interp1(x,y,1.5, ,"extrap"), cat(4, [1 1],[1.5 2.5]));  // OWrong size
+Ref = cat(4, [-1 1 3 5 ; 0 1 2 3]',[0.5:3.5 ; 1.5:4.5]');
+assert_checkequal(interp1(x,y, 0.5:4,   ,"extrap"), Ref);
+assert_checkequal(interp1(x,y,(0.5:4)', ,"extrap"), Ref);
+
+Ref = matrix([-1 5 1 7 3 9 0 3 1 4 2 5 0.5 3.5 1.5 4.5 2.5 5.5 1.5 4.5 2.5 5.5 3.5 6.5],..
+             [2 3 2 1 2]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], ,"extrap"), Ref);
+Ref = matrix([-1 1 3 5 7 9 0 1 2 3 4 5 0.5 1.5 2.5 3.5 4.5 5.5 1.5 2.5 3.5 4.5 5.5 6.5],..
+             [1 3 2 2 1 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), ,"extrap"), Ref);
+
+// extrapolation by padding
+assert_checkequal(interp1(x,y, 0, , %pi), cat(4,[%pi %pi],[%pi %pi]));   // OWrong size and values
+assert_checkequal(interp1(x,y,1.5, , %pi), cat(4,[1 1],[1.5 2.5]));
+Ref = cat(4,[%pi 1 3 5 ; %pi 1 2 3]', [%pi 1.5:3.5 ; %pi 2.5:4.5]');
+assert_checkequal(interp1(x,y, 0.5:4, , %pi), Ref);
+assert_checkequal(interp1(x,y,(0.5:4)', , %pi), Ref);
+Ref = matrix([%pi 5 1 %pi 3 %pi %pi 3 1 %pi 2 %pi %pi 3.5 1.5 %pi ..
+   2.5 %pi %pi 4.5 2.5 %pi 3.5 %pi],[2 3 2 1 2]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5;3.5 4.5 5.5], , %pi), Ref);
+Ref = matrix([%pi 1 3 5 %pi %pi %pi 1 2 3 %pi %pi %pi 1.5 2.5 3.5 ..
+              %pi %pi %pi 2.5 3.5 4.5 %pi %pi],[1 3 2 2 1 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), , %pi), Ref);
+
+// periodic extrapolation
+y = cat(4,[0 2 4 0 ; 0.5 1.5 2.5 0.5]', [1 2 3 1 ; 2 3 4 2]');
+assert_checkequal(interp1(x,y, 0, , "periodic"), cat(4,[4 2.5],[3 4]));
+assert_checkequal(interp1(x,y,1.5, , "periodic"), cat(4,[1 1],[1.5 2.5]));
+Ref = cat(4,[2 1 3 2 ; 1.5 1 2 1.5]', [2 1.5 2.5 2 ; 3 2.5 3.5 3]');
+assert_checkequal(interp1(x,y, 0.5:4, , "periodic"), Ref);
+assert_checkequal(interp1(x,y,(0.5:4)', , "periodic"), Ref);
+Ref = matrix([2 2 1 1 3 3 1.5 1.5 1 1 2 2 2 2 1.5 1.5 ..
+   2.5 2.5 3 3 2.5 2.5 3.5 3.5],[2 3 2 1 2]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], , "periodic"), Ref, 2*%eps);
+Ref = matrix([2 1 3 2 1 3 1.5 1 2 1.5 1 2 2 1.5 2.5 2 ..
+              1.5 2.5 3 2.5 3.5 3 2.5 3.5],[1 3 2 2 1 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), , "periodic"), Ref, 2*%eps);
+
+// extrapolation by edgevalue
+e = "edgevalue";
+y = cat(4,[0 2 4 6 ; 0.5 1.5 2.5 3.5]', [1:4 ; 2:5]');
+assert_checkequal(interp1(x,y, 0, , e) , cat(4,[0 0.5],[1 2]));
+assert_checkequal(interp1(x,y,1.5, , e), cat(4,[1 1],[1.5 2.5]));
+Ref = cat(4,[0 1 3 5 ; 0.5 1 2 3]', [1 1.5:3.5 ; 2 2.5:4.5]');
+assert_checkequal(interp1(x,y, 0.5:4, , e), Ref);
+assert_checkequal(interp1(x,y,(0.5:4)', , e), Ref);
+Ref = matrix([0 5 1 6 3 6 0.5 3 1 3.5 2 3.5 1 3.5 1.5 4 ..
+   2.5 4 2 4.5 2.5 5 3.5 5],[2 3 2 1 2]);
+assert_checkequal(interp1(x, y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], , e), Ref);
+Ref = matrix([0 1 3 5 6 6 0.5 1 2 3 3.5 3.5 1 1.5 2.5 3.5 ..
+              4 4 2 2.5 3.5 4.5 5 5],[1 3 2 2 1 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), , e), Ref);
+
+// method = "nearest"
+// ------------------
+m = "nearest";
+// extrapolation "by_nan" (default)
+assert_checkequal(interp1(x,y, [0 5.4], m), %nan*ones(2,2,1,2));
+assert_checkequal(interp1(x,y,1.5, m), cat(4,[2 1.5],[2 3]));
+Ref = cat(4, [%nan %nan ; 2 1.5 ; 4 2.5 ; 6 3.5], [%nan %nan ; 2 3 ; 3 4 ; 4 5]);
+assert_checkequal(interp1(x,y,[0.5:4], m), Ref);
+assert_checkequal(interp1(x,y,[0.5:4]', m), Ref);
+Ref = matrix([%nan 6 2 %nan 4 %nan %nan 3.5 1.5 %nan 2.5 %nan %nan 4 2 %nan 3 %nan %nan 5 3 %nan 4 %nan],[2 3 2 1 2]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m), Ref);           // OErr
+Ref = matrix([%nan 2 4 6 %nan %nan %nan 1.5 2.5 3.5 %nan %nan %nan 2 3 4 %nan %nan %nan 3 4 5 %nan%nan],[1 3 2 2 1 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m), Ref);  // OErr
+
+// extrapolation "extrap" = "nearest"
+for e = ["extrap" "nearest" "edgevalue"]
+    assert_checkequal(interp1(x,y, [0 5.4], m, e), cat(4,[0 0.5 ; 6 3.5],[1 2 ; 4 5]));
+    assert_checkequal(interp1(x,y,1.5, m, e), cat(4,[2 1.5],[2 3]));         // OWrong values
+    assert_checkequal(interp1(x,y,[0.5:4], m, e), y);
+    assert_checkequal(interp1(x,y,[0.5:4]', m, e), y);
+    Ref = matrix([0 6 2 6 4 6 0.5 3.5 1.5 3.5 2.5 3.5 1 4 2 4 3 4 2 5 3 5 4 5], [2 3 2 1 2]);
+    assert_checkequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m, e), Ref);
+    Ref = matrix([0 2 4 6 6 6 0.5 1.5 2.5 3.5 3.5 3.5 1 2 3 4 4 4 2 3 4 5 5 5],[1 3 2 2 1 2]);
+    assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, e), Ref);
+end
+
+// extrapolation by padding
+assert_checkequal(interp1(x,y, [0 5.4], m, %pi), %pi*ones(2,2,1,2));
+assert_checkequal(interp1(x,y,1.5, m, %pi), cat(4,[2 1.5],[2 3]));
+Ref = cat(4, [%pi %pi ; 2 1.5 ; 4 2.5 ; 6 3.5], [%pi %pi ; 2 3 ; 3 4 ; 4 5]);
+assert_checkequal(interp1(x,y,[0.5:4], m, %pi), Ref);
+assert_checkequal(interp1(x,y,[0.5:4]', m, %pi), Ref);
+Ref = matrix([%pi 6 2 %pi 4 %pi %pi 3.5 1.5 %pi 2.5 %pi %pi 4 2 %pi 3 %pi %pi 5 3 %pi 4 %pi],[2 3 2 1 2]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m, %pi), Ref);          // OErr
+Ref = matrix([%pi 2 4 6 %pi %pi %pi 1.5 2.5 3.5 %pi %pi %pi 2 3 4 %pi %pi %pi 3 4 5 %pi%pi],[1 3 2 2 1 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, %pi), Ref);  // OErr
+
+// method = "spline"
+// -----------------
+m = "spline";
+// extrapolation default == "spline" (=="natural")
+Ref = cat(4, [-2 -0.5 ; 8.8 4.9], [0 1 ; 5.4 6.4]);
+assert_checkalmostequal(interp1(x,y, [0 5.4], m), Ref, 5*%eps);
+assert_checkequal(interp1(x,y,1.5, m), cat(4,[1 1],[1.5 2.5]));              // OWrong values
+Ref = cat(4,[-1 1 3 5 ; 0 1 2 3]',[0.5:3.5 ; 1.5:4.5]');
+assert_checkequal(interp1(x,y,[0.5:4], m), Ref);
+assert_checkequal(interp1(x,y,[0.5:4]', m), Ref);
+Ref = matrix([-1 5 1 7 3 9 0 3 1 4 2 5 0.5 3.5 1.5 4.5 2.5 5.5 1.5 4.5 2.5 5.5 3.5 6.5],..
+             [2 3 2 1 2]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m), Ref, 5*%eps);
+Ref = matrix([-1 1 3 5 7 9 0 1 2 3 4 5 0.5 1.5 2.5 3.5 4.5 5.5 1.5 2.5 3.5 4.5 5.5 6.5],..
+             [1 3 2 2 1 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m), Ref, 5*%eps);
+
+// extrapolation "extrap" == "spline"
+Ref = cat(4, [-2 -0.5 ; 8.8 4.9], [0 1 ; 5.4 6.4]);
+assert_checkalmostequal(interp1(x,y, [0 5.4], m, "extrap"), Ref, 5*%eps);    // OWrong size
+assert_checkequal(interp1(x,y,1.5, m, "extrap"), cat(4, [1 1], [1.5 2.5]));  // OWrong size
+Ref = cat(4, [-1 1 3 5; 0:3]', [0.5:3.5 ; 1.5:4.5]');
+assert_checkequal(interp1(x,y,[0.5:4], m, "extrap"), Ref);
+assert_checkequal(interp1(x,y,[0.5:4]', m, "extrap"), Ref);
+Ref = matrix([-1 5 1 7 3 9 0 3 1 4 2 5 0.5 3.5 1.5 4.5 2.5 5.5 1.5 4.5 2.5 5.5 3.5 6.5],..
+             [2 3 2 1 2]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m, "extrap"), Ref, 5*%eps); // OWrong size
+Ref = matrix([-1 1 3 5 7 9 0 1 2 3 4 5 0.5 1.5 2.5 3.5 4.5 5.5 1.5 2.5 3.5 4.5 5.5 6.5],..
+             [1 3 2 2 1 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, "extrap"), Ref, 5*%eps);
+
+// extrapolation by padding
+assert_checkalmostequal(interp1(x,y, [0 5.4], m, %pi), %pi*ones(2,2,1,2), 5*%eps);
+assert_checkequal(interp1(x,y,1.5, m, %pi), cat(4,[1 1],[1.5 2.5]));         // OWrong values
+Ref = cat(4,[%pi 1 3 5 ; %pi 1 2 3]',[%pi 1.5 2.5 3.5 ; %pi 2.5 3.5 4.5]');
+assert_checkequal(interp1(x,y,[0.5:4], m, %pi), Ref);
+assert_checkequal(interp1(x,y,[0.5:4]', m, %pi), Ref);
+Ref = matrix([%pi 5 1 %pi 3 %pi %pi 3 1 %pi 2 %pi %pi 3.5 1.5 %pi 2.5 %pi %pi 4.5 2.5 %pi 3.5 %pi], [2 3 2 1 2]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m, %pi), Ref, 5*%eps);
+Ref = matrix([%pi 1 3 5 %pi %pi %pi 1 2 3 %pi %pi %pi 1.5 2.5 3.5 %pi %pi %pi 2.5 3.5 4.5 %pi %pi], [1 3 2 2 1 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, %pi), Ref, 5*%eps);
+
+// extrapolation by edgevalue
+e = "edgevalue";
+assert_checkequal(interp1(x,y, [0 5.4], m, e), cat(4, [0 0.5 ; 6 3.5],[1 2 ; 4 5]));
+assert_checkequal(interp1(x,y,1.5, m, e), cat(4,[1 1],[1.5 2.5]));
+Ref = cat(4, [0 1 3 5 ; 0.5 1 2 3]', [1 1.5 2.5 3.5 ; 2 2.5 3.5 4.5]');
+assert_checkequal(interp1(x,y,[0.5:4], m, e), Ref);
+assert_checkequal(interp1(x,y,[0.5:4]', m, e), Ref);
+Ref = matrix([0 5 1 6 3 6 0.5 3 1 3.5 2 3.5 1 3.5 1.5 4 2.5 4 2 4.5 2.5 5 3.5 5], [2 3 2 1 2]);
+assert_checkequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m, e), Ref);
+Ref = matrix([0 1 3 5 6 6 0.5 1 2 3 3.5 3.5 1 1.5 2.5 3.5 4 4 2 2.5 3.5 4.5 5 5], [1 3 2 2 1 2]);
+assert_checkequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, e), Ref);
+
+// periodic extrapolation
+y = cat(4,[0 2 4 0 ; 0.5 1.5 2.5 0.5]', [1 2 3 1 ; 2 3 4 2]');
+Ref = cat(4, [4 2.5 ; 3.136 2.068], [3 4 ; 2.568 3.568]);
+assert_checkalmostequal(interp1(x,y, [0 5.4], m, "periodic"), Ref, 5*%eps);
+assert_checkequal(interp1(x,y,1.5, m, "periodic"), cat(4,[0.625 0.8125],[1.3125 2.3125]));
+Ref = cat(4, [3.125 2.0625 ; 0.625 0.8125 ; 3.375 2.1875 ; 3.125 2.0625], ..
+             [2.5625 3.5625; 1.3125 2.3125 ; 2.6875 3.6875 ; 2.5625 3.5625]);
+assert_checkequal(interp1(x,y,[0.5:4], m, "periodic"), Ref);
+assert_checkequal(interp1(x,y,[0.5:4]', m, "periodic"), Ref);
+Ref = matrix([3.125 3.125 0.625 0.625 3.375 3.375 2.0625 2.0625 0.8125 0.8125 ..
+      2.1875 2.1875 2.5625 2.5625 1.3125 1.3125 2.6875 2.6875 3.5625 3.5625 ..
+      2.3125 2.3125 3.6875 3.6875], [2 3 2 1 2]);
+assert_checkalmostequal(interp1(x,y,[0.5 1.5 2.5 ; 3.5 4.5 5.5], m, "periodic"), Ref, 5*%eps);
+Ref = matrix([3.125 0.625 3.375 3.125 0.625 3.375 2.0625 0.8125 2.1875 2.0625 ..
+      0.8125 2.1875 2.5625 1.3125 2.6875 2.5625 1.3125 2.6875 3.5625 2.3125 ..
+      3.6875 3.5625 2.3125  3.6875], [1 3 2 2 1 2]);
+assert_checkalmostequal(interp1(x, y, cat(3,[0.5 1.5 2.5],[3.5 4.5 5.5]), m, "periodic"), Ref, 5*%eps);
+
+
+// y is complex
+// ============
+// method = "linear" (default)
+// ---------------------------
+assert_checkequal(interp1(0:3, (1:4)+%i, 0.5:2.5), (1.5:3.5)+%i);
+// extrapolation linear
+assert_checkequal(interp1(0:3,(1:4)+%i, 0.5:2:5, , "extrap"), [1.5 3.5 5.5]+%i);
+// Default padding with %nan
+assert_checkequal(interp1((1:4)+%i, 0.5:2.5), [%nan, 1.5+%i, 2.5+%i]);
+// Padding real
+assert_checkequal(interp1(1:3, (2:4)+%i, 0.5:2.5,,7), [7, 2.5+%i, 3.5+%i]);
+// Padding complex
+assert_checkequal(interp1(1:3, (2:4)+%i, [0.5 2.5 4.5],,7-%i), [7-%i, 3.5+%i, 7-%i]);
+
+// method = "nearest"
+// ------------------
+assert_checkequal(interp1(0:3, (1:4)+%i, 0.5:2.5, "nearest"), (2:4)+%i);
+// Default padding with %nan
+assert_checkequal(interp1((1:4)+%i, 0.5:2.5, "nearest"), [%nan, 2+%i, 3+%i]);
+// extrapolation with nearest
+assert_checkequal(interp1((1:4)+%i, 0.5:2:5, "nearest", "extrap"), [1 3 4]+%i);
+// Padding real
+assert_checkequal(interp1(1:3,(2:4)+%i, 0.5:2.5, "nearest", 7), [7, 3+%i, 4+%i]);
+// Padding complex
+assert_checkequal(interp1(1:3,(2:4)+%i, [0.5 2.5 4.5], "nearest",7-%i), [7-%i, 4+%i, 7-%i]);
+
+// method = "spline"
+// -----------------
+assert_checkequal(interp1(0:3, (1:4)+%i, 0.5:2.5, "spline"), (1.5:3.5)+%i);
+// Default padding with spline
+assert_checkequal(interp1((1:4)+%i, 0.5:2.5, "spline"), (0.5:2.5)+%i);
+assert_checkequal(interp1((1:4)+%i, 0.5:2.5, "spline", "extrap"), (0.5:2.5)+%i);
+// Padding real
+assert_checkequal(interp1(1:3,(2:4)+%i, 0.5:2.5, "spline", 7), [7, 2.5+%i, 3.5+%i]);
+// Padding complex
+assert_checkequal(interp1(1:3,(2:4)+%i, [0.5 2.5 4.5], "spline",7-%i), [7-%i, 3.5+%i, 7-%i]);
+
+
+// ==============
+// Error messages
+// ==============
+msg = _("%s: Wrong number of input arguments: Must be between %d and %d.\n")
+msg = msprintf(msg, "interp1", 2, 5);
+assert_checkerror("interp1()", msg);
+assert_checkerror("interp1(1)", msg);
+assert_checkerror("interp1(1,2,3,4,5,6)", msg);
+// x
+msg = msprintf(_("%s: Argument #%d: Real numbers expected.\n"),"interp1",1);
+assert_checkerror("interp1([1 %i],[1 2],1.5)", msg);
+msg = msprintf(_("%s: Argument #%d: Vector expected.\n"), "interp1", 1);
+assert_checkerror("interp1([1 2; 3 4],[1 2],1.5)", msg);
+assert_checkerror("interp1(1,[1 2],1.5)", msg);
+msg = msprintf(_("%s: Argument #%d: Nan value forbidden.\n"), "interp1", 1);
+assert_checkerror("interp1([1 2 %nan],[1 2],1.5)", msg);
+// y
+msg = msprintf(_("%s: Argument #%d: vector or at least 2 rows expected.\n"), "interp1", 2);
+assert_checkerror("interp1(1,1.5)", msg);
+msg = msprintf(_("%s: Argument #%d: Decimal or complex numbers expected.\n"), "interp1",1);
+assert_checkerror("interp1([1 %z],1.5)", msg);
+assert_checkerror("interp1([%t %f],1.5)", msg);
+//
+msg = msprintf(_("%s: Arguments #%d and #%d: Same numbers of elements expected.\n"),"interp1", 1, 2);
+assert_checkerror("interp1([1 2],1:3,3)", msg);
+assert_checkerror("interp1(1:3,1:2,3)", msg);
+msg = msprintf(_("%s: Arguments #%d and #%d: Same numbers of rows expected.\n"), "interp1", 1, 2);
+assert_checkerror("interp1(1:3,rand(2,3),3)", msg);
+    // periodicity
+msg = msprintf(_("%s: Argument #%d: periodicity y($)==y(1) expected.\n"), "interp1", 2);
+assert_checkerror("interp1(1:4,1:4,3,,""periodic"")", msg);
+msg = msprintf(_("%s: Argument #%d: periodicity y($)==y(1) expected.\n"), "interp1", 1);
+assert_checkerror("interp1(1:4,3,""linear"",""periodic"")", msg);
+msg = msprintf(_("%s: Argument #%d: periodicity y($,:)==y(1,:) expected.\n"), "interp1", 2);
+assert_checkerror("interp1(1:3,rand(3,2),3,,""periodic"")", msg);
+msg = msprintf(_("%s: Argument #%d: periodicity y($,:)==y(1,:) expected.\n"), "interp1", 1);
+assert_checkerror("interp1(rand(3,2),3,""linear"",""periodic"")", msg);
+
+// xp
+msg = msprintf(_("%s: Argument #%d: Decimal numbers expected.\n"), "interp1",2);
+assert_checkerror("interp1(1:3,[%i 1 1])", msg);
+msg = msprintf(_("%s: Argument #%d: Decimal numbers expected.\n"), "interp1",3);
+assert_checkerror("interp1(1:3,rand(3,2),%i)", msg);
+assert_checkerror("interp1(1:3,2:4,,""linear"")", msg);
+// method
+msg = msprintf(_("%s: Argument #%d: Scalar (1 element) expected.\n"), "interp1", 3);
+assert_checkerror("interp1(1:3,2,[""linear"" ""spline""])", msg);
+msg = msprintf(_("%s: Argument #%d: Scalar (1 element) expected.\n"), "interp1", 4);
+assert_checkerror("interp1(1:3,2:4,5,[""linear"" ""spline""])", msg);
+Set = """linear"",""spline"",""nearest""";
+msg = msprintf(_("%s: Argument #%d: Must be in the set {%s}.\n"), "interp1", 3, Set);
+assert_checkerror("interp1(1:3,5,""splines"")", msg);
+msg = msprintf(_("%s: Argument #%d: Must be in the set {%s}.\n"), "interp1", 4, Set);
+assert_checkerror("interp1(2:4,1:3,5,""splines"")", msg);
+assert_checkerror("interp1(2:4,1:3,5,%z)", msg);
+// extrapolation
+msg0 = _("%s: Argument #%d: scalar number or one of {%s} expected.\n");
+msg = msprintf(msg0, "interp1", 5, """extrap"",""linear"",""periodic"",""edgevalue""");
+assert_checkerror("interp1(1:3,2:4,5,,""abc"")", msg);
+assert_checkerror("interp1(1:3,2:4,5,,%z)", msg);
+assert_checkerror("interp1(1:3,2:4,5,,1:2)", msg);
+assert_checkerror("interp1(1:3,2:4,5,,int8(3))", msg);
+msg = msprintf(msg0, "interp1", 5, """extrap"",""edgevalue"",""nearest""");
+assert_checkerror("interp1(1:3,2:4,5,""nearest"",""linear"")", msg);
+assert_checkerror("interp1(1:3,2:4,5,""nearest"",""periodic"")", msg);
+msg = msprintf(msg0, "interp1", 4, """extrap"",""linear"",""periodic"",""edgevalue""");
+assert_checkerror("interp1(1:3,2:4,""linear"",""abc"")", msg);
+assert_checkerror("interp1(1:3,2:4,""linear"",%z)", msg);
+assert_checkerror("interp1(1:3,2:4,""linear"",1:2)", msg);
+assert_checkerror("interp1(1:3,2:4,""linear"",int8(3))", msg);
+msg = msprintf(msg0, "interp1", 4, """extrap"",""edgevalue"",""nearest""");
+assert_checkerror("interp1(1:3,2:4,""nearest"",""linear"")", msg);
+assert_checkerror("interp1(1:3,2:4,""nearest"",""periodic"")", msg);