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 */
017package org.apache.commons.math3.analysis.interpolation;
018
019import org.apache.commons.math3.analysis.TrivariateFunction;
020import org.apache.commons.math3.exception.DimensionMismatchException;
021import org.apache.commons.math3.exception.NoDataException;
022import org.apache.commons.math3.exception.NonMonotonicSequenceException;
023import org.apache.commons.math3.exception.OutOfRangeException;
024import 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 * <blockquote>
031 *  Tricubic interpolation in three dimensions,
032 *  F. Lekien and J. Marsden,
033 *  <em>Int. J. Numer. Meth. Engng</em> 2005; <b>63</b>:455-471
034 * </blockquote>
035 *
036 * @since 2.2
037 * @deprecated To be removed in 4.0 (see MATH-1166).
038 */
039@Deprecated
040public class TricubicSplineInterpolatingFunction
041    implements TrivariateFunction {
042    /**
043     * Matrix to compute the spline coefficients from the function values
044     * and function derivatives values
045     */
046    private static final double[][] AINV = {
047        { 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 },
048        { 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 },
049        { -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 },
050        { 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 },
051        { 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 },
052        { 0,0,0,0,0,0,0,0,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 },
053        { 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 },
054        { 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 },
055        { -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 },
056        { 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 },
057        { 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 },
058        { -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 },
059        { 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 },
060        { 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 },
061        { -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 },
062        { 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 },
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,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 },
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,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 },
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,-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 },
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,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 },
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,1,0,0,0,0,0,0,0,0,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,0,0,0,0,0,0,0,0,1,0,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,-3,3,0,0,0,0,0,0,-2,-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,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 },
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,-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 },
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,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 },
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,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 },
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,-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 },
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,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 },
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,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 },
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,-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 },
078        { 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 },
079        {-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 },
080        { 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 },
081        { 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 },
082        { -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 },
083        { 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 },
084        { 0,0,0,0,0,0,0,0,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 },
085        { 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 },
086        { 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 },
087        { 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 },
088        { 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 },
089        { -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 },
090        { 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 },
091        { -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 },
092        { 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 },
093        { 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 },
094        { -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 },
095        { 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 },
096        { 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 },
097        { -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 },
098        { 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 },
099        { 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 },
100        { 0,0,0,0,0,0,0,0,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 },
101        { 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 },
102        { 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 },
103        { -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 },
104        { 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 },
105        { 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 },
106        { -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 },
107        { 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 },
108        { 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 },
109        { -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 },
110        { 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 }
111    };
112
113    /** Samples x-coordinates */
114    private final double[] xval;
115    /** Samples y-coordinates */
116    private final double[] yval;
117    /** Samples z-coordinates */
118    private final double[] zval;
119    /** Set of cubic splines pacthing the whole data grid */
120    private final TricubicSplineFunction[][][] splines;
121
122    /**
123     * @param x Sample values of the x-coordinate, in increasing order.
124     * @param y Sample values of the y-coordinate, in increasing order.
125     * @param z Sample values of the y-coordinate, in increasing order.
126     * @param f Values of the function on every grid point.
127     * @param dFdX Values of the partial derivative of function with respect to x on every grid point.
128     * @param dFdY Values of the partial derivative of function with respect to y on every grid point.
129     * @param dFdZ Values of the partial derivative of function with respect to z on every grid point.
130     * @param d2FdXdY Values of the cross partial derivative of function on every grid point.
131     * @param d2FdXdZ Values of the cross partial derivative of function on every grid point.
132     * @param d2FdYdZ Values of the cross partial derivative of function on every grid point.
133     * @param d3FdXdYdZ Values of the cross partial derivative of function on every grid point.
134     * @throws NoDataException if any of the arrays has zero length.
135     * @throws DimensionMismatchException if the various arrays do not contain the expected number of elements.
136     * @throws NonMonotonicSequenceException if {@code x}, {@code y} or {@code z} are not strictly increasing.
137     */
138    public TricubicSplineInterpolatingFunction(double[] x,
139                                               double[] y,
140                                               double[] z,
141                                               double[][][] f,
142                                               double[][][] dFdX,
143                                               double[][][] dFdY,
144                                               double[][][] dFdZ,
145                                               double[][][] d2FdXdY,
146                                               double[][][] d2FdXdZ,
147                                               double[][][] d2FdYdZ,
148                                               double[][][] d3FdXdYdZ)
149        throws NoDataException,
150               DimensionMismatchException,
151               NonMonotonicSequenceException {
152        final int xLen = x.length;
153        final int yLen = y.length;
154        final int zLen = z.length;
155
156        if (xLen == 0 || yLen == 0 || z.length == 0 || f.length == 0 || f[0].length == 0) {
157            throw new NoDataException();
158        }
159        if (xLen != f.length) {
160            throw new DimensionMismatchException(xLen, f.length);
161        }
162        if (xLen != dFdX.length) {
163            throw new DimensionMismatchException(xLen, dFdX.length);
164        }
165        if (xLen != dFdY.length) {
166            throw new DimensionMismatchException(xLen, dFdY.length);
167        }
168        if (xLen != dFdZ.length) {
169            throw new DimensionMismatchException(xLen, dFdZ.length);
170        }
171        if (xLen != d2FdXdY.length) {
172            throw new DimensionMismatchException(xLen, d2FdXdY.length);
173        }
174        if (xLen != d2FdXdZ.length) {
175            throw new DimensionMismatchException(xLen, d2FdXdZ.length);
176        }
177        if (xLen != d2FdYdZ.length) {
178            throw new DimensionMismatchException(xLen, d2FdYdZ.length);
179        }
180        if (xLen != d3FdXdYdZ.length) {
181            throw new DimensionMismatchException(xLen, d3FdXdYdZ.length);
182        }
183
184        MathArrays.checkOrder(x);
185        MathArrays.checkOrder(y);
186        MathArrays.checkOrder(z);
187
188        xval = x.clone();
189        yval = y.clone();
190        zval = z.clone();
191
192        final int lastI = xLen - 1;
193        final int lastJ = yLen - 1;
194        final int lastK = zLen - 1;
195        splines = new TricubicSplineFunction[lastI][lastJ][lastK];
196
197        for (int i = 0; i < lastI; i++) {
198            if (f[i].length != yLen) {
199                throw new DimensionMismatchException(f[i].length, yLen);
200            }
201            if (dFdX[i].length != yLen) {
202                throw new DimensionMismatchException(dFdX[i].length, yLen);
203            }
204            if (dFdY[i].length != yLen) {
205                throw new DimensionMismatchException(dFdY[i].length, yLen);
206            }
207            if (dFdZ[i].length != yLen) {
208                throw new DimensionMismatchException(dFdZ[i].length, yLen);
209            }
210            if (d2FdXdY[i].length != yLen) {
211                throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
212            }
213            if (d2FdXdZ[i].length != yLen) {
214                throw new DimensionMismatchException(d2FdXdZ[i].length, yLen);
215            }
216            if (d2FdYdZ[i].length != yLen) {
217                throw new DimensionMismatchException(d2FdYdZ[i].length, yLen);
218            }
219            if (d3FdXdYdZ[i].length != yLen) {
220                throw new DimensionMismatchException(d3FdXdYdZ[i].length, yLen);
221            }
222
223            final int ip1 = i + 1;
224            for (int j = 0; j < lastJ; j++) {
225                if (f[i][j].length != zLen) {
226                    throw new DimensionMismatchException(f[i][j].length, zLen);
227                }
228                if (dFdX[i][j].length != zLen) {
229                    throw new DimensionMismatchException(dFdX[i][j].length, zLen);
230                }
231                if (dFdY[i][j].length != zLen) {
232                    throw new DimensionMismatchException(dFdY[i][j].length, zLen);
233                }
234                if (dFdZ[i][j].length != zLen) {
235                    throw new DimensionMismatchException(dFdZ[i][j].length, zLen);
236                }
237                if (d2FdXdY[i][j].length != zLen) {
238                    throw new DimensionMismatchException(d2FdXdY[i][j].length, zLen);
239                }
240                if (d2FdXdZ[i][j].length != zLen) {
241                    throw new DimensionMismatchException(d2FdXdZ[i][j].length, zLen);
242                }
243                if (d2FdYdZ[i][j].length != zLen) {
244                    throw new DimensionMismatchException(d2FdYdZ[i][j].length, zLen);
245                }
246                if (d3FdXdYdZ[i][j].length != zLen) {
247                    throw new DimensionMismatchException(d3FdXdYdZ[i][j].length, zLen);
248                }
249
250                final int jp1 = j + 1;
251                for (int k = 0; k < lastK; k++) {
252                    final int kp1 = k + 1;
253
254                    final double[] beta = new double[] {
255                        f[i][j][k], f[ip1][j][k],
256                        f[i][jp1][k], f[ip1][jp1][k],
257                        f[i][j][kp1], f[ip1][j][kp1],
258                        f[i][jp1][kp1], f[ip1][jp1][kp1],
259
260                        dFdX[i][j][k], dFdX[ip1][j][k],
261                        dFdX[i][jp1][k], dFdX[ip1][jp1][k],
262                        dFdX[i][j][kp1], dFdX[ip1][j][kp1],
263                        dFdX[i][jp1][kp1], dFdX[ip1][jp1][kp1],
264
265                        dFdY[i][j][k], dFdY[ip1][j][k],
266                        dFdY[i][jp1][k], dFdY[ip1][jp1][k],
267                        dFdY[i][j][kp1], dFdY[ip1][j][kp1],
268                        dFdY[i][jp1][kp1], dFdY[ip1][jp1][kp1],
269
270                        dFdZ[i][j][k], dFdZ[ip1][j][k],
271                        dFdZ[i][jp1][k], dFdZ[ip1][jp1][k],
272                        dFdZ[i][j][kp1], dFdZ[ip1][j][kp1],
273                        dFdZ[i][jp1][kp1], dFdZ[ip1][jp1][kp1],
274
275                        d2FdXdY[i][j][k], d2FdXdY[ip1][j][k],
276                        d2FdXdY[i][jp1][k], d2FdXdY[ip1][jp1][k],
277                        d2FdXdY[i][j][kp1], d2FdXdY[ip1][j][kp1],
278                        d2FdXdY[i][jp1][kp1], d2FdXdY[ip1][jp1][kp1],
279
280                        d2FdXdZ[i][j][k], d2FdXdZ[ip1][j][k],
281                        d2FdXdZ[i][jp1][k], d2FdXdZ[ip1][jp1][k],
282                        d2FdXdZ[i][j][kp1], d2FdXdZ[ip1][j][kp1],
283                        d2FdXdZ[i][jp1][kp1], d2FdXdZ[ip1][jp1][kp1],
284
285                        d2FdYdZ[i][j][k], d2FdYdZ[ip1][j][k],
286                        d2FdYdZ[i][jp1][k], d2FdYdZ[ip1][jp1][k],
287                        d2FdYdZ[i][j][kp1], d2FdYdZ[ip1][j][kp1],
288                        d2FdYdZ[i][jp1][kp1], d2FdYdZ[ip1][jp1][kp1],
289
290                        d3FdXdYdZ[i][j][k], d3FdXdYdZ[ip1][j][k],
291                        d3FdXdYdZ[i][jp1][k], d3FdXdYdZ[ip1][jp1][k],
292                        d3FdXdYdZ[i][j][kp1], d3FdXdYdZ[ip1][j][kp1],
293                        d3FdXdYdZ[i][jp1][kp1], d3FdXdYdZ[ip1][jp1][kp1],
294                    };
295
296                    splines[i][j][k] = new TricubicSplineFunction(computeSplineCoefficients(beta));
297                }
298            }
299        }
300    }
301
302    /**
303     * {@inheritDoc}
304     *
305     * @throws OutOfRangeException if any of the variables is outside its interpolation range.
306     */
307    public double value(double x, double y, double z)
308        throws OutOfRangeException {
309        final int i = searchIndex(x, xval);
310        if (i == -1) {
311            throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]);
312        }
313        final int j = searchIndex(y, yval);
314        if (j == -1) {
315            throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]);
316        }
317        final int k = searchIndex(z, zval);
318        if (k == -1) {
319            throw new OutOfRangeException(z, zval[0], zval[zval.length - 1]);
320        }
321
322        final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
323        final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
324        final double zN = (z - zval[k]) / (zval[k + 1] - zval[k]);
325
326        return splines[i][j][k].value(xN, yN, zN);
327    }
328
329    /**
330     * @param c Coordinate.
331     * @param val Coordinate samples.
332     * @return the index in {@code val} corresponding to the interval containing {@code c}, or {@code -1}
333     *   if {@code c} is out of the range defined by the end values of {@code val}.
334     */
335    private int searchIndex(double c, double[] val) {
336        if (c < val[0]) {
337            return -1;
338        }
339
340        final int max = val.length;
341        for (int i = 1; i < max; i++) {
342            if (c <= val[i]) {
343                return i - 1;
344            }
345        }
346
347        return -1;
348    }
349
350    /**
351     * Compute the spline coefficients from the list of function values and
352     * function partial derivatives values at the four corners of a grid
353     * element. They must be specified in the following order:
354     * <ul>
355     *  <li>f(0,0,0)</li>
356     *  <li>f(1,0,0)</li>
357     *  <li>f(0,1,0)</li>
358     *  <li>f(1,1,0)</li>
359     *  <li>f(0,0,1)</li>
360     *  <li>f(1,0,1)</li>
361     *  <li>f(0,1,1)</li>
362     *  <li>f(1,1,1)</li>
363     *
364     *  <li>f<sub>x</sub>(0,0,0)</li>
365     *  <li>... <em>(same order as above)</em></li>
366     *  <li>f<sub>x</sub>(1,1,1)</li>
367     *
368     *  <li>f<sub>y</sub>(0,0,0)</li>
369     *  <li>... <em>(same order as above)</em></li>
370     *  <li>f<sub>y</sub>(1,1,1)</li>
371     *
372     *  <li>f<sub>z</sub>(0,0,0)</li>
373     *  <li>... <em>(same order as above)</em></li>
374     *  <li>f<sub>z</sub>(1,1,1)</li>
375     *
376     *  <li>f<sub>xy</sub>(0,0,0)</li>
377     *  <li>... <em>(same order as above)</em></li>
378     *  <li>f<sub>xy</sub>(1,1,1)</li>
379     *
380     *  <li>f<sub>xz</sub>(0,0,0)</li>
381     *  <li>... <em>(same order as above)</em></li>
382     *  <li>f<sub>xz</sub>(1,1,1)</li>
383     *
384     *  <li>f<sub>yz</sub>(0,0,0)</li>
385     *  <li>... <em>(same order as above)</em></li>
386     *  <li>f<sub>yz</sub>(1,1,1)</li>
387     *
388     *  <li>f<sub>xyz</sub>(0,0,0)</li>
389     *  <li>... <em>(same order as above)</em></li>
390     *  <li>f<sub>xyz</sub>(1,1,1)</li>
391     * </ul>
392     * where the subscripts indicate the partial derivative with respect to
393     * the corresponding variable(s).
394     *
395     * @param beta List of function values and function partial derivatives values.
396     * @return the spline coefficients.
397     */
398    private double[] computeSplineCoefficients(double[] beta) {
399        final int sz = 64;
400        final double[] a = new double[sz];
401
402        for (int i = 0; i < sz; i++) {
403            double result = 0;
404            final double[] row = AINV[i];
405            for (int j = 0; j < sz; j++) {
406                result += row[j] * beta[j];
407            }
408            a[i] = result;
409        }
410
411        return a;
412    }
413}
414
415/**
416 * 3D-spline function.
417 *
418 */
419class 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    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}