1 /**********************************************************************
3 FILENAME: uiuc_3Dinterpolation.cpp
5 ----------------------------------------------------------------------
7 DESCRIPTION: A 3D interpolator. Does a linear interpolation between
8 two values that were found from using the 2D
9 interpolator (3Dinterpolation()), or uses 3Dinterp_quick()
10 to perform a 3D linear interpolation on "nice" data
12 ----------------------------------------------------------------------
16 ----------------------------------------------------------------------
20 ----------------------------------------------------------------------
22 HISTORY: 11/07/2001 initial release
23 02/18/2002 (RD) Created uiuc_3Dinterp_quick() to take
24 advantage of the "nice" format of the
25 nonlinear Twin Otter data. Performs a
26 quicker 3D interpolation. Modified
27 uiuc_3Dinterpolation() to handle new input
29 ----------------------------------------------------------------------
31 AUTHOR(S): Robert Deters <rdeters@uiuc.edu>
33 ----------------------------------------------------------------------
37 ----------------------------------------------------------------------
41 ----------------------------------------------------------------------
45 ----------------------------------------------------------------------
47 CALLED BY: uiuc_coef_drag
54 ----------------------------------------------------------------------
56 CALLS TO: 2Dinterpolation
58 ----------------------------------------------------------------------
60 COPYRIGHT: (C) 2001 by Michael Selig
62 This program is free software; you can redistribute it and/or
63 modify it under the terms of the GNU General Public License
64 as published by the Free Software Foundation.
66 This program is distributed in the hope that it will be useful,
67 but WITHOUT ANY WARRANTY; without even the implied warranty of
68 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
69 GNU General Public License for more details.
71 You should have received a copy of the GNU General Public License
72 along with this program; if not, write to the Free Software
73 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
75 **********************************************************************/
76 #include <simgear/compiler.h> // MSVC: to disable C4244 d to f warning
78 #include "uiuc_3Dinterpolation.h"
80 void uiuc_1DtoSingle( int temp1Darray[30],
84 single_value = temp1Darray[filenumber];
87 void uiuc_2Dto1D( double temp2Darray[30][100],
93 for (count1=0; count1<=99; count1++)
95 array1D[count1] = temp2Darray[filenumber][count1];
99 void uiuc_2Dto1D_int( int temp2Darray[30][100],
105 for (count1=0; count1<=99; count1++)
107 array1D[count1] = temp2Darray[filenumber][count1];
111 void uiuc_3Dto2D( double temp3Darray[30][100][100],
113 double array2D[100][100])
117 for (count1=0; count1<=99; count1++)
119 for (count2=0; count2<=99; count2++)
121 array2D[count1][count2] = temp3Darray[filenumber][count1][count2];
126 double uiuc_3Dinterpolation( double third_Array[30],
127 double full_xArray[30][100][100],
128 double full_yArray[30][100],
129 double full_zArray[30][100][100],
130 int full_nxArray[30][100],
137 double reduced_xArray[100][100], reduced_yArray[100];
138 double reduced_zArray[100][100];
139 int reduced_nxArray[100], reduced_ny;
141 double interpmin, interpmax, third_u, third_l;
145 bool third_same=false;
147 if (third_bet <= third_Array[1])
152 else if (third_bet >= third_Array[third_max])
154 third_min = third_max;
159 while (third_Array[k] <= third_bet)
169 uiuc_3Dto2D(full_xArray, third_min, reduced_xArray);
170 uiuc_2Dto1D(full_yArray, third_min, reduced_yArray);
171 uiuc_3Dto2D(full_zArray, third_min, reduced_zArray);
172 uiuc_2Dto1D_int(full_nxArray, third_min, reduced_nxArray);
173 uiuc_1DtoSingle(full_ny, third_min, reduced_ny);
175 interpI = uiuc_2Dinterpolation(reduced_xArray,
185 uiuc_3Dto2D(full_xArray, third_min, reduced_xArray);
186 uiuc_2Dto1D(full_yArray, third_min, reduced_yArray);
187 uiuc_3Dto2D(full_zArray, third_min, reduced_zArray);
188 uiuc_2Dto1D_int(full_nxArray, third_min, reduced_nxArray);
189 uiuc_1DtoSingle(full_ny, third_min, reduced_ny);
191 interpmin = uiuc_2Dinterpolation(reduced_xArray,
199 uiuc_3Dto2D(full_xArray, third_max, reduced_xArray);
200 uiuc_2Dto1D(full_yArray, third_max, reduced_yArray);
201 uiuc_3Dto2D(full_zArray, third_max, reduced_zArray);
202 uiuc_2Dto1D_int(full_nxArray, third_max, reduced_nxArray);
203 uiuc_1DtoSingle(full_ny, third_max, reduced_ny);
205 interpmax = uiuc_2Dinterpolation(reduced_xArray,
213 third_u = third_Array[third_max];
214 third_l = third_Array[third_min];
216 interpI=interpmax - (third_u-third_bet)*(interpmax-interpmin)/(third_u-third_l);
224 double uiuc_3Dinterp_quick( double z[30],
227 double fxyz[30][100][100],
236 int xnuml, xnumu, ynuml, ynumu, znuml, znumu;
237 double xl, xu, yl, yu, zl, zu;
238 double ptxl, ptxu, ptyl, ptyu, ptylxl, ptylxu, ptyuxl, ptyuxu;
239 double ptzl, ptzu, ptzlxl, ptzlxu, ptzuxl, ptzuxu;
240 double ptzlyl, ptzlyu, ptzuyl, ptzuyu;
241 double ptzlylxl, ptzlylxu, ptzlyuxl, ptzlyuxu;
242 double ptzuylxl, ptzuylxu, ptzuyuxl, ptzuyuxu, data_point;
259 else if (zp >= z[zmax])
280 else if (yp >= y[ymax])
302 else if (xp >= x[xmax])
321 data_point = fxyz[znuml][ynuml][xnuml];
325 ptxl = fxyz[znuml][ynuml][xnuml];
326 ptxu = fxyz[znuml][ynuml][xnumu];
327 data_point = ptxu - (xu-xp)*(ptxu-ptxl)/(xu-xl);
331 ptyl = fxyz[znuml][ynuml][xnuml];
332 ptyu = fxyz[znuml][ynumu][xnuml];
333 data_point = ptyu - (yu-yp)*(ptyu-ptyl)/(yu-yl);
337 ptylxl = fxyz[znuml][ynuml][xnuml];
338 ptylxu = fxyz[znuml][ynuml][xnumu];
339 ptyuxl = fxyz[znuml][ynumu][xnuml];
340 ptyuxu = fxyz[znuml][ynumu][xnumu];
341 ptyl = ptylxu - (xu-xp)*(ptylxu-ptylxl)/(xu-xl);
342 ptyu = ptyuxu - (xu-xp)*(ptyuxu-ptyuxl)/(xu-xl);
343 data_point = ptyu - (yu-yp)*(ptyu-ptyl)/(yu-yl);
350 ptzl = fxyz[znuml][ynuml][xnuml];
351 ptzu = fxyz[znumu][ynuml][xnuml];
352 data_point = ptzu - (zu-zp)*(ptzu-ptzl)/(zu-zl);
356 ptzlxl = fxyz[znuml][ynuml][xnuml];
357 ptzlxu = fxyz[znuml][ynuml][xnumu];
358 ptzuxl = fxyz[znumu][ynuml][xnuml];
359 ptzuxu = fxyz[znumu][ynuml][xnumu];
360 ptzl = ptzlxu - (xu-xp)*(ptzlxu-ptzlxl)/(xu-xl);
361 ptzu = ptzuxu - (xu-xp)*(ptzuxu-ptzuxl)/(xu-xl);
362 data_point = ptzu - (zu-zp)*(ptzu-ptzl)/(zu-zl);
366 ptzlyl = fxyz[znuml][ynuml][xnuml];
367 ptzlyu = fxyz[znuml][ynumu][xnuml];
368 ptzuyl = fxyz[znumu][ynuml][xnuml];
369 ptzuyu = fxyz[znumu][ynumu][xnuml];
370 ptzl = ptzlyu - (yu-yp)*(ptzlyu-ptzlyl)/(yu-yl);
371 ptzu = ptzuyu - (yu-yp)*(ptzuyu-ptzuyl)/(yu-yl);
372 data_point = ptzu - (zu-zp)*(ptzu-ptzl)/(zu-zl);
376 ptzlylxl = fxyz[znuml][ynuml][xnuml];
377 ptzlylxu = fxyz[znuml][ynuml][xnumu];
378 ptzlyuxl = fxyz[znuml][ynumu][xnuml];
379 ptzlyuxu = fxyz[znuml][ynumu][xnumu];
380 ptzuylxl = fxyz[znumu][ynuml][xnuml];
381 ptzuylxu = fxyz[znumu][ynuml][xnumu];
382 ptzuyuxl = fxyz[znumu][ynumu][xnuml];
383 ptzuyuxu = fxyz[znumu][ynumu][xnumu];
384 ptzlyl = ptzlylxu - (xu-xp)*(ptzlylxu-ptzlylxl)/(xu-xl);
385 ptzlyu = ptzlyuxu - (xu-xp)*(ptzlyuxu-ptzlyuxl)/(xu-xl);
386 ptzuyl = ptzuylxu - (xu-xp)*(ptzuylxu-ptzuylxl)/(xu-xl);
387 ptzuyu = ptzuyuxu - (xu-xp)*(ptzuyuxu-ptzuyuxl)/(xu-xl);
388 ptzl = ptzlyu - (yu-yp)*(ptzlyu-ptzlyl)/(yu-yl);
389 ptzu = ptzuyu - (yu-yp)*(ptzuyu-ptzuyl)/(yu-yl);
390 data_point = ptzu - (zu-zp)*(ptzu-ptzl)/(zu-zl);