src/FDM/UIUCModel/uiuc_2Dinterpolation.cpp

   1 /**********************************************************************
   2
   3  FILENAME:     uiuc_2Dinterpolation.cpp
   4
   5 ----------------------------------------------------------------------
   6
   7  DESCRIPTION:  reads in the zData, yData, and xData arrays and the
   8                values of x and y to be interpolated on; performs 2D
   9                interpolation, i.e. z=f(x,y)
  10
  11 ----------------------------------------------------------------------
  12
  13  STATUS:       alpha version
  14
  15 ----------------------------------------------------------------------
  16
  17  REFERENCES:   similar to 2D interpolation in Selig's propid code
  18                (see alcl.f)
  19                mathematics based on linear interpolation functions
  20                (see 1Dinterpolation.cpp for references)
  21
  22 ----------------------------------------------------------------------
  23
  24  HISTORY:      02/06/2000   initial release
  25
  26 ----------------------------------------------------------------------
  27
  28  AUTHOR(S):    Jeff Scott         <jscott@mail.com>
  29
  30 ----------------------------------------------------------------------
  31
  32  VARIABLES:
  33
  34 ----------------------------------------------------------------------
  35
  36  INPUTS:       -2D array of x data for each y case
  37                -1D array of y data
  38                -2D array of z data for each y case
  39                -1D array of max number of x-z data sets for each y case
  40                -max number of y data
  41                -x value to be interpolated on
  42                -y value to be interpolated on
  43
  44 ----------------------------------------------------------------------
  45
  46  OUTPUTS:      -z as function of x and y
  47
  48 ----------------------------------------------------------------------
  49
  50  CALLED BY:    uiuc_coefficients.cpp
  51
  52 ----------------------------------------------------------------------
  53
  54  CALLS TO:     none
  55
  56 ----------------------------------------------------------------------
  57
  58  COPYRIGHT:    (C) 2000 by Michael Selig
  59
  60  This program is free software; you can redistribute it and/or
  61  modify it under the terms of the GNU General Public License
  62  as published by the Free Software Foundation.
  63
  64  This program is distributed in the hope that it will be useful,
  65  but WITHOUT ANY WARRANTY; without even the implied warranty of
  66  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  67  GNU General Public License for more details.
  68
  69  You should have received a copy of the GNU General Public License
  70  along with this program; if not, write to the Free Software
  71  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
  72  USA or view http://www.gnu.org/copyleft/gpl.html.
  73
  74 **********************************************************************/
  75 #include <simgear/compiler.h>    // MSVC: to disable C4244 d to f warning
  76
  77 #include "uiuc_2Dinterpolation.h"
  78
  79
  80 double uiuc_2Dinterpolation( double xData[100][100],
  81                              double yData[100],
  82                              double zData[100][100],
  83                              int xmax[100],
  84                              int ymax,
  85                              double x,
  86                              double y )
  87 {
  88   int xmaxl=0, xmaxu=0;
  89   double x11=0, x12=0, x21=0, x22=0, y1=0, y2=0;
  90   double z11=0, z12=0, z21=0, z22=0, z1=0, z2=0;
  91   double L11=0, L12=0, L21=0, L22=0, L1=0, L2=0;
  92   int i=2, j=2;
  93   float zfxy=0;
  94
  95   bool luby=false;       //upper bound on y (deflection)
  96   bool llby=false;       //lower bound on y
  97   bool lubx1=false;      //upper x bound on lower y
  98   bool llbx1=false;      //lower x bound on lower y
  99   bool lubx2=false;      //upper x bound on upper y
 100   bool llbx2=false;      //lower x bound on upper y
 101
 102   // check bounds on y (control deflection) to see if data range encloses it
 103   // NOTE: [1] is first element of all arrays, [0] not used
 104   if (y <= yData[1])          //if y less than lowest y
 105     {
 106       llby = true;
 107       y1 = yData[1];
 108       xmaxl = xmax[1];
 109       xmaxu = xmax[ymax];
 110     }
 111   else if (y >= yData[ymax])  //if y greater than greatest y
 112     {
 113       luby = true;
 114       y2 = yData[ymax];
 115       j = ymax;
 116       xmaxl = xmax[1];
 117       xmaxu = xmax[ymax];
 118     }
 119   else                        //y between ymax and ymin
 120     {
 121       /* loop increases j until y is less than a known y,
 122          e.g. elevator from LaRCsim less than de given in
 123          tabulated data; once this value is found, j becomes
 124          the upper bound and j-1 the lower bound on y */
 125       while (yData[j] <= y)    //bracket upper bound
 126         {
 127           j++;
 128         }
 129       y2 = yData[j];      //set upper bound on y
 130       y1 = yData[j-1];    //set lower bound on y
 131
 132       xmaxu = xmax[j];    //set max x on upper y
 133       xmaxl = xmax[j-1];  //set max x on lower y
 134     }
 135
 136   //check bounds on x (alpha) to see if data range encloses it
 137   //x less than lowest x on lower y curve:
 138   if (x <= xData[j-1][1])
 139     {
 140       llbx1 = true;
 141       z11 = zData[j-1][1];
 142       z1 = z2 = z11;
 143     }
 144   //x greater than greatest x on lower y curve:
 145   if (x >= xData[j-1][xmaxl])
 146     {
 147       lubx1 = true;
 148       z12 = zData[j-1][xmaxl];
 149       z1 = z2 = z12;
 150     }
 151   //x less than lowest x on upper y curve:
 152   if (x <= xData[j][1])
 153     {
 154       llbx2 = true;
 155       z21 = zData[j][1];
 156       z1 = z2 = z21;
 157     }
 158   //x greater than greatest x on upper y curve:
 159   if (x >= xData[j][xmaxu])
 160     {
 161       lubx2 = true;
 162       z22 = zData[j][xmaxu];
 163       z1 = z2 = z22;
 164     }
 165
 166   //x between xmax and x min
 167   //interpolate on lower y-curve
 168   if (llbx1 == false && lubx1 == false)
 169     {
 170       /* loop increases i until x is less than a known x,
 171          e.g. Alpha from LaRCsim less than Alpha given in
 172          tabulated data; once this value is found, i becomes
 173          the upper bound and i-1 the lower bound */
 174       //bracket x bounds on lower y curve
 175       while (xData[j-1][i] <= x)
 176         {
 177           i++;
 178         }
 179       x12 = xData[j-1][i];        //set upper x and z on lower y
 180       z12 = zData[j-1][i];
 181       x11 = xData[j-1][i-1];      //set lower x and z on lower y
 182       z11 = zData[j-1][i-1];
 183
 184       //do linear interpolation on x1 terms (lower y-curve)
 185       L11 = (x - x12) / (x11 - x12);
 186       L12 = (x - x11) / (x12 - x11);
 187       z1 = L11 * z11 + L12 * z12;
 188     }
 189   //interpolate on upper y-curve
 190   if (llbx2 == false && lubx2 == false)
 191     {
 192       //bracket x bounds on upper y curve
 193       i = 1;
 194       while (xData[j][i] <= x)
 195         {
 196           i++;
 197         }
 198       x22 = xData[j][i];          //set upper x and z on upper y
 199       z22 = zData[j][i];
 200       x21 = xData[j][i-1];        //set lower x and z on upper y
 201       z21 = zData[j][i-1];
 202
 203       //do linear interpolation on x2 terms (upper y-curve)
 204       L21 = (x - x22) / (x21 - x22);
 205       L22 = (x - x21) / (x22 - x21);
 206       z2 = L21 * z21 + L22 * z22;
 207     }
 208
 209   //now have all data needed to find coefficient, check cases:
 210   if (llby == true)
 211     {
 212       if (llbx1 == true || llbx2 == true)
 213         {
 214           z1 = z11;
 215           z2 = z21;
 216         }
 217       if (lubx1 == true || lubx2 == true)
 218         {
 219           z1 = z12;
 220           z2 = z22;
 221         }
 222       else
 223         {
 224           zfxy = z1;
 225         }
 226     }
 227   else if (luby == true)
 228     {
 229       if (llbx1 == true || llbx2 == true)
 230         {
 231           z1 = z11;
 232           z2 = z21;
 233         }
 234       if (lubx1 == true || lubx2 == true)
 235         {
 236           z1 = z12;
 237           z2 = z22;
 238         }
 239       else
 240         {
 241           zfxy = z2;
 242         }
 243     }
 244
 245   //do linear interpolation on y terms
 246   L1 = (y - y2) / (y1 - y2);
 247   L2 = (y - y1) / (y2 - y1);
 248
 249   //solve for z=f(x,y)
 250   zfxy = L1 * z1 + L2 * z2;
 251
 252   return zfxy;
 253 }
 254
 255 // end uiuc_2Dinterpolation.cpp