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., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
74 USA or view http://www.gnu.org/copyleft/gpl.html.
76 **********************************************************************/
77 #include <simgear/compiler.h> // MSVC: to disable C4244 d to f warning
79 #include "uiuc_3Dinterpolation.h"
81 void uiuc_1DtoSingle( int temp1Darray[30],
85 single_value = temp1Darray[filenumber];
88 void uiuc_2Dto1D( double temp2Darray[30][100],
94 for (count1=0; count1<=99; count1++)
96 array1D[count1] = temp2Darray[filenumber][count1];
100 void uiuc_2Dto1D_int( int temp2Darray[30][100],
106 for (count1=0; count1<=99; count1++)
108 array1D[count1] = temp2Darray[filenumber][count1];
112 void uiuc_3Dto2D( double temp3Darray[30][100][100],
114 double array2D[100][100])
118 for (count1=0; count1<=99; count1++)
120 for (count2=0; count2<=99; count2++)
122 array2D[count1][count2] = temp3Darray[filenumber][count1][count2];
127 double uiuc_3Dinterpolation( double third_Array[30],
128 double full_xArray[30][100][100],
129 double full_yArray[30][100],
130 double full_zArray[30][100][100],
131 int full_nxArray[30][100],
138 double reduced_xArray[100][100], reduced_yArray[100];
139 double reduced_zArray[100][100];
140 int reduced_nxArray[100], reduced_ny;
142 double interpmin, interpmax, third_u, third_l;
146 bool third_same=false;
148 if (third_bet <= third_Array[1])
153 else if (third_bet >= third_Array[third_max])
155 third_min = third_max;
160 while (third_Array[k] <= third_bet)
170 uiuc_3Dto2D(full_xArray, third_min, reduced_xArray);
171 uiuc_2Dto1D(full_yArray, third_min, reduced_yArray);
172 uiuc_3Dto2D(full_zArray, third_min, reduced_zArray);
173 uiuc_2Dto1D_int(full_nxArray, third_min, reduced_nxArray);
174 uiuc_1DtoSingle(full_ny, third_min, reduced_ny);
176 interpI = uiuc_2Dinterpolation(reduced_xArray,
186 uiuc_3Dto2D(full_xArray, third_min, reduced_xArray);
187 uiuc_2Dto1D(full_yArray, third_min, reduced_yArray);
188 uiuc_3Dto2D(full_zArray, third_min, reduced_zArray);
189 uiuc_2Dto1D_int(full_nxArray, third_min, reduced_nxArray);
190 uiuc_1DtoSingle(full_ny, third_min, reduced_ny);
192 interpmin = uiuc_2Dinterpolation(reduced_xArray,
200 uiuc_3Dto2D(full_xArray, third_max, reduced_xArray);
201 uiuc_2Dto1D(full_yArray, third_max, reduced_yArray);
202 uiuc_3Dto2D(full_zArray, third_max, reduced_zArray);
203 uiuc_2Dto1D_int(full_nxArray, third_max, reduced_nxArray);
204 uiuc_1DtoSingle(full_ny, third_max, reduced_ny);
206 interpmax = uiuc_2Dinterpolation(reduced_xArray,
214 third_u = third_Array[third_max];
215 third_l = third_Array[third_min];
217 interpI=interpmax - (third_u-third_bet)*(interpmax-interpmin)/(third_u-third_l);
225 double uiuc_3Dinterp_quick( double z[30],
228 double fxyz[30][100][100],
237 int xnuml, xnumu, ynuml, ynumu, znuml, znumu;
238 double xl, xu, yl, yu, zl, zu;
239 double ptxl, ptxu, ptyl, ptyu, ptylxl, ptylxu, ptyuxl, ptyuxu;
240 double ptzl, ptzu, ptzlxl, ptzlxu, ptzuxl, ptzuxu;
241 double ptzlyl, ptzlyu, ptzuyl, ptzuyu;
242 double ptzlylxl, ptzlylxu, ptzlyuxl, ptzlyuxu;
243 double ptzuylxl, ptzuylxu, ptzuyuxl, ptzuyuxu, data_point;
260 else if (zp >= z[zmax])
281 else if (yp >= y[ymax])
303 else if (xp >= x[xmax])
322 data_point = fxyz[znuml][ynuml][xnuml];
326 ptxl = fxyz[znuml][ynuml][xnuml];
327 ptxu = fxyz[znuml][ynuml][xnumu];
328 data_point = ptxu - (xu-xp)*(ptxu-ptxl)/(xu-xl);
332 ptyl = fxyz[znuml][ynuml][xnuml];
333 ptyu = fxyz[znuml][ynumu][xnuml];
334 data_point = ptyu - (yu-yp)*(ptyu-ptyl)/(yu-yl);
338 ptylxl = fxyz[znuml][ynuml][xnuml];
339 ptylxu = fxyz[znuml][ynuml][xnumu];
340 ptyuxl = fxyz[znuml][ynumu][xnuml];
341 ptyuxu = fxyz[znuml][ynumu][xnumu];
342 ptyl = ptylxu - (xu-xp)*(ptylxu-ptylxl)/(xu-xl);
343 ptyu = ptyuxu - (xu-xp)*(ptyuxu-ptyuxl)/(xu-xl);
344 data_point = ptyu - (yu-yp)*(ptyu-ptyl)/(yu-yl);
351 ptzl = fxyz[znuml][ynuml][xnuml];
352 ptzu = fxyz[znumu][ynuml][xnuml];
353 data_point = ptzu - (zu-zp)*(ptzu-ptzl)/(zu-zl);
357 ptzlxl = fxyz[znuml][ynuml][xnuml];
358 ptzlxu = fxyz[znuml][ynuml][xnumu];
359 ptzuxl = fxyz[znumu][ynuml][xnuml];
360 ptzuxu = fxyz[znumu][ynuml][xnumu];
361 ptzl = ptzlxu - (xu-xp)*(ptzlxu-ptzlxl)/(xu-xl);
362 ptzu = ptzuxu - (xu-xp)*(ptzuxu-ptzuxl)/(xu-xl);
363 data_point = ptzu - (zu-zp)*(ptzu-ptzl)/(zu-zl);
367 ptzlyl = fxyz[znuml][ynuml][xnuml];
368 ptzlyu = fxyz[znuml][ynumu][xnuml];
369 ptzuyl = fxyz[znumu][ynuml][xnuml];
370 ptzuyu = fxyz[znumu][ynumu][xnuml];
371 ptzl = ptzlyu - (yu-yp)*(ptzlyu-ptzlyl)/(yu-yl);
372 ptzu = ptzuyu - (yu-yp)*(ptzuyu-ptzuyl)/(yu-yl);
373 data_point = ptzu - (zu-zp)*(ptzu-ptzl)/(zu-zl);
377 ptzlylxl = fxyz[znuml][ynuml][xnuml];
378 ptzlylxu = fxyz[znuml][ynuml][xnumu];
379 ptzlyuxl = fxyz[znuml][ynumu][xnuml];
380 ptzlyuxu = fxyz[znuml][ynumu][xnumu];
381 ptzuylxl = fxyz[znumu][ynuml][xnuml];
382 ptzuylxu = fxyz[znumu][ynuml][xnumu];
383 ptzuyuxl = fxyz[znumu][ynumu][xnuml];
384 ptzuyuxu = fxyz[znumu][ynumu][xnumu];
385 ptzlyl = ptzlylxu - (xu-xp)*(ptzlylxu-ptzlylxl)/(xu-xl);
386 ptzlyu = ptzlyuxu - (xu-xp)*(ptzlyuxu-ptzlyuxl)/(xu-xl);
387 ptzuyl = ptzuylxu - (xu-xp)*(ptzuylxu-ptzuylxl)/(xu-xl);
388 ptzuyu = ptzuyuxu - (xu-xp)*(ptzuyuxu-ptzuyuxl)/(xu-xl);
389 ptzl = ptzlyu - (yu-yp)*(ptzlyu-ptzlyl)/(yu-yl);
390 ptzu = ptzuyu - (yu-yp)*(ptzuyu-ptzuyl)/(yu-yl);
391 data_point = ptzu - (zu-zp)*(ptzu-ptzl)/(zu-zl);