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);
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
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();
{
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)
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();
// 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)
//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;
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;
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);
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;
}
/****************************************************************************/
// 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
#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) /
#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"
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;
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 */
/************************************************************************/
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();
return AirPressure.getValue() * exp(FF);
}
-
-
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
--- /dev/null
+/*****************************************************************************
+
+ 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;
+}
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);
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);
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)
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)
}
/****************************************************************************/
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);
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 \
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
--- /dev/null
+/*
+ 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];
+}
+//---------------------------------------------------------------------------
+
#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 ();
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
int numPoints;
double** point;
double** tmppoint;
- T* f;
+ unsigned int* f;
double xmin, xmax, ymin, ymax;
int InTriangle (Vertex& v0, Vertex& v1, Vertex& v2, Vertex& test);
};
-#include "linintp2.inl"
-
#endif
--- /dev/null
+/*
+ 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);
+}
+//---------------------------------------------------------------------------
#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;
}
int numPoints;
double* theta;
double* phi;
- T* func;
- mgcLinInterp2D<T>* pInterp;
+ unsigned int* func;
+ mgcLinInterp2D* pInterp;
};
-#include "sphrintp.inl"
-
#endif
projects / flightgear.git / commitdiff