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    *      http://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  package org.apache.commons.math3.analysis.interpolation;
18  
19  import org.apache.commons.math3.analysis.TrivariateFunction;
20  import org.apache.commons.math3.exception.DimensionMismatchException;
21  import org.apache.commons.math3.exception.NoDataException;
22  import org.apache.commons.math3.exception.OutOfRangeException;
23  import org.apache.commons.math3.exception.NonMonotonicSequenceException;
24  import org.apache.commons.math3.util.MathArrays;
25  
26  /**
27   * Function that implements the
28   * <a href="http://en.wikipedia.org/wiki/Tricubic_interpolation">
29   * tricubic spline interpolation</a>, as proposed in
30   * <quote>
31   *  Tricubic interpolation in three dimensions<br/>
32   *  F. Lekien and J. Marsden<br/>
33   *  <em>Int. J. Numer. Meth. Engng</em> 2005; <b>63</b>:455-471
34   * </quote>
35   *
36   * @since 2.2
37   */
38  public class TricubicSplineInterpolatingFunction
39      implements TrivariateFunction {
40      /**
41       * Matrix to compute the spline coefficients from the function values
42       * and function derivatives values
43       */
44      private static final double[][] AINV = {
45          { 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 },
46          { 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 },
47          { -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 },
48          { 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 },
49          { 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 },
50          { 0,0,0,0,0,0,0,0,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 },
51          { 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 },
52          { 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 },
53          { -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 },
54          { 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 },
55          { 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 },
56          { -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 },
57          { 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 },
58          { 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 },
59          { -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 },
60          { 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 },
61          { 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 },
62          { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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 },
63          { 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 },
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,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 },
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,0,0,0,0,0,0,0,0,1,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,0,0,0,0,0,0,0,0,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 },
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,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 },
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,2,-2,0,0,0,0,0,0,1,1,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,-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 },
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,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,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,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 },
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,-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 },
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,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 },
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,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 },
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,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,-2,-1,-2,-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,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 },
77          {-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 },
78          { 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 },
79          { 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 },
80          { -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 },
81          { 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 },
82          { 0,0,0,0,0,0,0,0,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 },
83          { 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 },
84          { 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 },
85          { 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 },
86          { 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 },
87          { -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 },
88          { 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 },
89          { -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 },
90          { 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 },
91          { 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 },
92          { -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 },
93          { 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 },
94          { 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 },
95          { -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 },
96          { 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 },
97          { 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 },
98          { 0,0,0,0,0,0,0,0,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 },
99          { 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 },
100         { 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 },
101         { -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 },
102         { 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 },
103         { 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 },
104         { -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 },
105         { 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 },
106         { 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 },
107         { -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 },
108         { 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 }
109     };
110 
111     /** Samples x-coordinates */
112     private final double[] xval;
113     /** Samples y-coordinates */
114     private final double[] yval;
115     /** Samples z-coordinates */
116     private final double[] zval;
117     /** Set of cubic splines pacthing the whole data grid */
118     private final TricubicSplineFunction[][][] splines;
119 
120     /**
121      * @param x Sample values of the x-coordinate, in increasing order.
122      * @param y Sample values of the y-coordinate, in increasing order.
123      * @param z Sample values of the y-coordinate, in increasing order.
124      * @param f Values of the function on every grid point.
125      * @param dFdX Values of the partial derivative of function with respect to x on every grid point.
126      * @param dFdY Values of the partial derivative of function with respect to y on every grid point.
127      * @param dFdZ Values of the partial derivative of function with respect to z on every grid point.
128      * @param d2FdXdY Values of the cross partial derivative of function on every grid point.
129      * @param d2FdXdZ Values of the cross partial derivative of function on every grid point.
130      * @param d2FdYdZ Values of the cross partial derivative of function on every grid point.
131      * @param d3FdXdYdZ Values of the cross partial derivative of function on every grid point.
132      * @throws NoDataException if any of the arrays has zero length.
133      * @throws DimensionMismatchException if the various arrays do not contain the expected number of elements.
134      * @throws NonMonotonicSequenceException if {@code x}, {@code y} or {@code z} are not strictly increasing.
135      */
136     public TricubicSplineInterpolatingFunction(double[] x,
137                                                double[] y,
138                                                double[] z,
139                                                double[][][] f,
140                                                double[][][] dFdX,
141                                                double[][][] dFdY,
142                                                double[][][] dFdZ,
143                                                double[][][] d2FdXdY,
144                                                double[][][] d2FdXdZ,
145                                                double[][][] d2FdYdZ,
146                                                double[][][] d3FdXdYdZ)
147         throws NoDataException,
148                DimensionMismatchException,
149                NonMonotonicSequenceException {
150         final int xLen = x.length;
151         final int yLen = y.length;
152         final int zLen = z.length;
153 
154         if (xLen == 0 || yLen == 0 || z.length == 0 || f.length == 0 || f[0].length == 0) {
155             throw new NoDataException();
156         }
157         if (xLen != f.length) {
158             throw new DimensionMismatchException(xLen, f.length);
159         }
160         if (xLen != dFdX.length) {
161             throw new DimensionMismatchException(xLen, dFdX.length);
162         }
163         if (xLen != dFdY.length) {
164             throw new DimensionMismatchException(xLen, dFdY.length);
165         }
166         if (xLen != dFdZ.length) {
167             throw new DimensionMismatchException(xLen, dFdZ.length);
168         }
169         if (xLen != d2FdXdY.length) {
170             throw new DimensionMismatchException(xLen, d2FdXdY.length);
171         }
172         if (xLen != d2FdXdZ.length) {
173             throw new DimensionMismatchException(xLen, d2FdXdZ.length);
174         }
175         if (xLen != d2FdYdZ.length) {
176             throw new DimensionMismatchException(xLen, d2FdYdZ.length);
177         }
178         if (xLen != d3FdXdYdZ.length) {
179             throw new DimensionMismatchException(xLen, d3FdXdYdZ.length);
180         }
181 
182         MathArrays.checkOrder(x);
183         MathArrays.checkOrder(y);
184         MathArrays.checkOrder(z);
185 
186         xval = x.clone();
187         yval = y.clone();
188         zval = z.clone();
189 
190         final int lastI = xLen - 1;
191         final int lastJ = yLen - 1;
192         final int lastK = zLen - 1;
193         splines = new TricubicSplineFunction[lastI][lastJ][lastK];
194 
195         for (int i = 0; i < lastI; i++) {
196             if (f[i].length != yLen) {
197                 throw new DimensionMismatchException(f[i].length, yLen);
198             }
199             if (dFdX[i].length != yLen) {
200                 throw new DimensionMismatchException(dFdX[i].length, yLen);
201             }
202             if (dFdY[i].length != yLen) {
203                 throw new DimensionMismatchException(dFdY[i].length, yLen);
204             }
205             if (dFdZ[i].length != yLen) {
206                 throw new DimensionMismatchException(dFdZ[i].length, yLen);
207             }
208             if (d2FdXdY[i].length != yLen) {
209                 throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
210             }
211             if (d2FdXdZ[i].length != yLen) {
212                 throw new DimensionMismatchException(d2FdXdZ[i].length, yLen);
213             }
214             if (d2FdYdZ[i].length != yLen) {
215                 throw new DimensionMismatchException(d2FdYdZ[i].length, yLen);
216             }
217             if (d3FdXdYdZ[i].length != yLen) {
218                 throw new DimensionMismatchException(d3FdXdYdZ[i].length, yLen);
219             }
220 
221             final int ip1 = i + 1;
222             for (int j = 0; j < lastJ; j++) {
223                 if (f[i][j].length != zLen) {
224                     throw new DimensionMismatchException(f[i][j].length, zLen);
225                 }
226                 if (dFdX[i][j].length != zLen) {
227                     throw new DimensionMismatchException(dFdX[i][j].length, zLen);
228                 }
229                 if (dFdY[i][j].length != zLen) {
230                     throw new DimensionMismatchException(dFdY[i][j].length, zLen);
231                 }
232                 if (dFdZ[i][j].length != zLen) {
233                     throw new DimensionMismatchException(dFdZ[i][j].length, zLen);
234                 }
235                 if (d2FdXdY[i][j].length != zLen) {
236                     throw new DimensionMismatchException(d2FdXdY[i][j].length, zLen);
237                 }
238                 if (d2FdXdZ[i][j].length != zLen) {
239                     throw new DimensionMismatchException(d2FdXdZ[i][j].length, zLen);
240                 }
241                 if (d2FdYdZ[i][j].length != zLen) {
242                     throw new DimensionMismatchException(d2FdYdZ[i][j].length, zLen);
243                 }
244                 if (d3FdXdYdZ[i][j].length != zLen) {
245                     throw new DimensionMismatchException(d3FdXdYdZ[i][j].length, zLen);
246                 }
247 
248                 final int jp1 = j + 1;
249                 for (int k = 0; k < lastK; k++) {
250                     final int kp1 = k + 1;
251 
252                     final double[] beta = new double[] {
253                         f[i][j][k], f[ip1][j][k],
254                         f[i][jp1][k], f[ip1][jp1][k],
255                         f[i][j][kp1], f[ip1][j][kp1],
256                         f[i][jp1][kp1], f[ip1][jp1][kp1],
257 
258                         dFdX[i][j][k], dFdX[ip1][j][k],
259                         dFdX[i][jp1][k], dFdX[ip1][jp1][k],
260                         dFdX[i][j][kp1], dFdX[ip1][j][kp1],
261                         dFdX[i][jp1][kp1], dFdX[ip1][jp1][kp1],
262 
263                         dFdY[i][j][k], dFdY[ip1][j][k],
264                         dFdY[i][jp1][k], dFdY[ip1][jp1][k],
265                         dFdY[i][j][kp1], dFdY[ip1][j][kp1],
266                         dFdY[i][jp1][kp1], dFdY[ip1][jp1][kp1],
267 
268                         dFdZ[i][j][k], dFdZ[ip1][j][k],
269                         dFdZ[i][jp1][k], dFdZ[ip1][jp1][k],
270                         dFdZ[i][j][kp1], dFdZ[ip1][j][kp1],
271                         dFdZ[i][jp1][kp1], dFdZ[ip1][jp1][kp1],
272 
273                         d2FdXdY[i][j][k], d2FdXdY[ip1][j][k],
274                         d2FdXdY[i][jp1][k], d2FdXdY[ip1][jp1][k],
275                         d2FdXdY[i][j][kp1], d2FdXdY[ip1][j][kp1],
276                         d2FdXdY[i][jp1][kp1], d2FdXdY[ip1][jp1][kp1],
277 
278                         d2FdXdZ[i][j][k], d2FdXdZ[ip1][j][k],
279                         d2FdXdZ[i][jp1][k], d2FdXdZ[ip1][jp1][k],
280                         d2FdXdZ[i][j][kp1], d2FdXdZ[ip1][j][kp1],
281                         d2FdXdZ[i][jp1][kp1], d2FdXdZ[ip1][jp1][kp1],
282 
283                         d2FdYdZ[i][j][k], d2FdYdZ[ip1][j][k],
284                         d2FdYdZ[i][jp1][k], d2FdYdZ[ip1][jp1][k],
285                         d2FdYdZ[i][j][kp1], d2FdYdZ[ip1][j][kp1],
286                         d2FdYdZ[i][jp1][kp1], d2FdYdZ[ip1][jp1][kp1],
287 
288                         d3FdXdYdZ[i][j][k], d3FdXdYdZ[ip1][j][k],
289                         d3FdXdYdZ[i][jp1][k], d3FdXdYdZ[ip1][jp1][k],
290                         d3FdXdYdZ[i][j][kp1], d3FdXdYdZ[ip1][j][kp1],
291                         d3FdXdYdZ[i][jp1][kp1], d3FdXdYdZ[ip1][jp1][kp1],
292                     };
293 
294                     splines[i][j][k] = new TricubicSplineFunction(computeSplineCoefficients(beta));
295                 }
296             }
297         }
298     }
299 
300     /**
301      * {@inheritDoc}
302      *
303      * @throws OutOfRangeException if any of the variables is outside its interpolation range.
304      */
305     public double value(double x, double y, double z)
306         throws OutOfRangeException {
307         final int i = searchIndex(x, xval);
308         if (i == -1) {
309             throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]);
310         }
311         final int j = searchIndex(y, yval);
312         if (j == -1) {
313             throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]);
314         }
315         final int k = searchIndex(z, zval);
316         if (k == -1) {
317             throw new OutOfRangeException(z, zval[0], zval[zval.length - 1]);
318         }
319 
320         final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
321         final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
322         final double zN = (z - zval[k]) / (zval[k + 1] - zval[k]);
323 
324         return splines[i][j][k].value(xN, yN, zN);
325     }
326 
327     /**
328      * @param c Coordinate.
329      * @param val Coordinate samples.
330      * @return the index in {@code val} corresponding to the interval containing {@code c}, or {@code -1}
331      *   if {@code c} is out of the range defined by the end values of {@code val}.
332      */
333     private int searchIndex(double c, double[] val) {
334         if (c < val[0]) {
335             return -1;
336         }
337 
338         final int max = val.length;
339         for (int i = 1; i < max; i++) {
340             if (c <= val[i]) {
341                 return i - 1;
342             }
343         }
344 
345         return -1;
346     }
347 
348     /**
349      * Compute the spline coefficients from the list of function values and
350      * function partial derivatives values at the four corners of a grid
351      * element. They must be specified in the following order:
352      * <ul>
353      *  <li>f(0,0,0)</li>
354      *  <li>f(1,0,0)</li>
355      *  <li>f(0,1,0)</li>
356      *  <li>f(1,1,0)</li>
357      *  <li>f(0,0,1)</li>
358      *  <li>f(1,0,1)</li>
359      *  <li>f(0,1,1)</li>
360      *  <li>f(1,1,1)</li>
361      *
362      *  <li>f<sub>x</sub>(0,0,0)</li>
363      *  <li>... <em>(same order as above)</em></li>
364      *  <li>f<sub>x</sub>(1,1,1)</li>
365      *
366      *  <li>f<sub>y</sub>(0,0,0)</li>
367      *  <li>... <em>(same order as above)</em></li>
368      *  <li>f<sub>y</sub>(1,1,1)</li>
369      *
370      *  <li>f<sub>z</sub>(0,0,0)</li>
371      *  <li>... <em>(same order as above)</em></li>
372      *  <li>f<sub>z</sub>(1,1,1)</li>
373      *
374      *  <li>f<sub>xy</sub>(0,0,0)</li>
375      *  <li>... <em>(same order as above)</em></li>
376      *  <li>f<sub>xy</sub>(1,1,1)</li>
377      *
378      *  <li>f<sub>xz</sub>(0,0,0)</li>
379      *  <li>... <em>(same order as above)</em></li>
380      *  <li>f<sub>xz</sub>(1,1,1)</li>
381      *
382      *  <li>f<sub>yz</sub>(0,0,0)</li>
383      *  <li>... <em>(same order as above)</em></li>
384      *  <li>f<sub>yz</sub>(1,1,1)</li>
385      *
386      *  <li>f<sub>xyz</sub>(0,0,0)</li>
387      *  <li>... <em>(same order as above)</em></li>
388      *  <li>f<sub>xyz</sub>(1,1,1)</li>
389      * </ul>
390      * where the subscripts indicate the partial derivative with respect to
391      * the corresponding variable(s).
392      *
393      * @param beta List of function values and function partial derivatives values.
394      * @return the spline coefficients.
395      */
396     private double[] computeSplineCoefficients(double[] beta) {
397         final int sz = 64;
398         final double[] a = new double[sz];
399 
400         for (int i = 0; i < sz; i++) {
401             double result = 0;
402             final double[] row = AINV[i];
403             for (int j = 0; j < sz; j++) {
404                 result += row[j] * beta[j];
405             }
406             a[i] = result;
407         }
408 
409         return a;
410     }
411 }
412 
413 /**
414  * 3D-spline function.
415  *
416  */
417 class TricubicSplineFunction
418     implements TrivariateFunction {
419     /** Number of points. */
420     private static final short N = 4;
421     /** Coefficients */
422     private final double[][][] a = new double[N][N][N];
423 
424     /**
425      * @param aV List of spline coefficients.
426      */
427     public TricubicSplineFunction(double[] aV) {
428         for (int i = 0; i < N; i++) {
429             for (int j = 0; j < N; j++) {
430                 for (int k = 0; k < N; k++) {
431                     a[i][j][k] = aV[i + N * (j + N * k)];
432                 }
433             }
434         }
435     }
436 
437     /**
438      * @param x x-coordinate of the interpolation point.
439      * @param y y-coordinate of the interpolation point.
440      * @param z z-coordinate of the interpolation point.
441      * @return the interpolated value.
442      * @throws OutOfRangeException if {@code x}, {@code y} or
443      * {@code z} are not in the interval {@code [0, 1]}.
444      */
445     public double value(double x, double y, double z)
446         throws OutOfRangeException {
447         if (x < 0 || x > 1) {
448             throw new OutOfRangeException(x, 0, 1);
449         }
450         if (y < 0 || y > 1) {
451             throw new OutOfRangeException(y, 0, 1);
452         }
453         if (z < 0 || z > 1) {
454             throw new OutOfRangeException(z, 0, 1);
455         }
456 
457         final double x2 = x * x;
458         final double x3 = x2 * x;
459         final double[] pX = { 1, x, x2, x3 };
460 
461         final double y2 = y * y;
462         final double y3 = y2 * y;
463         final double[] pY = { 1, y, y2, y3 };
464 
465         final double z2 = z * z;
466         final double z3 = z2 * z;
467         final double[] pZ = { 1, z, z2, z3 };
468 
469         double result = 0;
470         for (int i = 0; i < N; i++) {
471             for (int j = 0; j < N; j++) {
472                 for (int k = 0; k < N; k++) {
473                     result += a[i][j][k] * pX[i] * pY[j] * pZ[k];
474                 }
475             }
476         }
477 
478         return result;
479     }
480 }