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