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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  72
  73 **********************************************************************/
  74 #include <simgear/compiler.h>    // MSVC: to disable C4244 d to f warning
  75
  76 #include "uiuc_2Dinterpolation.h"
  77
  78
  79 double uiuc_2Dinterpolation( double xData[100][100],
  80                              double yData[100],
  81                              double zData[100][100],
  82                              int xmax[100],
  83                              int ymax,
  84                              double x,
  85                              double y )
  86 {
  87   int xmaxl=0, xmaxu=0;
  88   double x11=0, x12=0, x21=0, x22=0, y1=0, y2=0;
  89   double z11=0, z12=0, z21=0, z22=0, z1=0, z2=0;
  90   double L11=0, L12=0, L21=0, L22=0, L1=0, L2=0;
  91   int i=2, j=2;
  92   float zfxy=0;
  93
  94   bool luby=false;       //upper bound on y (deflection)
  95   bool llby=false;       //lower bound on y
  96   bool lubx1=false;      //upper x bound on lower y
  97   bool llbx1=false;      //lower x bound on lower y
  98   bool lubx2=false;      //upper x bound on upper y
  99   bool llbx2=false;      //lower x bound on upper y
 100
 101   // check bounds on y (control deflection) to see if data range encloses it
 102   // NOTE: [1] is first element of all arrays, [0] not used
 103   if (y <= yData[1])          //if y less than lowest y
 104     {
 105       llby = true;
 106       y1 = yData[1];
 107       xmaxl = xmax[1];
 108       xmaxu = xmax[ymax];
 109     }
 110   else if (y >= yData[ymax])  //if y greater than greatest y
 111     {
 112       luby = true;
 113       y2 = yData[ymax];
 114       j = ymax;
 115       xmaxl = xmax[1];
 116       xmaxu = xmax[ymax];
 117     }
 118   else                        //y between ymax and ymin
 119     {
 120       /* loop increases j until y is less than a known y,
 121          e.g. elevator from LaRCsim less than de given in
 122          tabulated data; once this value is found, j becomes
 123          the upper bound and j-1 the lower bound on y */
 124       while (yData[j] <= y)    //bracket upper bound
 125         {
 126           j++;
 127         }
 128       y2 = yData[j];      //set upper bound on y
 129       y1 = yData[j-1];    //set lower bound on y
 130
 131       xmaxu = xmax[j];    //set max x on upper y
 132       xmaxl = xmax[j-1];  //set max x on lower y
 133     }
 134
 135   //check bounds on x (alpha) to see if data range encloses it
 136   //x less than lowest x on lower y curve:
 137   if (x <= xData[j-1][1])
 138     {
 139       llbx1 = true;
 140       z11 = zData[j-1][1];
 141       z1 = z2 = z11;
 142     }
 143   //x greater than greatest x on lower y curve:
 144   if (x >= xData[j-1][xmaxl])
 145     {
 146       lubx1 = true;
 147       z12 = zData[j-1][xmaxl];
 148       z1 = z2 = z12;
 149     }
 150   //x less than lowest x on upper y curve:
 151   if (x <= xData[j][1])
 152     {
 153       llbx2 = true;
 154       z21 = zData[j][1];
 155       z1 = z2 = z21;
 156     }
 157   //x greater than greatest x on upper y curve:
 158   if (x >= xData[j][xmaxu])
 159     {
 160       lubx2 = true;
 161       z22 = zData[j][xmaxu];
 162       z1 = z2 = z22;
 163     }
 164
 165   //x between xmax and x min
 166   //interpolate on lower y-curve
 167   if (llbx1 == false && lubx1 == false)
 168     {
 169       /* loop increases i until x is less than a known x,
 170          e.g. Alpha from LaRCsim less than Alpha given in
 171          tabulated data; once this value is found, i becomes
 172          the upper bound and i-1 the lower bound */
 173       //bracket x bounds on lower y curve
 174       while (xData[j-1][i] <= x)
 175         {
 176           i++;
 177         }
 178       x12 = xData[j-1][i];        //set upper x and z on lower y
 179       z12 = zData[j-1][i];
 180       x11 = xData[j-1][i-1];      //set lower x and z on lower y
 181       z11 = zData[j-1][i-1];
 182
 183       //do linear interpolation on x1 terms (lower y-curve)
 184       L11 = (x - x12) / (x11 - x12);
 185       L12 = (x - x11) / (x12 - x11);
 186       z1 = L11 * z11 + L12 * z12;
 187     }
 188   //interpolate on upper y-curve
 189   if (llbx2 == false && lubx2 == false)
 190     {
 191       //bracket x bounds on upper y curve
 192       i = 1;
 193       while (xData[j][i] <= x)
 194         {
 195           i++;
 196         }
 197       x22 = xData[j][i];          //set upper x and z on upper y
 198       z22 = zData[j][i];
 199       x21 = xData[j][i-1];        //set lower x and z on upper y
 200       z21 = zData[j][i-1];
 201
 202       //do linear interpolation on x2 terms (upper y-curve)
 203       L21 = (x - x22) / (x21 - x22);
 204       L22 = (x - x21) / (x22 - x21);
 205       z2 = L21 * z21 + L22 * z22;
 206     }
 207
 208   //now have all data needed to find coefficient, check cases:
 209   if (llby == true)
 210     {
 211       if (llbx1 == true || llbx2 == true)
 212         {
 213           z1 = z11;
 214           z2 = z21;
 215         }
 216       if (lubx1 == true || lubx2 == true)
 217         {
 218           z1 = z12;
 219           z2 = z22;
 220         }
 221       else
 222         {
 223           zfxy = z1;
 224         }
 225     }
 226   else if (luby == true)
 227     {
 228       if (llbx1 == true || llbx2 == true)
 229         {
 230           z1 = z11;
 231           z2 = z21;
 232         }
 233       if (lubx1 == true || lubx2 == true)
 234         {
 235           z1 = z12;
 236           z2 = z22;
 237         }
 238       else
 239         {
 240           zfxy = z2;
 241         }
 242     }
 243
 244   //do linear interpolation on y terms
 245   L1 = (y - y2) / (y1 - y2);
 246   L2 = (y - y1) / (y2 - y1);
 247
 248   //solve for z=f(x,y)
 249   zfxy = L1 * z1 + L2 * z2;
 250
 251   return zfxy;
 252 }
 253
 254 // end uiuc_2Dinterpolation.cpp