FGFS-4.1/simgear
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Modified Files:

Browse Source 
projects/VC7.1/SimGear.vcproj projects/VC8/SimGear.vcproj
	simgear/math/Makefile.am simgear/math/SGGeoc.hxx
	simgear/math/SGGeodesy.cxx simgear/math/SGGeodesy.hxx
	simgear/math/polar3d.hxx simgear/math/sg_geodesy.hxx
	simgear/math/sg_types.hxx
Removed Files:
	simgear/math/polar3d.cxx simgear/math/sg_geodesy.cxx
	simgear/math/sg_memory.h:
	Remove sg_memory.h It is unused anyway and should not be required
	in a c++ world. Move distance course functions to the SG* type
	system. Move the implementation into SGGeodesy.cxx. Remove some of
	the old Point3D Based sg* functions that are already unused.
This commit is contained in:
frohlich
2007-08-07 05:26:21 +00:00
parent a8a02d3a2e
commit 30529ccdf5
12 changed files with 423 additions and 712 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
projects/VC7.1/SimGear.vcproj
View File

@@ -307,21 +307,12 @@
<File
RelativePath="..\..\simgear\math\point3d.hxx">
</File>
<File
RelativePath="..\..\simgear\math\polar3d.cxx">
</File>
<File
RelativePath="..\..\simgear\math\polar3d.hxx">
</File>
<File
RelativePath="..\..\simgear\math\sg_geodesy.cxx">
</File>
<File
RelativePath="..\..\simgear\math\sg_geodesy.hxx">
</File>
<File
RelativePath="..\..\simgear\math\sg_memory.h">
</File>
<File
RelativePath="..\..\simgear\math\sg_random.c">
</File>

12
projects/VC8/SimGear.vcproj
View File

@@ -467,10 +467,6 @@
RelativePath="..\..\simgear\sg_inlines.h"
>
</File>
<File
RelativePath="..\..\simgear\math\sg_memory.h"
>
</File>
<File
RelativePath="..\..\simgear\misc\sg_path.hxx"
>
@@ -1039,10 +1035,6 @@
RelativePath="..\..\simgear\scene\model\placementtrans.cxx"
>
</File>
<File
RelativePath="..\..\simgear\math\polar3d.cxx"
>
</File>
<File
RelativePath="..\..\simgear\props\props.cxx"
>
@@ -1087,10 +1079,6 @@
RelativePath="..\..\simgear\io\sg_file.cxx"
>
</File>
<File
RelativePath="..\..\simgear\math\sg_geodesy.cxx"
>
</File>
<File
RelativePath="..\..\simgear\misc\sg_path.cxx"
>

3
simgear/math/Makefile.am
View File

@@ -17,7 +17,6 @@ include_HEADERS = \
point3d.hxx \
polar3d.hxx \
sg_geodesy.hxx \
sg_memory.h \
sg_random.h \
sg_types.hxx \
vector.hxx \
@@ -47,8 +46,6 @@ include_HEADERS = \
libsgmath_a_SOURCES = \
interpolater.cxx \
leastsqs.cxx \
polar3d.cxx \
sg_geodesy.cxx \
sg_random.c \
vector.cxx \
SGGeodesy.cxx

35
simgear/math/SGGeoc.hxx
View File

@@ -74,6 +74,11 @@ public:
/// Set the geocentric radius from the argument given in feet
void setRadiusFt(double radius);
SGGeoc advanceRadM(double course, double distance) const;
static double courseRad(const SGGeoc& from, const SGGeoc& to);
static double courseDeg(const SGGeoc& from, const SGGeoc& to);
static double distanceM(const SGGeoc& from, const SGGeoc& to);
private:
/// This one is private since construction is not unique if you do
/// not know the units of the arguments, use the factory methods for
@@ -288,6 +293,36 @@ SGGeoc::setRadiusFt(double radius)
_radius = radius*SG_FEET_TO_METER;
}
inline
SGGeoc
SGGeoc::advanceRadM(double course, double distance) const
{
SGGeoc result;
SGGeodesy::advanceRadM(*this, course, distance, result);
return result;
}
inline
double
SGGeoc::courseRad(const SGGeoc& from, const SGGeoc& to)
{
return SGGeodesy::courseRad(from, to);
}
inline
double
SGGeoc::courseDeg(const SGGeoc& from, const SGGeoc& to)
{
return SGMiscd::rad2deg(courseRad(from, to));
}
inline
double
SGGeoc::distanceM(const SGGeoc& from, const SGGeoc& to)
{
return SGGeodesy::distanceM(from, to);
}
/// Output to an ostream
template<typename char_type, typename traits_type>
inline

343
simgear/math/SGGeodesy.cxx
View File

