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
*****************************************************************************/
/****************************************************************************/
/* 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);
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.
*****************************************************************************/
/****************************************************************************/
#include <Aircraft/aircraft.hxx>
#include "FGLocalWeatherDatabase.h"
-#include "FGVoronoi.h"
#include "FGWeatherParse.h"
/********************************** 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;
if (theFGLocalWeatherDatabase)
{
- //FG_LOG( FG_GENERAL, FG_ALERT, "Error: only one local weather allowed" );
cerr << "Error: only one local weather allowed";
exit(-1);
}
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
{
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" );
+ unsigned int n = parsed_data->stored_stations();
- parsed_data.input( "weather/current.gz" );
- global = new FGGlobalWeatherDatabase; //initialize GlobalDatabase
- global->setDatabaseStatus(FGGlobalWeatherDatabase_working);
+ 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<FGPhysicalProperties>(n, p, f);
+
+ //and free my allocations:
+ delete[] p;
+ delete[] f;
cerr << "Finished weather init.\n";
}
case manual:
case default_mode:
{
-
- vector<sgVec2Wrap> 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<FGPhysicalProperties>(2,x,y,z,f);
}
break;
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;
}
/****************************************************************************/
/****************************************************************************/
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";
}
/****************************************************************************/
/****************************************************************************/
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);
}
/****************************************************************************/
/****************************************************************************/
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(),
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.
*****************************************************************************/
/****************************************************************************/
/****************************************************************************/
/* 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 <windows.h>
-#endif
-#include <GL/glut.h>
-#include <XGL/xgl.h>
-
#include <vector>
-#include "sg.h"
+#include <sg.h>
+
+#include <Math/sphrintp.h>
#include "FGPhysicalProperties.h"
-#include "FGGlobalWeatherDatabase.h"
-#include "FGMicroWeather.h"
+#include "FGPhysicalProperty.h"
+
+
#include "FGWeatherFeature.h"
#include "FGWeatherDefs.h"
#include "FGThunderstorm.h"
{
private:
protected:
- FGGlobalWeatherDatabase *global; //point to the global database
-
- typedef vector<FGMicroWeather> FGMicroWeatherList;
- typedef FGMicroWeatherList::iterator FGMicroWeatherListIt;
+ SphereInterpolate<FGPhysicalProperties> *database;
typedef vector<sgVec2> pointVector;
typedef vector<pointVector> tileVector;
/************************************************************************/
/* 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;
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
/************************************************************************/
/* 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
/************************************************************************/
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;
}
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
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: ";
}
+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<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);
+}
+
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
*****************************************************************************/
/****************************************************************************/
#include <vector>
#include <map>
-#include "sg.h"
+#include <sg.h>
#include "FGWeatherDefs.h"
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:
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();
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();
return *this;
}
-
inline void FGPhysicalProperties::WindAt(sgVec3 ret, const WeatherPrecision a) const
{
typedef map<FGPhysicalProperties::Altitude, FGWindItem>::const_iterator vector_iterator;
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
{
#include <vector>
-#include "sg.h"
+#include <sg.h>
#include "FGWeatherDefs.h"
#include "FGPhysicalProperties.h"
# include <windows.h>
#endif
-#include "sg.h"
+#include <sg.h>
#include "FGWeatherDefs.h"
/****************************************************************************/
#include <Include/compiler.h>
-#include "sg.h"
+#include <sg.h>
#include "FGWeatherDefs.h"
};
/****************************************************************************/
-#endif /*FGWeatherFeature_H*/
\ No newline at end of file
+#endif /*FGWeatherFeature_H*/
HISTORY
------------------------------------------------------------------------------
18.10.1999 Christian Mayer Created
+14.12.1999 Christian Mayer minor internal changes
*****************************************************************************/
/****************************************************************************/
/* INCLUDES */
/****************************************************************************/
-#include "Include/fg_constants.h"
+#include <Include/fg_constants.h>
#include "FGWeatherParse.h"
#include "FGWeatherUtils.h"
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
+
+ // 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";
}
-
- 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 );
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 );
+}
HISTORY
------------------------------------------------------------------------------
18.10.1999 Christian Mayer Created
+14.12.1999 Christian Mayer minor internal changes
*****************************************************************************/
/****************************************************************************/
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;
};
/****************************************************************************/
/****************************************************************************/
/* INCLUDES */
/****************************************************************************/
-#include "sg.h"
+#include <sg.h>
/****************************************************************************/
/* DEFINES */
# include <windows.h>
#endif
-#include "sg.h"
+#include <sg.h>
#include "FGWeatherDefs.h"
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 \
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 \
--- /dev/null
+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.
+
+
+
+
projects / flightgear.git / commitdiff