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. Eng</em> 2005; <b>63</b>:455-471 034 * </blockquote> 035 * 036 * @since 3.4. 037 */ 038public class TricubicInterpolatingFunction 039 implements TrivariateFunction { 040 /** 041 * Matrix to compute the spline coefficients from the function values 042 * and function derivatives values 043 */ 044 private static final double[][] AINV = { 045 { 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 }, 046 { 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 }, 047 { -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 }, 048 { 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 }, 049 { 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 }, 050 { 0,0,0,0,0,0,0,0,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 }, 051 { 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 }, 052 { 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 }, 053 { -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 }, 054 { 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 }, 055 { 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 }, 056 { -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 }, 057 { 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 }, 058 { 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 }, 059 { -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 }, 060 { 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 }, 061 { 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 }, 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,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 }, 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,-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 }, 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,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 }, 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,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 }, 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,0,0,0,0,0,0,0,0,1,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,-3,3,0,0,0,0,0,0,-2,-1,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,2,-2,0,0,0,0,0,0,1,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,-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 }, 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,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,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,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 }, 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,-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 }, 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,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 }, 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,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 }, 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,-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 }, 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,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 }, 077 {-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 }, 078 { 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 }, 079 { 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 }, 080 { -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 }, 081 { 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 }, 082 { 0,0,0,0,0,0,0,0,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 }, 083 { 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 }, 084 { 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 }, 085 { 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 }, 086 { 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 }, 087 { -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 }, 088 { 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 }, 089 { -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 }, 090 { 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 }, 091 { 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 }, 092 { -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 }, 093 { 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 }, 094 { 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 }, 095 { -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 }, 096 { 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 }, 097 { 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 }, 098 { 0,0,0,0,0,0,0,0,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 }, 099 { 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 patching the whole data grid */ 118 private final TricubicFunction[][][] 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 to x on every grid point. 126 * @param dFdY Values of the partial derivative of function with respect to y on every grid point. 127 * @param dFdZ Values of the partial derivative of function with respect to z on every grid point. 128 * @param d2FdXdY Values of the cross partial derivative of function on every grid point. 129 * @param d2FdXdZ Values of the cross partial derivative of function on every grid point. 130 * @param d2FdYdZ Values of the cross partial derivative of function on every grid point. 131 * @param d3FdXdYdZ Values of the cross partial derivative of function on every grid point. 132 * @throws NoDataException if any of the arrays has zero length. 133 * @throws DimensionMismatchException if the various arrays do not contain the expected number of elements. 134 * @throws NonMonotonicSequenceException if {@code x}, {@code y} or {@code z} are not strictly increasing. 135 */ 136 public TricubicInterpolatingFunction(double[] x, 137 double[] y, 138 double[] z, 139 double[][][] f, 140 double[][][] dFdX, 141 double[][][] dFdY, 142 double[][][] dFdZ, 143 double[][][] d2FdXdY, 144 double[][][] d2FdXdZ, 145 double[][][] d2FdYdZ, 146 double[][][] d3FdXdYdZ) 147 throws NoDataException, 148 DimensionMismatchException, 149 NonMonotonicSequenceException { 150 final int xLen = x.length; 151 final int yLen = y.length; 152 final int zLen = z.length; 153 154 if (xLen == 0 || yLen == 0 || z.length == 0 || f.length == 0 || f[0].length == 0) { 155 throw new NoDataException(); 156 } 157 if (xLen != f.length) { 158 throw new DimensionMismatchException(xLen, f.length); 159 } 160 if (xLen != dFdX.length) { 161 throw new DimensionMismatchException(xLen, dFdX.length); 162 } 163 if (xLen != dFdY.length) { 164 throw new DimensionMismatchException(xLen, dFdY.length); 165 } 166 if (xLen != dFdZ.length) { 167 throw new DimensionMismatchException(xLen, dFdZ.length); 168 } 169 if (xLen != d2FdXdY.length) { 170 throw new DimensionMismatchException(xLen, d2FdXdY.length); 171 } 172 if (xLen != d2FdXdZ.length) { 173 throw new DimensionMismatchException(xLen, d2FdXdZ.length); 174 } 175 if (xLen != d2FdYdZ.length) { 176 throw new DimensionMismatchException(xLen, d2FdYdZ.length); 177 } 178 if (xLen != d3FdXdYdZ.length) { 179 throw new DimensionMismatchException(xLen, d3FdXdYdZ.length); 180 } 181 182 MathArrays.checkOrder(x); 183 MathArrays.checkOrder(y); 184 MathArrays.checkOrder(z); 185 186 xval = x.clone(); 187 yval = y.clone(); 188 zval = z.clone(); 189 190 final int lastI = xLen - 1; 191 final int lastJ = yLen - 1; 192 final int lastK = zLen - 1; 193 splines = new TricubicFunction[lastI][lastJ][lastK]; 194 195 for (int i = 0; i < lastI; i++) { 196 if (f[i].length != yLen) { 197 throw new DimensionMismatchException(f[i].length, yLen); 198 } 199 if (dFdX[i].length != yLen) { 200 throw new DimensionMismatchException(dFdX[i].length, yLen); 201 } 202 if (dFdY[i].length != yLen) { 203 throw new DimensionMismatchException(dFdY[i].length, yLen); 204 } 205 if (dFdZ[i].length != yLen) { 206 throw new DimensionMismatchException(dFdZ[i].length, yLen); 207 } 208 if (d2FdXdY[i].length != yLen) { 209 throw new DimensionMismatchException(d2FdXdY[i].length, yLen); 210 } 211 if (d2FdXdZ[i].length != yLen) { 212 throw new DimensionMismatchException(d2FdXdZ[i].length, yLen); 213 } 214 if (d2FdYdZ[i].length != yLen) { 215 throw new DimensionMismatchException(d2FdYdZ[i].length, yLen); 216 } 217 if (d3FdXdYdZ[i].length != yLen) { 218 throw new DimensionMismatchException(d3FdXdYdZ[i].length, yLen); 219 } 220 221 final int ip1 = i + 1; 222 final double xR = xval[ip1] - xval[i]; 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 final double yR = yval[jp1] - yval[j]; 251 final double xRyR = xR * yR; 252 for (int k = 0; k < lastK; k++) { 253 final int kp1 = k + 1; 254 final double zR = zval[kp1] - zval[k]; 255 final double xRzR = xR * zR; 256 final double yRzR = yR * zR; 257 final double xRyRzR = xR * yRzR; 258 259 final double[] beta = new double[] { 260 f[i][j][k], f[ip1][j][k], 261 f[i][jp1][k], f[ip1][jp1][k], 262 f[i][j][kp1], f[ip1][j][kp1], 263 f[i][jp1][kp1], f[ip1][jp1][kp1], 264 265 dFdX[i][j][k] * xR, dFdX[ip1][j][k] * xR, 266 dFdX[i][jp1][k] * xR, dFdX[ip1][jp1][k] * xR, 267 dFdX[i][j][kp1] * xR, dFdX[ip1][j][kp1] * xR, 268 dFdX[i][jp1][kp1] * xR, dFdX[ip1][jp1][kp1] * xR, 269 270 dFdY[i][j][k] * yR, dFdY[ip1][j][k] * yR, 271 dFdY[i][jp1][k] * yR, dFdY[ip1][jp1][k] * yR, 272 dFdY[i][j][kp1] * yR, dFdY[ip1][j][kp1] * yR, 273 dFdY[i][jp1][kp1] * yR, dFdY[ip1][jp1][kp1] * yR, 274 275 dFdZ[i][j][k] * zR, dFdZ[ip1][j][k] * zR, 276 dFdZ[i][jp1][k] * zR, dFdZ[ip1][jp1][k] * zR, 277 dFdZ[i][j][kp1] * zR, dFdZ[ip1][j][kp1] * zR, 278 dFdZ[i][jp1][kp1] * zR, dFdZ[ip1][jp1][kp1] * zR, 279 280 d2FdXdY[i][j][k] * xRyR, d2FdXdY[ip1][j][k] * xRyR, 281 d2FdXdY[i][jp1][k] * xRyR, d2FdXdY[ip1][jp1][k] * xRyR, 282 d2FdXdY[i][j][kp1] * xRyR, d2FdXdY[ip1][j][kp1] * xRyR, 283 d2FdXdY[i][jp1][kp1] * xRyR, d2FdXdY[ip1][jp1][kp1] * xRyR, 284 285 d2FdXdZ[i][j][k] * xRzR, d2FdXdZ[ip1][j][k] * xRzR, 286 d2FdXdZ[i][jp1][k] * xRzR, d2FdXdZ[ip1][jp1][k] * xRzR, 287 d2FdXdZ[i][j][kp1] * xRzR, d2FdXdZ[ip1][j][kp1] * xRzR, 288 d2FdXdZ[i][jp1][kp1] * xRzR, d2FdXdZ[ip1][jp1][kp1] * xRzR, 289 290 d2FdYdZ[i][j][k] * yRzR, d2FdYdZ[ip1][j][k] * yRzR, 291 d2FdYdZ[i][jp1][k] * yRzR, d2FdYdZ[ip1][jp1][k] * yRzR, 292 d2FdYdZ[i][j][kp1] * yRzR, d2FdYdZ[ip1][j][kp1] * yRzR, 293 d2FdYdZ[i][jp1][kp1] * yRzR, d2FdYdZ[ip1][jp1][kp1] * yRzR, 294 295 d3FdXdYdZ[i][j][k] * xRyRzR, d3FdXdYdZ[ip1][j][k] * xRyRzR, 296 d3FdXdYdZ[i][jp1][k] * xRyRzR, d3FdXdYdZ[ip1][jp1][k] * xRyRzR, 297 d3FdXdYdZ[i][j][kp1] * xRyRzR, d3FdXdYdZ[ip1][j][kp1] * xRyRzR, 298 d3FdXdYdZ[i][jp1][kp1] * xRyRzR, d3FdXdYdZ[ip1][jp1][kp1] * xRyRzR, 299 }; 300 301 splines[i][j][k] = new TricubicFunction(computeCoefficients(beta)); 302 } 303 } 304 } 305 } 306 307 /** 308 * {@inheritDoc} 309 * 310 * @throws OutOfRangeException if any of the variables is outside its interpolation range. 311 */ 312 public double value(double x, double y, double z) 313 throws OutOfRangeException { 314 final int i = searchIndex(x, xval); 315 if (i == -1) { 316 throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]); 317 } 318 final int j = searchIndex(y, yval); 319 if (j == -1) { 320 throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]); 321 } 322 final int k = searchIndex(z, zval); 323 if (k == -1) { 324 throw new OutOfRangeException(z, zval[0], zval[zval.length - 1]); 325 } 326 327 final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]); 328 final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]); 329 final double zN = (z - zval[k]) / (zval[k + 1] - zval[k]); 330 331 return splines[i][j][k].value(xN, yN, zN); 332 } 333 334 /** 335 * Indicates whether a point is within the interpolation range. 336 * 337 * @param x First coordinate. 338 * @param y Second coordinate. 339 * @param z Third coordinate. 340 * @return {@code true} if (x, y, z) is a valid point. 341 */ 342 public boolean isValidPoint(double x, double y, double z) { 343 if (x < xval[0] || 344 x > xval[xval.length - 1] || 345 y < yval[0] || 346 y > yval[yval.length - 1] || 347 z < zval[0] || 348 z > zval[zval.length - 1]) { 349 return false; 350 } else { 351 return true; 352 } 353 } 354 355 /** 356 * @param c Coordinate. 357 * @param val Coordinate samples. 358 * @return the index in {@code val} corresponding to the interval containing {@code c}, or {@code -1} 359 * if {@code c} is out of the range defined by the end values of {@code val}. 360 */ 361 private int searchIndex(double c, double[] val) { 362 if (c < val[0]) { 363 return -1; 364 } 365 366 final int max = val.length; 367 for (int i = 1; i < max; i++) { 368 if (c <= val[i]) { 369 return i - 1; 370 } 371 } 372 373 return -1; 374 } 375 376 /** 377 * Compute the spline coefficients from the list of function values and 378 * function partial derivatives values at the four corners of a grid 379 * element. They must be specified in the following order: 380 * <ul> 381 * <li>f(0,0,0)</li> 382 * <li>f(1,0,0)</li> 383 * <li>f(0,1,0)</li> 384 * <li>f(1,1,0)</li> 385 * <li>f(0,0,1)</li> 386 * <li>f(1,0,1)</li> 387 * <li>f(0,1,1)</li> 388 * <li>f(1,1,1)</li> 389 * 390 * <li>f<sub>x</sub>(0,0,0)</li> 391 * <li>... <em>(same order as above)</em></li> 392 * <li>f<sub>x</sub>(1,1,1)</li> 393 * 394 * <li>f<sub>y</sub>(0,0,0)</li> 395 * <li>... <em>(same order as above)</em></li> 396 * <li>f<sub>y</sub>(1,1,1)</li> 397 * 398 * <li>f<sub>z</sub>(0,0,0)</li> 399 * <li>... <em>(same order as above)</em></li> 400 * <li>f<sub>z</sub>(1,1,1)</li> 401 * 402 * <li>f<sub>xy</sub>(0,0,0)</li> 403 * <li>... <em>(same order as above)</em></li> 404 * <li>f<sub>xy</sub>(1,1,1)</li> 405 * 406 * <li>f<sub>xz</sub>(0,0,0)</li> 407 * <li>... <em>(same order as above)</em></li> 408 * <li>f<sub>xz</sub>(1,1,1)</li> 409 * 410 * <li>f<sub>yz</sub>(0,0,0)</li> 411 * <li>... <em>(same order as above)</em></li> 412 * <li>f<sub>yz</sub>(1,1,1)</li> 413 * 414 * <li>f<sub>xyz</sub>(0,0,0)</li> 415 * <li>... <em>(same order as above)</em></li> 416 * <li>f<sub>xyz</sub>(1,1,1)</li> 417 * </ul> 418 * where the subscripts indicate the partial derivative with respect to 419 * the corresponding variable(s). 420 * 421 * @param beta List of function values and function partial derivatives values. 422 * @return the spline coefficients. 423 */ 424 private double[] computeCoefficients(double[] beta) { 425 final int sz = 64; 426 final double[] a = new double[sz]; 427 428 for (int i = 0; i < sz; i++) { 429 double result = 0; 430 final double[] row = AINV[i]; 431 for (int j = 0; j < sz; j++) { 432 result += row[j] * beta[j]; 433 } 434 a[i] = result; 435 } 436 437 return a; 438 } 439} 440 441/** 442 * 3D-spline function. 443 * 444 */ 445class TricubicFunction 446 implements TrivariateFunction { 447 /** Number of points. */ 448 private static final short N = 4; 449 /** Coefficients */ 450 private final double[][][] a = new double[N][N][N]; 451 452 /** 453 * @param aV List of spline coefficients. 454 */ 455 TricubicFunction(double[] aV) { 456 for (int i = 0; i < N; i++) { 457 for (int j = 0; j < N; j++) { 458 for (int k = 0; k < N; k++) { 459 a[i][j][k] = aV[i + N * (j + N * k)]; 460 } 461 } 462 } 463 } 464 465 /** 466 * @param x x-coordinate of the interpolation point. 467 * @param y y-coordinate of the interpolation point. 468 * @param z z-coordinate of the interpolation point. 469 * @return the interpolated value. 470 * @throws OutOfRangeException if {@code x}, {@code y} or 471 * {@code z} are not in the interval {@code [0, 1]}. 472 */ 473 public double value(double x, double y, double z) 474 throws OutOfRangeException { 475 if (x < 0 || x > 1) { 476 throw new OutOfRangeException(x, 0, 1); 477 } 478 if (y < 0 || y > 1) { 479 throw new OutOfRangeException(y, 0, 1); 480 } 481 if (z < 0 || z > 1) { 482 throw new OutOfRangeException(z, 0, 1); 483 } 484 485 final double x2 = x * x; 486 final double x3 = x2 * x; 487 final double[] pX = { 1, x, x2, x3 }; 488 489 final double y2 = y * y; 490 final double y3 = y2 * y; 491 final double[] pY = { 1, y, y2, y3 }; 492 493 final double z2 = z * z; 494 final double z3 = z2 * z; 495 final double[] pZ = { 1, z, z2, z3 }; 496 497 double result = 0; 498 for (int i = 0; i < N; i++) { 499 for (int j = 0; j < N; j++) { 500 for (int k = 0; k < N; k++) { 501 result += a[i][j][k] * pX[i] * pY[j] * pZ[k]; 502 } 503 } 504 } 505 506 return result; 507 } 508}