src/WeatherCM/FGPhysicalProperties.cpp

   1 /*****************************************************************************
   2
   3  Module:       FGPhysicalProperties.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 Initialice the FGPhysicalProperties struct to something sensible(?)
  30
  31 HISTORY
  32 ------------------------------------------------------------------------------
  33 29.05.1999 Christian Mayer      Created
  34 16.06.1999 Durk Talsma          Portability for Linux
  35 20.06.1999 Christian Mayer      added lots of consts
  36 11.10.1999 Christian Mayer      changed set<> to map<> on Bernie Bright's
  37                                 suggestion
  38 19.10.1999 Christian Mayer      change to use PLIB's sg instead of Point[2/3]D
  39                                 and lots of wee code cleaning
  40 *****************************************************************************/
  41
  42 /****************************************************************************/
  43 /* INCLUDES                                                                 */
  44 /****************************************************************************/
  45 #include "FGPhysicalProperties.h"
  46 #include "FGWeatherDefs.h"
  47
  48 #include "FGWeatherUtils.h"
  49
  50 /****************************************************************************/
  51 /********************************** CODE ************************************/
  52 /****************************************************************************/
  53 FGPhysicalProperties::FGPhysicalProperties()
  54 {
  55     sgVec3 zero;
  56     sgZeroVec3( zero );
  57     /************************************************************************/
  58     /* This standart constructor fills the class with a standard weather    */
  59     /************************************************************************/
  60
  61     Wind[-1000.0] = FGWindItem(zero);               //no Wind by default
  62     Wind[10000.0] = FGWindItem(zero);               //no Wind by default
  63
  64     Turbulence[-1000.0] = FGTurbulenceItem(zero);   //no Turbulence by default
  65     Turbulence[10000.0] = FGTurbulenceItem(zero);   //no Turbulence by default
  66
  67     //Initialice with the CINA atmosphere
  68     Temperature[    0.0] = +15.0 + 273.16;
  69     Temperature[11000.0] = -56.5 + 273.16;
  70     Temperature[20000.0] = -56.5 + 273.16;
  71
  72     AirPressure = FGAirPressureItem(101325.0);
  73
  74     VaporPressure[-1000.0] = FG_WEATHER_DEFAULT_VAPORPRESSURE;   //in Pa (I *only* accept SI!)
  75     VaporPressure[10000.0] = FG_WEATHER_DEFAULT_VAPORPRESSURE;   //in Pa (I *only* accept SI!)
  76
  77     //Clouds.insert(FGCloudItem())    => none
  78     SnowRainIntensity = 0.0;
  79     snowRainType = Rain;
  80     LightningProbability = 0.0;
  81 }
  82
  83 /*
  84 The Methods:
  85
  86     WeatherPrecision getWind_x( int number ) const;
  87     WeatherPrecision getWind_y( int number ) const;
  88     WeatherPrecision getWind_z( int number ) const;
  89     WeatherPrecision getWind_a( int number ) const;
  90     void setWind_x( int number, WeatherPrecision x);
  91     void setWind_y( int number, WeatherPrecision y);
  92     void setWind_z( int number, WeatherPrecision z);
  93     void setWind_a( int number, WeatherPrecision a);
  94     WeatherPrecision getTurbulence_x( int number ) const;
  95     WeatherPrecision getTurbulence_y( int number ) const;
  96     WeatherPrecision getTurbulence_z( int number ) const;
  97     WeatherPrecision getTurbulence_a( int number ) const;
  98     void setTurbulence_x( int number, WeatherPrecision x);
  99     void setTurbulence_y( int number, WeatherPrecision y);
 100     void setTurbulence_z( int number, WeatherPrecision z);
 101     void setTurbulence_a( int number, WeatherPrecision a);
 102     WeatherPrecision getTemperature_x( int number ) const;
 103     WeatherPrecision getTemperature_a( int number ) const;
 104     void setTemperature_x( int number, WeatherPrecision x);
 105     void setTemperature_a( int number, WeatherPrecision a);
 106     WeatherPrecision getVaporPressure_x( int number ) const;
 107     WeatherPrecision getVaporPressure_a( int number ) const;
 108     void setVaporPressure_x( int number, WeatherPrecision x);
 109     void setVaporPressure_a( int number, WeatherPrecision a);
 110
 111 are in the extra file FGPhysicalProperties_bind.cpp
 112 */
 113
 114 unsigned int FGPhysicalProperties::getNumberOfCloudLayers(void) const
 115 {
 116     return Clouds.size();
 117 }
 118
 119 FGCloudItem FGPhysicalProperties::getCloudLayer(unsigned int nr) const
 120 {
 121     map<WeatherPrecision,FGCloudItem>::const_iterator CloudsIt = Clouds.begin();
 122
 123     //set the iterator to the 'nr'th entry
 124     for (; nr > 0; nr--)
 125         CloudsIt++;
 126
 127     return CloudsIt->second;
 128 }
 129
 130 ostream& operator<< ( ostream& out, const FGPhysicalProperties2D& p )
 131 {
 132     typedef map<FGPhysicalProperties::Altitude, FGWindItem      >::const_iterator wind_iterator;
 133     typedef map<FGPhysicalProperties::Altitude, FGTurbulenceItem>::const_iterator turbulence_iterator;
 134     typedef map<FGPhysicalProperties::Altitude, WeatherPrecision>::const_iterator scalar_iterator;
 135
 136     out << "Position: (" << p.p[0] << ", " << p.p[1] << ", " << p.p[2] << ")\n";
 137
 138     out << "Stored Wind: ";
 139     for (wind_iterator       WindIt = p.Wind.begin();
 140                              WindIt != p.Wind.end();
 141                              WindIt++)
 142         out << "(" << WindIt->second.x() << ", " << WindIt->second.y() << ", " << WindIt->second.z() << ") m/s at (" << WindIt->first << ") m; ";
 143     out << "\n";
 144
 145     out << "Stored Turbulence: ";
 146     for (turbulence_iterator TurbulenceIt = p.Turbulence.begin();
 147                              TurbulenceIt != p.Turbulence.end();
 148                              TurbulenceIt++)
 149         out << "(" << TurbulenceIt->second.x() << ", " << TurbulenceIt->second.y() << ", " << TurbulenceIt->second.z() << ") m/s at (" << TurbulenceIt->first << ") m; ";
 150     out << "\n";
 151
 152     out << "Stored Temperature: ";
 153     for (scalar_iterator     TemperatureIt = p.Temperature.begin();
 154                              TemperatureIt != p.Temperature.end();
 155                              TemperatureIt++)
 156         out << Kelvin2Celsius(TemperatureIt->second) << " degC at " << TemperatureIt->first << "m; ";
 157     out << "\n";
 158
 159     out << "Stored AirPressure: ";
 160     out << p.AirPressure.getValue()/100.0 << " hPa at " << 0.0 << "m; ";
 161     out << "\n";
 162
 163     out << "Stored VaporPressure: ";
 164     for (scalar_iterator     VaporPressureIt = p.VaporPressure.begin();
 165                              VaporPressureIt != p.VaporPressure.end();
 166                              VaporPressureIt++)
 167         out << VaporPressureIt->second/100.0 << " hPa at " << VaporPressureIt->first << "m; ";
 168     out << "\n";
 169
 170     return out << "\n";
 171 }
 172
 173
 174 inline double F(const WeatherPrecision factor, const WeatherPrecision a, const WeatherPrecision b, const WeatherPrecision r, const WeatherPrecision x)
 175 {
 176     const double c = 1.0 / (-b + a * r);
 177     return factor * c * ( 1.0 / (r + x) + a * c * log(fabs((r + x) * (b + a * x))) );
 178 }
 179
 180 WeatherPrecision FGPhysicalProperties::AirPressureAt(const WeatherPrecision x) const
 181 {
 182     const double rho0 = (AirPressure.getValue()*FG_WEATHER_DEFAULT_AIRDENSITY*FG_WEATHER_DEFAULT_TEMPERATURE)/(TemperatureAt(0)*FG_WEATHER_DEFAULT_AIRPRESSURE);
 183     const double G = 6.673e-11;    //Gravity; in m^3 kg^-1 s^-2
 184     const double m = 5.977e24;      //mass of the earth in kg
 185     const double r = 6368e3;        //radius of the earth in metres
 186     const double factor = -(rho0 * TemperatureAt(0) * G * m) / AirPressure.getValue();
 187
 188     double a, b, FF = 0.0;
 189
 190     //ok, integrate from 0 to a now.
 191     if (Temperature.size() < 2)
 192     {   //take care of the case that there aren't enough points
 193         //actually this should be impossible...
 194
 195         if (Temperature.size() == 0)
 196         {
 197             cerr << "ERROR in FGPhysicalProperties: Air pressure at " << x << " metres altiude requested,\n";
 198             cerr << "      but there isn't enough data stored! No temperature is aviable!\n";
 199             return FG_WEATHER_DEFAULT_AIRPRESSURE;
 200         }
 201
 202         //ok, I've got only one point. So I'm assuming that that temperature is
 203         //the same for all altitudes.
 204         a = 1;
 205         b = TemperatureAt(0);
 206         FF += F(factor, a, b, r, x  );
 207         FF -= F(factor, a, b, r, 0.0);
 208     }
 209     else
 210     {   //I've got at least two entries now
 211
 212         //integrate 'backwards' by integrating the strip ]n,x] first, then ]n-1,n] ... to [0,n-m]
 213
 214         if (x>=0.0)
 215         {
 216             map<WeatherPrecision, WeatherPrecision>::const_iterator temp2 = Temperature.upper_bound(x);
 217             map<WeatherPrecision, WeatherPrecision>::const_iterator temp1 = temp2; temp1--;
 218
 219             if (temp1->first == x)
 220             {   //ignore that interval
 221                 temp1--; temp2--;
 222             }
 223
 224             bool first_pass = true;
 225             while(true)
 226             {
 227                 if (temp2 == Temperature.end())
 228                 {
 229                     //temp2 doesn't exist. So cheat by assuming that the slope is the
 230                     //same as between the two earlier temperatures
 231                     temp1--; temp2--;
 232                     a = (temp2->second - temp1->second)/(temp2->first - temp1->first);
 233                     b = temp1->second - a * temp1->first;
 234                     temp1++; temp2++;
 235                 }
 236                 else
 237                 {
 238                     a = (temp2->second - temp1->second)/(temp2->first - temp1->first);
 239                     b = temp1->second - a * temp1->first;
 240                 }
 241
 242                 if (first_pass)
 243                 {
 244
 245                     FF += F(factor, a, b, r, x);
 246                     first_pass = false;
 247                 }
 248                 else
 249                 {
 250                     FF += F(factor, a, b, r, temp2->first);
 251                 }
 252
 253                 if (temp1->first>0.0)
 254                 {
 255                     FF -= F(factor, a, b, r, temp1->first);
 256                     temp1--; temp2--;
 257                 }
 258                 else
 259                 {
 260                     FF -= F(factor, a, b, r, 0.0);
 261                     return AirPressure.getValue() * exp(FF);
 262                 }
 263             }
 264         }
 265         else
 266         {   //ok x is smaller than 0.0, so do everything in reverse
 267             map<WeatherPrecision, WeatherPrecision>::const_iterator temp2 = Temperature.upper_bound(x);
 268             map<WeatherPrecision, WeatherPrecision>::const_iterator temp1 = temp2; temp1--;
 269
 270             bool first_pass = true;
 271             while(true)
 272             {
 273                 if (temp2 == Temperature.begin())
 274                 {
 275                     //temp1 doesn't exist. So cheat by assuming that the slope is the
 276                     //same as between the two earlier temperatures
 277                     temp1 = Temperature.begin(); temp2++;
 278                     a = (temp2->second - temp1->second)/(temp2->first - temp1->first);
 279                     b = temp1->second - a * temp1->first;
 280                     temp2--;
 281                 }
 282                 else
 283                 {
 284                     a = (temp2->second - temp1->second)/(temp2->first - temp1->first);
 285                     b = temp1->second - a * temp1->first;
 286                 }
 287
 288                 if (first_pass)
 289                 {
 290
 291                     FF += F(factor, a, b, r, x);
 292                     first_pass = false;
 293                 }
 294                 else
 295                 {
 296                     FF += F(factor, a, b, r, temp2->first);
 297                 }
 298
 299                 if (temp2->first<0.0)
 300                 {
 301                     FF -= F(factor, a, b, r, temp1->first);
 302
 303                     if (temp2 == Temperature.begin())
 304                     {
 305                         temp1 = Temperature.begin(); temp2++;
 306                     }
 307                     else
 308                     {
 309                         temp1++; temp2++;
 310                     }
 311                 }
 312                 else
 313                 {
 314                     FF -= F(factor, a, b, r, 0.0);
 315                     return AirPressure.getValue() * exp(FF);
 316                 }
 317             }
 318         }
 319     }
 320
 321     return AirPressure.getValue() * exp(FF);
 322 }