@@ -155,3 +155,346 @@ SGGeodesy::SGGeocToCart(const SGGeoc& geoc, SGVec3<double>& cart)
double clon = cos(lon);
cart = geoc.getRadiusM()*SGVec3<double>(clat*clon, clat*slon, slat);
}
// Notes:
//
// The XYZ/cartesian coordinate system in use puts the X axis through
// zero lat/lon (off west Africa), the Z axis through the north pole,
// and the Y axis through 90 degrees longitude (in the Indian Ocean).
//
// All latitude and longitude values are in radians. Altitude is in
// meters, with zero on the WGS84 ellipsoid.
//
// The code below makes use of the notion of "squashed" space. This
// is a 2D cylindrical coordinate system where the radius from the Z
// axis is multiplied by SQUASH; the earth in this space is a perfect
// circle with a radius of POLRAD.
////////////////////////////////////////////////////////////////////////
//
// Direct and inverse distance functions
//
// Proceedings of the 7th International Symposium on Geodetic
// Computations, 1985
//
// "The Nested Coefficient Method for Accurate Solutions of Direct and
// Inverse Geodetic Problems With Any Length"
//
// Zhang Xue-Lian
// pp 747-763
//
// modified for FlightGear to use WGS84 only -- Norman Vine
static inline double M0( double e2 ) {
//double e4 = e2*e2;
return SGMiscd::pi()*0.5*(1.0 - e2*( 1.0/4.0 + e2*( 3.0/64.0 +
e2*(5.0/256.0) )));
}
// given, lat1, lon1, az1 and distance (s), calculate lat2, lon2
// and az2. Lat, lon, and azimuth are in degrees. distance in meters
static int _geo_direct_wgs_84 ( double lat1, double lon1, double az1,
double s, double *lat2, double *lon2,
double *az2 )
{
double a = SGGeodesy::EQURAD, rf = SGGeodesy::iFLATTENING;
double testv = 1.0E-10;
double f = ( rf > 0.0 ? 1.0/rf : 0.0 );
double b = a*(1.0-f);
double e2 = f*(2.0-f);
double phi1 = SGMiscd::deg2rad(lat1), lam1 = SGMiscd::deg2rad(lon1);
double sinphi1 = sin(phi1), cosphi1 = cos(phi1);
double azm1 = SGMiscd::deg2rad(az1);
double sinaz1 = sin(azm1), cosaz1 = cos(azm1);
if( fabs(s) < 0.01 ) { // distance < centimeter => congruency
*lat2 = lat1;
*lon2 = lon1;
*az2 = 180.0 + az1;
if( *az2 > 360.0 ) *az2 -= 360.0;
return 0;
} else if( SGLimitsd::min() < fabs(cosphi1) ) { // non-polar origin
// u1 is reduced latitude
double tanu1 = sqrt(1.0-e2)*sinphi1/cosphi1;
double sig1 = atan2(tanu1,cosaz1);
double cosu1 = 1.0/sqrt( 1.0 + tanu1*tanu1 ), sinu1 = tanu1*cosu1;
double sinaz = cosu1*sinaz1, cos2saz = 1.0-sinaz*sinaz;
double us = cos2saz*e2/(1.0-e2);
// Terms
double ta = 1.0+us*(4096.0+us*(-768.0+us*(320.0-175.0*us)))/16384.0,
tb = us*(256.0+us*(-128.0+us*(74.0-47.0*us)))/1024.0,
tc = 0;
// FIRST ESTIMATE OF SIGMA (SIG)
double first = s/(b*ta); // !!
double sig = first;
double c2sigm, sinsig,cossig, temp,denom,rnumer, dlams, dlam;
do {
c2sigm = cos(2.0*sig1+sig);
sinsig = sin(sig); cossig = cos(sig);
temp = sig;
sig = first +
tb*sinsig*(c2sigm+tb*(cossig*(-1.0+2.0*c2sigm*c2sigm) -
tb*c2sigm*(-3.0+4.0*sinsig*sinsig)
*(-3.0+4.0*c2sigm*c2sigm)/6.0)
/4.0);
} while( fabs(sig-temp) > testv);
// LATITUDE OF POINT 2
// DENOMINATOR IN 2 PARTS (TEMP ALSO USED LATER)
temp = sinu1*sinsig-cosu1*cossig*cosaz1;
denom = (1.0-f)*sqrt(sinaz*sinaz+temp*temp);
// NUMERATOR
rnumer = sinu1*cossig+cosu1*sinsig*cosaz1;
*lat2 = SGMiscd::rad2deg(atan2(rnumer,denom));
// DIFFERENCE IN LONGITUDE ON AUXILARY SPHERE (DLAMS )
rnumer = sinsig*sinaz1;
denom = cosu1*cossig-sinu1*sinsig*cosaz1;
dlams = atan2(rnumer,denom);
// TERM C
tc = f*cos2saz*(4.0+f*(4.0-3.0*cos2saz))/16.0;
// DIFFERENCE IN LONGITUDE
dlam = dlams-(1.0-tc)*f*sinaz*(sig+tc*sinsig*
(c2sigm+
tc*cossig*(-1.0+2.0*
c2sigm*c2sigm)));
*lon2 = SGMiscd::rad2deg(lam1+dlam);
if (*lon2 > 180.0 ) *lon2 -= 360.0;
if (*lon2 < -180.0 ) *lon2 += 360.0;
// AZIMUTH - FROM NORTH
*az2 = SGMiscd::rad2deg(atan2(-sinaz,temp));
if ( fabs(*az2) < testv ) *az2 = 0.0;
if( *az2 < 0.0) *az2 += 360.0;
return 0;
} else { // phi1 == 90 degrees, polar origin
double dM = a*M0(e2) - s;
double paz = ( phi1 < 0.0 ? 180.0 : 0.0 );
double zero = 0.0f;
return _geo_direct_wgs_84( zero, lon1, paz, dM, lat2, lon2, az2 );
}
}
bool
SGGeodesy::direct(const SGGeod& p1, double course1,
double distance, SGGeod& p2, double& course2)
{
double lat2, lon2;
int ret = _geo_direct_wgs_84(p1.getLatitudeDeg(), p1.getLongitudeDeg(),
course1, distance, &lat2, &lon2, &course2);
p2.setLatitudeDeg(lat2);
p2.setLongitudeDeg(lon2);
p2.setElevationM(0);
return ret == 0;
}
// given lat1, lon1, lat2, lon2, calculate starting and ending
// az1, az2 and distance (s). Lat, lon, and azimuth are in degrees.
// distance in meters
static int _geo_inverse_wgs_84( double lat1, double lon1, double lat2,
double lon2, double *az1, double *az2,
double *s )
{
double a = SGGeodesy::EQURAD, rf = SGGeodesy::iFLATTENING;
int iter=0;
double testv = 1.0E-10;
double f = ( rf > 0.0 ? 1.0/rf : 0.0 );
double b = a*(1.0-f);
// double e2 = f*(2.0-f); // unused in this routine
double phi1 = SGMiscd::deg2rad(lat1), lam1 = SGMiscd::deg2rad(lon1);
double sinphi1 = sin(phi1), cosphi1 = cos(phi1);
double phi2 = SGMiscd::deg2rad(lat2), lam2 = SGMiscd::deg2rad(lon2);
double sinphi2 = sin(phi2), cosphi2 = cos(phi2);
if( (fabs(lat1-lat2) < testv &&
( fabs(lon1-lon2) < testv) || fabs(lat1-90.0) < testv ) )
{
// TWO STATIONS ARE IDENTICAL : SET DISTANCE & AZIMUTHS TO ZERO */
*az1 = 0.0; *az2 = 0.0; *s = 0.0;
return 0;
} else if( fabs(cosphi1) < testv ) {
// initial point is polar
int k = _geo_inverse_wgs_84( lat2,lon2,lat1,lon1, az1,az2,s );
k = k; // avoid compiler error since return result is unused
b = *az1; *az1 = *az2; *az2 = b;
return 0;
} else if( fabs(cosphi2) < testv ) {
// terminal point is polar
double _lon1 = lon1 + 180.0f;
int k = _geo_inverse_wgs_84( lat1, lon1, lat1, _lon1,
az1, az2, s );
k = k; // avoid compiler error since return result is unused
*s /= 2.0;
*az2 = *az1 + 180.0;
if( *az2 > 360.0 ) *az2 -= 360.0;
return 0;
} else if( (fabs( fabs(lon1-lon2) - 180 ) < testv) &&
(fabs(lat1+lat2) < testv) )
{
// Geodesic passes through the pole (antipodal)
double s1,s2;
_geo_inverse_wgs_84( lat1,lon1, lat1,lon2, az1,az2, &s1 );
_geo_inverse_wgs_84( lat2,lon2, lat1,lon2, az1,az2, &s2 );
*az2 = *az1;
*s = s1 + s2;
return 0;
} else {
// antipodal and polar points don't get here
double dlam = lam2 - lam1, dlams = dlam;
double sdlams,cdlams, sig,sinsig,cossig, sinaz,
cos2saz, c2sigm;
double tc,temp, us,rnumer,denom, ta,tb;
double cosu1,sinu1, sinu2,cosu2;
// Reduced latitudes
temp = (1.0-f)*sinphi1/cosphi1;
cosu1 = 1.0/sqrt(1.0+temp*temp);
sinu1 = temp*cosu1;
temp = (1.0-f)*sinphi2/cosphi2;
cosu2 = 1.0/sqrt(1.0+temp*temp);
sinu2 = temp*cosu2;
do {
sdlams = sin(dlams), cdlams = cos(dlams);
sinsig = sqrt(cosu2*cosu2*sdlams*sdlams+
(cosu1*sinu2-sinu1*cosu2*cdlams)*
(cosu1*sinu2-sinu1*cosu2*cdlams));
cossig = sinu1*sinu2+cosu1*cosu2*cdlams;
sig = atan2(sinsig,cossig);
sinaz = cosu1*cosu2*sdlams/sinsig;
cos2saz = 1.0-sinaz*sinaz;
c2sigm = (sinu1 == 0.0 || sinu2 == 0.0 ? cossig :
cossig-2.0*sinu1*sinu2/cos2saz);
tc = f*cos2saz*(4.0+f*(4.0-3.0*cos2saz))/16.0;
temp = dlams;
dlams = dlam+(1.0-tc)*f*sinaz*
(sig+tc*sinsig*
(c2sigm+tc*cossig*(-1.0+2.0*c2sigm*c2sigm)));
if (fabs(dlams) > SGMiscd::pi() && iter++ > 50) {
return iter;
}
} while ( fabs(temp-dlams) > testv);
us = cos2saz*(a*a-b*b)/(b*b); // !!
// BACK AZIMUTH FROM NORTH
rnumer = -(cosu1*sdlams);
denom = sinu1*cosu2-cosu1*sinu2*cdlams;
*az2 = SGMiscd::rad2deg(atan2(rnumer,denom));
if( fabs(*az2) < testv ) *az2 = 0.0;
if(*az2 < 0.0) *az2 += 360.0;
// FORWARD AZIMUTH FROM NORTH
rnumer = cosu2*sdlams;
denom = cosu1*sinu2-sinu1*cosu2*cdlams;
*az1 = SGMiscd::rad2deg(atan2(rnumer,denom));
if( fabs(*az1) < testv ) *az1 = 0.0;
if(*az1 < 0.0) *az1 += 360.0;
// Terms a & b
ta = 1.0+us*(4096.0+us*(-768.0+us*(320.0-175.0*us)))/
16384.0;
tb = us*(256.0+us*(-128.0+us*(74.0-47.0*us)))/1024.0;
// GEODETIC DISTANCE
*s = b*ta*(sig-tb*sinsig*
(c2sigm+tb*(cossig*(-1.0+2.0*c2sigm*c2sigm)-tb*
c2sigm*(-3.0+4.0*sinsig*sinsig)*
(-3.0+4.0*c2sigm*c2sigm)/6.0)/
4.0));
return 0;
}
}
bool
SGGeodesy::inverse(const SGGeod& p1, const SGGeod& p2, double& course1,
double& course2, double& distance)
{
int ret = _geo_inverse_wgs_84(p1.getLatitudeDeg(), p1.getLongitudeDeg(),
p2.getLatitudeDeg(), p2.getLongitudeDeg(),
&course1, &course2, &distance);
return ret == 0;
}
/// Geocentric routines
void
SGGeodesy::advanceRadM(const SGGeoc& geoc, double course, double distance,
SGGeoc& result)
{
result.setRadiusM(geoc.getRadiusM());
// lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
// IF (cos(lat)=0)
// lon=lon1 // endpoint a pole
// ELSE
// lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi
// ENDIF
distance *= SG_METER_TO_NM * SG_NM_TO_RAD;
double sinDistance = sin(distance);
double cosDistance = cos(distance);
double sinLat = sin(geoc.getLatitudeRad()) * cosDistance +
cos(geoc.getLatitudeRad()) * sinDistance * cos(course);
sinLat = SGMiscd::clip(sinLat, -1, 1);
result.setLatitudeRad(asin(sinLat));
double cosLat = cos(result.getLatitudeRad());
if (cosLat <= SGLimitsd::min()) {
// endpoint a pole
result.setLongitudeRad(geoc.getLongitudeRad());
} else {
double tmp = SGMiscd::clip(sin(course) * sinDistance / cosLat, -1, 1);
double lon = SGMiscd::normalizeAngle(geoc.getLongitudeRad() - asin( tmp ));
result.setLongitudeRad(lon);
}
}
double
SGGeodesy::courseRad(const SGGeoc& from, const SGGeoc& to)
{
double diffLon = to.getLongitudeRad() - from.getLongitudeRad();
double sinLatFrom = sin(from.getLatitudeRad());
double cosLatFrom = cos(from.getLatitudeRad());
double sinLatTo = sin(to.getLatitudeRad());
double cosLatTo = cos(to.getLatitudeRad());
double x = cosLatTo*sin(diffLon);
double y = cosLatFrom*sinLatTo - sinLatFrom*cosLatTo*cos(diffLon);
// guard atan2 returning NaN's
if (fabs(x) <= SGLimitsd::min() && fabs(y) <= SGLimitsd::min())
return 0;
double c = atan2(x, y);
if (c >= 0)
return SGMiscd::twopi() - c;
else
return -c;
}
double
SGGeodesy::distanceM(const SGGeoc& from, const SGGeoc& to)
{
// d = 2*asin(sqrt((sin((lat1-lat2)/2))^2 +
// cos(lat1)*cos(lat2)*(sin((lon1-lon2)/2))^2))
double cosLatFrom = cos(from.getLatitudeRad());
double cosLatTo = cos(to.getLatitudeRad());
double tmp1 = sin(0.5*(from.getLatitudeRad() - to.getLatitudeRad()));
double tmp2 = sin(0.5*(from.getLongitudeRad() - to.getLongitudeRad()));
double square = tmp1*tmp1 + cosLatFrom*cosLatTo*tmp2*tmp2;
double s = SGMiscd::min(sqrt(SGMiscd::max(square, 0)), 1);
return 2 * asin(s) * SG_RAD_TO_NM * SG_NM_TO_METER;
}

