1 /*
2 * Copyright (C) 2023 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.numerical.interpolation;
17
18 import com.irurueta.algebra.Matrix;
19 import com.irurueta.algebra.WrongSizeException;
20
21 /**
22 * Interpolation in two dimensions.
23 * This implementation uses a higher order than {@link BilinearInterpolator} for accuracy reasons.
24 */
25 public class Polynomial2DInterpolator {
26
27 /**
28 * Length of x1v array.
29 */
30 private final int m;
31
32 /**
33 * Length of x2v array.
34 */
35 private final int n;
36
37 /**
38 * Number of rows of sub-block of ym values to be processed.
39 */
40 private final int mm;
41
42 /**
43 * Number of columns of sub-block of ym values to be processed.
44 */
45 private final int nn;
46
47 /**
48 * Matrix of tabulated function values yij.
49 */
50 private final Matrix y;
51
52 /**
53 * Temporary array containing interpolated values in one direction.
54 */
55 private final double[] yv;
56
57 /**
58 * One dimensional interpolator for x1v.
59 */
60 private final PolynomialInterpolator x1terp;
61
62 /**
63 * One dimensional interpolator for x2v.
64 */
65 private final PolynomialInterpolator x2terp;
66
67 /**
68 * Constructor.
69 *
70 * @param x1v array of x1v.
71 * @param x2v array of x2v.
72 * @param ym matrix of tabulated function values yij.
73 * @param mp defines number of rows of sub-block of ym values to be processed.
74 * @param np defined number of columns of sub-block of ym values to be processed.
75 */
76 public Polynomial2DInterpolator(final double[] x1v, final double[] x2v, final Matrix ym, final int mp,
77 final int np) {
78 m = x1v.length;
79 n = x2v.length;
80 mm = mp;
81 nn = np;
82 y = ym;
83 yv = new double[m];
84 // Dummy 1-dim interpolations for their locate and hunt functions
85 x1terp = new PolynomialInterpolator(x1v, yv, mm);
86 x2terp = new PolynomialInterpolator(x2v, new double[n], nn);
87 }
88
89 /**
90 * Constructor.
91 *
92 * @param x1v array of x1v.
93 * @param x2v array of x2v.
94 * @param ym matrix of tabulated function values yij.
95 */
96 public Polynomial2DInterpolator(final double[] x1v, final double[] x2v, final Matrix ym) {
97 this(x1v, x2v, ym, x1v.length, x2v.length);
98 }
99
100 /**
101 * Gets length of x1v array.
102 *
103 * @return length of x1v array.
104 */
105 public int getM() {
106 return m;
107 }
108
109 /**
110 * Gets length of x2v array.
111 *
112 * @return length of x2v array.
113 */
114 public int getN() {
115 return n;
116 }
117
118 /**
119 * Gets number of rows of sub-block of ym values to be processed.
120 *
121 * @return number of rows of sub-block of ym values to be processed.
122 */
123 public int getMm() {
124 return mm;
125 }
126
127 /**
128 * Gets number of columns of sub-block of ym values to be processed.
129 *
130 * @return number of columns of sub-block of ym values to be processed.
131 */
132 public int getNn() {
133 return nn;
134 }
135
136 /**
137 * Given values x1p an x2p, returns an interpolated value.
138 *
139 * @param x1p x1p value where interpolation is estimated.
140 * @param x2p x2p value where interpolation is estimated.
141 * @return interpolated value.
142 * @throws InterpolationException if interpolation fails.
143 */
144 public double interpolate(final double x1p, final double x2p) throws InterpolationException {
145 try {
146 final var i = x1terp.cor != 0 ? x1terp.hunt(x1p) : x1terp.locate(x1p);
147 final var j = x2terp.cor != 0 ? x2terp.hunt(x2p) : x2terp.locate(x2p);
148 int k;
149
150 // Find grid block
151 for (k = i; k < i + mm; k++) {
152 // "mm" interpolations in the x2 direction.
153 // copy k-row of matrix y
154 y.getSubmatrixAsArray(k, 0, k, y.getColumns() - 1, x2terp.yy);
155 yv[k] = x2terp.rawinterp(j, x2p);
156 }
157
158 // A final interpolation in the x1 direction.
159 return x1terp.rawinterp(i, x1p);
160 } catch (final WrongSizeException e) {
161 throw new InterpolationException(e);
162 }
163 }
164 }