001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    package org.apache.commons.math3.analysis.interpolation;
018    
019    import org.apache.commons.math3.analysis.TrivariateFunction;
020    import org.apache.commons.math3.exception.DimensionMismatchException;
021    import org.apache.commons.math3.exception.NoDataException;
022    import org.apache.commons.math3.exception.OutOfRangeException;
023    import org.apache.commons.math3.exception.NonMonotonicSequenceException;
024    import org.apache.commons.math3.util.MathArrays;
025    
026    /**
027     * Function that implements the
028     * <a href="http://en.wikipedia.org/wiki/Tricubic_interpolation">
029     * tricubic spline interpolation</a>, as proposed in
030     * <quote>
031     *  Tricubic interpolation in three dimensions<br/>
032     *  F. Lekien and J. Marsden<br/>
033     *  <em>Int. J. Numer. Meth. Engng</em> 2005; <b>63</b>:455-471
034     * </quote>
035     *
036     * @since 2.2
037     * @version $Id: TricubicSplineInterpolatingFunction.java 1385314 2012-09-16 16:35:49Z tn $
038     */
039    public class TricubicSplineInterpolatingFunction
040        implements TrivariateFunction {
041        /**
042         * Matrix to compute the spline coefficients from the function values
043         * and function derivatives values
044         */
045        private static final double[][] AINV = {
046            { 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 },
047            { 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 },
048            { -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 },
049            { 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 },
050            { 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 },
051            { 0,0,0,0,0,0,0,0,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 },
052            { 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 },
053            { 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 },
054            { -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 },
055            { 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 },
056            { 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 },
057            { -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 },
058            { 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 },
059            { 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 },
060            { -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 },
061            { 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 },
062            { 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 },
063            { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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 },
064            { 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 },
065            { 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 },
066            { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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 },
067            { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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 },
068            { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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 },
069            { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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 },
070            { 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 },
071            { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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 },
072            { 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 },
073            { 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 },
074            { 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 },
075            { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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 },
076            { 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 },
077            { 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 },
078            {-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 },
079            { 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 },
080            { 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 },
081            { -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 },
082            { 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 },
083            { 0,0,0,0,0,0,0,0,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 },
084            { 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 },
085            { 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 },
086            { 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 },
087            { 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 },
088            { -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 },
089            { 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 },
090            { -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 },
091            { 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 },
092            { 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 },
093            { -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 },
094            { 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 },
095            { 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 },
096            { -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 },
097            { 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 },
098            { 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 },
099            { 0,0,0,0,0,0,0,0,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    }