13
simgear/math/SGGeodesy.hxx
View File

@@ -45,6 +45,19 @@ public:
/// Takes a geocentric coordinate data and returns the cartesian
/// coordinates.
static void SGGeocToCart(const SGGeoc& geoc, SGVec3<double>& cart);
// Geodetic course/distance computation
static bool direct(const SGGeod& p1, double course1,
double distance, SGGeod& p2, double& course2);
static bool inverse(const SGGeod& p1, const SGGeod& p2, double& course1,
double& course2, double& distance);
// Geocentric course/distance computation
static void advanceRadM(const SGGeoc& geoc, double course, double distance,
SGGeoc& result);
static double courseRad(const SGGeoc& from, const SGGeoc& to);
static double distanceM(const SGGeoc& from, const SGGeoc& to);
};
#endif

225
simgear/math/polar3d.cxx
View File

@@ -1,225 +0,0 @@
// polar.cxx -- routines to deal with polar math and transformations
//
// Written by Curtis Olson, started June 1997.
//
// Copyright (C) 1997 Curtis L. Olson - http://www.flightgear.org/~curt
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Library General Public
// License as published by the Free Software Foundation; either
// version 2 of the License, or (at your option) any later version.
//
// This library 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
// Library 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
//
// $Id$
#ifdef HAVE_CONFIG_H
# include <simgear_config.h>
#endif
#include <math.h>
#include <simgear/constants.h>
#include "polar3d.hxx"
/**
* Calculate new lon/lat given starting lon/lat, and offset radial, and
* distance. NOTE: starting point is specifed in radians, distance is
* specified in meters (and converted internally to radians)
* ... assumes a spherical world.
* @param orig specified in polar coordinates
* @param course offset radial
* @param dist offset distance
* @return destination point in polar coordinates
*/
Point3D calc_gc_lon_lat( const Point3D& orig, double course,
double dist ) {
Point3D result;
// lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
// IF (cos(lat)=0)
// lon=lon1 // endpoint a pole
// ELSE
// lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi
// ENDIF
// printf("calc_lon_lat() offset.theta = %.2f offset.dist = %.2f\n",
// offset.theta, offset.dist);
dist *= SG_METER_TO_NM * SG_NM_TO_RAD;
result.sety( asin( sin(orig.y()) * cos(dist) +
cos(orig.y()) * sin(dist) * cos(course) ) );
if ( cos(result.y()) < SG_EPSILON ) {
result.setx( orig.x() ); // endpoint a pole
} else {
result.setx(
fmod(orig.x() - asin( sin(course) * sin(dist) /
cos(result.y()) )
+ SGD_PI, SGD_2PI) - SGD_PI );
}
return result;
}
/**
* Calculate course/dist given two spherical points.
* @param start starting point
* @param dest ending point
* @param course resulting course
* @param dist resulting distance
*/
void calc_gc_course_dist( const Point3D& start, const Point3D& dest,
double *course, double *dist )
{
if ( start == dest) {
*dist=0;
*course=0;
return;
}
// d = 2*asin(sqrt((sin((lat1-lat2)/2))^2 +
// cos(lat1)*cos(lat2)*(sin((lon1-lon2)/2))^2))
double cos_start_y = cos( start.y() );
double tmp1 = sin( (start.y() - dest.y()) * 0.5 );
double tmp2 = sin( (start.x() - dest.x()) * 0.5 );
double d = 2.0 * asin( sqrt( tmp1 * tmp1 +
cos_start_y * cos(dest.y()) * tmp2 * tmp2));
*dist = d * SG_RAD_TO_NM * SG_NM_TO_METER;
#if 1
double c1 = atan2(
cos(dest.y())*sin(dest.x()-start.x()),
cos(start.y())*sin(dest.y())-
sin(start.y())*cos(dest.y())*cos(dest.x()-start.x()));
if (c1 >= 0)
*course = SGD_2PI-c1;
else
*course = -c1;
#else
// We obtain the initial course, tc1, (at point 1) from point 1 to
// point 2 by the following. The formula fails if the initial
// point is a pole. We can special case this with:
//
// IF (cos(lat1) < EPS) // EPS a small number ~ machine precision
// IF (lat1 > 0)
// tc1= pi // starting from N pole
// ELSE
// tc1= 0 // starting from S pole
// ENDIF
// ENDIF
//
// For starting points other than the poles:
//
// IF sin(lon2-lon1)<0
// tc1=acos((sin(lat2)-sin(lat1)*cos(d))/(sin(d)*cos(lat1)))
// ELSE
// tc1=2*pi-acos((sin(lat2)-sin(lat1)*cos(d))/(sin(d)*cos(lat1)))
// ENDIF
// if ( cos(start.y()) < SG_EPSILON ) {
// doing it this way saves a transcendental call
double sin_start_y = sin( start.y() );
if ( fabs(1.0-sin_start_y) < SG_EPSILON ) {
// EPS a small number ~ machine precision
if ( start.y() > 0 ) {
*course = SGD_PI; // starting from N pole
} else {
*course = 0; // starting from S pole
}
} else {
// For starting points other than the poles:
// double tmp3 = sin(d)*cos_start_y);
// double tmp4 = sin(dest.y())-sin(start.y())*cos(d);
// double tmp5 = acos(tmp4/tmp3);
// Doing this way gaurentees that the temps are
// not stored into memory
double tmp5 = acos( (sin(dest.y()) - sin_start_y * cos(d)) /
(sin(d) * cos_start_y) );
// if ( sin( dest.x() - start.x() ) < 0 ) {
// the sin of the negative angle is just the opposite sign
// of the sin of the angle so tmp2 will have the opposite
// sign of sin( dest.x() - start.x() )
if ( tmp2 >= 0 ) {
*course = tmp5;
} else {
*course = SGD_2PI - tmp5;
}
}
#endif
}
#if 0
/**
* Calculate course/dist given two spherical points.
* @param start starting point
* @param dest ending point
* @param course resulting course
* @param dist resulting distance
*/
void calc_gc_course_dist( const Point3D& start, const Point3D& dest,
double *course, double *dist ) {
// d = 2*asin(sqrt((sin((lat1-lat2)/2))^2 +
// cos(lat1)*cos(lat2)*(sin((lon1-lon2)/2))^2))
double tmp1 = sin( (start.y() - dest.y()) / 2 );
double tmp2 = sin( (start.x() - dest.x()) / 2 );
double d = 2.0 * asin( sqrt( tmp1 * tmp1 +
cos(start.y()) * cos(dest.y()) * tmp2 * tmp2));
// We obtain the initial course, tc1, (at point 1) from point 1 to
// point 2 by the following. The formula fails if the initial
// point is a pole. We can special case this with:
//
// IF (cos(lat1) < EPS) // EPS a small number ~ machine precision
// IF (lat1 > 0)
// tc1= pi // starting from N pole
// ELSE
// tc1= 0 // starting from S pole
// ENDIF
// ENDIF
//
// For starting points other than the poles:
//
// IF sin(lon2-lon1)<0
// tc1=acos((sin(lat2)-sin(lat1)*cos(d))/(sin(d)*cos(lat1)))
// ELSE
// tc1=2*pi-acos((sin(lat2)-sin(lat1)*cos(d))/(sin(d)*cos(lat1)))
// ENDIF
double tc1;
if ( cos(start.y()) < SG_EPSILON ) {
// EPS a small number ~ machine precision
if ( start.y() > 0 ) {
tc1 = SGD_PI; // starting from N pole
} else {
tc1 = 0; // starting from S pole
}
}
// For starting points other than the poles:
double tmp3 = sin(d)*cos(start.y());
double tmp4 = sin(dest.y())-sin(start.y())*cos(d);
double tmp5 = acos(tmp4/tmp3);
if ( sin( dest.x() - start.x() ) < 0 ) {
tc1 = tmp5;
} else {
tc1 = SGD_2PI - tmp5;
}
*course = tc1;
*dist = d * SG_RAD_TO_NM * SG_NM_TO_METER;
}
#endif // 0

