source: src/main/java/agents/anac/y2019/harddealer/math3/analysis/interpolation/TricubicInterpolatingFunction.java

Last change on this file was 204, checked in by Katsuhide Fujita, 5 years ago

Fixed errors of ANAC2019 agents
  • Property svn:executable set to *
File size: 26.4 KB
Line 
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 agents.anac.y2019.harddealer.math3.analysis.interpolation;
18
19import agents.anac.y2019.harddealer.math3.analysis.TrivariateFunction;
20import agents.anac.y2019.harddealer.math3.exception.DimensionMismatchException;
21import agents.anac.y2019.harddealer.math3.exception.NoDataException;
22import agents.anac.y2019.harddealer.math3.exception.NonMonotonicSequenceException;
23import agents.anac.y2019.harddealer.math3.exception.OutOfRangeException;
24import agents.anac.y2019.harddealer.math3.util.MathArrays;
25
26/**
27 * Function that implements the
28 * <a href="http://en.wikipedia.org/wiki/Tricubic_interpolation">
29 * tricubic spline interpolation</a>, as proposed in
30 * <blockquote>
31 * Tricubic interpolation in three dimensions,
32 * F. Lekien and J. Marsden,
33 * <em>Int. J. Numer. Meth. Eng</em> 2005; <b>63</b>:455-471
34 * </blockquote>
35 *
36 * @since 3.4.
37 */
38public class TricubicInterpolatingFunction
39 implements TrivariateFunction {
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 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}
Note: See TracBrowser for help on using the repository browser.