fly/flightgear
1
0
Fork 0
Code Activity

Major weather update from Christian Mayer, tying the weather code into

Browse source 
the property system, among other things.  A separate integration into
the FDMs will follow shortly.

This code will be used only if the --with-new-environment option is
*not* passed to configure.
This commit is contained in:
david 2002-03-16 20:31:27 +00:00
parent ac99a82b98
commit d4c49d65ac
15 changed files with 1475 additions and 49 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

9
src/WeatherCM/FGAirPressureItem.h
View file

@ -80,11 +80,16 @@ public:
FGAirPressureItem(const WeatherPrecision v) {value = v; }
FGAirPressureItem() {value = FG_WEATHER_DEFAULT_AIRPRESSURE;}
WeatherPrecision getValue(void) const
WeatherPrecision getValue() const
{
return value;
};
void setValue(WeatherPrecision p)
{
value = p;
};
FGAirPressureItem& operator*=(const WeatherPrecision arg);
FGAirPressureItem& operator+=(const FGAirPressureItem& arg);
FGAirPressureItem& operator-=(const FGAirPressureItem& arg);

139
src/WeatherCM/FGLocalWeatherDatabase.cpp
View file

@ -88,7 +88,7 @@ void FGLocalWeatherDatabase::init( const WeatherPrecision visibility,
setWeatherVisibility(visibility);
DatabaseStatus = type;
database = 0; //just get sure...
database_logic = 0; //just get sure...
Thunderstorm = false;
//I don't need to set theThunderstorm as Thunderstorm == false
@ -105,8 +105,8 @@ void FGLocalWeatherDatabase::init( const WeatherPrecision visibility,
case use_internet:
{
FGWeatherParse *parsed_data = new FGWeatherParse();
sgVec2 *p;
FGPhysicalProperties *f;
sgVec2 *p;
unsigned int *f;
string path_to_weather = root + "/weather/current.txt.gz";
parsed_data->input( path_to_weather.c_str() );
unsigned int n = parsed_data->stored_stations();
@ -116,13 +116,14 @@ void FGLocalWeatherDatabase::init( const WeatherPrecision visibility,
{
m = n;
p = new sgVec2 [n];
f = new FGPhysicalProperties[n];
p = new sgVec2[n];
f = new unsigned int[n];
// fill the database
for (unsigned int i = 0; i < n; i++)
{
f[i] = parsed_data->getFGPhysicalProperties(i);
f[i] = i;
database_data[i] = parsed_data->getFGPhysicalProperties(i);
parsed_data->getPosition(i, p[i]);
if ( (i%100) == 0)
@ -150,8 +151,8 @@ void FGLocalWeatherDatabase::init( const WeatherPrecision visibility,
squared_distance[sgDistanceSquaredVec2(cur_pos, pos)] = i;
}
p = new sgVec2 [m];
f = new FGPhysicalProperties[m];
p = new sgVec2 [m];
f = new unsigned int[m];
map<float, unsigned int>::const_iterator ci;
ci = squared_distance.begin();
@ -159,7 +160,8 @@ void FGLocalWeatherDatabase::init( const WeatherPrecision visibility,
// fill the database
for ( i = 0; i < m; i++ )
{
f[i] = parsed_data->getFGPhysicalProperties(ci->second);
f[i] = i;
database_data.push_back( parsed_data->getFGPhysicalProperties(ci->second) );
parsed_data->getPosition(ci->second, p[i]);
if ( (i%100) == 0)
@ -175,7 +177,7 @@ void FGLocalWeatherDatabase::init( const WeatherPrecision visibility,
//and finally init the interpolation
cerr << "\nInitialiating Interpolation. (2-3 minutes on a PII-350 for ca. 3500 stations)\n";
database = new SphereInterpolate<FGPhysicalProperties>(m, p, f);
database_logic = new SphereInterpolate(m, p, f);
//and free my allocations:
delete[] p;
@ -195,9 +197,9 @@ void FGLocalWeatherDatabase::init( const WeatherPrecision visibility,
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 = new FGPhysicalProperties[2]; //make an standard weather that's the same at the whole world
database = new SphereInterpolate<FGPhysicalProperties>(2,x,y,z,f);
delete[] f;
unsigned int f[2] = {0, 0};
database_data.push_back( FGPhysicalProperties() ); // == database_date[0]
database_logic = new SphereInterpolate(2,x,y,z,f);
}
break;
@ -209,6 +211,10 @@ void FGLocalWeatherDatabase::init( const WeatherPrecision visibility,
cache->longitude_deg = fgGetNode("/position/longitude-deg");
cache->altitude_ft = fgGetNode("/position/altitude-ft" );
}
void FGLocalWeatherDatabase::bind()
{
fgTie("/environment/weather/wind-north-mps", this, &FGLocalWeatherDatabase::get_wind_north);
fgTie("/environment/weather/wind-east-mps", this, &FGLocalWeatherDatabase::get_wind_east);
fgTie("/environment/weather/wind-up-mps", this, &FGLocalWeatherDatabase::get_wind_up);
@ -216,12 +222,109 @@ void FGLocalWeatherDatabase::init( const WeatherPrecision visibility,
fgTie("/environment/weather/air-pressure-Pa", this, &FGLocalWeatherDatabase::get_air_pressure);
fgTie("/environment/weather/vapor-pressure-Pa", this, &FGLocalWeatherDatabase::get_vapor_pressure);
fgTie("/environment/weather/air-density", this, &FGLocalWeatherDatabase::get_air_density);
SGPropertyNode * station_nodes = fgGetNode("/environment/weather");
if (station_nodes == 0) {
cerr << "No weatherstations (/environment/weather)!!";
return;
}
int index = 0;
for(vector<FGPhysicalProperties>::iterator it = database_data.begin(); it != database_data.end(); it++)
{
SGPropertyNode * station = station_nodes->getNode("station", index, true);
station -> tie("air-pressure-Pa",
SGRawValueMethods<FGAirPressureItem,WeatherPrecision>(
database_data[0].AirPressure,
&FGAirPressureItem::getValue,
&FGAirPressureItem::setValue)
,false);
int i;
for( i = 0; i < database_data[index].Wind.size(); i++)
{
SGPropertyNode * wind = station->getNode("wind", i, true);
wind -> tie("north-mps",
SGRawValueMethodsIndexed<FGPhysicalProperties,WeatherPrecision>(
database_data[index], i,
&FGPhysicalProperties::getWind_x,
&FGPhysicalProperties::setWind_x)
,false);
wind -> tie("east-mps",
SGRawValueMethodsIndexed<FGPhysicalProperties,WeatherPrecision>(
database_data[index], i,
&FGPhysicalProperties::getWind_y,
&FGPhysicalProperties::setWind_y)
,false);
wind -> tie("up-mps",
SGRawValueMethodsIndexed<FGPhysicalProperties,WeatherPrecision>(
database_data[index], i,
&FGPhysicalProperties::getWind_z,
&FGPhysicalProperties::setWind_z)
,false);
wind -> tie("altitude-m",
SGRawValueMethodsIndexed<FGPhysicalProperties,WeatherPrecision>(
database_data[index], i,
&FGPhysicalProperties::getWind_a,
&FGPhysicalProperties::setWind_a)
,false);
}
for( i = 0; i < database_data[index].Temperature.size(); i++)
{
SGPropertyNode * temperature = station->getNode("temperature", i, true);
temperature -> tie("value-K",
SGRawValueMethodsIndexed<FGPhysicalProperties,WeatherPrecision>(
database_data[index], i,
&FGPhysicalProperties::getTemperature_x,
&FGPhysicalProperties::setTemperature_x)
,false);
temperature -> tie("altitude-m",
SGRawValueMethodsIndexed<FGPhysicalProperties,WeatherPrecision>(
database_data[index], i,
&FGPhysicalProperties::getTemperature_a,
&FGPhysicalProperties::setTemperature_a)
,false);
}
for( i = 0; i < database_data[index].VaporPressure.size(); i++)
{
SGPropertyNode * vaporpressure = station->getNode("vapor-pressure", i, true);
vaporpressure -> tie("value-Pa",
SGRawValueMethodsIndexed<FGPhysicalProperties,WeatherPrecision>(
database_data[index], i,
&FGPhysicalProperties::getVaporPressure_x,
&FGPhysicalProperties::setVaporPressure_x)
,false);
vaporpressure -> tie("altitude-m",
SGRawValueMethodsIndexed<FGPhysicalProperties,WeatherPrecision>(
database_data[index], i,
&FGPhysicalProperties::getVaporPressure_a,
&FGPhysicalProperties::setVaporPressure_a)
,false);
}
index++;
}
}
void FGLocalWeatherDatabase::unbind()
{
fgUntie("/environment/weather/wind-north-mps");
fgUntie("/environment/weather/wind-east-mps");
fgUntie("/environment/weather/wind-up-mps");
fgUntie("/environment/weather/temperature-K");
fgUntie("/environment/weather/air-pressure-Pa");
fgUntie("/environment/weather/vapor-pressure-Pa");
fgUntie("/environment/weather/air-density");
}
FGLocalWeatherDatabase::~FGLocalWeatherDatabase()
{
//Tidying up:
delete database;
delete database_logic;
}
/****************************************************************************/
@ -265,9 +368,9 @@ FGPhysicalProperty FGLocalWeatherDatabase::get(const sgVec3& p) const
// check for bogous altitudes. Dunno why, but FGFS want's to know the
// weather at an altitude of roughly -3000 meters...
if (p[2] < -500.0f)
return FGPhysicalProperty(database->Evaluate(p), -500.0f);
return FGPhysicalProperty(DatabaseEvaluate(p), -500.0f);
return FGPhysicalProperty(database->Evaluate(p), p[2]);
return FGPhysicalProperty(DatabaseEvaluate(p), p[2]);
}
#ifdef macintosh
@ -279,13 +382,13 @@ FGPhysicalProperty FGLocalWeatherDatabase::get(const sgVec3& p) const
#else
FGPhysicalProperties FGLocalWeatherDatabase::get(const sgVec2& p) const
{
return database->Evaluate(p);
return DatabaseEvaluate(p);
}
#endif
WeatherPrecision FGLocalWeatherDatabase::getAirDensity(const sgVec3& p) const
{
FGPhysicalProperty dummy(database->Evaluate(p), p[2]);
FGPhysicalProperty dummy(DatabaseEvaluate(p), p[2]);
return
(dummy.AirPressure*FG_WEATHER_DEFAULT_AIRDENSITY*FG_WEATHER_DEFAULT_TEMPERATURE) /

38
src/WeatherCM/FGLocalWeatherDatabase.h
View file

@ -60,8 +60,10 @@ HISTORY
#include <plib/sg.h>
#include <simgear/misc/props.hxx>
#include <Main/fgfs.hxx>
#include <simgear/constants.h>
#include <simgear/misc/props.hxx>
#include "sphrintp.h"
@ -92,11 +94,33 @@ struct _FGLocalWeatherDatabaseCache
FGPhysicalProperty last_known_property;
};
class FGLocalWeatherDatabase
class FGLocalWeatherDatabase : public FGSubsystem
{
private:
protected:
SphereInterpolate<FGPhysicalProperties> *database;
SphereInterpolate *database_logic;
vector<FGPhysicalProperties> database_data;
#ifndef macintosh
FGPhysicalProperties DatabaseEvaluate(const sgVec2& p) const
{
sgVec2 p_converted = {p[0]*(SGD_2PI/360.0),
p[1]*(SGD_2PI/360.0)};
EvaluateData d = database_logic->Evaluate(p_converted);
return database_data[d.index[0]]*d.percentage[0] +
database_data[d.index[1]]*d.percentage[1] +
database_data[d.index[2]]*d.percentage[2] ;
}
#endif
FGPhysicalProperties DatabaseEvaluate(const sgVec3& p) const
{
sgVec3 p_converted = {p[0]*(SGD_2PI/360.0),
p[1]*(SGD_2PI/360.0),
p[2] };
EvaluateData d = database_logic->Evaluate(p_converted);
return database_data[d.index[0]]*d.percentage[0] +
database_data[d.index[1]]*d.percentage[1] +
database_data[d.index[2]]*d.percentage[2] ;
}
typedef vector<sgVec2> pointVector;
typedef vector<pointVector> tileVector;
@ -190,6 +214,14 @@ public:
void update(const sgVec3& p); //position has changed
void update(const sgVec3& p, const WeatherPrecision dt); //time and/or position has changed
/************************************************************************/
/* define methods requited for FGSubsystem */
/************************************************************************/
virtual void init () { /* do nothing; that's done in the constructor */ };
virtual void bind ();
virtual void unbind ();
virtual void update (int dt) { update((float) dt); };
/************************************************************************/
/* Get the physical properties on the specified point p */
/************************************************************************/

33
src/WeatherCM/FGPhysicalProperties.cpp
View file

@ -80,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();
@ -289,5 +320,3 @@ WeatherPrecision FGPhysicalProperties::AirPressureAt(const WeatherPrecision x) c
return AirPressure.getValue() * exp(FF);
}

26
src/WeatherCM/FGPhysicalProperties.h
View file

@ -122,6 +122,32 @@ public:
FGPhysicalProperties& operator *= ( const WeatherPrecision d );
FGPhysicalProperties& operator += ( const FGPhysicalProperties& p );
FGPhysicalProperties& operator -= ( const FGPhysicalProperties& p );
//for easy binding to the property system
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);
};
class FGPhysicalProperties2D : public FGPhysicalProperties

480
src/WeatherCM/FGPhysicalProperties_bind.cpp Normal file
View file

@ -0,0 +1,480 @@
/*****************************************************************************
Module: FGPhysicalProperties_bind.cpp
Author: Christian Mayer
Date started: 15.03.02
Called by: main program
-------- Copyright (C) 2002 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
Foundation; either version 2 of the License, or (at your option) any later
version.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
details.
You should have received a copy of the GNU General Public License along with
this program; if not, write to the Free Software Foundation, Inc., 59 Temple
Place - Suite 330, Boston, MA 02111-1307, USA.
Further information about the GNU General Public License can also be found on
the world wide web at http://www.gnu.org.
FUNCTIONAL DESCRIPTION
------------------------------------------------------------------------------
Initialice the FGPhysicalProperties struct to something sensible(?)
HISTORY
------------------------------------------------------------------------------
15.03.2002 Christian Mayer Created
*****************************************************************************/
#include "FGPhysicalProperties.h"
WeatherPrecision FGPhysicalProperties::getWind_x( int number ) const
{
if (number >= Wind.size())
{
cerr << "Error: There are " << Wind.size() << " entries for wind, but number " << number << " got requested!\n";
return 0.0;
}
map<Altitude,FGWindItem>::const_iterator it = Wind.begin();
while( number > 0)
{
it++;
number--;
}
return it->second.x();
}
void FGPhysicalProperties::setWind_x( int number, WeatherPrecision x)
{
if (number >= Wind.size())
{
cerr << "Error: There are " << Wind.size() << " entries for wind, but number " << number << " got requested!\n";
}
map<Altitude,FGWindItem>::iterator it = Wind.begin();
while( number > 0)
{
it++;
number--;
}
it->second.x( x );
}
WeatherPrecision FGPhysicalProperties::getWind_y( int number ) const
{
if (number >= Wind.size())
{
cerr << "Error: There are " << Wind.size() << " entries for wind, but number " << number << " got requested!\n";
return 0.0;
}
map<Altitude,FGWindItem>::const_iterator it = Wind.begin();
while( number > 0)
{
it++;
number--;
}
return it->second.y();
}
void FGPhysicalProperties::setWind_y( int number, WeatherPrecision y)
{
if (number >= Wind.size())
{
cerr << "Error: There are " << Wind.size() << " entries for wind, but number " << number << " got requested!\n";
}
map<Altitude,FGWindItem>::iterator it = Wind.begin();
while( number > 0)
{
it++;
number--;
}
it->second.y( y );
}
WeatherPrecision FGPhysicalProperties::getWind_z( int number ) const
{
if (number >= Wind.size())
{
cerr << "Error: There are " << Wind.size() << " entries for wind, but number " << number << " got requested!\n";
return 0.0;
}
map<Altitude,FGWindItem>::const_iterator it = Wind.begin();
while( number > 0)
{
it++;
number--;
}
return it->second.z();
}
void FGPhysicalProperties::setWind_z( int number, WeatherPrecision z)
{
if (number >= Wind.size())
{
cerr << "Error: There are " << Wind.size() << " entries for wind, but number " << number << " got requested!\n";
}
map<Altitude,FGWindItem>::iterator it = Wind.begin();
while( number > 0)
{
it++;
number--;
}
it->second.z( z );
}
WeatherPrecision FGPhysicalProperties::getWind_a( int number ) const
{
if (number >= Wind.size())
{
cerr << "Error: There are " << Wind.size() << " entries for wind, but number " << number << " got requested!\n";
return 0.0;
}
map<Altitude,FGWindItem>::const_iterator it = Wind.begin();
while( number > 0)
{
it++;
number--;
}
return it->first;
}
void FGPhysicalProperties::setWind_a( int number, WeatherPrecision a)
{
if (number >= Wind.size())
{
cerr << "Error: There are " << Wind.size() << " entries for wind, but number " << number << " got requested!\n";
}
map<Altitude,FGWindItem>::iterator it = Wind.begin();
while( number > 0)
{
it++;
number--;
}
it->first;
}
WeatherPrecision FGPhysicalProperties::getTurbulence_x( int number ) const
{
if (number >= Turbulence.size())
{
cerr << "Error: There are " << Turbulence.size() << " entries for Turbulence, but number " << number << " got requested!\n";
return 0.0;
}
map<Altitude,FGTurbulenceItem>::const_iterator it = Turbulence.begin();
while( number > 0)
{
it++;
number--;
}
return it->second.x();
}
void FGPhysicalProperties::setTurbulence_x( int number, WeatherPrecision x)
{
if (number >= Turbulence.size())
{
cerr << "Error: There are " << Turbulence.size() << " entries for Turbulence, but number " << number << " got requested!\n";
}
map<Altitude,FGTurbulenceItem>::iterator it = Turbulence.begin();
while( number > 0)
{
it++;
number--;
}
it->second.x( x );
}
WeatherPrecision FGPhysicalProperties::getTurbulence_y( int number ) const
{
if (number >= Turbulence.size())
{
cerr << "Error: There are " << Turbulence.size() << " entries for Turbulence, but number " << number << " got requested!\n";
return 0.0;
}
map<Altitude,FGTurbulenceItem>::const_iterator it = Turbulence.begin();
while( number > 0)
{
it++;
number--;
}
return it->second.y();
}
void FGPhysicalProperties::setTurbulence_y( int number, WeatherPrecision y)
{
if (number >= Turbulence.size())
{
cerr << "Error: There are " << Turbulence.size() << " entries for Turbulence, but number " << number << " got requested!\n";
}
map<Altitude,FGTurbulenceItem>::iterator it = Turbulence.begin();
while( number > 0)
{
it++;
number--;
}
it->second.y( y );
}
WeatherPrecision FGPhysicalProperties::getTurbulence_z( int number ) const
{
if (number >= Turbulence.size())
{
cerr << "Error: There are " << Turbulence.size() << " entries for Turbulence, but number " << number << " got requested!\n";
return 0.0;
}
map<Altitude,FGTurbulenceItem>::const_iterator it = Turbulence.begin();
while( number > 0)
{
it++;
number--;
}
return it->second.z();
}
void FGPhysicalProperties::setTurbulence_z( int number, WeatherPrecision z)
{
if (number >= Turbulence.size())
{
cerr << "Error: There are " << Turbulence.size() << " entries for Turbulence, but number " << number << " got requested!\n";
}
map<Altitude,FGTurbulenceItem>::iterator it = Turbulence.begin();
while( number > 0)
{
it++;
number--;
}
it->second.z( z );
}
WeatherPrecision FGPhysicalProperties::getTurbulence_a( int number ) const
{
if (number >= Turbulence.size())
{
cerr << "Error: There are " << Turbulence.size() << " entries for Turbulence, but number " << number << " got requested!\n";
return 0.0;
}
map<Altitude,FGTurbulenceItem>::const_iterator it = Turbulence.begin();
while( number > 0)
{
it++;
number--;
}
return it->first;
}
void FGPhysicalProperties::setTurbulence_a( int number, WeatherPrecision a)
{
if (number >= Turbulence.size())
{
cerr << "Error: There are " << Turbulence.size() << " entries for Turbulence, but number " << number << " got requested!\n";
}
map<Altitude,FGTurbulenceItem>::iterator it = Turbulence.begin();
while( number > 0)
{
it++;
number--;
}
it->first;
}
WeatherPrecision FGPhysicalProperties::getTemperature_x( int number ) const
{
if (number >= Temperature.size())
{
cerr << "Error: There are " << Temperature.size() << " entries for Temperature, but number " << number << " got requested!\n";
return 0.0;
}
map<Altitude,WeatherPrecision>::const_iterator it = Temperature.begin();
while( number > 0)
{
it++;
number--;
}
return it->second;
}
void FGPhysicalProperties::setTemperature_x( int number, WeatherPrecision x)
{
if (number >= Temperature.size())
{
cerr << "Error: There are " << Temperature.size() << " entries for Temperature, but number " << number << " got requested!\n";
}
map<Altitude,WeatherPrecision>::iterator it = Temperature.begin();
while( number > 0)
{
it++;
number--;
}
it->second = x;
}
WeatherPrecision FGPhysicalProperties::getTemperature_a( int number ) const
{
if (number >= Temperature.size())
{
cerr << "Error: There are " << Temperature.size() << " entries for Temperature, but number " << number << " got requested!\n";
return 0.0;
}
map<Altitude,WeatherPrecision>::const_iterator it = Temperature.begin();
while( number > 0)
{
it++;
number--;
}
return it->first;
}
void FGPhysicalProperties::setTemperature_a( int number, WeatherPrecision a)
{
if (number >= Temperature.size())
{
cerr << "Error: There are " << Temperature.size() << " entries for Temperature, but number " << number << " got requested!\n";
}
map<Altitude,WeatherPrecision>::iterator it = Temperature.begin();
while( number > 0)
{
it++;
number--;
}
it->first;
}
WeatherPrecision FGPhysicalProperties::getVaporPressure_x( int number ) const
{
if (number >= VaporPressure.size())
{
cerr << "Error: There are " << VaporPressure.size() << " entries for VaporPressure, but number " << number << " got requested!\n";
return 0.0;
}
map<Altitude,WeatherPrecision>::const_iterator it = VaporPressure.begin();
while( number > 0)
{
it++;
number--;
}
return it->second;
}
void FGPhysicalProperties::setVaporPressure_x( int number, WeatherPrecision x)
{
if (number >= VaporPressure.size())
{
cerr << "Error: There are " << Temperature.size() << " entries for VaporPressure, but number " << number << " got requested!\n";
}
map<Altitude,WeatherPrecision>::iterator it = VaporPressure.begin();
while( number > 0)
{
it++;
number--;
}
it->second = x;
}
WeatherPrecision FGPhysicalProperties::getVaporPressure_a( int number ) const
{
if (number >= VaporPressure.size())
{
cerr << "Error: There are " << VaporPressure.size() << " entries for VaporPressure, but number " << number << " got requested!\n";
return 0.0;
}
map<Altitude,WeatherPrecision>::const_iterator it = VaporPressure.begin();
while( number > 0)
{
it++;
number--;
}
return it->first;
}
void FGPhysicalProperties::setVaporPressure_a( int number, WeatherPrecision a)
{
if (number >= VaporPressure.size())
{
cerr << "Error: There are " << VaporPressure.size() << " entries for VaporPressure, but number " << number << " got requested!\n";
}
map<Altitude,WeatherPrecision>::iterator it = VaporPressure.begin();
while( number > 0)
{
it++;
number--;
}
it->first;
}

3
src/WeatherCM/FGTemperatureItem.h
View file

@ -76,6 +76,9 @@ public:
WeatherPrecision getValue() const { return value; };
WeatherPrecision getAlt() const { return alt; };
void setAlt (WeatherPrecision x) { alt = x; }
void setValue(WeatherPrecision x) { value = x; }
FGTemperatureItem& operator*= (const WeatherPrecision& arg);
FGTemperatureItem& operator+= (const FGTemperatureItem& arg);
FGTemperatureItem& operator-= (const FGTemperatureItem& arg);

4
src/WeatherCM/FGTurbulenceItem.h
View file

@ -96,6 +96,10 @@ public:
WeatherPrecision y(void) const { return value[1]; };
WeatherPrecision z(void) const { return value[2]; };
void x(const WeatherPrecision x) { value[0] = x; };
void y(const WeatherPrecision y) { value[1] = y; };
void z(const WeatherPrecision z) { value[2] = z; };
FGTurbulenceItem& operator*= (const WeatherPrecision arg);
FGTurbulenceItem& operator+= (const FGTurbulenceItem& arg);
FGTurbulenceItem& operator-= (const FGTurbulenceItem& arg);

4
src/WeatherCM/FGWeatherUtils.h
View file

@ -198,7 +198,7 @@ inline WeatherPrecision JSBsim2SIdensity (const WeatherPrecision JSBsim)
inline WeatherPrecision psf2Pascal (const WeatherPrecision psf)
{
return psf / 0.0020885434; //lbs / square foot (used in JSBsim)
return psf / 0.020885434; //lbs / square foot (used in JSBsim)
}
inline WeatherPrecision Kelvin2Rankine (const WeatherPrecision kelvin)
@ -213,7 +213,7 @@ inline WeatherPrecision SIdensity2JSBsim (const WeatherPrecision SIdensity)
inline WeatherPrecision Pascal2psf (const WeatherPrecision Pascal)
{
return 0.0020885434 * Pascal; //lbs / square feet (used in JSBsim)
return 0.020885434 * Pascal; //lbs / square feet (used in JSBsim)
}
/****************************************************************************/

4
src/WeatherCM/FGWindItem.h
View file

@ -99,6 +99,10 @@ public:
WeatherPrecision y(void) const { return value[1]; };
WeatherPrecision z(void) const { return value[2]; };
void x(const WeatherPrecision x) { value[0] = x; };
void y(const WeatherPrecision y) { value[1] = y; };
void z(const WeatherPrecision z) { value[2] = z; };
FGWindItem& operator*= (const WeatherPrecision arg);
FGWindItem& operator+= (const FGWindItem& arg);
FGWindItem& operator-= (const FGWindItem& arg);

7
src/WeatherCM/Makefile.am
View file

@ -1,12 +1,13 @@
noinst_LIBRARIES = libWeatherCM.a
EXTRA_DIST = linintp2.h linintp2.inl sphrintp.h sphrintp.inl
EXTRA_DIST = linintp2.inl sphrintp.inl
libWeatherCM_a_SOURCES = \
FGAirPressureItem.cpp FGAirPressureItem.h \
FGCloud.h FGCloudItem.cpp FGCloudItem.h \
FGLocalWeatherDatabase.cpp FGLocalWeatherDatabase.h \
FGPhysicalProperties.cpp FGPhysicalProperties.h \
FGPhysicalProperties_bind.cpp \
FGPhysicalProperty.cpp FGPhysicalProperty.h \
FGSnowRain.h \
FGTemperatureItem.cpp FGTemperatureItem.h \
@ -16,7 +17,9 @@ libWeatherCM_a_SOURCES = \
FGWeatherDefs.h FGWeatherFeature.h FGWeatherUtils.h \
FGWeatherParse.cpp FGWeatherParse.h \
FGWeatherVectorWrap.h \
FGWindItem.cpp FGWindItem.h
FGWindItem.cpp FGWindItem.h \
linintp2.cpp linintp2.h \
sphrintp.cpp sphrint.h
if OLD_AUTOMAKE
INCLUDES += -I$(top_srcdir) -I$(top_srcdir)/src

548
src/WeatherCM/linintp2.cpp Normal file
View file

@ -0,0 +1,548 @@
/*
WARNING - Do not remove this header.
This code is a templated version of the 'magic-software' spherical
interpolation code by Dave Eberly. The original (un-hacked) code can be
obtained from here: http://www.magic-software.com/gr_appr.htm
This code is derived from linintp2.h/cpp and sphrintp.h/cpp.
Dave Eberly says that the conditions for use are:
* You may distribute the original source code to others at no charge.
* You may modify the original source code and distribute it to others at
no charge. The modified code must be documented to indicate that it is
not part of the original package.
* You may use this code for non-commercial purposes. You may also
incorporate this code into commercial packages. However, you may not
sell any of your source code which contains my original and/or modified
source code. In such a case, you need to factor out my code and freely
distribute it.
* The original code comes with absolutely no warranty and no guarantee is
made that the code is bug-free.
This does not seem incompatible with GPL - so this modified version
is hereby placed under GPL along with the rest of FlightGear.
Christian Mayer
*/
#include <simgear/compiler.h>
#include <iostream>
#include <float.h>
#include <math.h>
#include <stdlib.h>
#include "linintp2.h"
SG_USING_NAMESPACE(std);
SG_USING_STD(cout);
//---------------------------------------------------------------------------
mgcLinInterp2D::mgcLinInterp2D (int _numPoints, double* x, double* y,
unsigned int* _f)
{
if ( (numPoints = _numPoints) < 3 )
{
point = 0;
edge = 0;
triangle = 0;
numTriangles = 0;
return;
}
cout << "[ 20%] allocating memory \r";
point = new double*[numPoints];
tmppoint = new double*[numPoints+3];
f = new unsigned int[numPoints];
int i;
for (i = 0; i < numPoints; i++)
point[i] = new double[2];
for (i = 0; i < numPoints+3; i++)
tmppoint[i] = new double[2];
for (i = 0; i < numPoints; i++)
{
point[i][0] = tmppoint[i][0] = x[i];
point[i][1] = tmppoint[i][1] = y[i];
f[i] = _f[i];
}
cout << "[ 30%] creating delaunay diagram \r";
Delaunay2D();
}
//---------------------------------------------------------------------------
mgcLinInterp2D::~mgcLinInterp2D ()
{
if ( numPoints < 3 )
return;
int i;
if ( point )
{
for (i = 0; i < numPoints; i++)
delete[] point[i];
delete[] point;
}
if ( tmppoint )
{
for (i = 0; i < numPoints+3; i++)
delete[] tmppoint[i];
delete[] tmppoint;
}
delete[] f;
delete[] edge;
delete[] triangle;
}
//---------------------------------------------------------------------------
void mgcLinInterp2D::ComputeBarycenter (Vertex& v0, Vertex& v1, Vertex& v2,
Vertex& ver, double c[3])
{
double a0 = v0.x-v2.x;
double b0 = v0.y-v2.y;
double a1 = v1.x-v2.x;
double b1 = v1.y-v2.y;
double a2 = ver.x-v2.x;
double b2 = ver.y-v2.y;
double m00 = a0*a0+b0*b0;
double m01 = a0*a1+b0*b1;
double m11 = a1*a1+b1*b1;
double r0 = a2*a0+b2*b0;
double r1 = a2*a1+b2*b1;
double det = m00*m11-m01*m01;
c[0] = (m11*r0-m01*r1)/det;
c[1] = (m00*r1-m01*r0)/det;
c[2] = 1-c[0]-c[1];
}
//---------------------------------------------------------------------------
int mgcLinInterp2D::InTriangle (Vertex& v0, Vertex& v1, Vertex& v2,
Vertex& test)
{
const double eps = 1e-08;
double tx, ty, nx, ny;
// test against normal to first edge
tx = test.x - v0.x;
ty = test.y - v0.y;
nx = v0.y - v1.y;
ny = v1.x - v0.x;
if ( tx*nx + ty*ny < -eps )
return 0;
// test against normal to second edge
tx = test.x - v1.x;
ty = test.y - v1.y;
nx = v1.y - v2.y;
ny = v2.x - v1.x;
if ( tx*nx + ty*ny < -eps )
return 0;
// test against normal to third edge
tx = test.x - v2.x;
ty = test.y - v2.y;
nx = v2.y - v0.y;
ny = v0.x - v2.x;
if ( tx*nx + ty*ny < -eps )
return 0;
return 1;
}
//---------------------------------------------------------------------------
int mgcLinInterp2D::Evaluate (double x, double y, EvaluateData& F)
{
Vertex ver = { x, y };
// determine which triangle contains the target point
int i;
Vertex v0, v1, v2;
for (i = 0; i < numTriangles; i++)
{
Triangle& t = triangle[i];
v0.x = point[t.vertex[0]][0];
v0.y = point[t.vertex[0]][1];
v1.x = point[t.vertex[1]][0];
v1.y = point[t.vertex[1]][1];
v2.x = point[t.vertex[2]][0];
v2.y = point[t.vertex[2]][1];
if ( InTriangle(v0,v1,v2,ver) )
break;
}
if ( i == numTriangles ) // point is outside interpolation region
{
return 0;
}
Triangle& t = triangle[i]; // (x,y) is in this triangle
// compute barycentric coordinates with respect to subtriangle
double bary[3];
ComputeBarycenter(v0,v1,v2,ver,bary);
// compute barycentric combination of function values at vertices
F.index[0] = f[t.vertex[0]];
F.index[1] = f[t.vertex[1]];
F.index[2] = f[t.vertex[2]];
F.percentage[0] = bary[0];
F.percentage[1] = bary[1];
F.percentage[2] = bary[2];
return 1;
}
//---------------------------------------------------------------------------
int mgcLinInterp2D::Delaunay2D ()
{
int result;
const double EPSILON = 1e-12;
const int TSIZE = 75;
const double RANGE = 10.0;
xmin = tmppoint[0][0];
xmax = xmin;
ymin = tmppoint[0][1];
ymax = ymin;
int i;
for (i = 0; i < numPoints; i++)
{
double value = tmppoint[i][0];
if ( xmax < value )
xmax = value;
if ( xmin > value )
xmin = value;
value = tmppoint[i][1];
if ( ymax < value )
ymax = value;
if ( ymin > value )
ymin = value;
}
double xrange = xmax-xmin, yrange = ymax-ymin;
double maxrange = xrange;
if ( maxrange < yrange )
maxrange = yrange;
// need to scale the data later to do a correct triangle count
double maxrange2 = maxrange*maxrange;
// tweak the points by very small random numbers
double bgs = EPSILON*maxrange;
srand(367);
for (i = 0; i < numPoints; i++)
{
tmppoint[i][0] += bgs*(0.5 - rand()/double(RAND_MAX));
tmppoint[i][1] += bgs*(0.5 - rand()/double(RAND_MAX));
}
double wrk[2][3] =
{
{ 5*RANGE, -RANGE, -RANGE },
{ -RANGE, 5*RANGE, -RANGE }
};
for (i = 0; i < 3; i++)
{
tmppoint[numPoints+i][0] = xmin+xrange*wrk[0][i];
tmppoint[numPoints+i][1] = ymin+yrange*wrk[1][i];
}
int i0, i1, i2, i3, i4, i5, i6, i7, i8, i9, i11;
int nts, ii[3];
double xx;
int tsz = 2*TSIZE;
int** tmp = new int*[tsz+1];
tmp[0] = new int[2*(tsz+1)];
for (i0 = 1; i0 < tsz+1; i0++)
tmp[i0] = tmp[0] + 2*i0;
i1 = 2*(numPoints + 2);
int* id = new int[i1];
for (i0 = 0; i0 < i1; i0++)
id[i0] = i0;
int** a3s = new int*[i1];
a3s[0] = new int[3*i1];
for (i0 = 1; i0 < i1; i0++)
a3s[i0] = a3s[0] + 3*i0;
a3s[0][0] = numPoints;
a3s[0][1] = numPoints+1;
a3s[0][2] = numPoints+2;
double** ccr = new double*[i1]; // circumscribed centers and radii
ccr[0] = new double[3*i1];
for (i0 = 1; i0 < i1; i0++)
ccr[i0] = ccr[0] + 3*i0;
ccr[0][0] = 0.0;
ccr[0][1] = 0.0;
ccr[0][2] = FLT_MAX;
nts = 1; // number of triangles
i4 = 1;
cout << "[ 40%] create triangulation \r";
// compute triangulation
for (i0 = 0; i0 < numPoints; i0++)
{
i1 = i7 = -1;
i9 = 0;
for (i11 = 0; i11 < nts; i11++)
{
i1++;
while ( a3s[i1][0] < 0 )
i1++;
xx = ccr[i1][2];
for (i2 = 0; i2 < 2; i2++)
{
double z = tmppoint[i0][i2]-ccr[i1][i2];
xx -= z*z;
if ( xx < 0 )
goto Corner3;
}
i9--;
i4--;
id[i4] = i1;
for (i2 = 0; i2 < 3; i2++)
{
ii[0] = 0;
if (ii[0] == i2)
ii[0]++;
for (i3 = 1; i3 < 2; i3++)
{
ii[i3] = ii[i3-1] + 1;
if (ii[i3] == i2)
ii[i3]++;
}
if ( i7 > 1 )
{
i8 = i7;
for (i3 = 0; i3 <= i8; i3++)
{
for (i5 = 0; i5 < 2; i5++)
if ( a3s[i1][ii[i5]] != tmp[i3][i5] )
goto Corner1;
for (i6 = 0; i6 < 2; i6++)
tmp[i3][i6] = tmp[i8][i6];
i7--;
goto Corner2;
Corner1:;
}
}
if ( ++i7 > tsz )
{
// temporary storage exceeded, increase TSIZE
result = 0;
goto ExitDelaunay;
}
for (i3 = 0; i3 < 2; i3++)
tmp[i7][i3] = a3s[i1][ii[i3]];
Corner2:;
}
a3s[i1][0] = -1;
Corner3:;
}
for (i1 = 0; i1 <= i7; i1++)
{
for (i2 = 0; i2 < 2; i2++)
for (wrk[i2][2] = 0, i3 = 0; i3 < 2; i3++)
{
wrk[i2][i3] = tmppoint[tmp[i1][i2]][i3]-tmppoint[i0][i3];
wrk[i2][2] +=
0.5*wrk[i2][i3]*(tmppoint[tmp[i1][i2]][i3]+
tmppoint[i0][i3]);
}
xx = wrk[0][0]*wrk[1][1]-wrk[1][0]*wrk[0][1];
ccr[id[i4]][0] = (wrk[0][2]*wrk[1][1]-wrk[1][2]*wrk[0][1])/xx;
ccr[id[i4]][1] = (wrk[0][0]*wrk[1][2]-wrk[1][0]*wrk[0][2])/xx;
for (ccr[id[i4]][2] = 0, i2 = 0; i2 < 2; i2++)
{
double z = tmppoint[i0][i2]-ccr[id[i4]][i2];
ccr[id[i4]][2] += z*z;
a3s[id[i4]][i2] = tmp[i1][i2];
}
a3s[id[i4]][2] = i0;
i4++;
i9++;
}
nts += i9;
}
// count the number of triangles
cout << "[ 50%] count the number of triangles \r";
numTriangles = 0;
i0 = -1;
for (i11 = 0; i11 < nts; i11++)
{
i0++;
while ( a3s[i0][0] < 0 )
i0++;
if ( a3s[i0][0] < numPoints )
{
for (i1 = 0; i1 < 2; i1++)
for (i2 = 0; i2 < 2; i2++)
wrk[i1][i2] =
tmppoint[a3s[i0][i1]][i2]-tmppoint[a3s[i0][2]][i2];
if ( fabs(wrk[0][0]*wrk[1][1]-wrk[0][1]*wrk[1][0]) > EPSILON*maxrange2 )
numTriangles++;
}
}
// create the triangles
cout << "[ 60%] create the triangles \r";
triangle = new Triangle[numTriangles];
numTriangles = 0;
i0 = -1;
for (i11 = 0; i11 < nts; i11++)
{
i0++;
while ( a3s[i0][0] < 0 )
i0++;
if ( a3s[i0][0] < numPoints )
{
for (i1 = 0; i1 < 2; i1++)
for (i2 = 0; i2 < 2; i2++)
wrk[i1][i2] =
tmppoint[a3s[i0][i1]][i2]-tmppoint[a3s[i0][2]][i2];
xx = wrk[0][0]*wrk[1][1]-wrk[0][1]*wrk[1][0];
if ( fabs(xx) > EPSILON*maxrange2 )
{
int delta = xx < 0 ? 1 : 0;
Triangle& tri = triangle[numTriangles];
tri.vertex[0] = a3s[i0][0];
tri.vertex[1] = a3s[i0][1+delta];
tri.vertex[2] = a3s[i0][2-delta];
tri.adj[0] = -1;
tri.adj[1] = -1;
tri.adj[2] = -1;
numTriangles++;
}
}
}
// build edge table
cout << "[ 70%] build the edge table \r";
numEdges = 0;
edge = new Edge[3*numTriangles];
int j, j0, j1;
for (i = 0; i < numTriangles; i++)
{
if ( (i%500) == 0)
cout << "[ 7" << 10*i/numTriangles << "%] build the edge table \r";
Triangle& t = triangle[i];
for (j0 = 0, j1 = 1; j0 < 3; j0++, j1 = (j1+1)%3)
{
for (j = 0; j < numEdges; j++)
{
Edge& e = edge[j];
if ( (t.vertex[j0] == e.vertex[0]
&& t.vertex[j1] == e.vertex[1])
|| (t.vertex[j0] == e.vertex[1]
&& t.vertex[j1] == e.vertex[0]) )
break;
}
if ( j == numEdges ) // add edge to table
{
edge[j].vertex[0] = t.vertex[j0];
edge[j].vertex[1] = t.vertex[j1];
edge[j].triangle[0] = i;
edge[j].index[0] = j0;
edge[j].triangle[1] = -1;
numEdges++;
}
else // edge already exists, add triangle to table
{
edge[j].triangle[1] = i;
edge[j].index[1] = j0;
}
}
}
// establish links between adjacent triangles
cout << "[ 80%] establishing links between adjacent triangles \r";
for (i = 0; i < numEdges; i++)
{
if ( edge[i].triangle[1] != -1 )
{
j0 = edge[i].triangle[0];
j1 = edge[i].triangle[1];
triangle[j0].adj[edge[i].index[0]] = j1;
triangle[j1].adj[edge[i].index[1]] = j0;
}
}
result = 1;
ExitDelaunay:;
delete[] tmp[0];
delete[] tmp;
delete[] id;
delete[] a3s[0];
delete[] a3s;
delete[] ccr[0];
delete[] ccr;
cout << "[ 90%] finsishes delauney triangulation \r";
return result;
}
//---------------------------------------------------------------------------
void mgcLinInterp2D::GetPoint (int i, double& x, double& y)
{
// assumes i is valid [can use PointCount() before passing i]
x = point[i][0];
y = point[i][1];
}
//---------------------------------------------------------------------------
void mgcLinInterp2D::GetEdge (int i, double& x0, double& y0, double& x1,
double& y1)
{
// assumes i is valid [can use EdgeCount() before passing i]
int v0 = edge[i].vertex[0], v1 = edge[i].vertex[1];
x0 = point[v0][0];
y0 = point[v0][1];
x1 = point[v1][0];
y1 = point[v1][1];
}
//---------------------------------------------------------------------------
void mgcLinInterp2D::GetTriangle (int i, double& x0, double& y0, double& x1,
double& y1, double& x2, double& y2)
{
// assumes i is valid [can use TriangleCount() before passing i]
int v0 = triangle[i].vertex[0];
int v1 = triangle[i].vertex[1];
int v2 = triangle[i].vertex[2];
x0 = point[v0][0];
y0 = point[v0][1];
x1 = point[v1][0];
y1 = point[v1][1];
x2 = point[v2][0];
y2 = point[v2][1];
}
//---------------------------------------------------------------------------

15
src/WeatherCM/linintp2.h
View file

@ -32,11 +32,16 @@
#ifndef LININTP2_H
#define LININTP2_H
template<class T>
struct EvaluateData
{
unsigned int index[3];
double percentage[3];
};
class mgcLinInterp2D
{
public:
mgcLinInterp2D (int _numPoints, double* x, double* y, T* _f);
mgcLinInterp2D (int _numPoints, double* x, double* y, unsigned int* _f);
~mgcLinInterp2D ();
@ -57,7 +62,7 @@ public:
void GetTriangle (int i, double& x0, double& y0, double& x1, double& y1,
double& x2, double& y2);
int Evaluate (double x, double y, T& F);
int Evaluate (double x, double y, EvaluateData& F);
private:
typedef struct
@ -88,7 +93,7 @@ private:
int numPoints;
double** point;
double** tmppoint;
T* f;
unsigned int* f;
double xmin, xmax, ymin, ymax;
@ -105,6 +110,4 @@ private:
int InTriangle (Vertex& v0, Vertex& v1, Vertex& v2, Vertex& test);
};
#include "linintp2.inl"
#endif

172
src/WeatherCM/sphrintp.cpp Normal file
View file

@ -0,0 +1,172 @@
/*
WARNING - Do not remove this header.
This code is a templated version of the 'magic-software' spherical
interpolation code by Dave Eberly. The original (un-hacked) code can be
obtained from here: http://www.magic-software.com/gr_appr.htm
This code is derived from linintp2.h/cpp and sphrintp.h/cpp.
Dave Eberly says that the conditions for use are:
* You may distribute the original source code to others at no charge.
* You may modify the original source code and distribute it to others at
no charge. The modified code must be documented to indicate that it is
not part of the original package.
* You may use this code for non-commercial purposes. You may also
incorporate this code into commercial packages. However, you may not
sell any of your source code which contains my original and/or modified
source code. In such a case, you need to factor out my code and freely
distribute it.
* The original code comes with absolutely no warranty and no guarantee is
made that the code is bug-free.
This does not seem incompatible with GPL - so this modified version
is hereby placed under GPL along with the rest of FlightGear.
Christian Mayer
*/
#include <simgear/compiler.h>
#include <iostream>
#include <math.h>
#include "sphrintp.h"
SG_USING_NAMESPACE(std);
SG_USING_STD(cout);
static const double PI = 4.0*atan(1.0);
static const double TWOPI = 2.0*PI;
//---------------------------------------------------------------------------
SphereInterpolate::SphereInterpolate (int n, const double* x,
const double* y, const double* z,
const unsigned int* f)
{
// Assumes (x[i],y[i],z[i]) is unit length for all 0 <= i < n.
// For complete spherical coverage, include the two antipodal points
// (0,0,1,f(0,0,1)) and (0,0,-1,f(0,0,-1)) in the data set.
cout << "Initialising spherical interpolator.\n";
cout << "[ 0%] Allocating memory \r";
theta = new double[3*n];
phi = new double[3*n];
func = new unsigned int[3*n];
// convert data to spherical coordinates
int i;
for (i = 0; i < n; i++)
{
GetSphericalCoords(x[i],y[i],z[i],theta[i],phi[i]);
func[i] = f[i];
}
// use periodicity to get wrap-around in the Delaunay triangulation
cout << "[ 10%] copying vertices for wrap-around\r";
int j, k;
for (i = 0, j = n, k = 2*n; i < n; i++, j++, k++)
{
theta[j] = theta[i]+TWOPI;
theta[k] = theta[i]-TWOPI;
phi[j] = phi[i];
phi[k] = phi[i];
func[j] = func[i];
func[k] = func[i];
}
pInterp = new mgcLinInterp2D(3*n,theta,phi,func);
cout << "[100%] Finished initialising spherical interpolator. \n";
}
SphereInterpolate::SphereInterpolate (int n, const sgVec2* p, const unsigned int* f)
{
// Assumes (x[i],y[i],z[i]) is unit length for all 0 <= i < n.
// For complete spherical coverage, include the two antipodal points
// (0,0,1,f(0,0,1)) and (0,0,-1,f(0,0,-1)) in the data set.
cout << "Initialising spherical interpolator.\n";
cout << "[ 0%] Allocating memory \r";
theta = new double[3*n];
phi = new double[3*n];
func = new unsigned int[3*n];
// convert data to spherical coordinates
cout << "[ 10%] copying vertices for wrap-around \r";
int i, j, k;
for (i = 0, j = n, k = 2*n; i < n; i++, j++, k++)
{
phi[i] = p[i][0];
theta[i] = p[i][1];
func[i] = f[i];
// use periodicity to get wrap-around in the Delaunay triangulation
phi[j] = phi[i];
phi[k] = phi[i];
theta[j] = theta[i]+TWOPI;
theta[k] = theta[i]-TWOPI;
func[j] = func[i];
func[k] = func[i];
}
pInterp = new mgcLinInterp2D(3*n,theta,phi,func);
cout << "[100%] Finished initialising spherical interpolator. \n";
}
//---------------------------------------------------------------------------
SphereInterpolate::~SphereInterpolate ()
{
delete pInterp;
delete[] theta;
delete[] phi;
delete[] func;
}
//---------------------------------------------------------------------------
void SphereInterpolate::GetSphericalCoords (const double x, const double y, const double z,
double& thetaAngle,
double& phiAngle) const
{
// Assumes (x,y,z) is unit length. Returns -PI <= thetaAngle <= PI
// and 0 <= phiAngle <= PI.
if ( z < 1.0f )
{
if ( z > -1.0f )
{
thetaAngle = atan2(y,x);
phiAngle = acos(z);
}
else
{
thetaAngle = -PI;
phiAngle = PI;
}
}
else
{
thetaAngle = -PI;
phiAngle = 0.0f;
}
}
//---------------------------------------------------------------------------
int SphereInterpolate::Evaluate (const double x, const double y, const double z, EvaluateData& f) const
{
// assumes (x,y,z) is unit length
double thetaAngle, phiAngle;
GetSphericalCoords(x,y,z,thetaAngle,phiAngle);
return pInterp->Evaluate(thetaAngle,phiAngle,f);
}
//---------------------------------------------------------------------------
int SphereInterpolate::Evaluate (const double thetaAngle, const double phiAngle, EvaluateData& f) const
{
return pInterp->Evaluate(thetaAngle,phiAngle,f);
}
//---------------------------------------------------------------------------

42
src/WeatherCM/sphrintp.h
View file

@ -32,43 +32,59 @@
#ifndef SPHRINTP_H
#define SPHRINTP_H
#include <simgear/compiler.h>
#include <iostream>
#include "linintp2.h"
#include <plib/sg.h>
template<class T>
SG_USING_NAMESPACE(std);
SG_USING_STD(cout);
class SphereInterpolate
{
public:
SphereInterpolate (int n, const double* x, const double* y,
const double* z, const T* f);
SphereInterpolate (int n, const sgVec2* p, const T* f);
const double* z, const unsigned int* f);
SphereInterpolate (int n, const sgVec2* p, const unsigned int* f);
~SphereInterpolate ();
void GetSphericalCoords (const double x, const double y, const double z,
double& thetaAngle, double& phiAngle) const;
int Evaluate (const double x, const double y, const double z, T& f) const;
int Evaluate (const double thetaAngle, const double phiAngle, T& f) const;
int Evaluate (const double x, const double y, const double z, EvaluateData& f) const;
int Evaluate (const double thetaAngle, const double phiAngle, EvaluateData& f) const;
#ifndef macintosh
// CodeWarrior doesn't know the differece between sgVec2 and
// sgVec3, so I commented this out for Mac builds. This change is
// related to a similar change in FGLocalWeatherDatabase module.
T Evaluate(const sgVec2& p) const
EvaluateData Evaluate(const sgVec2& p) const
{
T retval;
EvaluateData retval;
Evaluate(p[1], p[0], retval);
return retval;
}
#endif
T Evaluate(const sgVec3& p) const
EvaluateData Evaluate(const sgVec3& p) const
{
T retval;
Evaluate(p[1], p[0], retval);
EvaluateData retval;
if (!Evaluate(p[1], p[0], retval))
{
cout << "Error during spherical interpolation. Point ("
<< p[0] << "/" << p[1] << "/" << p[2] << ") was in no triangle\n";
retval.index[0] = 0; //fake something
retval.index[1] = 0;
retval.index[2] = 0;
retval.percentage[0] = 1.0;
retval.percentage[1] = 0.0;
retval.percentage[2] = 0.0;
}
return retval;
}
@ -76,10 +92,8 @@ protected:
int numPoints;
double* theta;
double* phi;
T* func;
mgcLinInterp2D<T>* pInterp;
unsigned int* func;
mgcLinInterp2D* pInterp;
};
#include "sphrintp.inl"
#endif