72
simgear/math/polar3d.hxx
View File

@@ -32,56 +32,9 @@
# error This library requires C++
#endif
#include <simgear/constants.h>
#include <simgear/math/point3d.hxx>
#include "SGMath.hxx"
/**
* Find the Altitude above the Ellipsoid (WGS84) given the Earth
* Centered Cartesian coordinate vector Distances are specified in
* meters.
* @param cp point specified in cartesian coordinates
* @return altitude above the (wgs84) earth in meters
*/
inline double sgGeodAltFromCart(const Point3D& p)
{
SGGeod geod;
SGGeodesy::SGCartToGeod(SGVec3<double>(p.x(), p.y(), p.z()), geod);
return geod.getElevationM();
}
/**
* Convert a polar coordinate to a cartesian coordinate. Lon and Lat
* must be specified in radians. The SG convention is for distances
* to be specified in meters
* @param p point specified in polar coordinates
* @return the same point in cartesian coordinates
*/
inline Point3D sgPolarToCart3d(const Point3D& p)
{
SGVec3<double> cart;
SGGeodesy::SGGeocToCart(SGGeoc::fromRadM(p.lon(), p.lat(), p.radius()), cart);
return Point3D::fromSGVec3(cart);
}
/**
* Convert a cartesian coordinate to polar coordinates (lon/lat
* specified in radians. Distances are specified in meters.
* @param cp point specified in cartesian coordinates
* @return the same point in polar coordinates
*/
inline Point3D sgCartToPolar3d(const Point3D& p)
{
SGGeoc geoc;
SGGeodesy::SGCartToGeoc(SGVec3<double>(p.x(), p.y(), p.z()), geoc);
return Point3D::fromSGGeoc(geoc);
}
/**
* Calculate new lon/lat given starting lon/lat, and offset radial, and
* distance. NOTE: starting point is specifed in radians, distance is
@@ -92,7 +45,8 @@ inline Point3D sgCartToPolar3d(const Point3D& p)
* @param dist offset distance
* @return destination point in polar coordinates
*/
Point3D calc_gc_lon_lat( const Point3D& orig, double course, double dist );
inline Point3D calc_gc_lon_lat(const Point3D& orig, double course, double dist)
{ return Point3D::fromSGGeoc(orig.toSGGeoc().advanceRadM(course, dist)); }
/**
@@ -102,20 +56,14 @@ Point3D calc_gc_lon_lat( const Point3D& orig, double course, double dist );
* @param course resulting course
* @param dist resulting distance
*/
void calc_gc_course_dist( const Point3D& start, const Point3D& dest,
double *course, double *dist );
#if 0
/**
* Calculate course/dist given two spherical points.
* @param start starting point
* @param dest ending point
* @param course resulting course
* @param dist resulting distance
*/
void calc_gc_course_dist( const Point3D& start, const Point3D& dest,
double *course, double *dist );
#endif // 0
inline void calc_gc_course_dist( const Point3D& start, const Point3D& dest,
double *course, double *dist )
{
SGGeoc gs = start.toSGGeoc();
SGGeoc gd = dest.toSGGeoc();
*course = SGGeoc::courseRad(gs, gd);
*dist = SGGeoc::distanceM(gs, gd);
}
#endif // _POLAR3D_HXX

