RadLib
Loading...
Searching...
No Matches
multilinear_interpolation.h
Go to the documentation of this file.
1
16
17#include <vector>
18
19using namespace std;
20
22int get_location_dim(const int nx, const double *x, const double xP) {
23 if(xP <= x[0])
24 return 0;
25 else if(xP >= x[nx-1])
26 return nx-2;
27 else
28 return lower_bound(x, x+nx, xP) - x - 1;
29}
30
32
33double LI_1D(const int nx,
34 const double *x,
35 const double *f,
36 const double xP){
37
38 int i = get_location_dim(nx,x,xP);
39
40 int ip = i+1;
41
42 double xWta = (xP - x[i])/(x[i+1]-x[i]);
43
44 double xWtb = 1-xWta;
45
46 return (f[i ]*xWtb +
47 f[ip]*xWta);
48}
49
51
52double LI_2D(const int nx, // number of points in x
53 const int ny, // number of points in y
54 const double *x, // x grid values
55 const double *y, // y grid values
56 const double *f, // field of function values over the x, y grid (1-D row major)
57 const double xP, // interpolation point in x
58 const double yP){ // interpolation point in y
59
60 int i = get_location_dim(nx,x,xP);
61 int j = get_location_dim(ny,y,yP);
62
63 int ip = i+1;
64 int jp = j+1;
65
66 double xWta = (xP - x[i])/(x[i+1]-x[i]);
67 double yWta = (yP - y[j])/(y[j+1]-y[j]);
68
69 double xWtb = 1-xWta;
70 double yWtb = 1-yWta;
71
72 return (f[i *ny +j ]*xWtb + // return f[i ][j ]*xWtb*yWtb +
73 f[ip*ny +j ]*xWta) * yWtb + // f[ip][j ]*xWta*yWtb +
74 (f[i *ny +jp]*xWtb + // f[i ][jp]*xWtb*yWta +
75 f[ip*ny +jp]*xWta) * yWta; // f[ip][jp]*xWta*yWta;
76}
77
79
80double LI_3D(const int nx,
81 const int ny,
82 const int nz,
83 const double *x,
84 const double *y,
85 const double *z,
86 const double *f,
87 const double xP,
88 const double yP,
89 const double zP){
90
91 int i = get_location_dim(nx,x,xP);
92 int j = get_location_dim(ny,y,yP);
93 int k = get_location_dim(nz,z,zP);
94
95 int ip = i+1;
96 int jp = j+1;
97 int kp = k+1;
98
99 double xWta = (xP - x[i])/(x[i+1]-x[i]);
100 double yWta = (yP - y[j])/(y[j+1]-y[j]);
101 double zWta = (zP - z[k])/(z[k+1]-z[k]);
102
103 double xWtb = 1-xWta;
104 double yWtb = 1-yWta;
105 double zWtb = 1-zWta;
106
107 int nzy = nz*ny;
108
109 return ((f[i *nzy + j *nz + k ]*xWtb + // return f[i ][j ][k ]*xWtb*yWtb*zWtb +
110 f[ip*nzy + j *nz + k ]*xWta) * yWtb + // f[ip][j ][k ]*xWta*yWtb*zWtb +
111 (f[i *nzy + jp*nz + k ]*xWtb + // f[i ][jp][k ]*xWtb*yWta*zWtb +
112 f[ip*nzy + jp*nz + k ]*xWta) * yWta) * zWtb + // f[ip][jp][k ]*xWta*yWta*zWtb +
113 ((f[i *nzy + j *nz + kp]*xWtb + // f[i ][j ][kp]*xWtb*yWtb*zWta +
114 f[ip*nzy + j *nz + kp]*xWta) * yWtb + // f[ip][j ][kp]*xWta*yWtb*zWta +
115 (f[i *nzy + jp*nz + kp]*xWtb + // f[i ][jp][kp]*xWtb*yWta*zWta +
116 f[ip*nzy + jp*nz + kp]*xWta) * yWta) * zWta; // f[ip][jp][kp]*xWta*yWta*zWta;
117}
118
120
121double LI_4D(const int nx,
122 const int ny,
123 const int nz,
124 const int nw,
125 const double *x,
126 const double *y,
127 const double *z,
128 const double *w,
129 const double *f,
130 const double xP,
131 const double yP,
132 const double zP,
133 const double wP){
134
135 int i = get_location_dim(nx,x,xP);
136 int j = get_location_dim(ny,y,yP);
137 int k = get_location_dim(nz,z,zP);
138 int l = get_location_dim(nw,w,wP);
139
140 int ip = i+1;
141 int jp = j+1;
142 int kp = k+1;
143 int lp = l+1;
144
145 double xWta = (xP - x[i])/(x[i+1]-x[i]);
146 double yWta = (yP - y[j])/(y[j+1]-y[j]);
147 double zWta = (zP - z[k])/(z[k+1]-z[k]);
148 double wWta = (wP - w[l])/(w[l+1]-w[l]);
149
150 double xWtb = 1-xWta;
151 double yWtb = 1-yWta;
152 double zWtb = 1-zWta;
153 double wWtb = 1-wWta;
154
155 int nwzy = nw*nz*ny;
156 int nwz = nw*nz;
157
158 return (((f[i *nwzy + j *nwz + k *nw + l ]*xWtb + // return f[i ][j ][k ][l ]*xWtb*yWtb*zWtb*wWtb +
159 f[ip*nwzy + j *nwz + k *nw + l ]*xWta) * yWtb + // f[ip][j ][k ][l ]*xWta*yWtb*zWtb*wWtb +
160 (f[i *nwzy + jp*nwz + k *nw + l ]*xWtb + // f[i ][jp][k ][l ]*xWtb*yWta*zWtb*wWtb +
161 f[ip*nwzy + jp*nwz + k *nw + l ]*xWta) * yWta) * zWtb + // f[ip][jp][k ][l ]*xWta*yWta*zWtb*wWtb +
162 ((f[i *nwzy + j *nwz + kp*nw + l ]*xWtb + // f[i ][j ][kp][l ]*xWtb*yWtb*zWta*wWtb +
163 f[ip*nwzy + j *nwz + kp*nw + l ]*xWta) * yWtb + // f[ip][j ][kp][l ]*xWta*yWtb*zWta*wWtb +
164 (f[i *nwzy + jp*nwz + kp*nw + l ]*xWtb + // f[i ][jp][kp][l ]*xWtb*yWta*zWta*wWtb +
165 f[ip*nwzy + jp*nwz + kp*nw + l ]*xWta) * yWta) * zWta) * wWtb + // f[ip][jp][kp][l ]*xWta*yWta*zWta*wWtb +
166 (((f[i *nwzy + j *nwz + k *nw + lp]*xWtb + // f[i ][j ][k ][lp]*xWtb*yWtb*zWtb*wWta +
167 f[ip*nwzy + j *nwz + k *nw + lp]*xWta) * yWtb + // f[ip][j ][k ][lp]*xWta*yWtb*zWtb*wWta +
168 (f[i *nwzy + jp*nwz + k *nw + lp]*xWtb + // f[i ][jp][k ][lp]*xWtb*yWta*zWtb*wWta +
169 f[ip*nwzy + jp*nwz + k *nw + lp]*xWta) * yWta) * zWtb + // f[ip][jp][k ][lp]*xWta*yWta*zWtb*wWta +
170 ((f[i *nwzy + j *nwz + kp*nw + lp]*xWtb + // f[i ][j ][kp][lp]*xWtb*yWtb*zWta*wWta +
171 f[ip*nwzy + j *nwz + kp*nw + lp]*xWta) * yWtb + // f[ip][j ][kp][lp]*xWta*yWtb*zWta*wWta +
172 (f[i *nwzy + jp*nwz + kp*nw + lp]*xWtb + // f[i ][jp][kp][lp]*xWtb*yWta*zWta*wWta +
173 f[ip*nwzy + jp*nwz + kp*nw + lp]*xWta) * yWta) * zWta) * wWta; // f[ip][jp][kp][lp]*xWta*yWta*zWta*wWta;
174}
175
177
178double LI_5D(const int nx,
179 const int ny,
180 const int nz,
181 const int nw,
182 const int nv,
183 const double *x,
184 const double *y,
185 const double *z,
186 const double *w,
187 const double *v,
188 const double *f,
189 const double xP,
190 const double yP,
191 const double zP,
192 const double wP,
193 const double vP){
194
195 int i = get_location_dim(nx,x,xP);
196 int j = get_location_dim(ny,y,yP);
197 int k = get_location_dim(nz,z,zP);
198 int l = get_location_dim(nw,w,wP);
199 int m = get_location_dim(nv,v,vP);
200
201 int ip = i+1;
202 int jp = j+1;
203 int kp = k+1;
204 int lp = l+1;
205 int mp = m+1;
206
207 double xWta = (xP - x[i])/(x[i+1]-x[i]);
208 double yWta = (yP - y[j])/(y[j+1]-y[j]);
209 double zWta = (zP - z[k])/(z[k+1]-z[k]);
210 double wWta = (wP - w[l])/(w[l+1]-w[l]);
211 double vWta = (vP - v[m])/(v[m+1]-v[m]);
212
213 double xWtb = 1-xWta;
214 double yWtb = 1-yWta;
215 double zWtb = 1-zWta;
216 double wWtb = 1-wWta;
217 double vWtb = 1-vWta;
218
219 int nvwzy = nv*nw*nz*ny;
220 int nvwz = nv*nw*nz;
221 int nvw = nv*nw;
222
223 return ((((f[i *nvwzy + j *nvwz + k *nvw + l *nv + m ]*xWtb + // return f[i ][j ][k ][l ][m ]*xWtb*yWtb*zWtb*wWtb*vWtb +
224 f[ip*nvwzy + j *nvwz + k *nvw + l *nv + m ]*xWta) * yWtb + // f[ip][j ][k ][l ][m ]*xWta*yWtb*zWtb*wWtb*vWtb +
225 (f[i *nvwzy + jp*nvwz + k *nvw + l *nv + m ]*xWtb + // f[i ][jp][k ][l ][m ]*xWtb*yWta*zWtb*wWtb*vWtb +
226 f[ip*nvwzy + jp*nvwz + k *nvw + l *nv + m ]*xWta) * yWta) * zWtb + // f[ip][jp][k ][l ][m ]*xWta*yWta*zWtb*wWtb*vWtb +
227 ((f[i *nvwzy + j *nvwz + kp*nvw + l *nv + m ]*xWtb + // f[i ][j ][kp][l ][m ]*xWtb*yWtb*zWta*wWtb*vWtb +
228 f[ip*nvwzy + j *nvwz + kp*nvw + l *nv + m ]*xWta) * yWtb + // f[ip][j ][kp][l ][m ]*xWta*yWtb*zWta*wWtb*vWtb +
229 (f[i *nvwzy + jp*nvwz + kp*nvw + l *nv + m ]*xWtb + // f[i ][jp][kp][l ][m ]*xWtb*yWta*zWta*wWtb*vWtb +
230 f[ip*nvwzy + jp*nvwz + kp*nvw + l *nv + m ]*xWta) * yWta) * zWta) * wWtb + // f[ip][jp][kp][l ][m ]*xWta*yWta*zWta*wWtb*vWtb +
231 (((f[i *nvwzy + j *nvwz + k *nvw + lp*nv + m ]*xWtb + // f[i ][j ][k ][lp][m ]*xWtb*yWtb*zWtb*wWta*vWtb +
232 f[ip*nvwzy + j *nvwz + k *nvw + lp*nv + m ]*xWta) * yWtb + // f[ip][j ][k ][lp][m ]*xWta*yWtb*zWtb*wWta*vWtb +
233 (f[i *nvwzy + jp*nvwz + k *nvw + lp*nv + m ]*xWtb + // f[i ][jp][k ][lp][m ]*xWtb*yWta*zWtb*wWta*vWtb +
234 f[ip*nvwzy + jp*nvwz + k *nvw + lp*nv + m ]*xWta) * yWta) * zWtb + // f[ip][jp][k ][lp][m ]*xWta*yWta*zWtb*wWta*vWtb +
235 ((f[i *nvwzy + j *nvwz + kp*nvw + lp*nv + m ]*xWtb + // f[i ][j ][kp][lp][m ]*xWtb*yWtb*zWta*wWta*vWtb +
236 f[ip*nvwzy + j *nvwz + kp*nvw + lp*nv + m ]*xWta) * yWtb + // f[ip][j ][kp][lp][m ]*xWta*yWtb*zWta*wWta*vWtb +
237 (f[i *nvwzy + jp*nvwz + kp*nvw + lp*nv + m ]*xWtb + // f[i ][jp][kp][lp][m ]*xWtb*yWta*zWta*wWta*vWtb +
238 f[ip*nvwzy + jp*nvwz + kp*nvw + lp*nv + m ]*xWta) * yWta) * zWta) * wWta) * vWtb + // f[ip][jp][kp][lp][m ]*xWta*yWta*zWta*wWta*vWtb +
239 ((((f[i *nvwzy + j *nvwz + k *nvw + l *nv + mp]*xWtb + // f[i ][j ][k ][l ][mp]*xWtb*yWtb*zWtb*wWtb*vWta +
240 f[ip*nvwzy + j *nvwz + k *nvw + l *nv + mp]*xWta) * yWtb + // f[ip][j ][k ][l ][mp]*xWta*yWtb*zWtb*wWtb*vWta +
241 (f[i *nvwzy + jp*nvwz + k *nvw + l *nv + mp]*xWtb + // f[i ][jp][k ][l ][mp]*xWtb*yWta*zWtb*wWtb*vWta +
242 f[ip*nvwzy + jp*nvwz + k *nvw + l *nv + mp]*xWta) * yWta) * zWtb + // f[ip][jp][k ][l ][mp]*xWta*yWta*zWtb*wWtb*vWta +
243 ((f[i *nvwzy + j *nvwz + kp*nvw + l *nv + mp]*xWtb + // f[i ][j ][kp][l ][mp]*xWtb*yWtb*zWta*wWtb*vWta +
244 f[ip*nvwzy + j *nvwz + kp*nvw + l *nv + mp]*xWta) * yWtb + // f[ip][j ][kp][l ][mp]*xWta*yWtb*zWta*wWtb*vWta +
245 (f[i *nvwzy + jp*nvwz + kp*nvw + l *nv + mp]*xWtb + // f[i ][jp][kp][l ][mp]*xWtb*yWta*zWta*wWtb*vWta +
246 f[ip*nvwzy + jp*nvwz + kp*nvw + l *nv + mp]*xWta) * yWta) * zWta) * wWtb + // f[ip][jp][kp][l ][mp]*xWta*yWta*zWta*wWtb*vWta +
247 (((f[i *nvwzy + j *nvwz + k *nvw + lp*nv + mp]*xWtb + // f[i ][j ][k ][lp][mp]*xWtb*yWtb*zWtb*wWta*vWta +
248 f[ip*nvwzy + j *nvwz + k *nvw + lp*nv + mp]*xWta) * yWtb + // f[ip][j ][k ][lp][mp]*xWta*yWtb*zWtb*wWta*vWta +
249 (f[i *nvwzy + jp*nvwz + k *nvw + lp*nv + mp]*xWtb + // f[i ][jp][k ][lp][mp]*xWtb*yWta*zWtb*wWta*vWta +
250 f[ip*nvwzy + jp*nvwz + k *nvw + lp*nv + mp]*xWta) * yWta) * zWtb + // f[ip][jp][k ][lp][mp]*xWta*yWta*zWtb*wWta*vWta +
251 ((f[i *nvwzy + j *nvwz + kp*nvw + lp*nv + mp]*xWtb + // f[i ][j ][kp][lp][mp]*xWtb*yWtb*zWta*wWta*vWta +
252 f[ip*nvwzy + j *nvwz + kp*nvw + lp*nv + mp]*xWta) * yWtb + // f[ip][j ][kp][lp][mp]*xWta*yWtb*zWta*wWta*vWta +
253 (f[i *nvwzy + jp*nvwz + kp*nvw + lp*nv + mp]*xWtb + // f[i ][jp][kp][lp][mp]*xWtb*yWta*zWta*wWta*vWta +
254 f[ip*nvwzy + jp*nvwz + kp*nvw + lp*nv + mp]*xWta) * yWta) * zWta) * wWta) * vWtb; // f[ip][jp][kp][lp][mp]*xWta*yWta*zWta*wWta*vWta;
255}
LI_1D
double LI_1D(const int nx, const double *x, const double *f, const double xP)
Definition multilinear_interpolation.h:33
LI_4D
double LI_4D(const int nx, const int ny, const int nz, const int nw, const double *x, const double *y, const double *z, const double *w, const double *f, const double xP, const double yP, const double zP, const double wP)
Definition multilinear_interpolation.h:121
get_location_dim
int get_location_dim(const int nx, const double *x, const double xP)
Definition multilinear_interpolation.h:22
LI_3D
double LI_3D(const int nx, const int ny, const int nz, const double *x, const double *y, const double *z, const double *f, const double xP, const double yP, const double zP)
Definition multilinear_interpolation.h:80
LI_2D
double LI_2D(const int nx, const int ny, const double *x, const double *y, const double *f, const double xP, const double yP)
Definition multilinear_interpolation.h:52
LI_5D
double LI_5D(const int nx, const int ny, const int nz, const int nw, const int nv, const double *x, const double *y, const double *z, const double *w, const double *v, const double *f, const double xP, const double yP, const double zP, const double wP, const double vP)
Definition multilinear_interpolation.h:178