src/WeatherCM/FGGlobalWeatherDatabase.cpp

   1 /*****************************************************************************
   2
   3  Module:       FGGlobalWeatherDatabase.cpp
   4  Author:       Christian Mayer
   5  Date started: 28.05.99
   6  Called by:    main program
   7
   8  -------- Copyright (C) 1999 Christian Mayer (fgfs@christianmayer.de) --------
   9
  10  This program is free software; you can redistribute it and/or modify it under
  11  the terms of the GNU General Public License as published by the Free Software
  12  Foundation; either version 2 of the License, or (at your option) any later
  13  version.
  14
  15  This program is distributed in the hope that it will be useful, but WITHOUT
  16  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  17  FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
  18  details.
  19
  20  You should have received a copy of the GNU General Public License along with
  21  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
  22  Place - Suite 330, Boston, MA  02111-1307, USA.
  23
  24  Further information about the GNU General Public License can also be found on
  25  the world wide web at http://www.gnu.org.
  26
  27 FUNCTIONAL DESCRIPTION
  28 ------------------------------------------------------------------------------
  29 Database for the global weather
  30 This database is only called by the local database and by the weather
  31 simulator driving this database
  32
  33 HISTORY
  34 ------------------------------------------------------------------------------
  35 28.05.1999 Christian Mayer      Created
  36 16.06.1999 Durk Talsma          Portability for Linux
  37 20.06.1999 Christian Mayer      added lots of consts
  38 11.10.1999 Christian Mayer      changed set<> to map<> on Bernie Bright's
  39                                 suggestion
  40 19.10.1999 Christian Mayer      change to use PLIB's sg instead of Point[2/3]D
  41                                 and lots of wee code cleaning
  42 *****************************************************************************/
  43
  44 /****************************************************************************/
  45 /* INCLUDES                                                                 */
  46 /****************************************************************************/
  47 #include "FGGlobalWeatherDatabase.h"
  48
  49 /****************************************************************************/
  50 /********************************** CODE ************************************/
  51 /****************************************************************************/
  52
  53 template<class V>
  54 V triangle_interpolate(const sgVec2& x1, const V& v1, const sgVec2& x2, const V& v2, const sgVec2& x3, const V& v3, const sgVec2& p)
  55 {
  56     /************************************************************************/
  57     /* First I have to solve the two equations. Rewritten they look like:   */
  58     /*                                                                      */
  59     /*  a11 * x1 + a12 * x2 = b1                                            */
  60     /*  a21 * x1 + a22 * x2 = b2                                            */
  61     /*                                                                      */
  62     /* with                                                                 */
  63     /*                                                                      */
  64     /*  a11 = x2[0] - x1[0]     a12 = x3[0] - x1[0]     b1 = p[0] - x1[0]   */
  65     /*  a21 = x2[1] - x1[1]     a22 = x3[1] - x1[1]     b2 = p[1] - x1[1]   */
  66     /*                                                                      */
  67     /* So I can easily get the solution by saying:                          */
  68     /*                                                                      */
  69     /*       | a11 a12 |                                                    */
  70     /*  D  = |         |                                                    */
  71     /*       | a21 a22 |                                                    */
  72     /*                                                                      */
  73     /*       | b1  a12 |             | a11  b1 |                            */
  74     /*       |         |             |         |                            */
  75     /*       | b2  a22 |             | a21  b2 |                            */
  76     /*  x1 = -----------        x2 = -----------                            */
  77     /*            D                       D                                 */
  78     /*                                                                      */
  79     /*  I just need to take care then that D != 0 or I would get no         */
  80     /*  solution or an infinite amount. Both wouildn't be good...           */
  81     /************************************************************************/
  82
  83     float D = (x2[0] - x1[0]) * (x3[1] - x1[1]) - (x2[1] - x1[1]) * (x3[0] - x1[0]);
  84
  85     if (D == 0.0)
  86         return v1;  //BAD THING HAPPENED!!! I should throw an exeption
  87
  88     float x = ( (p [0] - x1[0]) * (x3[1] - x1[1]) - (p [1] - x1[1]) * (x3[0] - x1[0]) ) / D;
  89     float y = ( (x2[0] - x1[0]) * (p [1] - x1[1]) - (x2[1] - x1[1]) * (p [0] - x1[0]) ) / D;
  90
  91     return v1 + x * (v2 - v1) + y * (v3 - v1);
  92 }
  93
  94 /****************************************************************************/
  95 /* Constructor and Destructor                                               */
  96 /****************************************************************************/
  97 FGGlobalWeatherDatabase::FGGlobalWeatherDatabase(const FGGlobalWeatherDatabaseStatus s)
  98 {
  99     DatabaseStatus = s;
 100 }
 101
 102 FGGlobalWeatherDatabase::~FGGlobalWeatherDatabase()
 103 {
 104 }
 105
 106 /****************************************************************************/
 107 /* Get the physical properties on the specified point p                     */
 108 /* do this by interpolating between the 3 closest points                    */
 109 /****************************************************************************/
 110 FGPhysicalProperties FGGlobalWeatherDatabase::get(const sgVec2& p) const
 111 {
 112     WeatherPrecision distance[3];               //store the 3 closest distances
 113     FGPhysicalProperties2DVectorConstIt iterator[3];    //and the coresponding iterators
 114     WeatherPrecision d;
 115
 116     distance[0] = 9.46e15;      //init with a distance that every calculated
 117     distance[1] = 9.46e15;      //distance is guranteed to be shorter as
 118     distance[2] = 9.46e15;      //9.46e15 metres are 1 light year...
 119
 120     for (FGPhysicalProperties2DVectorConstIt it=database.begin(); it!=database.end(); it++)
 121     {   //go through the whole database
 122         d = sgDistanceVec2(it->p, p);
 123
 124         if (d<distance[0])
 125         {
 126             distance[2] = distance[1]; distance[1] = distance[0]; distance[0] = d;
 127             iterator[2] = iterator[1]; iterator[1] = iterator[0]; iterator[0] = it;
 128             //NOTE: The last line causes a warning that an unitialized variable
 129             //is used. You can ignore this warning here.
 130         }
 131         else if (d<distance[1])
 132         {
 133             distance[2] = distance[1]; distance[1] = d;
 134             iterator[2] = iterator[1]; iterator[1] = it;
 135         }
 136         else if (d<distance[2])
 137         {
 138             distance[2] = d;
 139             iterator[2] = it;
 140         }
 141     }
 142
 143     //now I've got the closest entry in xx[0], the 2nd closest in xx[1] and the
 144     //3rd in xx[2];
 145
 146     //interpolate now:
 147     return triangle_interpolate(
 148         iterator[0]->p, (FGPhysicalProperties)*iterator[0],
 149         iterator[1]->p, (FGPhysicalProperties)*iterator[1],
 150         iterator[2]->p, (FGPhysicalProperties)*iterator[2], p);
 151 }
 152
 153 /****************************************************************************/
 154 /* update the database. Since the last call we had dt seconds               */
 155 /****************************************************************************/
 156 void FGGlobalWeatherDatabase::update(const WeatherPrecision dt)
 157 {
 158     // I've got nothing to update here (yet...)
 159 }
 160
 161 /****************************************************************************/
 162 /* Add a physical property on the specified point p                         */
 163 /****************************************************************************/
 164 void FGGlobalWeatherDatabase::add(const sgVec2& p, const FGPhysicalProperties& x)
 165 {
 166     FGPhysicalProperties2D e;
 167
 168     sgCopyVec2(e.p, p);
 169
 170     e.Wind = x.Wind;
 171     e.Turbulence = x.Turbulence;
 172     e.Temperature = x.Temperature;
 173     e.AirPressure = x.AirPressure;
 174     e.VaporPressure = x.VaporPressure;
 175
 176     e.Clouds = x.Clouds;
 177     e.SnowRainIntensity = x.SnowRainIntensity;
 178     e.snowRainType = x.snowRainType;
 179     e.LightningProbability = x.LightningProbability;
 180
 181     database.push_back(e);
 182 }
 183
 184 /****************************************************************************/
 185 /* Change the closest physical property to p. If p is further away than     */
 186 /* tolerance I'm returning false otherwise true                             */
 187 /****************************************************************************/
 188 bool FGGlobalWeatherDatabase::change(const FGPhysicalProperties2D& p, const WeatherPrecision tolerance)
 189 {
 190     for (FGPhysicalProperties2DVectorIt it = database.begin(); it != database.end(); it++)
 191     {
 192         if (sgScalarProductVec2(it->p, p.p) < (tolerance*tolerance))
 193         {   //assume that's my point
 194             (*it) = p;
 195             return true;
 196         }
 197     }
 198
 199     return false;
 200 }
 201
 202 /****************************************************************************/
 203 /* Get all, but at least min, stored point in the circle around p with the  */
 204 /* radius r                                                                 */
 205 /****************************************************************************/
 206 FGPhysicalProperties2DVector FGGlobalWeatherDatabase::getAll(const sgVec2& p, const WeatherPrecision r, const unsigned int min)
 207 {
 208     FGPhysicalProperties2DVector ret_list;
 209
 210     if (  (DatabaseStatus == FGGlobalWeatherDatabase_only_static)
 211         ||(DatabaseStatus == FGGlobalWeatherDatabase_working    ) )
 212     {   //doest it make sense?
 213
 214         FGPhysicalProperties2DVectorIt *it;         //store the closest entries
 215         WeatherPrecision *d;
 216         unsigned int act_it = 0;
 217         int i;
 218
 219         it = new FGPhysicalProperties2DVectorIt[min+1];
 220         d = new WeatherPrecision[min+1];
 221
 222         for (it[0]=database.begin(); it[act_it]!=database.end(); it[act_it]++)
 223         {       //go through the whole database
 224             d[act_it] = sgScalarProductVec2(it[act_it]->p, p);
 225
 226             if (r >= d[act_it])
 227             {   //add it
 228                 ret_list.push_back(*it[act_it]);
 229             }
 230             else
 231             {
 232                 if (act_it>0)
 233                 {   //figure out if this distance belongs to the closest ones
 234                     WeatherPrecision dummy;
 235                     FGPhysicalProperties2DVectorIt dummyIt;
 236
 237                     for (i = act_it++; i >= 0;)
 238                     {
 239                         if (d[i] >= d[--i])
 240                         {
 241                             act_it--;
 242                             break;  //nope => stop
 243                         }
 244
 245                         //swap both
 246                         dummy =d[i]; d[i] = d[i+1]; d[i+1] = dummy;
 247                         dummyIt = it[i]; it[i] = it[i+1]; it[i+1] = dummyIt;
 248                     }
 249                 }
 250             }
 251         }
 252
 253         if (ret_list.size()<min)
 254         {
 255             for(i = 0; (i < (int)(min - ret_list.size())) && (ret_list.size() < database.size()); i++)
 256                 ret_list.push_back(*it[i]);
 257         }
 258
 259         delete d;
 260         delete it;
 261     }
 262
 263     return ret_list;
 264 }
 265