283
simgear/math/sg_geodesy.cxx
View File

@@ -1,283 +0,0 @@
#ifdef HAVE_CONFIG_H
# include <simgear_config.h>
#endif
#include <simgear/constants.h>
#include "SGMath.hxx"
#include "sg_geodesy.hxx"
// Notes:
//
// The XYZ/cartesian coordinate system in use puts the X axis through
// zero lat/lon (off west Africa), the Z axis through the north pole,
// and the Y axis through 90 degrees longitude (in the Indian Ocean).
//
// All latitude and longitude values are in radians. Altitude is in
// meters, with zero on the WGS84 ellipsoid.
//
// The code below makes use of the notion of "squashed" space. This
// is a 2D cylindrical coordinate system where the radius from the Z
// axis is multiplied by SQUASH; the earth in this space is a perfect
// circle with a radius of POLRAD.
//
// Performance: with full optimization, a transformation from
// lat/lon/alt to XYZ and back takes 5263 CPU cycles on my 2.2GHz
// Pentium 4. About 83% of this is spent in the iterative sgCartToGeod()
// algorithm.
// These are hard numbers from the WGS84 standard. DON'T MODIFY
// unless you want to change the datum.
static const double EQURAD = 6378137;
static const double iFLATTENING = 298.257223563;
// These are derived quantities more useful to the code:
#if 0
static const double SQUASH = 1 - 1/iFLATTENING;
static const double STRETCH = 1/SQUASH;
static const double POLRAD = EQURAD * SQUASH;
#else
// High-precision versions of the above produced with an arbitrary
// precision calculator (the compiler might lose a few bits in the FPU
// operations). These are specified to 81 bits of mantissa, which is
// higher than any FPU known to me:
static const double SQUASH = 0.9966471893352525192801545;
static const double STRETCH = 1.0033640898209764189003079;
static const double POLRAD = 6356752.3142451794975639668;
#endif
////////////////////////////////////////////////////////////////////////
//
// Direct and inverse distance functions
//
// Proceedings of the 7th International Symposium on Geodetic
// Computations, 1985
//
// "The Nested Coefficient Method for Accurate Solutions of Direct and
// Inverse Geodetic Problems With Any Length"
//
// Zhang Xue-Lian
// pp 747-763
//
// modified for FlightGear to use WGS84 only -- Norman Vine
static const double GEOD_INV_PI = SGD_PI;
// s == distance
// az = azimuth
static inline double M0( double e2 ) {
//double e4 = e2*e2;
return GEOD_INV_PI*(1.0 - e2*( 1.0/4.0 + e2*( 3.0/64.0 +
e2*(5.0/256.0) )))/2.0;
}
// given, lat1, lon1, az1 and distance (s), calculate lat2, lon2
// and az2. Lat, lon, and azimuth are in degrees. distance in meters
int geo_direct_wgs_84 ( double lat1, double lon1, double az1,
double s, double *lat2, double *lon2,
double *az2 )
{
double a = EQURAD, rf = iFLATTENING;
double RADDEG = (GEOD_INV_PI)/180.0, testv = 1.0E-10;
double f = ( rf > 0.0 ? 1.0/rf : 0.0 );
double b = a*(1.0-f);
double e2 = f*(2.0-f);
double phi1 = lat1*RADDEG, lam1 = lon1*RADDEG;
double sinphi1 = sin(phi1), cosphi1 = cos(phi1);
double azm1 = az1*RADDEG;
double sinaz1 = sin(azm1), cosaz1 = cos(azm1);
if( fabs(s) < 0.01 ) { // distance < centimeter => congruency
*lat2 = lat1;
*lon2 = lon1;
*az2 = 180.0 + az1;
if( *az2 > 360.0 ) *az2 -= 360.0;
return 0;
} else if( cosphi1 ) { // non-polar origin
// u1 is reduced latitude
double tanu1 = sqrt(1.0-e2)*sinphi1/cosphi1;
double sig1 = atan2(tanu1,cosaz1);
double cosu1 = 1.0/sqrt( 1.0 + tanu1*tanu1 ), sinu1 = tanu1*cosu1;
double sinaz = cosu1*sinaz1, cos2saz = 1.0-sinaz*sinaz;
double us = cos2saz*e2/(1.0-e2);
// Terms
double ta = 1.0+us*(4096.0+us*(-768.0+us*(320.0-175.0*us)))/16384.0,
tb = us*(256.0+us*(-128.0+us*(74.0-47.0*us)))/1024.0,
tc = 0;
// FIRST ESTIMATE OF SIGMA (SIG)
double first = s/(b*ta); // !!
double sig = first;
double c2sigm, sinsig,cossig, temp,denom,rnumer, dlams, dlam;
do {
c2sigm = cos(2.0*sig1+sig);
sinsig = sin(sig); cossig = cos(sig);
temp = sig;
sig = first +
tb*sinsig*(c2sigm+tb*(cossig*(-1.0+2.0*c2sigm*c2sigm) -
tb*c2sigm*(-3.0+4.0*sinsig*sinsig)
*(-3.0+4.0*c2sigm*c2sigm)/6.0)
/4.0);
} while( fabs(sig-temp) > testv);
// LATITUDE OF POINT 2
// DENOMINATOR IN 2 PARTS (TEMP ALSO USED LATER)
temp = sinu1*sinsig-cosu1*cossig*cosaz1;
denom = (1.0-f)*sqrt(sinaz*sinaz+temp*temp);
// NUMERATOR
rnumer = sinu1*cossig+cosu1*sinsig*cosaz1;
*lat2 = atan2(rnumer,denom)/RADDEG;
// DIFFERENCE IN LONGITUDE ON AUXILARY SPHERE (DLAMS )
rnumer = sinsig*sinaz1;
denom = cosu1*cossig-sinu1*sinsig*cosaz1;
dlams = atan2(rnumer,denom);
// TERM C
tc = f*cos2saz*(4.0+f*(4.0-3.0*cos2saz))/16.0;
// DIFFERENCE IN LONGITUDE
dlam = dlams-(1.0-tc)*f*sinaz*(sig+tc*sinsig*
(c2sigm+
tc*cossig*(-1.0+2.0*
c2sigm*c2sigm)));
*lon2 = (lam1+dlam)/RADDEG;
if (*lon2 > 180.0 ) *lon2 -= 360.0;
if (*lon2 < -180.0 ) *lon2 += 360.0;
// AZIMUTH - FROM NORTH
*az2 = atan2(-sinaz,temp)/RADDEG;
if ( fabs(*az2) < testv ) *az2 = 0.0;
if( *az2 < 0.0) *az2 += 360.0;
return 0;
} else { // phi1 == 90 degrees, polar origin
double dM = a*M0(e2) - s;
double paz = ( phi1 < 0.0 ? 180.0 : 0.0 );
double zero = 0.0f;
return geo_direct_wgs_84( zero, lon1, paz, dM, lat2, lon2, az2 );
}
}
// given lat1, lon1, lat2, lon2, calculate starting and ending
// az1, az2 and distance (s). Lat, lon, and azimuth are in degrees.
// distance in meters
int geo_inverse_wgs_84( double lat1, double lon1, double lat2,
double lon2, double *az1, double *az2,
double *s )
{
double a = EQURAD, rf = iFLATTENING;
int iter=0;
double RADDEG = (GEOD_INV_PI)/180.0, testv = 1.0E-10;
double f = ( rf > 0.0 ? 1.0/rf : 0.0 );
double b = a*(1.0-f);
// double e2 = f*(2.0-f); // unused in this routine
double phi1 = lat1*RADDEG, lam1 = lon1*RADDEG;
double sinphi1 = sin(phi1), cosphi1 = cos(phi1);
double phi2 = lat2*RADDEG, lam2 = lon2*RADDEG;
double sinphi2 = sin(phi2), cosphi2 = cos(phi2);
if( (fabs(lat1-lat2) < testv &&
( fabs(lon1-lon2) < testv) || fabs(lat1-90.0) < testv ) )
{
// TWO STATIONS ARE IDENTICAL : SET DISTANCE & AZIMUTHS TO ZERO */
*az1 = 0.0; *az2 = 0.0; *s = 0.0;
return 0;
} else if( fabs(cosphi1) < testv ) {
// initial point is polar
int k = geo_inverse_wgs_84( lat2,lon2,lat1,lon1, az1,az2,s );
k = k; // avoid compiler error since return result is unused
b = *az1; *az1 = *az2; *az2 = b;
return 0;
} else if( fabs(cosphi2) < testv ) {
// terminal point is polar
double _lon1 = lon1 + 180.0f;
int k = geo_inverse_wgs_84( lat1, lon1, lat1, _lon1,
az1, az2, s );
k = k; // avoid compiler error since return result is unused
*s /= 2.0;
*az2 = *az1 + 180.0;
if( *az2 > 360.0 ) *az2 -= 360.0;
return 0;
} else if( (fabs( fabs(lon1-lon2) - 180 ) < testv) &&
(fabs(lat1+lat2) < testv) )
{
// Geodesic passes through the pole (antipodal)
double s1,s2;
geo_inverse_wgs_84( lat1,lon1, lat1,lon2, az1,az2, &s1 );
geo_inverse_wgs_84( lat2,lon2, lat1,lon2, az1,az2, &s2 );
*az2 = *az1;
*s = s1 + s2;
return 0;
} else {
// antipodal and polar points don't get here
double dlam = lam2 - lam1, dlams = dlam;
double sdlams,cdlams, sig,sinsig,cossig, sinaz,
cos2saz, c2sigm;
double tc,temp, us,rnumer,denom, ta,tb;
double cosu1,sinu1, sinu2,cosu2;
// Reduced latitudes
temp = (1.0-f)*sinphi1/cosphi1;
cosu1 = 1.0/sqrt(1.0+temp*temp);
sinu1 = temp*cosu1;
temp = (1.0-f)*sinphi2/cosphi2;
cosu2 = 1.0/sqrt(1.0+temp*temp);
sinu2 = temp*cosu2;
do {
sdlams = sin(dlams), cdlams = cos(dlams);
sinsig = sqrt(cosu2*cosu2*sdlams*sdlams+
(cosu1*sinu2-sinu1*cosu2*cdlams)*
(cosu1*sinu2-sinu1*cosu2*cdlams));
cossig = sinu1*sinu2+cosu1*cosu2*cdlams;
sig = atan2(sinsig,cossig);
sinaz = cosu1*cosu2*sdlams/sinsig;
cos2saz = 1.0-sinaz*sinaz;
c2sigm = (sinu1 == 0.0 || sinu2 == 0.0 ? cossig :
cossig-2.0*sinu1*sinu2/cos2saz);
tc = f*cos2saz*(4.0+f*(4.0-3.0*cos2saz))/16.0;
temp = dlams;
dlams = dlam+(1.0-tc)*f*sinaz*
(sig+tc*sinsig*
(c2sigm+tc*cossig*(-1.0+2.0*c2sigm*c2sigm)));
if (fabs(dlams) > GEOD_INV_PI && iter++ > 50) {
return iter;
}
} while ( fabs(temp-dlams) > testv);
us = cos2saz*(a*a-b*b)/(b*b); // !!
// BACK AZIMUTH FROM NORTH
rnumer = -(cosu1*sdlams);
denom = sinu1*cosu2-cosu1*sinu2*cdlams;
*az2 = atan2(rnumer,denom)/RADDEG;
if( fabs(*az2) < testv ) *az2 = 0.0;
if(*az2 < 0.0) *az2 += 360.0;
// FORWARD AZIMUTH FROM NORTH
rnumer = cosu2*sdlams;
denom = cosu1*sinu2-sinu1*cosu2*cdlams;
*az1 = atan2(rnumer,denom)/RADDEG;
if( fabs(*az1) < testv ) *az1 = 0.0;
if(*az1 < 0.0) *az1 += 360.0;
// Terms a & b
ta = 1.0+us*(4096.0+us*(-768.0+us*(320.0-175.0*us)))/
16384.0;
tb = us*(256.0+us*(-128.0+us*(74.0-47.0*us)))/1024.0;
// GEODETIC DISTANCE
*s = b*ta*(sig-tb*sinsig*
(c2sigm+tb*(cossig*(-1.0+2.0*c2sigm*c2sigm)-tb*
c2sigm*(-3.0+4.0*sinsig*sinsig)*
(-3.0+4.0*c2sigm*c2sigm)/6.0)/
4.0));
return 0;
}
}

