From d4c49d65ac3a23203f4dc785bb6e3b1109afffc6 Mon Sep 17 00:00:00 2001 From: david Date: Sat, 16 Mar 2002 20:31:27 +0000 Subject: [PATCH] Major weather update from Christian Mayer, tying the weather code into 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. --- src/WeatherCM/FGAirPressureItem.h | 9 +- src/WeatherCM/FGLocalWeatherDatabase.cpp | 139 ++++- src/WeatherCM/FGLocalWeatherDatabase.h | 38 +- src/WeatherCM/FGPhysicalProperties.cpp | 33 +- src/WeatherCM/FGPhysicalProperties.h | 26 + src/WeatherCM/FGPhysicalProperties_bind.cpp | 480 +++++++++++++++++ src/WeatherCM/FGTemperatureItem.h | 3 + src/WeatherCM/FGTurbulenceItem.h | 4 + src/WeatherCM/FGWeatherUtils.h | 4 +- src/WeatherCM/FGWindItem.h | 4 + src/WeatherCM/Makefile.am | 7 +- src/WeatherCM/linintp2.cpp | 548 ++++++++++++++++++++ src/WeatherCM/linintp2.h | 15 +- src/WeatherCM/sphrintp.cpp | 172 ++++++ src/WeatherCM/sphrintp.h | 42 +- 15 files changed, 1475 insertions(+), 49 deletions(-) create mode 100644 src/WeatherCM/FGPhysicalProperties_bind.cpp create mode 100644 src/WeatherCM/linintp2.cpp create mode 100644 src/WeatherCM/sphrintp.cpp diff --git a/src/WeatherCM/FGAirPressureItem.h b/src/WeatherCM/FGAirPressureItem.h index 5c80eefc4..21192904b 100644 --- a/src/WeatherCM/FGAirPressureItem.h +++ b/src/WeatherCM/FGAirPressureItem.h @@ -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); diff --git a/src/WeatherCM/FGLocalWeatherDatabase.cpp b/src/WeatherCM/FGLocalWeatherDatabase.cpp index ee7d4e88f..0ad64b9fa 100644 --- a/src/WeatherCM/FGLocalWeatherDatabase.cpp +++ b/src/WeatherCM/FGLocalWeatherDatabase.cpp @@ -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::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(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(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::iterator it = database_data.begin(); it != database_data.end(); it++) + { + SGPropertyNode * station = station_nodes->getNode("station", index, true); + + station -> tie("air-pressure-Pa", + SGRawValueMethods( + 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( + database_data[index], i, + &FGPhysicalProperties::getWind_x, + &FGPhysicalProperties::setWind_x) + ,false); + wind -> tie("east-mps", + SGRawValueMethodsIndexed( + database_data[index], i, + &FGPhysicalProperties::getWind_y, + &FGPhysicalProperties::setWind_y) + ,false); + wind -> tie("up-mps", + SGRawValueMethodsIndexed( + database_data[index], i, + &FGPhysicalProperties::getWind_z, + &FGPhysicalProperties::setWind_z) + ,false); + wind -> tie("altitude-m", + SGRawValueMethodsIndexed( + 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( + database_data[index], i, + &FGPhysicalProperties::getTemperature_x, + &FGPhysicalProperties::setTemperature_x) + ,false); + temperature -> tie("altitude-m", + SGRawValueMethodsIndexed( + 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( + database_data[index], i, + &FGPhysicalProperties::getVaporPressure_x, + &FGPhysicalProperties::setVaporPressure_x) + ,false); + vaporpressure -> tie("altitude-m", + SGRawValueMethodsIndexed( + 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) / diff --git a/src/WeatherCM/FGLocalWeatherDatabase.h b/src/WeatherCM/FGLocalWeatherDatabase.h index 655356aeb..60063fa17 100644 --- a/src/WeatherCM/FGLocalWeatherDatabase.h +++ b/src/WeatherCM/FGLocalWeatherDatabase.h @@ -60,8 +60,10 @@ HISTORY #include -#include +#include
+ #include +#include #include "sphrintp.h" @@ -92,11 +94,33 @@ struct _FGLocalWeatherDatabaseCache FGPhysicalProperty last_known_property; }; -class FGLocalWeatherDatabase +class FGLocalWeatherDatabase : public FGSubsystem { private: protected: - SphereInterpolate *database; + SphereInterpolate *database_logic; + vector 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 pointVector; typedef vector 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 */ /************************************************************************/ diff --git a/src/WeatherCM/FGPhysicalProperties.cpp b/src/WeatherCM/FGPhysicalProperties.cpp index 2802ebbbf..636ff130e 100644 --- a/src/WeatherCM/FGPhysicalProperties.cpp +++ b/src/WeatherCM/FGPhysicalProperties.cpp @@ -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); } - - diff --git a/src/WeatherCM/FGPhysicalProperties.h b/src/WeatherCM/FGPhysicalProperties.h index 1e20e35ce..2e39d5b6b 100644 --- a/src/WeatherCM/FGPhysicalProperties.h +++ b/src/WeatherCM/FGPhysicalProperties.h @@ -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 diff --git a/src/WeatherCM/FGPhysicalProperties_bind.cpp b/src/WeatherCM/FGPhysicalProperties_bind.cpp new file mode 100644 index 000000000..cd351c69e --- /dev/null +++ b/src/WeatherCM/FGPhysicalProperties_bind.cpp @@ -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::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::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::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::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::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::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::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::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::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::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::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::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::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::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::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::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::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::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::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::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::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::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::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::iterator it = VaporPressure.begin(); + + while( number > 0) + { + it++; + number--; + } + + it->first; +} diff --git a/src/WeatherCM/FGTemperatureItem.h b/src/WeatherCM/FGTemperatureItem.h index 14f19934e..f6cc5e04e 100644 --- a/src/WeatherCM/FGTemperatureItem.h +++ b/src/WeatherCM/FGTemperatureItem.h @@ -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); diff --git a/src/WeatherCM/FGTurbulenceItem.h b/src/WeatherCM/FGTurbulenceItem.h index 379748d07..1afabd955 100644 --- a/src/WeatherCM/FGTurbulenceItem.h +++ b/src/WeatherCM/FGTurbulenceItem.h @@ -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); diff --git a/src/WeatherCM/FGWeatherUtils.h b/src/WeatherCM/FGWeatherUtils.h index 59283d3e4..505fcdc58 100644 --- a/src/WeatherCM/FGWeatherUtils.h +++ b/src/WeatherCM/FGWeatherUtils.h @@ -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) } /****************************************************************************/ diff --git a/src/WeatherCM/FGWindItem.h b/src/WeatherCM/FGWindItem.h index 3142de0f2..9aa6ab395 100644 --- a/src/WeatherCM/FGWindItem.h +++ b/src/WeatherCM/FGWindItem.h @@ -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); diff --git a/src/WeatherCM/Makefile.am b/src/WeatherCM/Makefile.am index 036aca2e1..2ac495940 100644 --- a/src/WeatherCM/Makefile.am +++ b/src/WeatherCM/Makefile.am @@ -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 diff --git a/src/WeatherCM/linintp2.cpp b/src/WeatherCM/linintp2.cpp new file mode 100644 index 000000000..0c321f23c --- /dev/null +++ b/src/WeatherCM/linintp2.cpp @@ -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 +#include + +#include +#include +#include +#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]; +} +//--------------------------------------------------------------------------- + diff --git a/src/WeatherCM/linintp2.h b/src/WeatherCM/linintp2.h index f0e8a1a9e..c0925b440 100644 --- a/src/WeatherCM/linintp2.h +++ b/src/WeatherCM/linintp2.h @@ -32,11 +32,16 @@ #ifndef LININTP2_H #define LININTP2_H -template +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 diff --git a/src/WeatherCM/sphrintp.cpp b/src/WeatherCM/sphrintp.cpp new file mode 100644 index 000000000..41ff38c83 --- /dev/null +++ b/src/WeatherCM/sphrintp.cpp @@ -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 +#include + +#include +#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); +} +//--------------------------------------------------------------------------- diff --git a/src/WeatherCM/sphrintp.h b/src/WeatherCM/sphrintp.h index e5debcab7..4f0e781c9 100644 --- a/src/WeatherCM/sphrintp.h +++ b/src/WeatherCM/sphrintp.h @@ -32,43 +32,59 @@ #ifndef SPHRINTP_H #define SPHRINTP_H +#include +#include + #include "linintp2.h" #include -template +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* pInterp; + unsigned int* func; + mgcLinInterp2D* pInterp; }; -#include "sphrintp.inl" - #endif