src/FDM/UIUCModel/uiuc_flapdata.cpp

   1 /*flapdata.cpp
   2 Implements the flapping data class
   3 Written by Theresa Robinson
   4 robinst@ecf.toronto.edu
   5 */
   6
   7 //#ifndef flapdata_cpp
   8 //#define flapdata_cpp
   9 #include "uiuc_flapdata.h"
  10 //#include <fstream>
  11 #include <cassert>
  12
  13 ///////////////////////////////////////////////////////////
  14 //Implementation of FlapStruct public methods
  15 //////////////////////////////////////////////////////////
  16
  17 flapStruct::flapStruct(){
  18         Lift=0;
  19         Thrust=0;
  20         Inertia=0;
  21         Moment=0;
  22 }
  23
  24 flapStruct::flapStruct(const flapStruct &rhs){
  25         Lift=rhs.getLift();
  26         Thrust=rhs.getThrust();
  27         Inertia=rhs.getInertia();
  28         Moment=rhs.getMoment();
  29 }
  30
  31 flapStruct::flapStruct(double newLift, double newThrust,  double newMoment, double newInertia){
  32         Lift=newLift;
  33         Thrust=newThrust;
  34         Inertia=newInertia;
  35         Moment=newMoment;
  36 }
  37
  38 double flapStruct::getLift() const{
  39         return Lift;
  40 }
  41
  42 double flapStruct::getThrust() const{
  43         return Thrust;
  44 }
  45
  46 double flapStruct::getInertia() const{
  47         return Inertia;
  48 }
  49
  50 double flapStruct::getMoment() const{
  51         return Moment;
  52 }
  53
  54
  55 /////////////////////////////////////////////////////////////////
  56 //Implementation of FlapData public methods
  57 ////////////////////////////////////////////////////////////////
  58
  59 FlapData::FlapData(){
  60         liftTable=NULL;
  61         thrustTable=NULL;
  62         momentTable=NULL;
  63         inertiaTable=NULL;
  64
  65         alphaArray=NULL;
  66         speedArray=NULL;
  67         freqArray=NULL;
  68         phiArray=NULL;
  69
  70         lastAlphaIndex=0;
  71         lastSpeedIndex=0;
  72         lastPhiIndex=0;
  73         lastFreqIndex=0;
  74 }
  75
  76 //A constructor that takes a file name:
  77 //Opens that file and fills all the arrays from it
  78 //sets the guesses to zero for the speed and halfway
  79 //along the array for the alpha and frequency
  80 //All it does is call init
  81
  82 FlapData::FlapData(const char* filename){
  83   //  printf("init flapdata\n");
  84         init(filename);
  85         lastAlphaIndex=0;
  86         lastSpeedIndex=0;
  87         lastPhiIndex=0;
  88         lastFreqIndex=0;
  89 }
  90
  91 //The destructor:
  92 //Frees all memory associated with this object
  93 FlapData::~FlapData(){
  94   //  printf("deleting flapdata\n");
  95         delete liftTable;
  96         delete thrustTable;
  97         delete momentTable;
  98         delete inertiaTable;
  99         delete alphaArray;
 100         delete speedArray;
 101         delete freqArray;
 102         delete phiArray;
 103 }
 104
 105 //An initialization function that does the same thing
 106 //as the second constructor
 107 //returns zero if it was successful
 108 int FlapData::init(const char* filename){
 109
 110         ifstream* f=new ifstream(filename);     //open file for reading in text (ascii) mode
 111         if (f==NULL) {  //file open error
 112                 return(1);
 113         }
 114         if(readIn(f)){ //read the file, if there's a problem
 115                 return(2);
 116         }
 117         delete f;
 118         return 0;
 119   //close the file, return the success of the file close
 120 }
 121
 122 //A function that returns the interpolated values
 123 //for all four associated numbers
 124 //given the angle of attack, speed, and flapping frequency
 125 flapStruct FlapData::flapper(double alpha, double speed, double freq, double phi){
 126
 127         flapStruct results;
 128         int i,j,k,l;
 129         double lift,thrust,moment,inertia;
 130         if(speed<speedArray[0]){
 131                 speed=speedArray[0];
 132         }
 133         if(speed>speedArray[speedLength-1]){
 134                 speed=speedArray[speedLength-1];
 135         }
 136         if(alpha<alphaArray[0]){
 137                 alpha=alphaArray[0];
 138         }
 139         if(alpha>alphaArray[alphaLength-1]){
 140                 alpha=alphaArray[alphaLength-1];
 141         }
 142         i=findIndex(alphaArray,alphaLength,alpha,lastAlphaIndex);
 143         j=findIndex(speedArray,speedLength,speed,lastSpeedIndex);
 144         k=findIndex(freqArray,freqLength,freq,lastFreqIndex);
 145         l=findIndex(phiArray,phiLength,phi,lastPhiIndex);
 146
 147         lift=interpolate(liftTable, i, j, k, l, alpha, speed, freq, phi);
 148         thrust=interpolate(thrustTable, i, j, k, l, alpha, speed, freq, phi);
 149         moment=interpolate(momentTable, i, j, k, l, alpha, speed, freq, phi);
 150         inertia=interpolate(inertiaTable, i, j, k, l, alpha, speed, freq, phi);
 151         results=flapStruct(lift,thrust,moment,inertia);
 152         return results;
 153 }
 154
 155 //////////////////////////////////////////////////////////////////
 156 //Implementation of private  FlapData methods
 157 //////////////////////////////////////////////////////////////////
 158
 159 //A function that returns an index i such that
 160 // array[i] < value < array[i+1]
 161 //The function returns -1 if
 162 // (value < array[0]) OR (value > array[n-1])
 163 //(i.e. the value is not within the bounds of the array)
 164 //It performs a linear search starting at guess
 165 int FlapData::findIndex(double array[], double n, double value, int i){
 166
 167         while(value<array[i]){  //less than the lower end of interval i
 168                 if(i==0){            //if we're at the start of the array
 169                         return(-1);                     //there's a problem
 170                 }
 171                 i--;                 //otherwise move to the next lower interval
 172         }
 173         while(value>array[i+1]){        //more than the higher end of interval i
 174                 if(i==n-1){             //if we're at the end of the array
 175                         return(-1);                             //there's a problem
 176                 }
 177                 i++;                    //otherwise move to the next higher interval
 178         }
 179 //        errmsg("In findIndex: array[" << i << "]= " << array[i] << "<=" << value << "<= array[" << (i+1) << "]=" << array[i+1]);
 180         return(i);
 181 }
 182
 183 //A function that performs a linear interpolation based on the
 184 //eight points surrounding the value required
 185 double FlapData::interpolate(double**** table, int i, int j, int k, int l, double alpha, double speed, double freq, double phi){
 186 //        errmsg("\t\t\t\t\t\t\t\tGetting Values");
 187         double f0000=table[i][j][k][l];
 188         double f0001=table[i][j][k][l+1];
 189         double f0010=table[i][j][k+1][l];
 190         double f0011=table[i][j][k+1][l+1];
 191         double f0100=table[i][j+1][k][l];
 192         double f0101=table[i][j+1][k][l+1];
 193         double f0110=table[i][j+1][k+1][l];
 194         double f0111=table[i][j+1][k+1][l+1];
 195         double f1000=table[i+1][j][k][l];
 196         double f1001=table[i+1][j][k][l+1];
 197         double f1010=table[i+1][j][k+1][l];
 198         double f1011=table[i+1][j][k+1][l+1];
 199         double f1100=table[i+1][j+1][k][l];
 200         double f1101=table[i+1][j+1][k][l+1];
 201         double f1110=table[i+1][j+1][k+1][l];
 202         double f1111=table[i+1][j+1][k+1][l+1];
 203
 204   //      errmsg("\t\t\t\t\t\t\t\t1st pass (3)");
 205     //    errmsg("phi[" << l << "]=" << phiArray[l] << "; phi[" << (l+1) <<"]=" << phiArray[l+1]);
 206       //  errmsg("Finding " << phi <<endl;
 207         double f000=interpolate(phiArray[l],f0000,phiArray[l+1],f0001,phi);
 208         double f001=interpolate(phiArray[l],f0010,phiArray[l+1],f0011,phi);
 209         double f010=interpolate(phiArray[l],f0100,phiArray[l+1],f0101,phi);
 210         double f011=interpolate(phiArray[l],f0110,phiArray[l+1],f0111,phi);
 211         double f100=interpolate(phiArray[l],f1000,phiArray[l+1],f1001,phi);
 212         double f101=interpolate(phiArray[l],f1010,phiArray[l+1],f1011,phi);
 213         double f110=interpolate(phiArray[l],f1100,phiArray[l+1],f1101,phi);
 214         double f111=interpolate(phiArray[l],f1110,phiArray[l+1],f1111,phi);
 215
 216 //        errmsg("\t\t\t\t\t\t\t\t2nd pass (2)");
 217         double f00=interpolate(freqArray[k],f000,freqArray[k+1],f001,freq);
 218         double f01=interpolate(freqArray[k],f010,freqArray[k+1],f011,freq);
 219         double f10=interpolate(freqArray[k],f100,freqArray[k+1],f101,freq);
 220         double f11=interpolate(freqArray[k],f110,freqArray[k+1],f111,freq);
 221
 222   //      errmsg("\t\t\t\t\t\t\t\t3rd pass (1)");
 223         double f0=interpolate(speedArray[j],f00,speedArray[j+1],f01,speed);
 224         double f1=interpolate(speedArray[j],f10,speedArray[j+1],f11,speed);
 225
 226     //    errmsg("\t\t\t\t\t\t\t\t4th pass (0)");
 227         double f=interpolate(alphaArray[i],f0,alphaArray[i+1],f1,alpha);
 228         return(f);
 229 }
 230
 231 //A function that performs a linear interpolation based
 232 //on the two nearest points
 233 double FlapData::interpolate(double x0, double y0, double x1, double y1, double x){
 234         double slope,y;
 235         assert(x1!=x0);
 236         slope=(y1-y0)/(x1-x0);
 237         y=y0+slope*(x-x0);
 238         return y;
 239 }
 240
 241 //A function called by init that reads in the file
 242 //of the correct format and stores it in the arrays and tables
 243 int FlapData::readIn (ifstream* f){
 244
 245         int i,j,k,l;
 246         int count=0;
 247         char numstr[200];
 248
 249         f->getline(numstr,200);
 250         sscanf(numstr,"%d,%d,%d,%d",&alphaLength,&speedLength,&freqLength,&phiLength);
 251
 252         //Check to see if the first line is 0 0 0 0
 253         //If so, tell user to download data file
 254         //Quits FlightGear
 255         if (alphaLength==0 && speedLength==0 && freqLength==0 && phiLength==0)
 256           uiuc_warnings_errors(7,"");
 257
 258         alphaArray=new double[alphaLength];
 259         speedArray=new double[speedLength];
 260         freqArray=new double[freqLength];
 261         phiArray=new double[phiLength];
 262
 263         for(i=0;i<alphaLength;i++){
 264                 f->get(numstr,20,',');
 265                 sscanf(numstr,"%lf",&alphaArray[i]);
 266                 f->get();
 267         }
 268         f->get();
 269         for(i=0;i<speedLength;i++){
 270                 f->get(numstr,20,',');
 271                 sscanf(numstr,"%lf",&speedArray[i]);
 272                 f->get();
 273         }
 274         f->get();
 275         for(i=0;i<freqLength;i++){
 276                 f->get(numstr,20,',');
 277                 sscanf(numstr,"%lf",&freqArray[i]);
 278                 f->get();
 279         }
 280         f->get();
 281         for(i=0;i<phiLength;i++){
 282                 f->get(numstr,20,',');
 283                 sscanf(numstr,"%lf",&phiArray[i]);
 284                 f->get();
 285         }
 286         f->get();
 287         liftTable=new double***[alphaLength];
 288         thrustTable=new double***[alphaLength];
 289         momentTable=new double***[alphaLength];
 290         inertiaTable=new double***[alphaLength];
 291         for(i=0;i<alphaLength;i++){
 292                 liftTable[i]=new double**[speedLength];
 293                 thrustTable[i]=new double**[speedLength];
 294                 momentTable[i]=new double**[speedLength];
 295                 inertiaTable[i]=new double**[speedLength];
 296                 for(j=0;j<speedLength;j++){
 297                         liftTable[i][j]=new double*[freqLength];
 298                         thrustTable[i][j]=new double*[freqLength];
 299                         momentTable[i][j]=new double*[freqLength];
 300                         inertiaTable[i][j]=new double*[freqLength];
 301                         for(k=0;k<freqLength;k++){
 302                                 assert((liftTable[i][j][k]=new double[phiLength])!=NULL);
 303                                 assert((thrustTable[i][j][k]=new double[phiLength])!=NULL);
 304                                 assert((momentTable[i][j][k]=new double[phiLength])!=NULL);
 305                                 assert((inertiaTable[i][j][k]=new double[phiLength])!=NULL);
 306                         }
 307                 }
 308         }
 309
 310         for(i=0;i<alphaLength;i++){
 311                 for(j=0;j<speedLength;j++){
 312                         for(k=0;k<freqLength;k++){
 313                                 for(l=0;l<phiLength;l++){
 314                                         f->getline(numstr,200);
 315                                         sscanf(numstr,"%lf %lf %lf %lf",&liftTable[i][j][k][l],&thrustTable[i][j][k][l],&momentTable[i][j][k][l],&inertiaTable[i][j][k][l]);
 316                                 }
 317                         }
 318                 }
 319         }
 320         return 0;
 321 };
 322
 323 //#endif