View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.analysis.interpolation;
23  
24  import org.hipparchus.analysis.TrivariateFunction;
25  import org.hipparchus.exception.LocalizedCoreFormats;
26  import org.hipparchus.exception.MathIllegalArgumentException;
27  import org.hipparchus.util.MathArrays;
28  import org.hipparchus.util.MathUtils;
29  
30  /**
31   * Function that implements the
32   * <a href="http://en.wikipedia.org/wiki/Tricubic_interpolation">
33   * tricubic spline interpolation</a>, as proposed in
34   * <blockquote>
35   *  Tricubic interpolation in three dimensions<br>
36   *  F. Lekien and J. Marsden<br>
37   *  <em>Int. J. Numer. Meth. Eng</em> 2005; <b>63</b>:455-471<br>
38   * </blockquote>
39   *
40   */
41  public class TricubicInterpolatingFunction
42      implements TrivariateFunction {
43      /**
44       * Matrix to compute the spline coefficients from the function values
45       * and function derivatives values
46       */
47      private static final double[][] AINV = {
48          { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
49          { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
50          { -3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
51          { 2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
52          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
53          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
54          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
55          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
56          { -3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
57          { 0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
58          { 9,-9,-9,9,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,0,0,0,0,0,0,0,0,4,2,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
59          { -6,6,6,-6,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,0,0,0,0,0,0,0,0,-2,-2,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
60          { 2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
61          { 0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
62          { -6,6,6,-6,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
63          { 4,-4,-4,4,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
64          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
65          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
66          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
67          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
68          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
69          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
70          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0 },
71          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0 },
72          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0 },
73          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
74          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,-9,9,0,0,0,0,0,0,0,0,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,4,2,2,1,0,0,0,0 },
75          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,-2,-2,-1,-1,0,0,0,0 },
76          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0 },
77          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0 },
78          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,-2,-1,-2,-1,0,0,0,0 },
79          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,-4,4,0,0,0,0,0,0,0,0,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,1,1,1,1,0,0,0,0 },
80          {-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
81          { 0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
82          { 9,-9,0,0,-9,9,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,0,0,0,0,0,0,0,0,4,2,0,0,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
83          { -6,6,0,0,6,-6,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,0,0,0,0,0,0,0,0,0,0,-2,-2,0,0,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
84          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0 },
85          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0 },
86          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,0,0,-9,9,0,0,0,0,0,0,0,0,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,4,2,0,0,2,1,0,0 },
87          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,0,0,-2,-2,0,0,-1,-1,0,0 },
88          { 9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0,0,0,0,0,0,0,0,0 },
89          { 0,0,0,0,0,0,0,0,9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0 },
90          { -27,27,27,-27,27,-27,-27,27,-18,-9,18,9,18,9,-18,-9,-18,18,-9,9,18,-18,9,-9,-18,18,18,-18,-9,9,9,-9,-12,-6,-6,-3,12,6,6,3,-12,-6,12,6,-6,-3,6,3,-12,12,-6,6,-6,6,-3,3,-8,-4,-4,-2,-4,-2,-2,-1 },
91          { 18,-18,-18,18,-18,18,18,-18,9,9,-9,-9,-9,-9,9,9,12,-12,6,-6,-12,12,-6,6,12,-12,-12,12,6,-6,-6,6,6,6,3,3,-6,-6,-3,-3,6,6,-6,-6,3,3,-3,-3,8,-8,4,-4,4,-4,2,-2,4,4,2,2,2,2,1,1 },
92          { -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0,0,0,0,0,0,0,0,0 },
93          { 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0 },
94          { 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,9,-9,9,-9,-9,9,-9,9,12,-12,-12,12,6,-6,-6,6,6,3,6,3,-6,-3,-6,-3,8,4,-8,-4,4,2,-4,-2,6,-6,6,-6,3,-3,3,-3,4,2,4,2,2,1,2,1 },
95          { -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-6,6,-6,6,6,-6,6,-6,-8,8,8,-8,-4,4,4,-4,-3,-3,-3,-3,3,3,3,3,-4,-4,4,4,-2,-2,2,2,-4,4,-4,4,-2,2,-2,2,-2,-2,-2,-2,-1,-1,-1,-1 },
96          { 2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
97          { 0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
98          { -6,6,0,0,6,-6,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
99          { 4,-4,0,0,-4,4,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
100         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
101         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0 },
102         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,-2,-1,0,0,-2,-1,0,0 },
103         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,0,0,-4,4,0,0,0,0,0,0,0,0,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,1,1,0,0,1,1,0,0 },
104         { -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0,0,0,0,0,0,0,0,0 },
105         { 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0 },
106         { 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,12,-12,6,-6,-12,12,-6,6,9,-9,-9,9,9,-9,-9,9,8,4,4,2,-8,-4,-4,-2,6,3,-6,-3,6,3,-6,-3,6,-6,3,-3,6,-6,3,-3,4,2,2,1,4,2,2,1 },
107         { -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-8,8,-4,4,8,-8,4,-4,-6,6,6,-6,-6,6,6,-6,-4,-4,-2,-2,4,4,2,2,-3,-3,3,3,-3,-3,3,3,-4,4,-2,2,-4,4,-2,2,-2,-2,-1,-1,-2,-2,-1,-1 },
108         { 4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0 },
109         { 0,0,0,0,0,0,0,0,4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0 },
110         { -12,12,12,-12,12,-12,-12,12,-8,-4,8,4,8,4,-8,-4,-6,6,-6,6,6,-6,6,-6,-6,6,6,-6,-6,6,6,-6,-4,-2,-4,-2,4,2,4,2,-4,-2,4,2,-4,-2,4,2,-3,3,-3,3,-3,3,-3,3,-2,-1,-2,-1,-2,-1,-2,-1 },
111         { 8,-8,-8,8,-8,8,8,-8,4,4,-4,-4,-4,-4,4,4,4,-4,4,-4,-4,4,-4,4,4,-4,-4,4,4,-4,-4,4,2,2,2,2,-2,-2,-2,-2,2,2,-2,-2,2,2,-2,-2,2,-2,2,-2,2,-2,2,-2,1,1,1,1,1,1,1,1 }
112     };
113 
114     /** Samples x-coordinates */
115     private final double[] xval;
116     /** Samples y-coordinates */
117     private final double[] yval;
118     /** Samples z-coordinates */
119     private final double[] zval;
120     /** Set of cubic splines pacthing the whole data grid */
121     private final TricubicFunction[][][] splines;
122 
123     /** Simple constructor.
124      * @param x Sample values of the x-coordinate, in increasing order.
125      * @param y Sample values of the y-coordinate, in increasing order.
126      * @param z Sample values of the y-coordinate, in increasing order.
127      * @param f Values of the function on every grid point.
128      * @param dFdX Values of the partial derivative of function with respect to x on every grid point.
129      * @param dFdY Values of the partial derivative of function with respect to y on every grid point.
130      * @param dFdZ Values of the partial derivative of function with respect to z on every grid point.
131      * @param d2FdXdY Values of the cross partial derivative of function on every grid point.
132      * @param d2FdXdZ Values of the cross partial derivative of function on every grid point.
133      * @param d2FdYdZ Values of the cross partial derivative of function on every grid point.
134      * @param d3FdXdYdZ Values of the cross partial derivative of function on every grid point.
135      * @throws MathIllegalArgumentException if any of the arrays has zero length.
136      * @throws MathIllegalArgumentException if the various arrays do not contain the expected number of elements.
137      * @throws MathIllegalArgumentException if {@code x}, {@code y} or {@code z} are not strictly increasing.
138      */
139     public TricubicInterpolatingFunction(double[] x,
140                                          double[] y,
141                                          double[] z,
142                                          double[][][] f,
143                                          double[][][] dFdX,
144                                          double[][][] dFdY,
145                                          double[][][] dFdZ,
146                                          double[][][] d2FdXdY,
147                                          double[][][] d2FdXdZ,
148                                          double[][][] d2FdYdZ,
149                                          double[][][] d3FdXdYdZ)
150         throws MathIllegalArgumentException {
151         final int xLen = x.length;
152         final int yLen = y.length;
153         final int zLen = z.length;
154 
155         if (xLen == 0 || yLen == 0 || z.length == 0 || f.length == 0 || f[0].length == 0) {
156             throw new MathIllegalArgumentException(LocalizedCoreFormats.NO_DATA);
157         }
158         if (xLen != f.length) {
159             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
160                                                    xLen, f.length);
161         }
162         if (xLen != dFdX.length) {
163             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
164                                                    xLen, dFdX.length);
165         }
166         if (xLen != dFdY.length) {
167             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
168                                                    xLen, dFdY.length);
169         }
170         if (xLen != dFdZ.length) {
171             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
172                                                    xLen, dFdZ.length);
173         }
174         if (xLen != d2FdXdY.length) {
175             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
176                                                    xLen, d2FdXdY.length);
177         }
178         if (xLen != d2FdXdZ.length) {
179             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
180                                                    xLen, d2FdXdZ.length);
181         }
182         if (xLen != d2FdYdZ.length) {
183             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
184                                                    xLen, d2FdYdZ.length);
185         }
186         if (xLen != d3FdXdYdZ.length) {
187             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
188                                                    xLen, d3FdXdYdZ.length);
189         }
190 
191         MathArrays.checkOrder(x);
192         MathArrays.checkOrder(y);
193         MathArrays.checkOrder(z);
194 
195         xval = x.clone();
196         yval = y.clone();
197         zval = z.clone();
198 
199         final int lastI = xLen - 1;
200         final int lastJ = yLen - 1;
201         final int lastK = zLen - 1;
202         splines = new TricubicFunction[lastI][lastJ][lastK];
203 
204         for (int i = 0; i < lastI; i++) {
205             if (f[i].length != yLen) {
206                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
207                                                        f[i].length, yLen);
208             }
209             if (dFdX[i].length != yLen) {
210                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
211                                                        dFdX[i].length, yLen);
212             }
213             if (dFdY[i].length != yLen) {
214                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
215                                                        dFdY[i].length, yLen);
216             }
217             if (dFdZ[i].length != yLen) {
218                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
219                                                        dFdZ[i].length, yLen);
220             }
221             if (d2FdXdY[i].length != yLen) {
222                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
223                                                        d2FdXdY[i].length, yLen);
224             }
225             if (d2FdXdZ[i].length != yLen) {
226                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
227                                                        d2FdXdZ[i].length, yLen);
228             }
229             if (d2FdYdZ[i].length != yLen) {
230                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
231                                                        d2FdYdZ[i].length, yLen);
232             }
233             if (d3FdXdYdZ[i].length != yLen) {
234                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
235                                                        d3FdXdYdZ[i].length, yLen);
236             }
237 
238             final int ip1 = i + 1;
239             final double xR = xval[ip1] - xval[i];
240             for (int j = 0; j < lastJ; j++) {
241                 if (f[i][j].length != zLen) {
242                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
243                                                            f[i][j].length, zLen);
244                 }
245                 if (dFdX[i][j].length != zLen) {
246                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
247                                                            dFdX[i][j].length, zLen);
248                 }
249                 if (dFdY[i][j].length != zLen) {
250                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
251                                                            dFdY[i][j].length, zLen);
252                 }
253                 if (dFdZ[i][j].length != zLen) {
254                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
255                                                            dFdZ[i][j].length, zLen);
256                 }
257                 if (d2FdXdY[i][j].length != zLen) {
258                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
259                                                            d2FdXdY[i][j].length, zLen);
260                 }
261                 if (d2FdXdZ[i][j].length != zLen) {
262                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
263                                                            d2FdXdZ[i][j].length, zLen);
264                 }
265                 if (d2FdYdZ[i][j].length != zLen) {
266                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
267                                                            d2FdYdZ[i][j].length, zLen);
268                 }
269                 if (d3FdXdYdZ[i][j].length != zLen) {
270                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
271                                                            d3FdXdYdZ[i][j].length, zLen);
272                 }
273 
274                 final int jp1 = j + 1;
275                 final double yR = yval[jp1] - yval[j];
276                 final double xRyR = xR * yR;
277                 for (int k = 0; k < lastK; k++) {
278                     final int kp1 = k + 1;
279                     final double zR = zval[kp1] - zval[k];
280                     final double xRzR = xR * zR;
281                     final double yRzR = yR * zR;
282                     final double xRyRzR = xR * yRzR;
283 
284                     final double[] beta = {
285                         f[i][j][k], f[ip1][j][k],
286                         f[i][jp1][k], f[ip1][jp1][k],
287                         f[i][j][kp1], f[ip1][j][kp1],
288                         f[i][jp1][kp1], f[ip1][jp1][kp1],
289 
290                         dFdX[i][j][k] * xR, dFdX[ip1][j][k] * xR,
291                         dFdX[i][jp1][k] * xR, dFdX[ip1][jp1][k] * xR,
292                         dFdX[i][j][kp1] * xR, dFdX[ip1][j][kp1] * xR,
293                         dFdX[i][jp1][kp1] * xR, dFdX[ip1][jp1][kp1] * xR,
294 
295                         dFdY[i][j][k] * yR, dFdY[ip1][j][k] * yR,
296                         dFdY[i][jp1][k] * yR, dFdY[ip1][jp1][k] * yR,
297                         dFdY[i][j][kp1] * yR, dFdY[ip1][j][kp1] * yR,
298                         dFdY[i][jp1][kp1] * yR, dFdY[ip1][jp1][kp1] * yR,
299 
300                         dFdZ[i][j][k] * zR, dFdZ[ip1][j][k] * zR,
301                         dFdZ[i][jp1][k] * zR, dFdZ[ip1][jp1][k] * zR,
302                         dFdZ[i][j][kp1] * zR, dFdZ[ip1][j][kp1] * zR,
303                         dFdZ[i][jp1][kp1] * zR, dFdZ[ip1][jp1][kp1] * zR,
304 
305                         d2FdXdY[i][j][k] * xRyR, d2FdXdY[ip1][j][k] * xRyR,
306                         d2FdXdY[i][jp1][k] * xRyR, d2FdXdY[ip1][jp1][k] * xRyR,
307                         d2FdXdY[i][j][kp1] * xRyR, d2FdXdY[ip1][j][kp1] * xRyR,
308                         d2FdXdY[i][jp1][kp1] * xRyR, d2FdXdY[ip1][jp1][kp1] * xRyR,
309 
310                         d2FdXdZ[i][j][k] * xRzR, d2FdXdZ[ip1][j][k] * xRzR,
311                         d2FdXdZ[i][jp1][k] * xRzR, d2FdXdZ[ip1][jp1][k] * xRzR,
312                         d2FdXdZ[i][j][kp1] * xRzR, d2FdXdZ[ip1][j][kp1] * xRzR,
313                         d2FdXdZ[i][jp1][kp1] * xRzR, d2FdXdZ[ip1][jp1][kp1] * xRzR,
314 
315                         d2FdYdZ[i][j][k] * yRzR, d2FdYdZ[ip1][j][k] * yRzR,
316                         d2FdYdZ[i][jp1][k] * yRzR, d2FdYdZ[ip1][jp1][k] * yRzR,
317                         d2FdYdZ[i][j][kp1] * yRzR, d2FdYdZ[ip1][j][kp1] * yRzR,
318                         d2FdYdZ[i][jp1][kp1] * yRzR, d2FdYdZ[ip1][jp1][kp1] * yRzR,
319 
320                         d3FdXdYdZ[i][j][k] * xRyRzR, d3FdXdYdZ[ip1][j][k] * xRyRzR,
321                         d3FdXdYdZ[i][jp1][k] * xRyRzR, d3FdXdYdZ[ip1][jp1][k] * xRyRzR,
322                         d3FdXdYdZ[i][j][kp1] * xRyRzR, d3FdXdYdZ[ip1][j][kp1] * xRyRzR,
323                         d3FdXdYdZ[i][jp1][kp1] * xRyRzR, d3FdXdYdZ[ip1][jp1][kp1] * xRyRzR,
324                     };
325 
326                     splines[i][j][k] = new TricubicFunction(computeCoefficients(beta));
327                 }
328             }
329         }
330     }
331 
332     /**
333      * {@inheritDoc}
334      *
335      * @throws MathIllegalArgumentException if any of the variables is outside its interpolation range.
336      */
337     @Override
338     public double value(double x, double y, double z)
339         throws MathIllegalArgumentException {
340         final int i = searchIndex(x, xval);
341         if (i == -1) {
342             throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE,
343                                                    x, xval[0], xval[xval.length - 1]);
344         }
345         final int j = searchIndex(y, yval);
346         if (j == -1) {
347             throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE,
348                                                    y, yval[0], yval[yval.length - 1]);
349         }
350         final int k = searchIndex(z, zval);
351         if (k == -1) {
352             throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE,
353                                                    z, zval[0], zval[zval.length - 1]);
354         }
355 
356         final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
357         final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
358         final double zN = (z - zval[k]) / (zval[k + 1] - zval[k]);
359 
360         return splines[i][j][k].value(xN, yN, zN);
361     }
362 
363     /**
364      * Indicates whether a point is within the interpolation range.
365      *
366      * @param x First coordinate.
367      * @param y Second coordinate.
368      * @param z Third coordinate.
369      * @return {@code true} if (x, y, z) is a valid point.
370      */
371     public boolean isValidPoint(double x, double y, double z) {
372         if (x < xval[0] ||
373             x > xval[xval.length - 1] ||
374             y < yval[0] ||
375             y > yval[yval.length - 1] ||
376             z < zval[0] ||
377             z > zval[zval.length - 1]) {
378             return false;
379         } else {
380             return true;
381         }
382     }
383 
384     /**
385      * @param c Coordinate.
386      * @param val Coordinate samples.
387      * @return the index in {@code val} corresponding to the interval containing {@code c}, or {@code -1}
388      *   if {@code c} is out of the range defined by the end values of {@code val}.
389      */
390     private int searchIndex(double c, double[] val) {
391         if (c < val[0]) {
392             return -1;
393         }
394 
395         final int max = val.length;
396         for (int i = 1; i < max; i++) {
397             if (c <= val[i]) {
398                 return i - 1;
399             }
400         }
401 
402         return -1;
403     }
404 
405     /**
406      * Compute the spline coefficients from the list of function values and
407      * function partial derivatives values at the four corners of a grid
408      * element. They must be specified in the following order:
409      * <ul>
410      *  <li>f(0,0,0)</li>
411      *  <li>f(1,0,0)</li>
412      *  <li>f(0,1,0)</li>
413      *  <li>f(1,1,0)</li>
414      *  <li>f(0,0,1)</li>
415      *  <li>f(1,0,1)</li>
416      *  <li>f(0,1,1)</li>
417      *  <li>f(1,1,1)</li>
418      *
419      *  <li>f<sub>x</sub>(0,0,0)</li>
420      *  <li>... <em>(same order as above)</em></li>
421      *  <li>f<sub>x</sub>(1,1,1)</li>
422      *
423      *  <li>f<sub>y</sub>(0,0,0)</li>
424      *  <li>... <em>(same order as above)</em></li>
425      *  <li>f<sub>y</sub>(1,1,1)</li>
426      *
427      *  <li>f<sub>z</sub>(0,0,0)</li>
428      *  <li>... <em>(same order as above)</em></li>
429      *  <li>f<sub>z</sub>(1,1,1)</li>
430      *
431      *  <li>f<sub>xy</sub>(0,0,0)</li>
432      *  <li>... <em>(same order as above)</em></li>
433      *  <li>f<sub>xy</sub>(1,1,1)</li>
434      *
435      *  <li>f<sub>xz</sub>(0,0,0)</li>
436      *  <li>... <em>(same order as above)</em></li>
437      *  <li>f<sub>xz</sub>(1,1,1)</li>
438      *
439      *  <li>f<sub>yz</sub>(0,0,0)</li>
440      *  <li>... <em>(same order as above)</em></li>
441      *  <li>f<sub>yz</sub>(1,1,1)</li>
442      *
443      *  <li>f<sub>xyz</sub>(0,0,0)</li>
444      *  <li>... <em>(same order as above)</em></li>
445      *  <li>f<sub>xyz</sub>(1,1,1)</li>
446      * </ul>
447      * where the subscripts indicate the partial derivative with respect to
448      * the corresponding variable(s).
449      *
450      * @param beta List of function values and function partial derivatives values.
451      * @return the spline coefficients.
452      */
453     private double[] computeCoefficients(double[] beta) {
454         final int sz = 64;
455         final double[] a = new double[sz];
456 
457         for (int i = 0; i < sz; i++) {
458             double result = 0;
459             final double[] row = AINV[i];
460             for (int j = 0; j < sz; j++) {
461                 result += row[j] * beta[j];
462             }
463             a[i] = result;
464         }
465 
466         return a;
467     }
468 }
469 
470 /**
471  * 3D-spline function.
472  *
473  */
474 class TricubicFunction
475     implements TrivariateFunction {
476     /** Number of points. */
477     private static final short N = 4;
478     /** Coefficients */
479     private final double[][][] a = new double[N][N][N];
480 
481     /**
482      * @param aV List of spline coefficients.
483      */
484     TricubicFunction(double[] aV) {
485         for (int i = 0; i < N; i++) {
486             for (int j = 0; j < N; j++) {
487                 for (int k = 0; k < N; k++) {
488                     a[i][j][k] = aV[i + N * (j + N * k)];
489                 }
490             }
491         }
492     }
493 
494     /**
495      * @param x x-coordinate of the interpolation point.
496      * @param y y-coordinate of the interpolation point.
497      * @param z z-coordinate of the interpolation point.
498      * @return the interpolated value.
499      * @throws MathIllegalArgumentException if {@code x}, {@code y} or
500      * {@code z} are not in the interval {@code [0, 1]}.
501      */
502     @Override
503     public double value(double x, double y, double z) throws MathIllegalArgumentException {
504         MathUtils.checkRangeInclusive(x, 0, 1);
505         MathUtils.checkRangeInclusive(y, 0, 1);
506         MathUtils.checkRangeInclusive(z, 0, 1);
507 
508         final double x2 = x * x;
509         final double x3 = x2 * x;
510         final double[] pX = { 1, x, x2, x3 };
511 
512         final double y2 = y * y;
513         final double y3 = y2 * y;
514         final double[] pY = { 1, y, y2, y3 };
515 
516         final double z2 = z * z;
517         final double z3 = z2 * z;
518         final double[] pZ = { 1, z, z2, z3 };
519 
520         double result = 0;
521         for (int i = 0; i < N; i++) {
522             for (int j = 0; j < N; j++) {
523                 for (int k = 0; k < N; k++) {
524                     result += a[i][j][k] * pX[i] * pY[j] * pZ[k];
525                 }
526             }
527         }
528 
529         return result;
530     }
531 }