From 210e87ec3a91b15d0edd34d255e8a3b3c654952f Mon Sep 17 00:00:00 2001 From: curt Date: Thu, 23 Dec 1999 16:54:54 +0000 Subject: [PATCH] The WeatherDatabase doesn't need the voronoi code anymore but uses Dave Eberly's spherical interpolation code (found in the Lib/Math directory). So it would be great if you could give him also a place in the thanks file. Changing the WeatherDatabse made actually a heavy internal redesign necessary but no code outside the database is affected (isn't code hiding great?). --- src/WeatherCM/FGAirPressureItem.h | 18 +- src/WeatherCM/FGLocalWeatherDatabase.cpp | 247 ++++----------------- src/WeatherCM/FGLocalWeatherDatabase.h | 77 ++----- src/WeatherCM/FGPhysicalProperties.cpp | 154 ++++++++++++- src/WeatherCM/FGPhysicalProperties.h | 19 +- src/WeatherCM/FGPhysicalProperty.h | 2 +- src/WeatherCM/FGTurbulenceItem.h | 2 +- src/WeatherCM/FGWeatherFeature.h | 4 +- src/WeatherCM/FGWeatherParse.cpp | 101 ++++++--- src/WeatherCM/FGWeatherParse.h | 4 +- src/WeatherCM/FGWeatherVectorWrap.h | 2 +- src/WeatherCM/FGWindItem.h | 2 +- src/WeatherCM/Makefile.am | 3 - src/WeatherCM/air-pressure-explanation.txt | 212 ++++++++++++++++++ 14 files changed, 517 insertions(+), 330 deletions(-) create mode 100644 src/WeatherCM/air-pressure-explanation.txt diff --git a/src/WeatherCM/FGAirPressureItem.h b/src/WeatherCM/FGAirPressureItem.h index 3b4d422f0..5c80eefc4 100644 --- a/src/WeatherCM/FGAirPressureItem.h +++ b/src/WeatherCM/FGAirPressureItem.h @@ -37,6 +37,10 @@ HISTORY suggestion 19.10.1999 Christian Mayer change to use PLIB's sg instead of Point[2/3]D and lots of wee code cleaning +15.12.1999 Christian Mayer changed the air pressure calculation to a much + more realistic formula. But as I need for that + the temperature I moved the code to + FGPhysicalProperties *****************************************************************************/ /****************************************************************************/ @@ -62,25 +66,23 @@ FGAirPressureItem operator-(const FGAirPressureItem& arg); /* CLASS DECLARATION */ /* NOTE: The value stored in 'value' is the air preasure that we'd have at */ /* an altitude of 0.0 The correct airpreasure at the stored altitude */ -/* gets calulated when getValue() is called. */ +/* gets calulated in FGPhyiscalProperties as I need to know the */ +/* temperatures at the different altitudes for that. */ /****************************************************************************/ class FGAirPressureItem { private: - WeatherPrecision value; + WeatherPrecision value; //that's the airpressure at 0 metres protected: public: - FGAirPressureItem(const WeatherPrecision v) {value = v;} + FGAirPressureItem(const WeatherPrecision v) {value = v; } FGAirPressureItem() {value = FG_WEATHER_DEFAULT_AIRPRESSURE;} - WeatherPrecision getValue(const WeatherPrecision& alt) const + WeatherPrecision getValue(void) const { - return (WeatherPrecision)((value / 101325.0) * - ( - 1.01325e5 + alt * (-1.19459535223623e1 + alt * (5.50461110007561e-4 + alt * (-1.13574703113648e-8 + alt * 8.61601726143988e-14))) - )); + return value; }; FGAirPressureItem& operator*=(const WeatherPrecision arg); diff --git a/src/WeatherCM/FGLocalWeatherDatabase.cpp b/src/WeatherCM/FGLocalWeatherDatabase.cpp index 6a91b4508..3e8fadbd1 100644 --- a/src/WeatherCM/FGLocalWeatherDatabase.cpp +++ b/src/WeatherCM/FGLocalWeatherDatabase.cpp @@ -38,6 +38,10 @@ HISTORY suggestion 19.10.1999 Christian Mayer change to use PLIB's sg instead of Point[2/3]D and lots of wee code cleaning +14.12.1999 Christian Mayer Changed the internal structure to use Dave + Eberly's spherical interpolation code. This + stops our dependancy on the (ugly) voronoi + code and simplyfies the code structure a lot. *****************************************************************************/ /****************************************************************************/ @@ -49,7 +53,6 @@ HISTORY #include #include "FGLocalWeatherDatabase.h" -#include "FGVoronoi.h" #include "FGWeatherParse.h" @@ -57,37 +60,6 @@ HISTORY /********************************** CODE ************************************/ /****************************************************************************/ -/****************************************************************************/ -/* return the index (better: ID) of the area with point p */ -/****************************************************************************/ -unsigned int FGLocalWeatherDatabase::AreaWith(const sgVec2& p) const -{ - - for (FGMicroWeatherList::size_type i = 0; i != WeatherAreas.size(); i++) - { - if (WeatherAreas[i].hasPoint(p) == true) - return i+1; - } - - return 0; //nothing found -} - -/****************************************************************************/ -/* make tiles out of points on a 2D plane */ -/****************************************************************************/ -void FGLocalWeatherDatabase::tileLocalWeather(const FGPhysicalProperties2DVector& EntryList) -{ - FGVoronoiInputList input; - - for (FGPhysicalProperties2DVectorConstIt it1 = EntryList.begin(); it1 != EntryList.end(); it1++) - input.push_back(FGVoronoiInput(it1->p, *it1)); - - FGVoronoiOutputList output = Voronoiate(input); - - for (FGVoronoiOutputList::iterator it2 = output.begin(); it2 != output.end(); it2++) - WeatherAreas.push_back(FGMicroWeather(it2->value, it2->boundary)); -} - FGLocalWeatherDatabase* FGLocalWeatherDatabase::theFGLocalWeatherDatabase = 0; FGLocalWeatherDatabase *WeatherDatabase; @@ -98,7 +70,6 @@ void FGLocalWeatherDatabase::init(const WeatherPrecision visibility, const Datab if (theFGLocalWeatherDatabase) { - //FG_LOG( FG_GENERAL, FG_ALERT, "Error: only one local weather allowed" ); cerr << "Error: only one local weather allowed"; exit(-1); } @@ -106,7 +77,7 @@ void FGLocalWeatherDatabase::init(const WeatherPrecision visibility, const Datab setWeatherVisibility(visibility); DatabaseStatus = type; - global = 0; //just get sure... + database = 0; //just get sure... Thunderstorm = false; //I don't need to set theThunderstorm as Thunderstorm == false @@ -115,33 +86,42 @@ void FGLocalWeatherDatabase::init(const WeatherPrecision visibility, const Datab { case use_global: { - global = new FGGlobalWeatherDatabase; //initialize GlobalDatabase - global->setDatabaseStatus(FGGlobalWeatherDatabase_working); - tileLocalWeather(global->getAll(last_known_position, WeatherVisibility, 3)); + cerr << "Error: there's no global database anymore!\n"; + exit(-1); } break; case use_internet: { - FGWeatherParse parsed_data; + FGWeatherParse *parsed_data = new FGWeatherParse(); - parsed_data.input( "weather/current.gz" ); - global = new FGGlobalWeatherDatabase; //initialize GlobalDatabase - global->setDatabaseStatus(FGGlobalWeatherDatabase_working); + parsed_data->input( "weather/current.gz" ); + unsigned int n = parsed_data->stored_stations(); + + sgVec2 *p = new sgVec2 [n]; + FGPhysicalProperties *f = new FGPhysicalProperties[n]; // fill the database - for (unsigned int i = 0; i != parsed_data.stored_stations(); i++) + for (unsigned int i = 0; i < n; i++) { - global->add( parsed_data.getFGPhysicalProperties2D(i) ); - //cerr << parsed_data.getFGPhysicalProperties2D(i); + f[i] = parsed_data->getFGPhysicalProperties(i); + parsed_data->getPosition(i, p[i]); if ( (i%100) == 0) cerr << "."; } - cerr << "\n"; - //and finally tile it - tileLocalWeather(global->getAll(last_known_position, WeatherVisibility, 3)); + // free the memory of the parsed data to ease the required memory + // for the very memory consuming spherical interpolation + delete parsed_data; + + //and finally init the interpolation + cerr << "\nInitialiating Interpolation. (2-3 minutes on a PII-350)\n"; + database = new SphereInterpolate(n, p, f); + + //and free my allocations: + delete[] p; + delete[] f; cerr << "Finished weather init.\n"; } @@ -154,9 +134,11 @@ void FGLocalWeatherDatabase::init(const WeatherPrecision visibility, const Datab case manual: case default_mode: { - - vector emptyList; - WeatherAreas.push_back(FGMicroWeather(FGPhysicalProperties2D(), emptyList)); //in these cases I've only got one tile + double x[2] = {0.0, 0.0}; //make an standard weather that's the same at the whole world + double y[2] = {0.0, 0.0}; //make an standard weather that's the same at the whole world + double z[2] = {1.0, -1.0}; //make an standard weather that's the same at the whole world + FGPhysicalProperties f[2]; //make an standard weather that's the same at the whole world + database = new SphereInterpolate(2,x,y,z,f); } break; @@ -168,13 +150,7 @@ void FGLocalWeatherDatabase::init(const WeatherPrecision visibility, const Datab FGLocalWeatherDatabase::~FGLocalWeatherDatabase() { //Tidying up: - - //delete every stored area - WeatherAreas.erase(WeatherAreas.begin(), WeatherAreas.end()); - - //delete global database if necessary - if (DatabaseStatus == use_global) - delete global; + delete database; } /****************************************************************************/ @@ -182,16 +158,7 @@ FGLocalWeatherDatabase::~FGLocalWeatherDatabase() /****************************************************************************/ void FGLocalWeatherDatabase::reset(const DatabaseWorkingType type) { - //delete global database if necessary - if ((DatabaseStatus == use_global) && (type != use_global)) - delete global; - - DatabaseStatus = type; - if (DatabaseStatus == use_global) - tileLocalWeather(global->getAll(last_known_position, WeatherVisibility, 3)); - - //delete every stored area - WeatherAreas.erase(WeatherAreas.begin(), WeatherAreas.end()); + cerr << "FGLocalWeatherDatabase::reset isn't supported yet\n"; } /****************************************************************************/ @@ -199,45 +166,29 @@ void FGLocalWeatherDatabase::reset(const DatabaseWorkingType type) /****************************************************************************/ void FGLocalWeatherDatabase::update(const WeatherPrecision dt) { - if (DatabaseStatus==use_global) - global->update(dt); + //if (DatabaseStatus==use_global) + // global->update(dt); } void FGLocalWeatherDatabase::update(const sgVec3& p) //position has changed { sgCopyVec3(last_known_position, p); - if ( AreaWith(p) == 0 ) - { //I have moved out of my local area... - - //now I should erase all my areas and get the new ones - //but that takes too long :-( I have to take care about that soon - } - //uncomment this when you are using the GlobalDatabase /* cerr << "****\nupdate(p) inside\n"; cerr << "Parameter: " << p[0] << "/" << p[1] << "/" << p[2] << "\n"; sgVec2 p2d; sgSetVec2( p2d, p[0], p[1] ); - cerr << FGPhysicalProperties2D(global->get(p2d), p2d); + cerr << FGPhysicalProperties2D(get(p2d), p2d); cerr << "****\n"; */ + } void FGLocalWeatherDatabase::update(const sgVec3& p, const WeatherPrecision dt) //time and/or position has changed { sgCopyVec3(last_known_position, p); - - if ( AreaWith(p) == 0 ) - { //I have moved out of my local area... - - //now I should erase all my areas and get the new ones - //but that takes too long :-( I have to take care about that soon - } - - if (DatabaseStatus==use_global) - global->update(dt); } /****************************************************************************/ @@ -245,154 +196,46 @@ void FGLocalWeatherDatabase::update(const sgVec3& p, const WeatherPrecision dt) /****************************************************************************/ FGPhysicalProperty FGLocalWeatherDatabase::get(const sgVec3& p) const { - unsigned int a = AreaWith(p); - if (a != 0) - return WeatherAreas[a-1].get(p[3]); - else //point is outside => ask GlobalWeatherDatabase - return global->get(p); + return FGPhysicalProperty(database->Evaluate(p), p[3]); } FGPhysicalProperties FGLocalWeatherDatabase::get(const sgVec2& p) const { - sgVec3 temp; - sgSetVec3(temp, p[0], p[1], 0.0); - - unsigned int a = AreaWith(temp); - if (a != 0) - return WeatherAreas[a-1].get(); - else //point is outside => ask GlobalWeatherDatabase - return global->get(p); + return database->Evaluate(p); } WeatherPrecision FGLocalWeatherDatabase::getAirDensity(const sgVec3& p) const { - FGPhysicalProperty dummy; - unsigned int a = AreaWith(p); - if (a != 0) - dummy = WeatherAreas[a-1].get(p[3]); - else //point is outside => ask GlobalWeatherDatabase - dummy = global->get(p); + FGPhysicalProperty dummy(database->Evaluate(p), p[3]); return (dummy.AirPressure*FG_WEATHER_DEFAULT_AIRDENSITY*FG_WEATHER_DEFAULT_TEMPERATURE) / (dummy.Temperature*FG_WEATHER_DEFAULT_AIRPRESSURE); } -/****************************************************************************/ -/* Add a weather feature at the point p and surrounding area */ -/****************************************************************************/ -void FGLocalWeatherDatabase::addWind(const WeatherPrecision alt, const sgVec3& x, const sgVec2& p) -{ - unsigned int a = AreaWith(p); - if (a != 0) - WeatherAreas[a-1].addWind(alt, x); -} - -void FGLocalWeatherDatabase::addTurbulence(const WeatherPrecision alt, const sgVec3& x, const sgVec2& p) -{ - unsigned int a = AreaWith(p); - if (a != 0) - WeatherAreas[a-1].addTurbulence(alt, x); -} - -void FGLocalWeatherDatabase::addTemperature(const WeatherPrecision alt, const WeatherPrecision x, const sgVec2& p) -{ - unsigned int a = AreaWith(p); - if (a != 0) - WeatherAreas[a-1].addTemperature(alt, x); -} - -void FGLocalWeatherDatabase::addAirPressure(const WeatherPrecision alt, const WeatherPrecision x, const sgVec2& p) -{ - unsigned int a = AreaWith(p); - if (a != 0) - WeatherAreas[a-1].addAirPressure(alt, x); -} - -void FGLocalWeatherDatabase::addVaporPressure(const WeatherPrecision alt, const WeatherPrecision x, const sgVec2& p) -{ - unsigned int a = AreaWith(p); - if (a != 0) - WeatherAreas[a-1].addVaporPressure(alt, x); -} - -void FGLocalWeatherDatabase::addCloud(const WeatherPrecision alt, const FGCloudItem& x, const sgVec2& p) -{ - unsigned int a = AreaWith(p); - if (a != 0) - WeatherAreas[a-1].addCloud(alt, x); -} void FGLocalWeatherDatabase::setSnowRainIntensity(const WeatherPrecision x, const sgVec2& p) { - unsigned int a = AreaWith(p); - if (a != 0) - WeatherAreas[a-1].setSnowRainIntensity(x); + /* not supported yet */ } void FGLocalWeatherDatabase::setSnowRainType(const SnowRainType x, const sgVec2& p) { - unsigned int a = AreaWith(p); - if (a != 0) - WeatherAreas[a-1].setSnowRainType(x); + /* not supported yet */ } void FGLocalWeatherDatabase::setLightningProbability(const WeatherPrecision x, const sgVec2& p) { - unsigned int a = AreaWith(p); - if (a != 0) - WeatherAreas[a-1].setLightningProbability(x); -} - -void FGLocalWeatherDatabase::addProperties(const FGPhysicalProperties2D& x) -{ - if (DatabaseStatus==use_global) - { - global->add(x); - - //BAD, BAD, BAD thing I'm doing here: I'm adding to the global database a point that - //changes my voronoi diagram but I don't update it! instead I'm changing one local value - //that could be anywhere!! - //This only *might* work when the plane moves so far so fast that the diagram gets newly - //calculated soon... - unsigned int a = AreaWith(x.p); - if (a != 0) - WeatherAreas[a-1].setStoredWeather(x); - } - else - { - unsigned int a = AreaWith(x.p); - if (a != 0) - WeatherAreas[a-1].setStoredWeather(x); - } + /* not supported yet */ } void FGLocalWeatherDatabase::setProperties(const FGPhysicalProperties2D& x) { - if (DatabaseStatus==use_global) - { - global->change(x); - - //BAD, BAD, BAD thing I'm doing here: I'm adding to the global database a point that - //changes my voronoi diagram but I don't update it! Instead I'm changing one local value - //that could be anywhere!! - //This only *might* work when the plane moves so far so fast that the diagram gets newly - //calculated soon... - unsigned int a = AreaWith(x.p); - if (a != 0) - WeatherAreas[a-1].setStoredWeather(x); - } - else - { - unsigned int a = AreaWith(x.p); - if (a != 0) - WeatherAreas[a-1].setStoredWeather(x); - } + /* not supported yet */ } void fgUpdateWeatherDatabase(void) { - //cerr << "FGLocalWeatherDatabase::update()\n"; sgVec3 position; sgSetVec3(position, current_aircraft.fdm_state->get_Latitude(), diff --git a/src/WeatherCM/FGLocalWeatherDatabase.h b/src/WeatherCM/FGLocalWeatherDatabase.h index b9a34bb4d..817ea34b7 100644 --- a/src/WeatherCM/FGLocalWeatherDatabase.h +++ b/src/WeatherCM/FGLocalWeatherDatabase.h @@ -38,6 +38,10 @@ HISTORY suggestion 19.10.1999 Christian Mayer change to use PLIB's sg instead of Point[2/3]D and lots of wee code cleaning +14.12.1999 Christian Mayer Changed the internal structure to use Dave + Eberly's spherical interpolation code. This + stops our dependancy on the (ugly) voronoi + code and simplyfies the code structure a lot. *****************************************************************************/ /****************************************************************************/ @@ -49,21 +53,16 @@ HISTORY /****************************************************************************/ /* INCLUDES */ /****************************************************************************/ -//This is only here for smoother code change. In the end the WD should be clean -//of *any* OpenGL: -#ifdef HAVE_WINDOWS_H -# include -#endif -#include -#include - #include -#include "sg.h" +#include + +#include #include "FGPhysicalProperties.h" -#include "FGGlobalWeatherDatabase.h" -#include "FGMicroWeather.h" +#include "FGPhysicalProperty.h" + + #include "FGWeatherFeature.h" #include "FGWeatherDefs.h" #include "FGThunderstorm.h" @@ -81,10 +80,7 @@ class FGLocalWeatherDatabase { private: protected: - FGGlobalWeatherDatabase *global; //point to the global database - - typedef vector FGMicroWeatherList; - typedef FGMicroWeatherList::iterator FGMicroWeatherListIt; + SphereInterpolate *database; typedef vector pointVector; typedef vector tileVector; @@ -92,10 +88,6 @@ protected: /************************************************************************/ /* make tiles out of points on a 2D plane */ /************************************************************************/ - void tileLocalWeather(const FGPhysicalProperties2DVector& EntryList); - - FGMicroWeatherList WeatherAreas; - WeatherPrecision WeatherVisibility; //how far do I need to simulate the //local weather? Unit: metres sgVec3 last_known_position; @@ -103,23 +95,11 @@ protected: bool Thunderstorm; //is there a thunderstorm near by? FGThunderstorm *theThunderstorm; //pointer to the thunderstorm. - /************************************************************************/ - /* return the index of the area with point p */ - /************************************************************************/ - unsigned int AreaWith(const sgVec2& p) const; - unsigned int AreaWith(const sgVec3& p) const - { - sgVec2 temp; - sgSetVec2(temp, p[0], p[1]); - - return AreaWith(temp); - } - public: static FGLocalWeatherDatabase *theFGLocalWeatherDatabase; enum DatabaseWorkingType { - use_global, //use global database for data + use_global, //use global database for data !!obsolete!! use_internet, //use the weather data that came from the internet manual, //use only user inputs distant, //use distant information, e.g. like LAN when used in @@ -186,19 +166,12 @@ public: /************************************************************************/ /* Add a weather feature at the point p and surrounding area */ /************************************************************************/ - - void addWind (const WeatherPrecision alt, const sgVec3& x, const sgVec2& p); - void addTurbulence (const WeatherPrecision alt, const sgVec3& x, const sgVec2& p); - void addTemperature (const WeatherPrecision alt, const WeatherPrecision x, const sgVec2& p); - void addAirPressure (const WeatherPrecision alt, const WeatherPrecision x, const sgVec2& p); - void addVaporPressure(const WeatherPrecision alt, const WeatherPrecision x, const sgVec2& p); - void addCloud (const WeatherPrecision alt, const FGCloudItem& x, const sgVec2& p); + // !! Adds aren't supported anymore !! void setSnowRainIntensity (const WeatherPrecision x, const sgVec2& p); void setSnowRainType (const SnowRainType x, const sgVec2& p); void setLightningProbability(const WeatherPrecision x, const sgVec2& p); - void addProperties(const FGPhysicalProperties2D& x); //add a property void setProperties(const FGPhysicalProperties2D& x); //change a property /************************************************************************/ @@ -231,34 +204,10 @@ void inline FGLocalWeatherDatabase::setWeatherVisibility(const WeatherPrecision WeatherVisibility = visibility; else WeatherVisibility = MINIMUM_WEATHER_VISIBILITY; - -#if 0 - //This code doesn't belong here as this is the optical visibility and not - //the visibility of the weather database (that should be bigger...). The - //optical visibility should be calculated from the vapor pressure e.g. - //But for the sake of a smoother change from the old way to the new one... - - GLfloat fog_exp_density; - GLfloat fog_exp2_density; - - // for GL_FOG_EXP - fog_exp_density = -log(0.01 / WeatherVisibility); - - // for GL_FOG_EXP2 - fog_exp2_density = sqrt( -log(0.01) ) / WeatherVisibility; - - // Set correct opengl fog density - xglFogf (GL_FOG_DENSITY, fog_exp2_density); - - // FG_LOG( FG_INPUT, FG_DEBUG, "Fog density = " << w->fog_density ); - //cerr << "FGLocalWeatherDatabase::setWeatherVisibility(" << visibility << "):\n"; - //cerr << "Fog density = " << fog_exp_density << "\n"; -#endif } WeatherPrecision inline FGLocalWeatherDatabase::getWeatherVisibility(void) const { - //cerr << "FGLocalWeatherDatabase::getWeatherVisibility() = " << WeatherVisibility << "\n"; return WeatherVisibility; } diff --git a/src/WeatherCM/FGPhysicalProperties.cpp b/src/WeatherCM/FGPhysicalProperties.cpp index 044e6805b..966b321c0 100644 --- a/src/WeatherCM/FGPhysicalProperties.cpp +++ b/src/WeatherCM/FGPhysicalProperties.cpp @@ -71,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 @@ -126,7 +126,7 @@ ostream& operator<< ( ostream& out, const FGPhysicalProperties2D& p ) out << "\n"; out << "Stored AirPressure: "; - out << p.AirPressure.getValue(0)/100.0 << " hPa at " << 0.0 << "m; "; + out << p.AirPressure.getValue()/100.0 << " hPa at " << 0.0 << "m; "; out << "\n"; out << "Stored VaporPressure: "; @@ -140,4 +140,154 @@ ostream& operator<< ( ostream& out, const FGPhysicalProperties2D& p ) } +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(abs((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); +} + diff --git a/src/WeatherCM/FGPhysicalProperties.h b/src/WeatherCM/FGPhysicalProperties.h index f6e8f23d9..6d88b6e81 100644 --- a/src/WeatherCM/FGPhysicalProperties.h +++ b/src/WeatherCM/FGPhysicalProperties.h @@ -38,6 +38,10 @@ HISTORY suggestion 19.10.1999 Christian Mayer change to use PLIB's sg instead of Point[2/3]D and lots of wee code cleaning +15.12.1999 Christian Mayer changed the air pressure calculation to a much + more realistic formula. But as I need for that + the temperature I moved the code to + FGPhysicalProperties *****************************************************************************/ /****************************************************************************/ @@ -63,7 +67,7 @@ HISTORY #include #include -#include "sg.h" +#include #include "FGWeatherDefs.h" @@ -107,7 +111,7 @@ public: void WindAt (sgVec3 ret, const WeatherPrecision a) const; void TurbulenceAt (sgVec3 ret, const WeatherPrecision a) const; WeatherPrecision TemperatureAt (const WeatherPrecision a) const; - WeatherPrecision AirPressureAt (const WeatherPrecision a) const; + WeatherPrecision AirPressureAt (const WeatherPrecision x) const; //x is used here instead of a on purpose WeatherPrecision VaporPressureAt(const WeatherPrecision a) const; //for easier access to the cloud stuff: @@ -214,7 +218,7 @@ inline FGPhysicalProperties& FGPhysicalProperties::operator += (const FGPhysical if (!Temperature.insert(*TemperatureIt).second) Temperature[TemperatureIt->first] += TemperatureIt->second; - AirPressure += p.AirPressure.getValue(0.0); + AirPressure += p.AirPressure.getValue(); for (scalar_iterator VaporPressureIt = p.VaporPressure.begin(); VaporPressureIt != p.VaporPressure.end(); @@ -249,7 +253,7 @@ inline FGPhysicalProperties& FGPhysicalProperties::operator -= (const FGPhysical if (!Temperature.insert( make_pair(TemperatureIt->first, -TemperatureIt->second) ).second) Temperature[TemperatureIt->first] -= TemperatureIt->second; - AirPressure -= p.AirPressure.getValue(0.0); + AirPressure -= p.AirPressure.getValue(); for (scalar_iterator VaporPressureIt = p.VaporPressure.begin(); VaporPressureIt != p.VaporPressure.end(); @@ -261,7 +265,6 @@ inline FGPhysicalProperties& FGPhysicalProperties::operator -= (const FGPhysical return *this; } - inline void FGPhysicalProperties::WindAt(sgVec3 ret, const WeatherPrecision a) const { typedef map::const_iterator vector_iterator; @@ -302,10 +305,8 @@ inline WeatherPrecision FGPhysicalProperties::TemperatureAt(const WeatherPrecisi return ( (it2->second - it->second)/(it2->first - it->first) ) * (a - it2->first) + it2->second; } -inline WeatherPrecision FGPhysicalProperties::AirPressureAt(const WeatherPrecision a) const -{ - return AirPressure.getValue(a); -} +//inline WeatherPrecision FGPhysicalProperties::AirPressureAt(const WeatherPrecision x) const +//moved to FGPhysicalProperties.cpp as it got too complex to inline inline WeatherPrecision FGPhysicalProperties::VaporPressureAt(const WeatherPrecision a) const { diff --git a/src/WeatherCM/FGPhysicalProperty.h b/src/WeatherCM/FGPhysicalProperty.h index ef8098cd5..fe6bca4a2 100644 --- a/src/WeatherCM/FGPhysicalProperty.h +++ b/src/WeatherCM/FGPhysicalProperty.h @@ -62,7 +62,7 @@ HISTORY #include -#include "sg.h" +#include #include "FGWeatherDefs.h" #include "FGPhysicalProperties.h" diff --git a/src/WeatherCM/FGTurbulenceItem.h b/src/WeatherCM/FGTurbulenceItem.h index 41117caaa..8818eb94c 100644 --- a/src/WeatherCM/FGTurbulenceItem.h +++ b/src/WeatherCM/FGTurbulenceItem.h @@ -58,7 +58,7 @@ HISTORY # include #endif -#include "sg.h" +#include #include "FGWeatherDefs.h" diff --git a/src/WeatherCM/FGWeatherFeature.h b/src/WeatherCM/FGWeatherFeature.h index 4f5836744..20031ae4b 100644 --- a/src/WeatherCM/FGWeatherFeature.h +++ b/src/WeatherCM/FGWeatherFeature.h @@ -51,7 +51,7 @@ HISTORY /****************************************************************************/ #include -#include "sg.h" +#include #include "FGWeatherDefs.h" @@ -86,4 +86,4 @@ public: }; /****************************************************************************/ -#endif /*FGWeatherFeature_H*/ \ No newline at end of file +#endif /*FGWeatherFeature_H*/ diff --git a/src/WeatherCM/FGWeatherParse.cpp b/src/WeatherCM/FGWeatherParse.cpp index 1f9a33f11..f83b4a1d0 100644 --- a/src/WeatherCM/FGWeatherParse.cpp +++ b/src/WeatherCM/FGWeatherParse.cpp @@ -50,12 +50,13 @@ You can also visit his homepage at http://www.wetterzentrale.de HISTORY ------------------------------------------------------------------------------ 18.10.1999 Christian Mayer Created +14.12.1999 Christian Mayer minor internal changes *****************************************************************************/ /****************************************************************************/ /* INCLUDES */ /****************************************************************************/ -#include "Include/fg_constants.h" +#include #include "FGWeatherParse.h" #include "FGWeatherUtils.h" @@ -81,48 +82,75 @@ void FGWeatherParse::input(const char *file) cerr << "Parsing \"" << file << "\" for weather datas:\n"; in.open( file ); - - while ( in ) + + if (! in.is_open() ) { - entry temp; - - in >> temp.year; - in >> temp.month; - in >> temp.day; - in >> temp.hour; - in >> temp.station_number; - in >> temp.lat; - in >> temp.lon; - in >> temp.x_wind; - in >> temp.y_wind; - in >> temp.temperature; - in >> temp.dewpoint; - in >> temp.airpressure; - in >> temp.airpressure_history[0]; - in >> temp.airpressure_history[1]; - in >> temp.airpressure_history[2]; - in >> temp.airpressure_history[3]; - - weather_station.push_back( temp ); - - // output a point to ease the waiting - if ( ((nr++)%100) == 0 ) - cerr << "."; + cerr << "Couldn't find that file!\nExiting...\n"; + exit(-1); } + else + { + bool skip = false; + while ( in ) + { + entry temp; + + in >> temp.year; + in >> temp.month; + in >> temp.day; + in >> temp.hour; + in >> temp.station_number; + in >> temp.lat; + in >> temp.lon; + in >> temp.x_wind; + in >> temp.y_wind; + in >> temp.temperature; + in >> temp.dewpoint; + in >> temp.airpressure; + in >> temp.airpressure_history[0]; + in >> temp.airpressure_history[1]; + in >> temp.airpressure_history[2]; + in >> temp.airpressure_history[3]; + + for (int i = 0; i < weather_station.size(); i++) + { + if ((weather_station[i].lat == temp.lat) && (weather_station[i].lon == temp.lon)) + { + // Two weatherstations are at the same positon + // => averageing both - cerr << "\n" << nr << " stations read\n"; + // just taking care of the stuff that metters for us + weather_station[i].x_wind += temp.x_wind; weather_station[i].x_wind *= 0.5; + weather_station[i].y_wind += temp.y_wind; weather_station[i].y_wind *= 0.5; + weather_station[i].temperature += temp.temperature; weather_station[i].temperature *= 0.5; + weather_station[i].dewpoint += temp.dewpoint; weather_station[i].dewpoint *= 0.5; + weather_station[i].airpressure += temp.airpressure; weather_station[i].airpressure *= 0.5; + + skip = true; + } + } + + if (skip == false) + weather_station.push_back( temp ); + + skip = false; + + // output a point to ease the waiting + if ( ((nr++)%100) == 0 ) + cerr << "."; + } + + cerr << "\n" << nr << " stations read\n"; + } } -FGPhysicalProperties2D FGWeatherParse::getFGPhysicalProperties2D(const unsigned int nr) const +FGPhysicalProperties FGWeatherParse::getFGPhysicalProperties(const unsigned int nr) const { - FGPhysicalProperties2D ret_val; + FGPhysicalProperties ret_val; //chache this entry entry this_entry = weather_station[nr]; - //set the position of the station - sgSetVec2( ret_val.p, this_entry.lat * DEG_TO_RAD, this_entry.lon * DEG_TO_RAD ); - ret_val.Wind[-1000.0] = FGWindItem(this_entry.x_wind, this_entry.y_wind, 0.0); ret_val.Wind[10000.0] = FGWindItem(this_entry.x_wind, this_entry.y_wind, 0.0); ret_val.Temperature[0.0] = Celsius( this_entry.temperature ); @@ -138,6 +166,9 @@ FGPhysicalProperties2D FGWeatherParse::getFGPhysicalProperties2D(const unsigned return ret_val; } - - +void FGWeatherParse::getPosition(const unsigned int nr, sgVec2 pos) const +{ + //set the position of the station + sgSetVec2( pos, weather_station[nr].lat * DEG_TO_RAD, weather_station[nr].lon * DEG_TO_RAD ); +} diff --git a/src/WeatherCM/FGWeatherParse.h b/src/WeatherCM/FGWeatherParse.h index 68077db27..ff6cd1756 100644 --- a/src/WeatherCM/FGWeatherParse.h +++ b/src/WeatherCM/FGWeatherParse.h @@ -49,6 +49,7 @@ You can also visit his homepage at http://www.wetterzentrale.de HISTORY ------------------------------------------------------------------------------ 18.10.1999 Christian Mayer Created +14.12.1999 Christian Mayer minor internal changes *****************************************************************************/ /****************************************************************************/ @@ -121,7 +122,8 @@ public: return weather_station[nr]; } - FGPhysicalProperties2D getFGPhysicalProperties2D(const unsigned int nr) const; + FGPhysicalProperties getFGPhysicalProperties(const unsigned int nr) const; + void getPosition(const unsigned int nr, sgVec2 pos) const; }; /****************************************************************************/ diff --git a/src/WeatherCM/FGWeatherVectorWrap.h b/src/WeatherCM/FGWeatherVectorWrap.h index e4226589e..76b3f95fd 100644 --- a/src/WeatherCM/FGWeatherVectorWrap.h +++ b/src/WeatherCM/FGWeatherVectorWrap.h @@ -43,7 +43,7 @@ HISTORY /****************************************************************************/ /* INCLUDES */ /****************************************************************************/ -#include "sg.h" +#include /****************************************************************************/ /* DEFINES */ diff --git a/src/WeatherCM/FGWindItem.h b/src/WeatherCM/FGWindItem.h index b416a0802..ceffd6d02 100644 --- a/src/WeatherCM/FGWindItem.h +++ b/src/WeatherCM/FGWindItem.h @@ -58,7 +58,7 @@ HISTORY # include #endif -#include "sg.h" +#include #include "FGWeatherDefs.h" diff --git a/src/WeatherCM/Makefile.am b/src/WeatherCM/Makefile.am index c3c801ac3..a2dcee186 100644 --- a/src/WeatherCM/Makefile.am +++ b/src/WeatherCM/Makefile.am @@ -3,9 +3,7 @@ noinst_LIBRARIES = libWeatherCM.a libWeatherCM_a_SOURCES = \ FGAirPressureItem.cpp FGAirPressureItem.h \ FGCloud.h FGCloudItem.cpp FGCloudItem.h \ - FGGlobalWeatherDatabase.cpp FGGlobalWeatherDatabase.h \ FGLocalWeatherDatabase.cpp FGLocalWeatherDatabase.h \ - FGMicroWeather.cpp FGMicroWeather.h \ FGPhysicalProperties.cpp FGPhysicalProperties.h \ FGPhysicalProperty.cpp FGPhysicalProperty.h \ FGSnowRain.h \ @@ -13,7 +11,6 @@ libWeatherCM_a_SOURCES = \ FGThunderstorm.cpp FGThunderstorm.h \ FGTurbulenceItem.cpp FGTurbulenceItem.h \ FGVaporPressureItem.cpp FGVaporPressureItem.h \ - FGVoronoi.cpp FGVoronoi.h \ FGWeatherDefs.h FGWeatherFeature.h FGWeatherUtils.h \ FGWeatherParse.cpp FGWeatherParse.h \ FGWeatherVectorWrap.h \ diff --git a/src/WeatherCM/air-pressure-explanation.txt b/src/WeatherCM/air-pressure-explanation.txt new file mode 100644 index 000000000..b26f118cc --- /dev/null +++ b/src/WeatherCM/air-pressure-explanation.txt @@ -0,0 +1,212 @@ +The formula p(x) for calculating the air pressure at a given altitude +--------------------------------------------------------------------- + +Well known is the baromertic(?) formula + + rho0 + ------ * g * x + p0 + p(x) = p0 * e + +with p0 being the airpressure and rho0 being the air density at an altitude of +0 metres above sea level and g being the gravity constant of 9.81 m/sq. sec + +This formula can easily be derivated when you know, that: + + * the pressure difference is + + dp = - rho * g * dx + + * Boyle-Mariotte says: + + p0 : p = rho0 : rho + +Combinig the terms and changing them around I get: + + dp [ rho0 ] + -- = - rho * g = - [ ------ * p(x) ] * g + dx [ p0 ] + + rho0 + p'(x) = - ------ * p(x) * g + p0 + +Solving that differential equation and knowing that p(0) = p0 I get: + + rho0 + - ------ * g * x + p0 + p(x) = p0 * e + +q.e.d. + +------------------------------------------------------------------------------- + +The problem with that equation is that it doesn't take different temperatures +at different altitudes into account. And the inaccuracies due to it are huge. +That's why this formula is only used in very low altitudes. + +So to get a usefull formula for FlightGear I need to extend it. And as I'm +already 'recreating' that formula I'm taking the change in g into account, too. +This doesn't make such a dramatic difference to the result as the inclusion of +temperature change does, but it doesn't complicate the final formula too much. + +So I get three formulas that I'm combining in the end: + + * the change of g with the altitude: + + G * m + g(x) = ----------- + (x + r)^2 + G is the universal gravity constant(?) and is 6.673e-11 m^3 kg^-1 s^-2 + m is the mass of the earth and is 5.977e24 kg + r is the radius of the earth and is 6368 km + + * The pressure difference stays the same: + + dp = - rho * g(x) * dx + + * If I combine Boyle-Mariotte with Gay-Lussac I get: + + rho0 * T0 p + rho = ----------- * --- + p0 T + +Combining the terms again I get this time: + + dp [ rho0 * T0 p(x) ] + -- = - rho * g(x) = - [ ----------- * ------ ] * g(x) + dx [ p0 T(x) ] + + rho0 * T0 p(x) * g(x) + p'(x) = - ----------- * ------------- + p0 T(x) + +This DE isn't that easy to solve as the one above, it by looking into the right +books you'll see the general solution for: + + y' + f(x)*y = 0 + +is + x + /\ + - | f(x) dx + \/ + n + y = m * e + +and P(m,n) will be a point on the graph. + +For q = n = 0 metres altitude we get y = m. As y is p(x) we know that m has to +be p0. + +So our final formuala is + + ho0 * T0 g(x) + f1(x) = ----------- * ------ + p0 T(x) + + + x x + /\ /\ + - | f1(x) dx | f(x) dx + \/ \/ + 0 0 F(x) - F(0) + p(x) = p0 * e = p0 * e = p0 * e + +The only disturbing thing we've got left is the integral. Luckily there is a +great service at http://integrals.wolfram.com/ that helps me doing it :-) + +But the f(x) is still too general so I'm substituting: + + rho0 * T0 * G * m + f(x) = - ----------------------- + p0 * (x + r)^2 * T(x) + +but even that isn't good enough. But as I'm linearily interpolating between +two different temperatures I can say that T(x) = a*x + b for the x inbetween +two different stored temperatures. So I just need to integrate every pice +independandly. But anyway, I get: + + rho0 * T0 * G * m + f(x) = - ------------------------------ + p0 * (x + r)^2 * (a * x + b) + +Integrating that I get: + + rho0 * T0 * G * m [ 1 + F(x) = - ------------------- * [ ------------------------ - + p0 [ (-b + a * r) * (r + x) + + + a * log|r + x| a * log|b + a * x| ] + ---------------- + -------------------- ] + (b - a * r)^2 (b - a * r)^2 ] + +To lower the computional cost I can transfere the equation. + + * I'm defining + + rho0 * T0 * G * m + factor = - ------------------- + p0 + 1 + c = -------------- + (-b + a * r) + + * now I can write + + [ c ] + F(x) = factor * [ --------- - a * c * c * [log|r + x| + log|b + a * x|] ] + [ (r + x) ] + + * and simplyfy it to + + [ 1 ] + F(x) = factor * c * [ --------- - a * c * log|(r + x) * (b + a * x)| ] + [ (r + x) ] + +------------------------------------------------------------------------------- +The following table shows quite nicely how accurate my formula is: + + Altitude[m] | Airpressure [hPa] | Error [%] + | Official | My formula | + ------------+---------------+---------------+--------------- + -200 | 1037.51 | 1037.24 | 0.0260 + -100 | 1025.32 | 1025.19 | 0.0127 + 0 | 1013.25 | 1013.25 | 0.0 + 500 | 954.59 | 955.224 | 0.0664 + 1000 | 898.70 | 899.912 | 0.1349 + 2000 | 794.88 | 797.042 | 0.2720 + 3000 | 700.99 | 703.885 | 0.4130 + 4000 | 616.28 | 619.727 | 0.5593 + 5000 | 540.07 | 543.89 | 0.7073 + 6000 | 471.67 | 475.731 | 0.8610 + 7000 | 410.46 | 414.643 | 1.0191 + 8000 | 355.84 | 360.054 | 1.1842 + 9000 | 307.27 | 311.422 | 1.3513 + 10000 | 264.21 | 268.238 | 1.5245 + 20000 | 54.670/55.3 | 55.7971 | 2.0616/0.8989 + 30000 | 11.8 | 11.3149 | 1.5441 + 40000 | 3.0 | 2.74665 | 18.9703 + 50000 | 0.88 | 0.753043 | 41.9183 + 60000 | 0.257 | 0.221907 | 57.9802 + 70000 | 0.0602 | 0.0530785 | 61.9153 + 80000 | 0.0101 | 0.00905461 | 51.5725 + 100000 | 2.14e-4 | 2.03984e-4 | 5.5131 + +The official values are from the CINA atmosphere which assumes a air pressure +of 1013.25 hPa and a temperature of 15 degC at sea level and a temperature +gradient of -6.5 deg/km. The CINA atmosphere gives only values for altiudes +up to 20 km. The values for 20 km and above are from the 1959 ARDC atmosphere. +That's why I've got two values at 20000 metres. +The temperature changes dramtically in the altitudes over 20 km which I didn't +take care of in my calculations which explains the huge errors at that altitude +range. But you can see nicely that the values are at least quite close to the +official values. +Using a better temperature model for the altitudes above 20 km should +dramatically increase the accuracy there. + + + +