TricubicInterpolatingFunction.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      https://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.analysis.interpolation;

  24. import org.hipparchus.analysis.TrivariateFunction;
  25. import org.hipparchus.exception.LocalizedCoreFormats;
  26. import org.hipparchus.exception.MathIllegalArgumentException;
  27. import org.hipparchus.util.MathArrays;
  28. import org.hipparchus.util.MathUtils;

  30. /**
  31.  * Function that implements the
  32.  * <a href="http://en.wikipedia.org/wiki/Tricubic_interpolation">
  33.  * tricubic spline interpolation</a>, as proposed in
  34.  * <blockquote>
  35.  *  Tricubic interpolation in three dimensions<br>
  36.  *  F. Lekien and J. Marsden<br>
  37.  *  <em>Int. J. Numer. Meth. Eng</em> 2005; <b>63</b>:455-471<br>
  38.  * </blockquote>
  39.  *
  40.  */
  41. public class TricubicInterpolatingFunction
  42.     implements TrivariateFunction {
  43.     /**
  44.      * Matrix to compute the spline coefficients from the function values
  45.      * and function derivatives values
  46.      */
  47.     private static final double[][] AINV = {
  48.         { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  49.         { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  50.         { -3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  51.         { 2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  52.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  53.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  54.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  55.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  56.         { -3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  57.         { 0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  58.         { 9,-9,-9,9,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,0,0,0,0,0,0,0,0,4,2,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  59.         { -6,6,6,-6,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,0,0,0,0,0,0,0,0,-2,-2,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  60.         { 2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  61.         { 0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  62.         { -6,6,6,-6,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  63.         { 4,-4,-4,4,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  64.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  65.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  66.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  67.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  68.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  69.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
  70.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0 },
  71.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0 },
  72.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  73.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
  74.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,-9,9,0,0,0,0,0,0,0,0,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,4,2,2,1,0,0,0,0 },
  75.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,-2,-2,-1,-1,0,0,0,0 },
  76.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  77.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0 },
  78.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,-2,-1,-2,-1,0,0,0,0 },
  79.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,-4,4,0,0,0,0,0,0,0,0,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,1,1,1,1,0,0,0,0 },
  80.         {-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  81.         { 0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  82.         { 9,-9,0,0,-9,9,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,0,0,0,0,0,0,0,0,4,2,0,0,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  83.         { -6,6,0,0,6,-6,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,0,0,0,0,0,0,0,0,0,0,-2,-2,0,0,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  84.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0 },
  85.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0 },
  86.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,0,0,-9,9,0,0,0,0,0,0,0,0,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,4,2,0,0,2,1,0,0 },
  87.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,0,0,-2,-2,0,0,-1,-1,0,0 },
  88.         { 9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0,0,0,0,0,0,0,0,0 },
  89.         { 0,0,0,0,0,0,0,0,9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0 },
  90.         { -27,27,27,-27,27,-27,-27,27,-18,-9,18,9,18,9,-18,-9,-18,18,-9,9,18,-18,9,-9,-18,18,18,-18,-9,9,9,-9,-12,-6,-6,-3,12,6,6,3,-12,-6,12,6,-6,-3,6,3,-12,12,-6,6,-6,6,-3,3,-8,-4,-4,-2,-4,-2,-2,-1 },
  91.         { 18,-18,-18,18,-18,18,18,-18,9,9,-9,-9,-9,-9,9,9,12,-12,6,-6,-12,12,-6,6,12,-12,-12,12,6,-6,-6,6,6,6,3,3,-6,-6,-3,-3,6,6,-6,-6,3,3,-3,-3,8,-8,4,-4,4,-4,2,-2,4,4,2,2,2,2,1,1 },
  92.         { -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0,0,0,0,0,0,0,0,0 },
  93.         { 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0 },
  94.         { 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,9,-9,9,-9,-9,9,-9,9,12,-12,-12,12,6,-6,-6,6,6,3,6,3,-6,-3,-6,-3,8,4,-8,-4,4,2,-4,-2,6,-6,6,-6,3,-3,3,-3,4,2,4,2,2,1,2,1 },
  95.         { -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-6,6,-6,6,6,-6,6,-6,-8,8,8,-8,-4,4,4,-4,-3,-3,-3,-3,3,3,3,3,-4,-4,4,4,-2,-2,2,2,-4,4,-4,4,-2,2,-2,2,-2,-2,-2,-2,-1,-1,-1,-1 },
  96.         { 2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  97.         { 0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  98.         { -6,6,0,0,6,-6,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  99.         { 4,-4,0,0,-4,4,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  100.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
  101.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0 },
  102.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,-2,-1,0,0,-2,-1,0,0 },
  103.         { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,0,0,-4,4,0,0,0,0,0,0,0,0,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,1,1,0,0,1,1,0,0 },
  104.         { -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0,0,0,0,0,0,0,0,0 },
  105.         { 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0 },
  106.         { 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,12,-12,6,-6,-12,12,-6,6,9,-9,-9,9,9,-9,-9,9,8,4,4,2,-8,-4,-4,-2,6,3,-6,-3,6,3,-6,-3,6,-6,3,-3,6,-6,3,-3,4,2,2,1,4,2,2,1 },
  107.         { -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-8,8,-4,4,8,-8,4,-4,-6,6,6,-6,-6,6,6,-6,-4,-4,-2,-2,4,4,2,2,-3,-3,3,3,-3,-3,3,3,-4,4,-2,2,-4,4,-2,2,-2,-2,-1,-1,-2,-2,-1,-1 },
  108.         { 4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0 },
  109.         { 0,0,0,0,0,0,0,0,4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0 },
  110.         { -12,12,12,-12,12,-12,-12,12,-8,-4,8,4,8,4,-8,-4,-6,6,-6,6,6,-6,6,-6,-6,6,6,-6,-6,6,6,-6,-4,-2,-4,-2,4,2,4,2,-4,-2,4,2,-4,-2,4,2,-3,3,-3,3,-3,3,-3,3,-2,-1,-2,-1,-2,-1,-2,-1 },
  111.         { 8,-8,-8,8,-8,8,8,-8,4,4,-4,-4,-4,-4,4,4,4,-4,4,-4,-4,4,-4,4,4,-4,-4,4,4,-4,-4,4,2,2,2,2,-2,-2,-2,-2,2,2,-2,-2,2,2,-2,-2,2,-2,2,-2,2,-2,2,-2,1,1,1,1,1,1,1,1 }
  112.     };

  114.     /** Samples x-coordinates */
  115.     private final double[] xval;
  116.     /** Samples y-coordinates */
  117.     private final double[] yval;
  118.     /** Samples z-coordinates */
  119.     private final double[] zval;
  120.     /** Set of cubic splines pacthing the whole data grid */
  121.     private final TricubicFunction[][][] splines;

  123.     /** Simple constructor.
  124.      * @param x Sample values of the x-coordinate, in increasing order.
  125.      * @param y Sample values of the y-coordinate, in increasing order.
  126.      * @param z Sample values of the y-coordinate, in increasing order.
  127.      * @param f Values of the function on every grid point.
  128.      * @param dFdX Values of the partial derivative of function with respect to x on every grid point.
  129.      * @param dFdY Values of the partial derivative of function with respect to y on every grid point.
  130.      * @param dFdZ Values of the partial derivative of function with respect to z on every grid point.
  131.      * @param d2FdXdY Values of the cross partial derivative of function on every grid point.
  132.      * @param d2FdXdZ Values of the cross partial derivative of function on every grid point.
  133.      * @param d2FdYdZ Values of the cross partial derivative of function on every grid point.
  134.      * @param d3FdXdYdZ Values of the cross partial derivative of function on every grid point.
  135.      * @throws MathIllegalArgumentException if any of the arrays has zero length.
  136.      * @throws MathIllegalArgumentException if the various arrays do not contain the expected number of elements.
  137.      * @throws MathIllegalArgumentException if {@code x}, {@code y} or {@code z} are not strictly increasing.
  138.      */
  139.     public TricubicInterpolatingFunction(double[] x,
  140.                                          double[] y,
  141.                                          double[] z,
  142.                                          double[][][] f,
  143.                                          double[][][] dFdX,
  144.                                          double[][][] dFdY,
  145.                                          double[][][] dFdZ,
  146.                                          double[][][] d2FdXdY,
  147.                                          double[][][] d2FdXdZ,
  148.                                          double[][][] d2FdYdZ,
  149.                                          double[][][] d3FdXdYdZ)
  150.         throws MathIllegalArgumentException {
  151.         final int xLen = x.length;
  152.         final int yLen = y.length;
  153.         final int zLen = z.length;

  155.         if (xLen == 0 || yLen == 0 || z.length == 0 || f.length == 0 || f[0].length == 0) {
  156.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NO_DATA);
  157.         }
  158.         if (xLen != f.length) {
  159.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  160.                                                    xLen, f.length);
  161.         }
  162.         if (xLen != dFdX.length) {
  163.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  164.                                                    xLen, dFdX.length);
  165.         }
  166.         if (xLen != dFdY.length) {
  167.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  168.                                                    xLen, dFdY.length);
  169.         }
  170.         if (xLen != dFdZ.length) {
  171.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  172.                                                    xLen, dFdZ.length);
  173.         }
  174.         if (xLen != d2FdXdY.length) {
  175.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  176.                                                    xLen, d2FdXdY.length);
  177.         }
  178.         if (xLen != d2FdXdZ.length) {
  179.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  180.                                                    xLen, d2FdXdZ.length);
  181.         }
  182.         if (xLen != d2FdYdZ.length) {
  183.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  184.                                                    xLen, d2FdYdZ.length);
  185.         }
  186.         if (xLen != d3FdXdYdZ.length) {
  187.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  188.                                                    xLen, d3FdXdYdZ.length);
  189.         }

  191.         MathArrays.checkOrder(x);
  192.         MathArrays.checkOrder(y);
  193.         MathArrays.checkOrder(z);

  195.         xval = x.clone();
  196.         yval = y.clone();
  197.         zval = z.clone();

  199.         final int lastI = xLen - 1;
  200.         final int lastJ = yLen - 1;
  201.         final int lastK = zLen - 1;
  202.         splines = new TricubicFunction[lastI][lastJ][lastK];

  204.         for (int i = 0; i < lastI; i++) {
  205.             if (f[i].length != yLen) {
  206.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  207.                                                        f[i].length, yLen);
  208.             }
  209.             if (dFdX[i].length != yLen) {
  210.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  211.                                                        dFdX[i].length, yLen);
  212.             }
  213.             if (dFdY[i].length != yLen) {
  214.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  215.                                                        dFdY[i].length, yLen);
  216.             }
  217.             if (dFdZ[i].length != yLen) {
  218.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  219.                                                        dFdZ[i].length, yLen);
  220.             }
  221.             if (d2FdXdY[i].length != yLen) {
  222.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  223.                                                        d2FdXdY[i].length, yLen);
  224.             }
  225.             if (d2FdXdZ[i].length != yLen) {
  226.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  227.                                                        d2FdXdZ[i].length, yLen);
  228.             }
  229.             if (d2FdYdZ[i].length != yLen) {
  230.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  231.                                                        d2FdYdZ[i].length, yLen);
  232.             }
  233.             if (d3FdXdYdZ[i].length != yLen) {
  234.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  235.                                                        d3FdXdYdZ[i].length, yLen);
  236.             }

  238.             final int ip1 = i + 1;
  239.             final double xR = xval[ip1] - xval[i];
  240.             for (int j = 0; j < lastJ; j++) {
  241.                 if (f[i][j].length != zLen) {
  242.                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  243.                                                            f[i][j].length, zLen);
  244.                 }
  245.                 if (dFdX[i][j].length != zLen) {
  246.                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  247.                                                            dFdX[i][j].length, zLen);
  248.                 }
  249.                 if (dFdY[i][j].length != zLen) {
  250.                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  251.                                                            dFdY[i][j].length, zLen);
  252.                 }
  253.                 if (dFdZ[i][j].length != zLen) {
  254.                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  255.                                                            dFdZ[i][j].length, zLen);
  256.                 }
  257.                 if (d2FdXdY[i][j].length != zLen) {
  258.                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  259.                                                            d2FdXdY[i][j].length, zLen);
  260.                 }
  261.                 if (d2FdXdZ[i][j].length != zLen) {
  262.                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  263.                                                            d2FdXdZ[i][j].length, zLen);
  264.                 }
  265.                 if (d2FdYdZ[i][j].length != zLen) {
  266.                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  267.                                                            d2FdYdZ[i][j].length, zLen);
  268.                 }
  269.                 if (d3FdXdYdZ[i][j].length != zLen) {
  270.                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  271.                                                            d3FdXdYdZ[i][j].length, zLen);
  272.                 }

  274.                 final int jp1 = j + 1;
  275.                 final double yR = yval[jp1] - yval[j];
  276.                 final double xRyR = xR * yR;
  277.                 for (int k = 0; k < lastK; k++) {
  278.                     final int kp1 = k + 1;
  279.                     final double zR = zval[kp1] - zval[k];
  280.                     final double xRzR = xR * zR;
  281.                     final double yRzR = yR * zR;
  282.                     final double xRyRzR = xR * yRzR;

  284.                     final double[] beta = {
  285.                         f[i][j][k], f[ip1][j][k],
  286.                         f[i][jp1][k], f[ip1][jp1][k],
  287.                         f[i][j][kp1], f[ip1][j][kp1],
  288.                         f[i][jp1][kp1], f[ip1][jp1][kp1],

  290.                         dFdX[i][j][k] * xR, dFdX[ip1][j][k] * xR,
  291.                         dFdX[i][jp1][k] * xR, dFdX[ip1][jp1][k] * xR,
  292.                         dFdX[i][j][kp1] * xR, dFdX[ip1][j][kp1] * xR,
  293.                         dFdX[i][jp1][kp1] * xR, dFdX[ip1][jp1][kp1] * xR,

  295.                         dFdY[i][j][k] * yR, dFdY[ip1][j][k] * yR,
  296.                         dFdY[i][jp1][k] * yR, dFdY[ip1][jp1][k] * yR,
  297.                         dFdY[i][j][kp1] * yR, dFdY[ip1][j][kp1] * yR,
  298.                         dFdY[i][jp1][kp1] * yR, dFdY[ip1][jp1][kp1] * yR,

  300.                         dFdZ[i][j][k] * zR, dFdZ[ip1][j][k] * zR,
  301.                         dFdZ[i][jp1][k] * zR, dFdZ[ip1][jp1][k] * zR,
  302.                         dFdZ[i][j][kp1] * zR, dFdZ[ip1][j][kp1] * zR,
  303.                         dFdZ[i][jp1][kp1] * zR, dFdZ[ip1][jp1][kp1] * zR,

  305.                         d2FdXdY[i][j][k] * xRyR, d2FdXdY[ip1][j][k] * xRyR,
  306.                         d2FdXdY[i][jp1][k] * xRyR, d2FdXdY[ip1][jp1][k] * xRyR,
  307.                         d2FdXdY[i][j][kp1] * xRyR, d2FdXdY[ip1][j][kp1] * xRyR,
  308.                         d2FdXdY[i][jp1][kp1] * xRyR, d2FdXdY[ip1][jp1][kp1] * xRyR,

  310.                         d2FdXdZ[i][j][k] * xRzR, d2FdXdZ[ip1][j][k] * xRzR,
  311.                         d2FdXdZ[i][jp1][k] * xRzR, d2FdXdZ[ip1][jp1][k] * xRzR,
  312.                         d2FdXdZ[i][j][kp1] * xRzR, d2FdXdZ[ip1][j][kp1] * xRzR,
  313.                         d2FdXdZ[i][jp1][kp1] * xRzR, d2FdXdZ[ip1][jp1][kp1] * xRzR,

  315.                         d2FdYdZ[i][j][k] * yRzR, d2FdYdZ[ip1][j][k] * yRzR,
  316.                         d2FdYdZ[i][jp1][k] * yRzR, d2FdYdZ[ip1][jp1][k] * yRzR,
  317.                         d2FdYdZ[i][j][kp1] * yRzR, d2FdYdZ[ip1][j][kp1] * yRzR,
  318.                         d2FdYdZ[i][jp1][kp1] * yRzR, d2FdYdZ[ip1][jp1][kp1] * yRzR,

  320.                         d3FdXdYdZ[i][j][k] * xRyRzR, d3FdXdYdZ[ip1][j][k] * xRyRzR,
  321.                         d3FdXdYdZ[i][jp1][k] * xRyRzR, d3FdXdYdZ[ip1][jp1][k] * xRyRzR,
  322.                         d3FdXdYdZ[i][j][kp1] * xRyRzR, d3FdXdYdZ[ip1][j][kp1] * xRyRzR,
  323.                         d3FdXdYdZ[i][jp1][kp1] * xRyRzR, d3FdXdYdZ[ip1][jp1][kp1] * xRyRzR,
  324.                     };

  326.                     splines[i][j][k] = new TricubicFunction(computeCoefficients(beta));
  327.                 }
  328.             }
  329.         }
  330.     }

  332.     /**
  333.      * {@inheritDoc}
  334.      *
  335.      * @throws MathIllegalArgumentException if any of the variables is outside its interpolation range.
  336.      */
  337.     @Override
  338.     public double value(double x, double y, double z)
  339.         throws MathIllegalArgumentException {
  340.         final int i = searchIndex(x, xval);
  341.         if (i == -1) {
  342.             throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE,
  343.                                                    x, xval[0], xval[xval.length - 1]);
  344.         }
  345.         final int j = searchIndex(y, yval);
  346.         if (j == -1) {
  347.             throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE,
  348.                                                    y, yval[0], yval[yval.length - 1]);
  349.         }
  350.         final int k = searchIndex(z, zval);
  351.         if (k == -1) {
  352.             throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE,
  353.                                                    z, zval[0], zval[zval.length - 1]);
  354.         }

  356.         final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
  357.         final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
  358.         final double zN = (z - zval[k]) / (zval[k + 1] - zval[k]);

  360.         return splines[i][j][k].value(xN, yN, zN);
  361.     }

  363.     /**
  364.      * Indicates whether a point is within the interpolation range.
  365.      *
  366.      * @param x First coordinate.
  367.      * @param y Second coordinate.
  368.      * @param z Third coordinate.
  369.      * @return {@code true} if (x, y, z) is a valid point.
  370.      */
  371.     public boolean isValidPoint(double x, double y, double z) {
  372.         if (x < xval[0] ||
  373.             x > xval[xval.length - 1] ||
  374.             y < yval[0] ||
  375.             y > yval[yval.length - 1] ||
  376.             z < zval[0] ||
  377.             z > zval[zval.length - 1]) {
  378.             return false;
  379.         } else {
  380.             return true;
  381.         }
  382.     }

  384.     /**
  385.      * @param c Coordinate.
  386.      * @param val Coordinate samples.
  387.      * @return the index in {@code val} corresponding to the interval containing {@code c}, or {@code -1}
  388.      *   if {@code c} is out of the range defined by the end values of {@code val}.
  389.      */
  390.     private int searchIndex(double c, double[] val) {
  391.         if (c < val[0]) {
  392.             return -1;
  393.         }

  395.         final int max = val.length;
  396.         for (int i = 1; i < max; i++) {
  397.             if (c <= val[i]) {
  398.                 return i - 1;
  399.             }
  400.         }

  402.         return -1;
  403.     }

  405.     /**
  406.      * Compute the spline coefficients from the list of function values and
  407.      * function partial derivatives values at the four corners of a grid
  408.      * element. They must be specified in the following order:
  409.      * <ul>
  410.      *  <li>f(0,0,0)</li>
  411.      *  <li>f(1,0,0)</li>
  412.      *  <li>f(0,1,0)</li>
  413.      *  <li>f(1,1,0)</li>
  414.      *  <li>f(0,0,1)</li>
  415.      *  <li>f(1,0,1)</li>
  416.      *  <li>f(0,1,1)</li>
  417.      *  <li>f(1,1,1)</li>
  418.      *
  419.      *  <li>f<sub>x</sub>(0,0,0)</li>
  420.      *  <li>... <em>(same order as above)</em></li>
  421.      *  <li>f<sub>x</sub>(1,1,1)</li>
  422.      *
  423.      *  <li>f<sub>y</sub>(0,0,0)</li>
  424.      *  <li>... <em>(same order as above)</em></li>
  425.      *  <li>f<sub>y</sub>(1,1,1)</li>
  426.      *
  427.      *  <li>f<sub>z</sub>(0,0,0)</li>
  428.      *  <li>... <em>(same order as above)</em></li>
  429.      *  <li>f<sub>z</sub>(1,1,1)</li>
  430.      *
  431.      *  <li>f<sub>xy</sub>(0,0,0)</li>
  432.      *  <li>... <em>(same order as above)</em></li>
  433.      *  <li>f<sub>xy</sub>(1,1,1)</li>
  434.      *
  435.      *  <li>f<sub>xz</sub>(0,0,0)</li>
  436.      *  <li>... <em>(same order as above)</em></li>
  437.      *  <li>f<sub>xz</sub>(1,1,1)</li>
  438.      *
  439.      *  <li>f<sub>yz</sub>(0,0,0)</li>
  440.      *  <li>... <em>(same order as above)</em></li>
  441.      *  <li>f<sub>yz</sub>(1,1,1)</li>
  442.      *
  443.      *  <li>f<sub>xyz</sub>(0,0,0)</li>
  444.      *  <li>... <em>(same order as above)</em></li>
  445.      *  <li>f<sub>xyz</sub>(1,1,1)</li>
  446.      * </ul>
  447.      * where the subscripts indicate the partial derivative with respect to
  448.      * the corresponding variable(s).
  449.      *
  450.      * @param beta List of function values and function partial derivatives values.
  451.      * @return the spline coefficients.
  452.      */
  453.     private double[] computeCoefficients(double[] beta) {
  454.         final int sz = 64;
  455.         final double[] a = new double[sz];

  457.         for (int i = 0; i < sz; i++) {
  458.             double result = 0;
  459.             final double[] row = AINV[i];
  460.             for (int j = 0; j < sz; j++) {
  461.                 result += row[j] * beta[j];
  462.             }
  463.             a[i] = result;
  464.         }

  466.         return a;
  467.     }
  468. }

  470. /**
  471.  * 3D-spline function.
  472.  *
  473.  */
  474. class TricubicFunction
  475.     implements TrivariateFunction {
  476.     /** Number of points. */
  477.     private static final short N = 4;
  478.     /** Coefficients */
  479.     private final double[][][] a = new double[N][N][N];

  481.     /**
  482.      * @param aV List of spline coefficients.
  483.      */
  484.     TricubicFunction(double[] aV) {
  485.         for (int i = 0; i < N; i++) {
  486.             for (int j = 0; j < N; j++) {
  487.                 for (int k = 0; k < N; k++) {
  488.                     a[i][j][k] = aV[i + N * (j + N * k)];
  489.                 }
  490.             }
  491.         }
  492.     }

  494.     /**
  495.      * @param x x-coordinate of the interpolation point.
  496.      * @param y y-coordinate of the interpolation point.
  497.      * @param z z-coordinate of the interpolation point.
  498.      * @return the interpolated value.
  499.      * @throws MathIllegalArgumentException if {@code x}, {@code y} or
  500.      * {@code z} are not in the interval {@code [0, 1]}.
  501.      */
  502.     @Override
  503.     public double value(double x, double y, double z) throws MathIllegalArgumentException {
  504.         MathUtils.checkRangeInclusive(x, 0, 1);
  505.         MathUtils.checkRangeInclusive(y, 0, 1);
  506.         MathUtils.checkRangeInclusive(z, 0, 1);

  508.         final double x2 = x * x;
  509.         final double x3 = x2 * x;
  510.         final double[] pX = { 1, x, x2, x3 };

  512.         final double y2 = y * y;
  513.         final double y3 = y2 * y;
  514.         final double[] pY = { 1, y, y2, y3 };

  516.         final double z2 = z * z;
  517.         final double z3 = z2 * z;
  518.         final double[] pZ = { 1, z, z2, z3 };

  520.         double result = 0;
  521.         for (int i = 0; i < N; i++) {
  522.             for (int j = 0; j < N; j++) {
  523.                 for (int k = 0; k < N; k++) {
  524.                     result += a[i][j][k] * pX[i] * pY[j] * pZ[k];
  525.                 }
  526.             }
  527.         }

  529.         return result;
  530.     }
  531. }
Created with JaCoCo 0.8.12.202403310830