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