Date started: 28.05.99
Called by: main program
- ---------- Copyright (C) 1999 Christian Mayer (vader@t-online.de) ----------
+ -------- Copyright (C) 1999 Christian Mayer (fgfs@christianmayer.de) --------
This program is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
29.05.1999 Christian Mayer Created
16.06.1999 Durk Talsma Portability for Linux
20.06.1999 Christian Mayer added lots of consts
+11.10.1999 Christian Mayer changed set<> to map<> on Bernie Bright's
+ suggestion
+19.10.1999 Christian Mayer change to use PLIB's sg instead of Point[2/3]D
+ and lots of wee code cleaning
*****************************************************************************/
/****************************************************************************/
#include "FGPhysicalProperties.h"
#include "FGWeatherDefs.h"
+#include "FGWeatherUtils.h"
+
/****************************************************************************/
/********************************** CODE ************************************/
/****************************************************************************/
FGPhysicalProperties::FGPhysicalProperties()
{
+ sgVec3 zero;
+ sgZeroVec3( zero );
/************************************************************************/
/* This standart constructor fills the class with a standard weather */
/************************************************************************/
- Wind.insert(FGWindItem(-1000.0, Point3D(0.0))); //no Wind by default
- Wind.insert(FGWindItem(10000.0, Point3D(0.0))); //no Wind by default
+ Wind[-1000.0] = FGWindItem(zero); //no Wind by default
+ Wind[10000.0] = FGWindItem(zero); //no Wind by default
- Turbulence.insert(FGTurbulenceItem(-1000.0, Point3D(0.0))); //no Turbulence by default
- Turbulence.insert(FGTurbulenceItem(10000.0, Point3D(0.0))); //no Turbulence by default
+ Turbulence[-1000.0] = FGTurbulenceItem(zero); //no Turbulence by default
+ Turbulence[10000.0] = FGTurbulenceItem(zero); //no Turbulence by default
//Initialice with the CINA atmosphere
- Temperature.insert(FGTemperatureItem( 0.0, (+15.0+273.16)));
- Temperature.insert(FGTemperatureItem(11000.0, (-56.5+273.16)));
- Temperature.insert(FGTemperatureItem(20000.0, (-56.5+273.16)));
+ Temperature[ 0.0] = +15.0 + 273.16;
+ Temperature[11000.0] = -56.5 + 273.16;
+ Temperature[20000.0] = -56.5 + 273.16;
AirPressure = FGAirPressureItem(101325.0);
- VaporPressure.insert(FGVaporPressureItem( 0.0, FG_WEATHER_DEFAULT_VAPORPRESSURE)); //in Pa (I *only* accept SI!)
- VaporPressure.insert(FGVaporPressureItem(10000.0, FG_WEATHER_DEFAULT_VAPORPRESSURE)); //in Pa (I *only* accept SI!)
+ VaporPressure[-1000.0] = FG_WEATHER_DEFAULT_VAPORPRESSURE; //in Pa (I *only* accept SI!)
+ VaporPressure[10000.0] = FG_WEATHER_DEFAULT_VAPORPRESSURE; //in Pa (I *only* accept SI!)
//Clouds.insert(FGCloudItem()) => none
SnowRainIntensity = 0.0;
LightningProbability = 0.0;
}
+/*
+The Methods:
+
+ WeatherPrecision getWind_x( int number ) const;
+ WeatherPrecision getWind_y( int number ) const;
+ WeatherPrecision getWind_z( int number ) const;
+ WeatherPrecision getWind_a( int number ) const;
+ void setWind_x( int number, WeatherPrecision x);
+ void setWind_y( int number, WeatherPrecision y);
+ void setWind_z( int number, WeatherPrecision z);
+ void setWind_a( int number, WeatherPrecision a);
+ WeatherPrecision getTurbulence_x( int number ) const;
+ WeatherPrecision getTurbulence_y( int number ) const;
+ WeatherPrecision getTurbulence_z( int number ) const;
+ WeatherPrecision getTurbulence_a( int number ) const;
+ void setTurbulence_x( int number, WeatherPrecision x);
+ void setTurbulence_y( int number, WeatherPrecision y);
+ void setTurbulence_z( int number, WeatherPrecision z);
+ void setTurbulence_a( int number, WeatherPrecision a);
+ WeatherPrecision getTemperature_x( int number ) const;
+ WeatherPrecision getTemperature_a( int number ) const;
+ void setTemperature_x( int number, WeatherPrecision x);
+ void setTemperature_a( int number, WeatherPrecision a);
+ WeatherPrecision getVaporPressure_x( int number ) const;
+ WeatherPrecision getVaporPressure_a( int number ) const;
+ void setVaporPressure_x( int number, WeatherPrecision x);
+ void setVaporPressure_a( int number, WeatherPrecision a);
+
+are in the extra file FGPhysicalProperties_bind.cpp
+*/
+
+unsigned int FGPhysicalProperties::getNumberOfCloudLayers(void) const
+{
+ return Clouds.size();
+}
+
+FGCloudItem FGPhysicalProperties::getCloudLayer(unsigned int nr) const
+{
+ map<WeatherPrecision,FGCloudItem>::const_iterator CloudsIt = Clouds.begin();
+
+ //set the iterator to the 'nr'th entry
+ for (; nr > 0; nr--)
+ CloudsIt++;
+
+ return CloudsIt->second;
+}
+
+ostream& operator<< ( ostream& out, const FGPhysicalProperties2D& p )
+{
+ typedef map<FGPhysicalProperties::Altitude, FGWindItem >::const_iterator wind_iterator;
+ typedef map<FGPhysicalProperties::Altitude, FGTurbulenceItem>::const_iterator turbulence_iterator;
+ typedef map<FGPhysicalProperties::Altitude, WeatherPrecision>::const_iterator scalar_iterator;
+
+ out << "Position: (" << p.p[0] << ", " << p.p[1] << ", " << p.p[2] << ")\n";
+
+ out << "Stored Wind: ";
+ for (wind_iterator WindIt = p.Wind.begin();
+ WindIt != p.Wind.end();
+ WindIt++)
+ out << "(" << WindIt->second.x() << ", " << WindIt->second.y() << ", " << WindIt->second.z() << ") m/s at (" << WindIt->first << ") m; ";
+ out << "\n";
+
+ out << "Stored Turbulence: ";
+ for (turbulence_iterator TurbulenceIt = p.Turbulence.begin();
+ TurbulenceIt != p.Turbulence.end();
+ TurbulenceIt++)
+ out << "(" << TurbulenceIt->second.x() << ", " << TurbulenceIt->second.y() << ", " << TurbulenceIt->second.z() << ") m/s at (" << TurbulenceIt->first << ") m; ";
+ out << "\n";
+
+ out << "Stored Temperature: ";
+ for (scalar_iterator TemperatureIt = p.Temperature.begin();
+ TemperatureIt != p.Temperature.end();
+ TemperatureIt++)
+ out << Kelvin2Celsius(TemperatureIt->second) << " degC at " << TemperatureIt->first << "m; ";
+ out << "\n";
+
+ out << "Stored AirPressure: ";
+ out << p.AirPressure.getValue()/100.0 << " hPa at " << 0.0 << "m; ";
+ out << "\n";
+
+ out << "Stored VaporPressure: ";
+ for (scalar_iterator VaporPressureIt = p.VaporPressure.begin();
+ VaporPressureIt != p.VaporPressure.end();
+ VaporPressureIt++)
+ out << VaporPressureIt->second/100.0 << " hPa at " << VaporPressureIt->first << "m; ";
+ out << "\n";
+
+ return out << "\n";
+}
+
+inline double F(const WeatherPrecision factor, const WeatherPrecision a, const WeatherPrecision b, const WeatherPrecision r, const WeatherPrecision x)
+{
+ const double c = 1.0 / (-b + a * r);
+ return factor * c * ( 1.0 / (r + x) + a * c * log(fabs((r + x) * (b + a * x))) );
+}
+WeatherPrecision FGPhysicalProperties::AirPressureAt(const WeatherPrecision x) const
+{
+ const double rho0 = (AirPressure.getValue()*FG_WEATHER_DEFAULT_AIRDENSITY*FG_WEATHER_DEFAULT_TEMPERATURE)/(TemperatureAt(0)*FG_WEATHER_DEFAULT_AIRPRESSURE);
+ const double G = 6.673e-11; //Gravity; in m^3 kg^-1 s^-2
+ const double m = 5.977e24; //mass of the earth in kg
+ const double r = 6368e3; //radius of the earth in metres
+ const double factor = -(rho0 * TemperatureAt(0) * G * m) / AirPressure.getValue();
+
+ double a, b, FF = 0.0;
+
+ //ok, integrate from 0 to a now.
+ if (Temperature.size() < 2)
+ { //take care of the case that there aren't enough points
+ //actually this should be impossible...
+
+ if (Temperature.size() == 0)
+ {
+ cerr << "ERROR in FGPhysicalProperties: Air pressure at " << x << " metres altiude requested,\n";
+ cerr << " but there isn't enough data stored! No temperature is aviable!\n";
+ return FG_WEATHER_DEFAULT_AIRPRESSURE;
+ }
+ //ok, I've got only one point. So I'm assuming that that temperature is
+ //the same for all altitudes.
+ a = 1;
+ b = TemperatureAt(0);
+ FF += F(factor, a, b, r, x );
+ FF -= F(factor, a, b, r, 0.0);
+ }
+ else
+ { //I've got at least two entries now
+
+ //integrate 'backwards' by integrating the strip ]n,x] first, then ]n-1,n] ... to [0,n-m]
+
+ if (x>=0.0)
+ {
+ map<WeatherPrecision, WeatherPrecision>::const_iterator temp2 = Temperature.upper_bound(x);
+ map<WeatherPrecision, WeatherPrecision>::const_iterator temp1 = temp2; temp1--;
+
+ if (temp1->first == x)
+ { //ignore that interval
+ temp1--; temp2--;
+ }
+
+ bool first_pass = true;
+ while(true)
+ {
+ if (temp2 == Temperature.end())
+ {
+ //temp2 doesn't exist. So cheat by assuming that the slope is the
+ //same as between the two earlier temperatures
+ temp1--; temp2--;
+ a = (temp2->second - temp1->second)/(temp2->first - temp1->first);
+ b = temp1->second - a * temp1->first;
+ temp1++; temp2++;
+ }
+ else
+ {
+ a = (temp2->second - temp1->second)/(temp2->first - temp1->first);
+ b = temp1->second - a * temp1->first;
+ }
+
+ if (first_pass)
+ {
+
+ FF += F(factor, a, b, r, x);
+ first_pass = false;
+ }
+ else
+ {
+ FF += F(factor, a, b, r, temp2->first);
+ }
+
+ if (temp1->first>0.0)
+ {
+ FF -= F(factor, a, b, r, temp1->first);
+ temp1--; temp2--;
+ }
+ else
+ {
+ FF -= F(factor, a, b, r, 0.0);
+ return AirPressure.getValue() * exp(FF);
+ }
+ }
+ }
+ else
+ { //ok x is smaller than 0.0, so do everything in reverse
+ map<WeatherPrecision, WeatherPrecision>::const_iterator temp2 = Temperature.upper_bound(x);
+ map<WeatherPrecision, WeatherPrecision>::const_iterator temp1 = temp2; temp1--;
+
+ bool first_pass = true;
+ while(true)
+ {
+ if (temp2 == Temperature.begin())
+ {
+ //temp1 doesn't exist. So cheat by assuming that the slope is the
+ //same as between the two earlier temperatures
+ temp1 = Temperature.begin(); temp2++;
+ a = (temp2->second - temp1->second)/(temp2->first - temp1->first);
+ b = temp1->second - a * temp1->first;
+ temp2--;
+ }
+ else
+ {
+ a = (temp2->second - temp1->second)/(temp2->first - temp1->first);
+ b = temp1->second - a * temp1->first;
+ }
+
+ if (first_pass)
+ {
+
+ FF += F(factor, a, b, r, x);
+ first_pass = false;
+ }
+ else
+ {
+ FF += F(factor, a, b, r, temp2->first);
+ }
+
+ if (temp2->first<0.0)
+ {
+ FF -= F(factor, a, b, r, temp1->first);
+
+ if (temp2 == Temperature.begin())
+ {
+ temp1 = Temperature.begin(); temp2++;
+ }
+ else
+ {
+ temp1++; temp2++;
+ }
+ }
+ else
+ {
+ FF -= F(factor, a, b, r, 0.0);
+ return AirPressure.getValue() * exp(FF);
+ }
+ }
+ }
+ }
+
+ return AirPressure.getValue() * exp(FF);
+}