57
simgear/math/sg_geodesy.hxx
View File

@@ -4,26 +4,8 @@
#include <simgear/math/point3d.hxx>
#include "SGMath.hxx"
/**
* Convert from geocentric coordinates to geodetic coordinates
* @param lat_geoc (in) Geocentric latitude, radians, + = North
* @param radius (in) C.G. radius to earth center (meters)
* @param lat_geod (out) Geodetic latitude, radians, + = North
* @param alt (out) C.G. altitude above mean sea level (meters)
* @param sea_level_r (out) radius from earth center to sea level at
* local vertical (surface normal) of C.G. (meters)
*/
inline void sgGeocToGeod(double lat_geoc, double radius,
double *lat_geod, double *alt, double *sea_level_r)
{
SGVec3<double> cart;
SGGeodesy::SGGeocToCart(SGGeoc::fromRadM(0, lat_geoc, radius), cart);
SGGeod geod;
SGGeodesy::SGCartToGeod(cart, geod);
*lat_geod = geod.getLatitudeRad();
*alt = geod.getElevationM();
*sea_level_r = SGGeodesy::SGGeodToSeaLevelRadius(geod);
}
// Compatibility header.
// Please use the SGGeodesy and SGMath functions directly.
/**
* Convert from geodetic coordinates to geocentric coordinates.
@@ -128,9 +110,17 @@ inline Point3D sgGeodToCart(const Point3D& geod)
* @param lon2 (out) degrees
* @param az2 (out) return course in degrees
*/
int geo_direct_wgs_84 ( double lat1, double lon1, double az1,
double s, double *lat2, double *lon2,
double *az2 );
inline int geo_direct_wgs_84 ( double lat1, double lon1, double az1,
double s, double *lat2, double *lon2,
double *az2 )
{
SGGeod p2;
if (!SGGeodesy::direct(SGGeod::fromDeg(lon1, lat1), az1, s, p2, *az2))
return 1;
*lat2 = p2.getLatitudeDeg();
*lon2 = p2.getLongitudeDeg();
return 0;
}
inline int geo_direct_wgs_84 ( double alt, double lat1,
double lon1, double az1,
double s, double *lat2, double *lon2,
@@ -149,12 +139,7 @@ inline int geo_direct_wgs_84 ( double alt, double lat1,
inline int geo_direct_wgs_84(const SGGeod& p1, double az1,
double s, SGGeod& p2, double *az2 )
{
double lat2, lon2;
int ret = geo_direct_wgs_84(p1.getLatitudeDeg(), p1.getLongitudeDeg(),
az1, s, &lat2, &lon2, az2);
p2.setLatitudeDeg(lat2);
p2.setLongitudeDeg(lon2);
return ret;
return !SGGeodesy::direct(p1, az1, s, p2, *az2);
}
/**
@@ -169,9 +154,13 @@ inline int geo_direct_wgs_84(const SGGeod& p1, double az1,
* @param az2 (out) end heading degrees
* @param s (out) distance meters
*/
int geo_inverse_wgs_84( double lat1, double lon1, double lat2,
double lon2, double *az1, double *az2,
double *s );
inline int geo_inverse_wgs_84( double lat1, double lon1, double lat2,
double lon2, double *az1, double *az2,
double *s )
{
return !SGGeodesy::inverse(SGGeod::fromDeg(lon1, lat1),
SGGeod::fromDeg(lon2, lat2), *az1, *az2, *s);
}
inline int geo_inverse_wgs_84( double alt, double lat1,
double lon1, double lat2,
double lon2, double *az1, double *az2,
@@ -191,9 +180,7 @@ inline int geo_inverse_wgs_84( double alt, double lat1,
inline int geo_inverse_wgs_84(const SGGeod& p1, const SGGeod& p2,
double *az1, double *az2, double *s )
{
return geo_inverse_wgs_84(p1.getLatitudeDeg(), p1.getLongitudeDeg(),
p2.getLatitudeDeg(), p2.getLongitudeDeg(),
az1, az2, s);
return !SGGeodesy::inverse(p1, p2, *az1, *az2, *s);
}
#endif // _SG_GEODESY_HXX

63
simgear/math/sg_memory.h
View File

@@ -1,63 +0,0 @@
/**
* \file sg_memory.h
* memcpy/bcopy portability declarations (not actually used by anything
* as best as I can tell.
*/
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Library General Public
// License as published by the Free Software Foundation; either
// version 2 of the License, or (at your option) any later version.
//
// This library 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
// Library 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
//
// $Id$
#ifndef _SG_MEMORY_H
#define _SG_MEMORY_H
#ifdef HAVE_CONFIG_H
# include <simgear_config.h>
#endif
#ifdef HAVE_MEMCPY
# ifdef HAVE_MEMORY_H
# include <memory.h>
# endif
# define sgmemcmp memcmp
# define sgmemcpy memcpy
# define sgmemzero(dest,len) memset(dest,0,len)
#elif defined(HAVE_BCOPY)
# define sgmemcmp bcmp
# define sgmemcpy(dest,src,n) bcopy(src,dest,n)
# define sgmemzero bzero
#else
/*
* Neither memcpy() or bcopy() available.
* Use substitutes provided be zlib.
*/
# include <zutil.h>
# define sgmemcmp zmemcmp
# define sgmemcpy zmemcpy
# define sgmemzero zmemzero
#endif
#endif // _SG_MEMORY_H

20
simgear/math/sg_types.hxx
View File

@@ -63,25 +63,5 @@ typedef vector < string > string_list;
typedef string_list::iterator string_list_iterator;
typedef string_list::const_iterator const_string_list_iterator;
/**
* Simple 2d point class where members can be accessed as x, dist, or lon
* and y, theta, or lat
*/
class point2d {
public:
union {
double x;
double dist;
double lon;
};
union {
double y;
double theta;
double lat;
};
};
#endif // _SG_TYPES_HXX