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 unsigned int FGPhysicalProperties::getNumberOfCloudLayers(void) const
  84 {
  85     return Clouds.size();
  86 }
  87
  88 FGCloudItem FGPhysicalProperties::getCloudLayer(unsigned int nr) const
  89 {
  90     map<WeatherPrecision,FGCloudItem>::const_iterator CloudsIt = Clouds.begin();
  91
  92     //set the iterator to the 'nr'th entry
  93     for (; nr > 0; nr--)
  94         CloudsIt++;
  95
  96     return CloudsIt->second;
  97 }
  98
  99 ostream& operator<< ( ostream& out, const FGPhysicalProperties2D& p )
 100 {
 101     typedef map<FGPhysicalProperties::Altitude, FGWindItem      >::const_iterator wind_iterator;
 102     typedef map<FGPhysicalProperties::Altitude, FGTurbulenceItem>::const_iterator turbulence_iterator;
 103     typedef map<FGPhysicalProperties::Altitude, WeatherPrecision>::const_iterator scalar_iterator;
 104
 105     out << "Position: (" << p.p[0] << ", " << p.p[1] << ", " << p.p[2] << ")\n";
 106
 107     out << "Stored Wind: ";
 108     for (wind_iterator       WindIt = p.Wind.begin();
 109                              WindIt != p.Wind.end();
 110                              WindIt++)
 111         out << "(" << WindIt->second.x() << ", " << WindIt->second.y() << ", " << WindIt->second.z() << ") m/s at (" << WindIt->first << ") m; ";
 112     out << "\n";
 113
 114     out << "Stored Turbulence: ";
 115     for (turbulence_iterator TurbulenceIt = p.Turbulence.begin();
 116                              TurbulenceIt != p.Turbulence.end();
 117                              TurbulenceIt++)
 118         out << "(" << TurbulenceIt->second.x() << ", " << TurbulenceIt->second.y() << ", " << TurbulenceIt->second.z() << ") m/s at (" << TurbulenceIt->first << ") m; ";
 119     out << "\n";
 120
 121     out << "Stored Temperature: ";
 122     for (scalar_iterator     TemperatureIt = p.Temperature.begin();
 123                              TemperatureIt != p.Temperature.end();
 124                              TemperatureIt++)
 125         out << Kelvin2Celsius(TemperatureIt->second) << " degC at " << TemperatureIt->first << "m; ";
 126     out << "\n";
 127
 128     out << "Stored AirPressure: ";
 129     out << p.AirPressure.getValue()/100.0 << " hPa at " << 0.0 << "m; ";
 130     out << "\n";
 131
 132     out << "Stored VaporPressure: ";
 133     for (scalar_iterator     VaporPressureIt = p.VaporPressure.begin();
 134                              VaporPressureIt != p.VaporPressure.end();
 135                              VaporPressureIt++)
 136         out << VaporPressureIt->second/100.0 << " hPa at " << VaporPressureIt->first << "m; ";
 137     out << "\n";
 138
 139     return out << "\n";
 140 }
 141
 142
 143 inline double F(const WeatherPrecision factor, const WeatherPrecision a, const WeatherPrecision b, const WeatherPrecision r, const WeatherPrecision x)
 144 {
 145     const double c = 1.0 / (-b + a * r);
 146     return factor * c * ( 1.0 / (r + x) + a * c * log(fabs((r + x) * (b + a * x))) );
 147 }
 148
 149 WeatherPrecision FGPhysicalProperties::AirPressureAt(const WeatherPrecision x) const
 150 {
 151     const double rho0 = (AirPressure.getValue()*FG_WEATHER_DEFAULT_AIRDENSITY*FG_WEATHER_DEFAULT_TEMPERATURE)/(TemperatureAt(0)*FG_WEATHER_DEFAULT_AIRPRESSURE);
 152     const double G = 6.673e-11;    //Gravity; in m^3 kg^-1 s^-2
 153     const double m = 5.977e24;      //mass of the earth in kg
 154     const double r = 6368e3;        //radius of the earth in metres
 155     const double factor = -(rho0 * TemperatureAt(0) * G * m) / AirPressure.getValue();
 156
 157     double a, b, FF = 0.0;
 158
 159     //ok, integrate from 0 to a now.
 160     if (Temperature.size() < 2)
 161     {   //take care of the case that there aren't enough points
 162         //actually this should be impossible...
 163
 164         if (Temperature.size() == 0)
 165         {
 166             cerr << "ERROR in FGPhysicalProperties: Air pressure at " << x << " metres altiude requested,\n";
 167             cerr << "      but there isn't enough data stored! No temperature is aviable!\n";
 168             return FG_WEATHER_DEFAULT_AIRPRESSURE;
 169         }
 170
 171         //ok, I've got only one point. So I'm assuming that that temperature is
 172         //the same for all altitudes.
 173         a = 1;
 174         b = TemperatureAt(0);
 175         FF += F(factor, a, b, r, x  );
 176         FF -= F(factor, a, b, r, 0.0);
 177     }
 178     else
 179     {   //I've got at least two entries now
 180
 181         //integrate 'backwards' by integrating the strip ]n,x] first, then ]n-1,n] ... to [0,n-m]
 182
 183         if (x>=0.0)
 184         {
 185             map<WeatherPrecision, WeatherPrecision>::const_iterator temp2 = Temperature.upper_bound(x);
 186             map<WeatherPrecision, WeatherPrecision>::const_iterator temp1 = temp2; temp1--;
 187
 188             if (temp1->first == x)
 189             {   //ignore that interval
 190                 temp1--; temp2--;
 191             }
 192
 193             bool first_pass = true;
 194             while(true)
 195             {
 196                 if (temp2 == Temperature.end())
 197                 {
 198                     //temp2 doesn't exist. So cheat by assuming that the slope is the
 199                     //same as between the two earlier temperatures
 200                     temp1--; temp2--;
 201                     a = (temp2->second - temp1->second)/(temp2->first - temp1->first);
 202                     b = temp1->second - a * temp1->first;
 203                     temp1++; temp2++;
 204                 }
 205                 else
 206                 {
 207                     a = (temp2->second - temp1->second)/(temp2->first - temp1->first);
 208                     b = temp1->second - a * temp1->first;
 209                 }
 210
 211                 if (first_pass)
 212                 {
 213
 214                     FF += F(factor, a, b, r, x);
 215                     first_pass = false;
 216                 }
 217                 else
 218                 {
 219                     FF += F(factor, a, b, r, temp2->first);
 220                 }
 221
 222                 if (temp1->first>0.0)
 223                 {
 224                     FF -= F(factor, a, b, r, temp1->first);
 225                     temp1--; temp2--;
 226                 }
 227                 else
 228                 {
 229                     FF -= F(factor, a, b, r, 0.0);
 230                     return AirPressure.getValue() * exp(FF);
 231                 }
 232             }
 233         }
 234         else
 235         {   //ok x is smaller than 0.0, so do everything in reverse
 236             map<WeatherPrecision, WeatherPrecision>::const_iterator temp2 = Temperature.upper_bound(x);
 237             map<WeatherPrecision, WeatherPrecision>::const_iterator temp1 = temp2; temp1--;
 238
 239             bool first_pass = true;
 240             while(true)
 241             {
 242                 if (temp2 == Temperature.begin())
 243                 {
 244                     //temp1 doesn't exist. So cheat by assuming that the slope is the
 245                     //same as between the two earlier temperatures
 246                     temp1 = Temperature.begin(); temp2++;
 247                     a = (temp2->second - temp1->second)/(temp2->first - temp1->first);
 248                     b = temp1->second - a * temp1->first;
 249                     temp2--;
 250                 }
 251                 else
 252                 {
 253                     a = (temp2->second - temp1->second)/(temp2->first - temp1->first);
 254                     b = temp1->second - a * temp1->first;
 255                 }
 256
 257                 if (first_pass)
 258                 {
 259
 260                     FF += F(factor, a, b, r, x);
 261                     first_pass = false;
 262                 }
 263                 else
 264                 {
 265                     FF += F(factor, a, b, r, temp2->first);
 266                 }
 267
 268                 if (temp2->first<0.0)
 269                 {
 270                     FF -= F(factor, a, b, r, temp1->first);
 271
 272                     if (temp2 == Temperature.begin())
 273                     {
 274                         temp1 = Temperature.begin(); temp2++;
 275                     }
 276                     else
 277                     {
 278                         temp1++; temp2++;
 279                     }
 280                 }
 281                 else
 282                 {
 283                     FF -= F(factor, a, b, r, 0.0);
 284                     return AirPressure.getValue() * exp(FF);
 285                 }
 286             }
 287         }
 288     }
 289
 290     return AirPressure.getValue() * exp(FF);
 291 }
 292
 293