X-Git-Url: https://git.mxchange.org/?a=blobdiff_plain;f=src%2FWeatherCM%2FFGPhysicalProperties.cpp;h=636ff130e9b1661f1c3cd5070101943493d120d2;hb=a4e81f4ff075e6a3c0c2ea1b5a29c0bcdfdbc671;hp=df3d1c7b254344b4139ffe9e036ed5a5ca473c39;hpb=275154500802cc9405b5905b3b055e5c4a17b33f;p=flightgear.git diff --git a/src/WeatherCM/FGPhysicalProperties.cpp b/src/WeatherCM/FGPhysicalProperties.cpp index df3d1c7b2..636ff130e 100644 --- a/src/WeatherCM/FGPhysicalProperties.cpp +++ b/src/WeatherCM/FGPhysicalProperties.cpp @@ -45,6 +45,8 @@ HISTORY #include "FGPhysicalProperties.h" #include "FGWeatherDefs.h" +#include "FGWeatherUtils.h" + /****************************************************************************/ /********************************** CODE ************************************/ /****************************************************************************/ @@ -69,7 +71,7 @@ FGPhysicalProperties::FGPhysicalProperties() AirPressure = FGAirPressureItem(101325.0); - VaporPressure[ 0.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 @@ -78,6 +80,37 @@ FGPhysicalProperties::FGPhysicalProperties() 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(); @@ -94,7 +127,196 @@ FGCloudItem FGPhysicalProperties::getCloudLayer(unsigned int nr) const return CloudsIt->second; } +ostream& operator<< ( ostream& out, const FGPhysicalProperties2D& p ) +{ + typedef map::const_iterator wind_iterator; + typedef map::const_iterator turbulence_iterator; + typedef map::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::const_iterator temp2 = Temperature.upper_bound(x); + map::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::const_iterator temp2 = Temperature.upper_bound(x); + map::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); +}