From 7832a026b123a02ac44b3e66a575b982308a3bcb Mon Sep 17 00:00:00 2001 From: Alexander Kmoch Date: Fri, 1 Nov 2024 01:55:51 +0200 Subject: [PATCH] add cmath and/or _USE_MATH_DEFINES for MSVC in DGGRID files --- src/lib/dglib/include/dglib/DgEllipsoidRF.h | 560 ++++++----- src/lib/dglib/lib/DgEllipsoidRF.cpp | 983 ++++++++++---------- src/lib/dglib/lib/DgProjGnomonicRF.cpp | 242 +++-- 3 files changed, 908 insertions(+), 877 deletions(-) diff --git a/src/lib/dglib/include/dglib/DgEllipsoidRF.h b/src/lib/dglib/include/dglib/DgEllipsoidRF.h index efb1a116..3ba8553f 100644 --- a/src/lib/dglib/include/dglib/DgEllipsoidRF.h +++ b/src/lib/dglib/include/dglib/DgEllipsoidRF.h @@ -26,19 +26,20 @@ #define DGELLIPSOIDRF_H #include +#include #include #include #include #include -//#include +// #include using namespace std; typedef struct GeoCoord { - long double lat; - long double lon; + long double lat; + long double lon; } GeoCoord; @@ -47,81 +48,85 @@ class DgPolygon; //////////////////////////////////////////////////////////////////////////////// class DgGeoCoord : public DgDVec2D { - public: - - static const DgGeoCoord undefGeoCoord; +public: + static const DgGeoCoord undefGeoCoord; - static const long double tolerance; // decimal degrees + static const long double tolerance; // decimal degrees - static long double gcDist (const DgGeoCoord& g1, const DgGeoCoord& g2, + static long double gcDist(const DgGeoCoord &g1, const DgGeoCoord &g2, bool rads = true); - static DgGeoCoord gcIntersect (const DgGeoCoord& g11, - const DgGeoCoord& g12, const DgGeoCoord& g21, - const DgGeoCoord& g22); + static DgGeoCoord gcIntersect(const DgGeoCoord &g11, const DgGeoCoord &g12, + const DgGeoCoord &g21, const DgGeoCoord &g22); - static long double geoTriArea (const DgGeoCoord& g1, const DgGeoCoord& g2, - const DgGeoCoord& g3); + static long double geoTriArea(const DgGeoCoord &g1, const DgGeoCoord &g2, + const DgGeoCoord &g3); - static long double geoPolyArea (const DgPolygon& poly, + static long double geoPolyArea(const DgPolygon &poly, const DgGeoCoord center); - DgGeoCoord (void) { } - - DgGeoCoord (const DgGeoCoord& coord) = default; - - DgGeoCoord (const DgDVec2D& coord, bool rads = true) - { if (rads) *this = coord; else *this = coord * M_PI_180; } - - DgGeoCoord (long double lon, long double lat, bool rads = true) - { if (rads) *this = DgDVec2D(lon, lat); - else *this = DgDVec2D(lon * M_PI_180, lat * M_PI_180); } - - DgGeoCoord (const GeoCoord& coord) - { *this = DgDVec2D(coord.lon, coord.lat); } - - DgGeoCoord& operator= (const DgGeoCoord& coord) - { DgDVec2D::operator=(coord); return *this; } - - DgGeoCoord& operator= (const DgDVec2D& coord) - { DgDVec2D::operator=(coord); return *this; } - - operator string (void) const - { return "(" + dgg::util::to_string(lonDegs()) + ", " + - dgg::util::to_string(latDegs()) + ")"; } - - long double lat (void) const { return y(); } // latitude in radians - long double lon (void) const { return x(); } // longitude in radians - - long double latDegs (void) const { return lat() * M_180_PI; } - long double lonDegs (void) const { return lon() * M_180_PI; } - - void setLat (long double rads) { setY(rads); } - void setLon (long double rads) { setX(rads); } - - void setLatDeg (long double degs) { setY(degs * M_PI_180); } - void setLonDeg (long double degs) { setX(degs * M_PI_180); } - - void normalize (void) - { - if (fabsl(latDegs()) >= (90.0L - tolerance)) - { - setLon(0.0L); - } - else - { - while (lonDegs() <= -180.0L) setLonDeg(lonDegs() + 360.0L); - while (lonDegs() > 180.0L) setLonDeg(lonDegs() - 360.0L); - } - } - + DgGeoCoord(void) {} + + DgGeoCoord(const DgGeoCoord &coord) = default; + + DgGeoCoord(const DgDVec2D &coord, bool rads = true) { + if (rads) + *this = coord; + else + *this = coord * M_PI_180; + } + + DgGeoCoord(long double lon, long double lat, bool rads = true) { + if (rads) + *this = DgDVec2D(lon, lat); + else + *this = DgDVec2D(lon * M_PI_180, lat * M_PI_180); + } + + DgGeoCoord(const GeoCoord &coord) { *this = DgDVec2D(coord.lon, coord.lat); } + + DgGeoCoord &operator=(const DgGeoCoord &coord) { + DgDVec2D::operator=(coord); + return *this; + } + + DgGeoCoord &operator=(const DgDVec2D &coord) { + DgDVec2D::operator=(coord); + return *this; + } + + operator string(void) const { + return "(" + dgg::util::to_string(lonDegs()) + ", " + + dgg::util::to_string(latDegs()) + ")"; + } + + long double lat(void) const { return y(); } // latitude in radians + long double lon(void) const { return x(); } // longitude in radians + + long double latDegs(void) const { return lat() * M_180_PI; } + long double lonDegs(void) const { return lon() * M_180_PI; } + + void setLat(long double rads) { setY(rads); } + void setLon(long double rads) { setX(rads); } + + void setLatDeg(long double degs) { setY(degs * M_PI_180); } + void setLonDeg(long double degs) { setX(degs * M_PI_180); } + + void normalize(void) { + if (fabsl(latDegs()) >= (90.0L - tolerance)) { + setLon(0.0L); + } else { + while (lonDegs() <= -180.0L) + setLonDeg(lonDegs() + 360.0L); + while (lonDegs() > 180.0L) + setLonDeg(lonDegs() - 360.0L); + } + } }; //////////////////////////////////////////////////////////////////////////////// -inline ostream& -operator<< (ostream& stream, const DgGeoCoord& pt) -{ - return stream << string(pt); +inline ostream &operator<<(ostream &stream, const DgGeoCoord &pt) { + return stream << string(pt); } // ostream& operator<< @@ -169,198 +174,193 @@ operator<< (ostream& stream, const DgGeoCoordData& pt) //////////////////////////////////////////////////////////////////////////////// class DgEllipsoidRF : public DgGeoDatumRF { - public: +public: + DgEllipsoidRF &operator=(const DgEllipsoidRF &ell) { + if (&ell != this) { + a_ = ell.a(); + b_ = ell.b(); + f_ = ell.f(); + e_ = ell.e(); + es_ = ell.es(); + as_ = ell.as(); + bs_ = ell.bs(); + ra_ = ell.ra(); + one_es_ = ell.one_es(); + rone_es_ = ell.rone_es(); + } - DgEllipsoidRF& operator= (const DgEllipsoidRF& ell) - { - if (&ell != this) - { - a_ = ell.a(); - b_ = ell.b(); - f_ = ell.f(); - e_ = ell.e(); - es_ = ell.es(); - as_ = ell.as(); - bs_ = ell.bs(); - ra_ = ell.ra(); - one_es_ = ell.one_es(); - rone_es_ = ell.rone_es(); - } - - return *this; - - } // DgEllipsoidRF::operator= + return *this; -/* - bool operator== (const DgEllipsoidRF& ell) const - { - return (a() == ell.a() && - b() == ell.b()); + } // DgEllipsoidRF::operator= - } // bool DgEllipsoidRF::operator== + /* + bool operator== (const DgEllipsoidRF& ell) const + { + return (a() == ell.a() && + b() == ell.b()); - bool operator!= (const DgEllipsoidRF& ell) const - { - return !operator==(ell); + } // bool DgEllipsoidRF::operator== - } // bool DgEllipsoidRF::operator!= -*/ + bool operator!= (const DgEllipsoidRF& ell) const + { + return !operator==(ell); - // get methods + } // bool DgEllipsoidRF::operator!= + */ - long double a (void) const { return a_; } + // get methods - long double b (void) const { return b_; } + long double a(void) const { return a_; } - long double f (void) const { return f_; } + long double b(void) const { return b_; } - long double e (void) const { return e_; } + long double f(void) const { return f_; } - long double es (void) const { return es_; } + long double e(void) const { return e_; } - long double as (void) const { return as_; } + long double es(void) const { return es_; } - long double bs (void) const { return bs_; } + long double as(void) const { return as_; } - long double ra (void) const { return ra_; } + long double bs(void) const { return bs_; } - long double one_es (void) const { return one_es_; } + long double ra(void) const { return ra_; } - long double rone_es (void) const { return rone_es_; } + long double one_es(void) const { return one_es_; } - // set methods + long double rone_es(void) const { return rone_es_; } - void setA (long double aIn) { a_ = aIn; } + // set methods - void setB (long double bIn) { b_ = bIn; } + void setA(long double aIn) { a_ = aIn; } - void setF (long double fIn) { f_ = fIn; } + void setB(long double bIn) { b_ = bIn; } - void setE (long double eIn) { e_ = eIn; } + void setF(long double fIn) { f_ = fIn; } - void setEs (long double esIn) { es_ = esIn; } + void setE(long double eIn) { e_ = eIn; } - void setAs (long double asIn) { as_ = asIn; } + void setEs(long double esIn) { es_ = esIn; } - void setBs (long double bsIn) { bs_ = bsIn; } + void setAs(long double asIn) { as_ = asIn; } - void setRa (long double raIn) { ra_ = raIn; } + void setBs(long double bsIn) { bs_ = bsIn; } - void setOne_es (long double one_esIn) { one_es_ = one_esIn; } + void setRa(long double raIn) { ra_ = raIn; } - void setRone_es (long double rone_esIn) { rone_es_ = rone_esIn; } + void setOne_es(long double one_esIn) { one_es_ = one_esIn; } - virtual const DgGeoCoord& undefAddress (void) const - { return DgGeoCoord::undefGeoCoord; } + void setRone_es(long double rone_esIn) { rone_es_ = rone_esIn; } - // the following routines are "back-doors" included for speed; - // use with care! + virtual const DgGeoCoord &undefAddress(void) const { + return DgGeoCoord::undefGeoCoord; + } - virtual DgAddressBase* vecAddress (const DgDVec2D& v) const - { return new DgAddress(DgGeoCoord(v, false)); } + // the following routines are "back-doors" included for speed; + // use with care! - virtual DgLocation* vecLocation (const DgDVec2D& v) const - { return makeLocation(DgGeoCoord(v, false)); } + virtual DgAddressBase *vecAddress(const DgDVec2D &v) const { + return new DgAddress(DgGeoCoord(v, false)); + } - virtual DgDVec2D getVecAddress (const DgAddressBase& add) const - { - const DgGeoCoord& g = - static_cast< const DgAddress &>(add).address(); - return DgDVec2D(g.lonDegs(), g.latDegs()); - } + virtual DgLocation *vecLocation(const DgDVec2D &v) const { + return makeLocation(DgGeoCoord(v, false)); + } - virtual DgDVec2D getVecLocation (const DgLocation& loc) const - { return DgDVec2D(getAddress(loc)->lonDegs(), - getAddress(loc)->latDegs()); } + virtual DgDVec2D getVecAddress(const DgAddressBase &add) const { + const DgGeoCoord &g = + static_cast &>(add).address(); + return DgDVec2D(g.lonDegs(), g.latDegs()); + } - // distance in km; not currently defined (except for DgGeoSphRF) + virtual DgDVec2D getVecLocation(const DgLocation &loc) const { + return DgDVec2D(getAddress(loc)->lonDegs(), getAddress(loc)->latDegs()); + } - virtual long double dist (const DgGeoCoord&, const DgGeoCoord&) const - { return 1.0L; } + // distance in km; not currently defined (except for DgGeoSphRF) - virtual string add2str (const DgGeoCoord& add) const - { return string(add); } + virtual long double dist(const DgGeoCoord &, const DgGeoCoord &) const { + return 1.0L; + } - virtual string add2str (const DgGeoCoord& add, char delimiter) const - { return dgg::util::to_string(add.lonDegs(), formatStr()) + delimiter + - dgg::util::to_string(add.latDegs(), formatStr()); } + virtual string add2str(const DgGeoCoord &add) const { return string(add); } - virtual const char* str2add (DgGeoCoord* add, const char* str, - char delimiter) const; + virtual string add2str(const DgGeoCoord &add, char delimiter) const { + return dgg::util::to_string(add.lonDegs(), formatStr()) + delimiter + + dgg::util::to_string(add.latDegs(), formatStr()); + } - virtual string dist2str (const long double& dist) const - { return dgg::util::to_string(dist, formatStr()); } + virtual const char *str2add(DgGeoCoord *add, const char *str, + char delimiter) const; - virtual long double dist2dbl (const long double& dist) const - { return dist; } + virtual string dist2str(const long double &dist) const { + return dgg::util::to_string(dist, formatStr()); + } - virtual unsigned long long int dist2int (const long double& dist) const - { return static_cast(dist); } + virtual long double dist2dbl(const long double &dist) const { return dist; } - protected: + virtual unsigned long long int dist2int(const long double &dist) const { + return static_cast(dist); + } - DgEllipsoidRF (DgRFNetwork& networkIn, const string& nameIn) - : DgGeoDatumRF (networkIn, nameIn), - a_ (-1.0L), b_ (-1.0L), f_ (-1.0L), e_ (-1.0L), es_ (-1.0L), - as_ (-1.0L), bs_ (-1.0L), ra_ (-1.0L), one_es_ (-1.0L), - rone_es_ (-1.0L) { } +protected: + DgEllipsoidRF(DgRFNetwork &networkIn, const string &nameIn) + : DgGeoDatumRF(networkIn, nameIn), a_(-1.0L), + b_(-1.0L), f_(-1.0L), e_(-1.0L), es_(-1.0L), as_(-1.0L), bs_(-1.0L), + ra_(-1.0L), one_es_(-1.0L), rone_es_(-1.0L) {} - DgEllipsoidRF (DgRFNetwork& networkIn, const string& nameIn, - const DgEllipsoidRF& ell) - : DgGeoDatumRF (networkIn, nameIn), - a_ (ell.a()), b_ (ell.b()), f_ (ell.f()), - e_ (ell.e()), es_ (ell.es()), as_ (ell.as()), bs_ (ell.bs()), - ra_ (ell.ra()), one_es_ (ell.one_es()), - rone_es_ (ell.rone_es()) - { undefLoc_ = makeLocation(undefAddress()); } + DgEllipsoidRF(DgRFNetwork &networkIn, const string &nameIn, + const DgEllipsoidRF &ell) + : DgGeoDatumRF(networkIn, nameIn), a_(ell.a()), + b_(ell.b()), f_(ell.f()), e_(ell.e()), es_(ell.es()), as_(ell.as()), + bs_(ell.bs()), ra_(ell.ra()), one_es_(ell.one_es()), + rone_es_(ell.rone_es()) { + undefLoc_ = makeLocation(undefAddress()); + } - DgEllipsoidRF (DgRFNetwork& networkIn, const string& nameIn, - long double aIn, long double bIn) - : DgGeoDatumRF (networkIn, nameIn), - a_ (aIn), b_ (bIn) - { - setUndefLoc(makeLocation(undefAddress())); + DgEllipsoidRF(DgRFNetwork &networkIn, const string &nameIn, long double aIn, + long double bIn) + : DgGeoDatumRF(networkIn, nameIn), a_(aIn), + b_(bIn) { + setUndefLoc(makeLocation(undefAddress())); - f_ = (a() - b()) / a(); - as_ = a() * a(); - bs_ = b() * b(); - es_ = 1.0L - bs() / as(); - e_ = sqrtl(es()); - ra_ = 1.0L / a(); - one_es_ = 1.0L - es(); - rone_es_ = 1.0L / one_es(); - } + f_ = (a() - b()) / a(); + as_ = a() * a(); + bs_ = b() * b(); + es_ = 1.0L - bs() / as(); + e_ = sqrtl(es()); + ra_ = 1.0L / a(); + one_es_ = 1.0L - es(); + rone_es_ = 1.0L / one_es(); + } - private: - - long double a_; // semi-major axis in meters +private: + long double a_; // semi-major axis in meters - long double b_; // semi-minor axis in meters + long double b_; // semi-minor axis in meters - long double f_; // flattening = (a - b) / a + long double f_; // flattening = (a - b) / a - long double e_; // eccentricity + long double e_; // eccentricity - long double es_; // e squared + long double es_; // e squared - long double as_; // a squared + long double as_; // a squared - long double bs_; // b squared + long double bs_; // b squared - long double ra_; // 1 / a + long double ra_; // 1 / a - long double one_es_; // 1 / es + long double one_es_; // 1 / es - long double rone_es_; // 1 / one_es + long double rone_es_; // 1 / one_es }; // class DgEllipsoidRF //////////////////////////////////////////////////////////////////////////////// -inline ostream& -operator<< (ostream& stream, const DgEllipsoidRF& ell) -{ - return stream << "a: " << ell.a() << " " << "b: " << ell.b() << " " - << "f: " << ell.f() << " " << "e: " << ell.e(); +inline ostream &operator<<(ostream &stream, const DgEllipsoidRF &ell) { + return stream << "a: " << ell.a() << " " << "b: " << ell.b() << " " + << "f: " << ell.f() << " " << "e: " << ell.e(); } // ostream& operator<< @@ -378,8 +378,8 @@ operator<< (ostream& stream, const DgEllipsoidRF& ell) /**** macros ****/ -#define UNDEFVAL LDBL_MAX -#define UNDEFINT INT_MAX +#define UNDEFVAL LDBL_MAX +#define UNDEFINT INT_MAX #define PRECISION 0.0000000000005L /* precise 12 digits after dot */ @@ -387,45 +387,41 @@ operator<< (ostream& stream, const DgEllipsoidRF& ell) typedef struct Vec2D { - long double x; - long double y; + long double x; + long double y; - Vec2D() - : x(0.0L), y(0.0L) - {} + Vec2D() : x(0.0L), y(0.0L) {} } Vec2D; typedef struct Vec3D { - long double x; - long double y; - long double z; + long double x; + long double y; + long double z; - Vec3D() - : x(0.0L), y(0.0L), z(0.0L) - {} + Vec3D() : x(0.0L), y(0.0L), z(0.0L) {} } Vec3D; typedef struct SphTri { - long long int code; - GeoCoord verts[3]; /* vertices in degrees */ - long double edges[3]; /* edges opposite to verts in central angle degrees */ - long double angles[3]; /* interior angles in degrees corresponding to verts */ - long double area; /* area in km */ - long double perimeter; /* perimeter in km */ - long double compactness; + long long int code; + GeoCoord verts[3]; /* vertices in degrees */ + long double edges[3]; /* edges opposite to verts in central angle degrees */ + long double angles[3]; /* interior angles in degrees corresponding to verts */ + long double area; /* area in km */ + long double perimeter; /* perimeter in km */ + long double compactness; } SphTri; typedef struct PlaneTri { - long long int code; - long long int direction; /* 0: up, 1: down */ - Vec2D points[3]; /* points in plane triangle in km */ - Vec2D cenpoint; /* central point of plane triangle in km */ + long long int code; + long long int direction; /* 0: up, 1: down */ + Vec2D points[3]; /* points in plane triangle in km */ + Vec2D cenpoint; /* central point of plane triangle in km */ } PlaneTri; @@ -438,11 +434,11 @@ typedef struct SCtri { long double area; long double perimeter; long double sclat[3]; - /* latitude from pole of small circle v1v2 in radius */ -/* long double darea[3]; - area surounded by great circle and small circle of two points. - positive means small circle is more convex than great circle - from the center of the triangle. */ + /* latitude from pole of small circle v1v2 in radius */ + /* long double darea[3]; + area surounded by great circle and small circle of two points. + positive means small circle is more convex than great circle + from the center of the triangle. */ } SCtri; @@ -450,11 +446,11 @@ typedef struct SCtri { typedef struct PreCompGeo { - GeoCoord pt; - long double sinLat; - long double sinLon; - long double cosLat; - long double cosLon; + GeoCoord pt; + long double sinLat; + long double sinLon; + long double cosLat; + long double cosLon; } PreCompGeo; @@ -463,45 +459,45 @@ typedef struct PreCompGeo { typedef struct PreCompInTri { - long double p0x0; - long double p0y0; - long double p0z0; + long double p0x0; + long double p0y0; + long double p0z0; - long double p0x1; - long double p0y1; - long double p0z1; + long double p0x1; + long double p0y1; + long double p0z1; - long double p0x2; - long double p0y2; - long double p0z2; + long double p0x2; + long double p0y2; + long double p0z2; - long double t00; - long double t01; - long double t02; + long double t00; + long double t01; + long double t02; } PreCompInTri; typedef struct SphIcosa { - /* sufficient icosa placement data */ + /* sufficient icosa placement data */ - GeoCoord pt; - long double azimuth; + GeoCoord pt; + long double azimuth; - /* pre-calculated projection values */ + /* pre-calculated projection values */ - GeoCoord icoverts[12]; - GeoCoord icotri[20][3]; - PreCompGeo triCen[20]; - PreCompInTri ptin[20]; - long double dazh[20]; + GeoCoord icoverts[12]; + GeoCoord icotri[20][3]; + PreCompGeo triCen[20]; + PreCompInTri ptin[20]; + long double dazh[20]; } SphIcosa; typedef struct IcosaGridPt { - Vec2D pt; /* the point itself in triangle coordinates */ - int triangle; /* index (0 - 19) of the triangle */ + Vec2D pt; /* the point itself in triangle coordinates */ + int triangle; /* index (0 - 19) of the triangle */ } IcosaGridPt; @@ -509,64 +505,64 @@ typedef struct IcosaGridPt { /* printing functions */ -void printVal (long double val); /* print val or 'UNDEFVAL' to stdout */ -void printGeoCoord (const GeoCoord& p); /* print p to stdout */ -void printVec2D (const Vec2D& p); /* print p to stdout */ -void printVec3D (const Vec3D& p); /* print p to stdout */ -void printSphTri (const SphTri& tri); /* print tri to stdout */ -void printPlaneTri (const PlaneTri& tri);/* print plane tri to stdout */ +void printVal(long double val); /* print val or 'UNDEFVAL' to stdout */ +void printGeoCoord(const GeoCoord &p); /* print p to stdout */ +void printVec2D(const Vec2D &p); /* print p to stdout */ +void printVec3D(const Vec3D &p); /* print p to stdout */ +void printSphTri(const SphTri &tri); /* print tri to stdout */ +void printPlaneTri(const PlaneTri &tri); /* print plane tri to stdout */ /* Vec3D functions */ -Vec3D vecAdd (const Vec3D& A, const Vec3D& B); /* return A + B */ -Vec3D vecSub (const Vec3D& A, const Vec3D& B); /* return A - B */ -Vec3D vecCross (const Vec3D& A, const Vec3D& B); /* return A x B */ -long double vecDot (const Vec3D& A, const Vec3D& B); /* return A . B */ -long double vecMag (const Vec3D& V); /* return |V| */ -Vec3D vecNormalize (const Vec3D& V); /* return V as a unit vector */ +Vec3D vecAdd(const Vec3D &A, const Vec3D &B); /* return A + B */ +Vec3D vecSub(const Vec3D &A, const Vec3D &B); /* return A - B */ +Vec3D vecCross(const Vec3D &A, const Vec3D &B); /* return A x B */ +long double vecDot(const Vec3D &A, const Vec3D &B); /* return A . B */ +long double vecMag(const Vec3D &V); /* return |V| */ +Vec3D vecNormalize(const Vec3D &V); /* return V as a unit vector */ /* spherical triangle functions */ -void sphTriInit (SphTri* tri); /* initialize with UNDEFVAL values */ +void sphTriInit(SphTri *tri); /* initialize with UNDEFVAL values */ -void planeTriInit (PlaneTri* tri); /* initialize with UNDEFVAL values */ +void planeTriInit(PlaneTri *tri); /* initialize with UNDEFVAL values */ /* solve for three sides,three angles and area */ -void sphTriSolve (SphTri* tri); +void sphTriSolve(SphTri *tri); /* calculate the center point of a sphere triangle */ GeoCoord sphTricenpoint(GeoCoord sp[3]); /* sphere functions */ -GeoCoord xyzll(const Vec3D& v); /* transform a point from xyz to latlon */ +GeoCoord xyzll(const Vec3D &v); /* transform a point from xyz to latlon */ -Vec3D llxyz(const GeoCoord& sv); /* transform a point from latlon to xyz */ +Vec3D llxyz(const GeoCoord &sv); /* transform a point from latlon to xyz */ /* calculate the chord distance */ -long double chorddist(const GeoCoord& ll1,const GeoCoord& ll2); +long double chorddist(const GeoCoord &ll1, const GeoCoord &ll2); /* calculate the sphere distance */ -long double spheredist(const GeoCoord& ll1,const GeoCoord& ll2); +long double spheredist(const GeoCoord &ll1, const GeoCoord &ll2); /* decide if pt is in sphere triangle */ -int ptinsphTri(const GeoCoord& pt, GeoCoord sTri[3]); +int ptinsphTri(const GeoCoord &pt, GeoCoord sTri[3]); /* return the middle point of the great circle */ -GeoCoord GCmidpoint(const GeoCoord& pp1, const GeoCoord& pp2); +GeoCoord GCmidpoint(const GeoCoord &pp1, const GeoCoord &pp2); /* return intersect point of two great circle */ -GeoCoord GCintersect(const GeoCoord& sv11, const GeoCoord& sv12, - const GeoCoord& sv21, const GeoCoord& sv22, int sign); +GeoCoord GCintersect(const GeoCoord &sv11, const GeoCoord &sv12, + const GeoCoord &sv21, const GeoCoord &sv22, int sign); /* return latitude of a point with known longtitude on great circle */ -long double GCptlat(long double lon, const GeoCoord& sv1, const GeoCoord& sv2); +long double GCptlat(long double lon, const GeoCoord &sv1, const GeoCoord &sv2); /* return azimuth from pt1 to pt2 */ -long double Azimuth(const GeoCoord& pt1, const GeoCoord& pt2); +long double Azimuth(const GeoCoord &pt1, const GeoCoord &pt2); /* return Geopoint that have great-circle-distance d along azimuth az from pt */ -GeoCoord GCdaz(const GeoCoord& pt, long double d, long double az); +GeoCoord GCdaz(const GeoCoord &pt, long double d, long double az); /******************************************************************************/ diff --git a/src/lib/dglib/lib/DgEllipsoidRF.cpp b/src/lib/dglib/lib/DgEllipsoidRF.cpp index 0c016ee1..21360de5 100644 --- a/src/lib/dglib/lib/DgEllipsoidRF.cpp +++ b/src/lib/dglib/lib/DgEllipsoidRF.cpp @@ -22,6 +22,7 @@ // //////////////////////////////////////////////////////////////////////////////// +#include #include #include #include @@ -30,9 +31,8 @@ const DgGeoCoord DgGeoCoord::undefGeoCoord(LDBL_MAX, LDBL_MAX); const long double DgGeoCoord::tolerance = 0.0000000005L; //////////////////////////////////////////////////////////////////////////////// -long double -DgGeoCoord::gcDist (const DgGeoCoord& g1, const DgGeoCoord& g2, - bool rads) +long double DgGeoCoord::gcDist(const DgGeoCoord &g1, const DgGeoCoord &g2, + bool rads) /* return great circle distance in radians between g1 and g2. @@ -41,43 +41,48 @@ DgGeoCoord::gcDist (const DgGeoCoord& g1, const DgGeoCoord& g2, Kevin Sahr, 3/10/99 */ { - // use spherical triangle with g1 at A, g2 at B, and north pole at C + // use spherical triangle with g1 at A, g2 at B, and north pole at C - long double bigC = fabsl(g2.lon() - g1.lon()); + long double bigC = fabsl(g2.lon() - g1.lon()); - if (bigC > M_PI) // assume we want the complement - { - // note that in this case they can't both be negative + if (bigC > M_PI) // assume we want the complement + { + // note that in this case they can't both be negative - long double lon1 = g1.lon(); - if (lon1 < 0.0L) lon1 += 2.0L * M_PI; - long double lon2 = g2.lon(); - if (lon2 < 0.0L) lon2 += 2.0L * M_PI; + long double lon1 = g1.lon(); + if (lon1 < 0.0L) + lon1 += 2.0L * M_PI; + long double lon2 = g2.lon(); + if (lon2 < 0.0L) + lon2 += 2.0L * M_PI; - bigC = fabsl(lon2 - lon1); - } + bigC = fabsl(lon2 - lon1); + } - long double b = M_PI_2 - g1.lat(); - long double a = M_PI_2 - g2.lat(); + long double b = M_PI_2 - g1.lat(); + long double a = M_PI_2 - g2.lat(); - // use law of cosines to find c + // use law of cosines to find c - long double cosc = cosl(a) * cosl(b) + sinl(a) * sinl(b) * cosl(bigC); - if (cosc > 1.0L) cosc = 1.0L; - if (cosc < -1.0L) cosc = -1.0L; + long double cosc = cosl(a) * cosl(b) + sinl(a) * sinl(b) * cosl(bigC); + if (cosc > 1.0L) + cosc = 1.0L; + if (cosc < -1.0L) + cosc = -1.0L; - long double retVal = acosl(cosc); + long double retVal = acosl(cosc); - if (!rads) retVal *= M_180_PI; + if (!rads) + retVal *= M_180_PI; - return retVal; + return retVal; } // long double DgGeoCoord::gcDist //////////////////////////////////////////////////////////////////////////////// -DgGeoCoord -DgGeoCoord::gcIntersect (const DgGeoCoord& g11, const DgGeoCoord& g12, - const DgGeoCoord& g21, const DgGeoCoord& /* g22 */) +DgGeoCoord DgGeoCoord::gcIntersect(const DgGeoCoord &g11, const DgGeoCoord &g12, + const DgGeoCoord &g21, + const DgGeoCoord & /* g22 */) /* Return point of intersection of the two great circle arc segments g11-g12 and g21-g22. @@ -85,22 +90,24 @@ DgGeoCoord::gcIntersect (const DgGeoCoord& g11, const DgGeoCoord& g12, Works by calling Lian Song's routine GCintersect. */ { - GeoCoord sv11, sv12, sv21; - sv11.lon = g11.lon(); sv11.lat = g11.lat(); - sv12.lon = g12.lon(); sv12.lat = g12.lat(); - sv21.lon = g21.lon(); sv21.lat = g21.lat(); - //sv22.lon = g22.lon(); sv22.lat = g22.lat(); + GeoCoord sv11, sv12, sv21; + sv11.lon = g11.lon(); + sv11.lat = g11.lat(); + sv12.lon = g12.lon(); + sv12.lat = g12.lat(); + sv21.lon = g21.lon(); + sv21.lat = g21.lat(); + // sv22.lon = g22.lon(); sv22.lat = g22.lat(); - GeoCoord ans = GCintersect(sv11, sv12, sv21, sv11, 1); + GeoCoord ans = GCintersect(sv11, sv12, sv21, sv11, 1); - return DgGeoCoord(ans.lon, ans.lat); + return DgGeoCoord(ans.lon, ans.lat); } // DgGeoCoord DgGeoCoord::gcIntersect //////////////////////////////////////////////////////////////////////////////// -long double -DgGeoCoord::geoTriArea (const DgGeoCoord& g1, const DgGeoCoord& g2, - const DgGeoCoord& g3) +long double DgGeoCoord::geoTriArea(const DgGeoCoord &g1, const DgGeoCoord &g2, + const DgGeoCoord &g3) /* return area in radians. @@ -109,47 +116,47 @@ DgGeoCoord::geoTriArea (const DgGeoCoord& g1, const DgGeoCoord& g2, Kevin Sahr, 3/10/99 */ { - // determine the edges + // determine the edges - long double a = DgGeoCoord::gcDist(g2, g3); - long double b = DgGeoCoord::gcDist(g1, g3); - long double c = DgGeoCoord::gcDist(g1, g2); + long double a = DgGeoCoord::gcDist(g2, g3); + long double b = DgGeoCoord::gcDist(g1, g3); + long double c = DgGeoCoord::gcDist(g1, g2); - // determine the angles using half-angle formulas + // determine the angles using half-angle formulas - long double s = (a + b + c) / 2.0L; + long double s = (a + b + c) / 2.0L; - long double sinsa = sinl(s - a); - long double sinsb = sinl(s - b); - long double sinsc = sinl(s - c); + long double sinsa = sinl(s - a); + long double sinsb = sinl(s - b); + long double sinsc = sinl(s - c); - long double k = sqrtl(sinsa * sinsb * sinsc / sinl(s)); + long double k = sqrtl(sinsa * sinsb * sinsc / sinl(s)); - long double bigA = 2.0L * atanl(k / sinsa); - long double bigB = 2.0L * atanl(k / sinsb); - long double bigC = 2.0L * atanl(k / sinsc); + long double bigA = 2.0L * atanl(k / sinsa); + long double bigB = 2.0L * atanl(k / sinsb); + long double bigC = 2.0L * atanl(k / sinsc); - long double E = bigA + bigB + bigC - M_PI; + long double E = bigA + bigB + bigC - M_PI; -/* - cout << "geoTriArea: " << g1 << " " << g2 << " " << g3 << endl - << " a: " << a - << " b: " << b - << " c: " << c << endl - << " A: " << bigA - << " B: " << bigB - << " C: " << bigC << endl - << " E: " << E << " " - << E / (4.0L * M_PI) << endl; -*/ + /* + cout << "geoTriArea: " << g1 << " " << g2 << " " << g3 << endl + << " a: " << a + << " b: " << b + << " c: " << c << endl + << " A: " << bigA + << " B: " << bigB + << " C: " << bigC << endl + << " E: " << E << " " + << E / (4.0L * M_PI) << endl; + */ - return E; + return E; } // long double DgGeoCoord::geoTriArea //////////////////////////////////////////////////////////////////////////////// -long double -DgGeoCoord::geoPolyArea (const DgPolygon& poly, const DgGeoCoord center) +long double DgGeoCoord::geoPolyArea(const DgPolygon &poly, + const DgGeoCoord center) // // returns area of spherical polygon in radians; // @@ -159,38 +166,37 @@ DgGeoCoord::geoPolyArea (const DgPolygon& poly, const DgGeoCoord center) // that center is an interior point. // { - long double totArea = 0.0L; + long double totArea = 0.0L; - const DgGeoSphRF* geoRF = dynamic_cast(&poly.rf()); - if (geoRF == 0) - report("DgGeoCoord::geoPolyArea() non-geo polygon", DgBase::Fatal); + const DgGeoSphRF *geoRF = dynamic_cast(&poly.rf()); + if (geoRF == 0) + report("DgGeoCoord::geoPolyArea() non-geo polygon", DgBase::Fatal); - // do each sub-triangle - for (int i = 0; i < poly.size(); i++) - { - const DgGeoCoord& v1 = *geoRF->getAddress(poly[i]); - const DgGeoCoord& v2 = *geoRF->getAddress(poly[(i + 1) % poly.size()]); + // do each sub-triangle + for (int i = 0; i < poly.size(); i++) { + const DgGeoCoord &v1 = *geoRF->getAddress(poly[i]); + const DgGeoCoord &v2 = *geoRF->getAddress(poly[(i + 1) % poly.size()]); - totArea += DgGeoCoord::geoTriArea(center, v1, v2); - } + totArea += DgGeoCoord::geoTriArea(center, v1, v2); + } - return totArea; + return totArea; } // long double DgGeoCoord::geoPolyArea //////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// -const char* -DgEllipsoidRF::str2add (DgGeoCoord* add, const char* str, char delimiter) const -{ - if (!add) add = new DgGeoCoord(); +const char *DgEllipsoidRF::str2add(DgGeoCoord *add, const char *str, + char delimiter) const { + if (!add) + add = new DgGeoCoord(); - DgDVec2D vec; - const char* tmp = vec.fromString(str, delimiter); + DgDVec2D vec; + const char *tmp = vec.fromString(str, delimiter); - *add = DgGeoCoord(vec.x(), vec.y(), false); + *add = DgGeoCoord(vec.x(), vec.y(), false); - return tmp; + return tmp; } // const char* DgEllipsoidRF::str2add @@ -213,252 +219,258 @@ DgEllipsoidRF::str2add (DgGeoCoord* add, const char* str, char delimiter) const #include /******************************************************************************/ -void sphTriInit (SphTri* tri) +void sphTriInit(SphTri *tri) /* Initialize tri with UNDEFVAL values. */ { - int i; + int i; - for (i = 0; i < 3; i++) - { - tri->verts[i].lon = tri->verts[i].lat = UNDEFVAL; - tri->edges[i] = tri->angles[i] = UNDEFVAL; - } - tri->code = UNDEFINT; - tri->area = tri->perimeter = UNDEFVAL; - tri->compactness = UNDEFVAL; + for (i = 0; i < 3; i++) { + tri->verts[i].lon = tri->verts[i].lat = UNDEFVAL; + tri->edges[i] = tri->angles[i] = UNDEFVAL; + } + tri->code = UNDEFINT; + tri->area = tri->perimeter = UNDEFVAL; + tri->compactness = UNDEFVAL; } /* void sphTriInit */ /******************************************************************************/ -void planeTriInit (PlaneTri* tri) +void planeTriInit(PlaneTri *tri) /* Initialize tri with UNDEFVAL values. */ { - int i; + int i; - for (i = 0; i < 3; i++) - tri->points[i].x = tri->points[i].y = UNDEFVAL; - tri->code = tri->direction = UNDEFINT; - tri->cenpoint.x = tri->cenpoint.y = UNDEFVAL; + for (i = 0; i < 3; i++) + tri->points[i].x = tri->points[i].y = UNDEFVAL; + tri->code = tri->direction = UNDEFINT; + tri->cenpoint.x = tri->cenpoint.y = UNDEFVAL; } /* void planeTriInit */ /******************************************************************************/ -void printInt (long long int val) +void printInt(long long int val) /* Print val or 'UNDEFVAL' to stdout as appropriate. */ { - if (val == UNDEFVAL) - dgcout << "UNDEFVAL"; - else - dgcout << val; + if (val == UNDEFVAL) + dgcout << "UNDEFVAL"; + else + dgcout << val; } /* void printInt */ /******************************************************************************/ -void printVal (long double val) +void printVal(long double val) /* Print val or 'UNDEFVAL' to stdout as appropriate. */ { - if (val == UNDEFVAL) - dgcout << "UNDEFVAL"; - else - dgcout << val; + if (val == UNDEFVAL) + dgcout << "UNDEFVAL"; + else + dgcout << val; } /* void printVal */ /******************************************************************************/ -void printGeoCoord (const GeoCoord& p) +void printGeoCoord(const GeoCoord &p) /* Print p to stdout. */ { - dgcout << "(" << p.lon * M_180_PI << ", " << p.lat * M_180_PI << ")"; + dgcout << "(" << p.lon * M_180_PI << ", " << p.lat * M_180_PI << ")"; } /* void printGeoCoord */ /******************************************************************************/ -void printVec2D (const Vec2D& p) +void printVec2D(const Vec2D &p) /* Print p to stdout. */ { - dgcout << "(" << p.x << ", " << p.y << ")"; + dgcout << "(" << p.x << ", " << p.y << ")"; } /* void printVec2D */ /******************************************************************************/ -void printVec3D (const Vec3D& p) +void printVec3D(const Vec3D &p) /* Print p to stdout. */ { - dgcout << "(" << p.x << ", " << p.y << ", " << p.z << ")"; + dgcout << "(" << p.x << ", " << p.y << ", " << p.z << ")"; } /* void printVec3D */ /******************************************************************************/ -void printSphTri (const SphTri& tri) +void printSphTri(const SphTri &tri) /* Print tri to stdout. */ { - int i; - - dgcout << "{\n code: "; - printInt(tri.code); - dgcout << "\n vertices: "; - for (i = 0; i < 3; i++) - { - dgcout << " "; - printGeoCoord(tri.verts[i]); - } - dgcout << "\n"; - - dgcout << " A: "; printVal(tri.edges[0]); - dgcout << " B: "; printVal(tri.edges[1]); - dgcout << " C: "; printVal(tri.edges[2]); - dgcout << "\n"; - - dgcout << " a: "; printVal(tri.angles[0]* M_180_PI); - dgcout << " b: "; printVal(tri.angles[1]* M_180_PI); - dgcout << " c: "; printVal(tri.angles[2]* M_180_PI); - dgcout << "\n"; - - dgcout << " area: "; printVal(tri.area); - dgcout << " perimeter: "; printVal(tri.perimeter); - dgcout << " compactness: "; printVal(tri.compactness); - dgcout << "\n"; + int i; - dgcout << "}\n"; + dgcout << "{\n code: "; + printInt(tri.code); + dgcout << "\n vertices: "; + for (i = 0; i < 3; i++) { + dgcout << " "; + printGeoCoord(tri.verts[i]); + } + dgcout << "\n"; + + dgcout << " A: "; + printVal(tri.edges[0]); + dgcout << " B: "; + printVal(tri.edges[1]); + dgcout << " C: "; + printVal(tri.edges[2]); + dgcout << "\n"; + + dgcout << " a: "; + printVal(tri.angles[0] * M_180_PI); + dgcout << " b: "; + printVal(tri.angles[1] * M_180_PI); + dgcout << " c: "; + printVal(tri.angles[2] * M_180_PI); + dgcout << "\n"; + + dgcout << " area: "; + printVal(tri.area); + dgcout << " perimeter: "; + printVal(tri.perimeter); + dgcout << " compactness: "; + printVal(tri.compactness); + dgcout << "\n"; + + dgcout << "}\n"; } /* void printSphTri */ /******************************************************************************/ -void printPlaneTri (const PlaneTri& tri) +void printPlaneTri(const PlaneTri &tri) /* Print tri to stdout. */ { - int i; + int i; - dgcout << "{\n code: "; printInt(tri.code); - dgcout << "\n vertices: "; - for (i = 0; i < 3; i++) - { - dgcout << " "; - printVec2D(tri.points[i]); - } - dgcout << "\n"; - printVec2D(tri.cenpoint); - dgcout << "\n"; + dgcout << "{\n code: "; + printInt(tri.code); + dgcout << "\n vertices: "; + for (i = 0; i < 3; i++) { + dgcout << " "; + printVec2D(tri.points[i]); + } + dgcout << "\n"; + printVec2D(tri.cenpoint); + dgcout << "\n"; - dgcout << "}\n"; + dgcout << "}\n"; } /* void PlaneTri */ /******************************************************************************/ -Vec3D vecAdd (const Vec3D& A, const Vec3D& B) +Vec3D vecAdd(const Vec3D &A, const Vec3D &B) /* Calculate and return A + B. */ { - Vec3D C; + Vec3D C; - C.x = A.x + B.x; - C.y = A.y + B.y; - C.z = A.z + B.z; + C.x = A.x + B.x; + C.y = A.y + B.y; + C.z = A.z + B.z; - return C; + return C; } /* Vec3D vecAdd */ /******************************************************************************/ -Vec3D vecSub (const Vec3D& A, const Vec3D& B) +Vec3D vecSub(const Vec3D &A, const Vec3D &B) /* Calculate and return A - B. */ { - Vec3D C; + Vec3D C; - C.x = A.x - B.x; - C.y = A.y - B.y; - C.z = A.z - B.z; + C.x = A.x - B.x; + C.y = A.y - B.y; + C.z = A.z - B.z; - return C; + return C; } /* Vec3D vecSub */ /******************************************************************************/ -Vec3D vecCross (const Vec3D& A, const Vec3D& B) +Vec3D vecCross(const Vec3D &A, const Vec3D &B) /* Calculate and return A x B. */ { - Vec3D C; + Vec3D C; - C.x = A.y * B.z - A.z * B.y; - C.y = A.z * B.x - A.x * B.z; - C.z = A.x * B.y - A.y * B.x; + C.x = A.y * B.z - A.z * B.y; + C.y = A.z * B.x - A.x * B.z; + C.z = A.x * B.y - A.y * B.x; - return C; + return C; } /* Vec3D vecCross */ /******************************************************************************/ -long double vecMag (const Vec3D& V) +long double vecMag(const Vec3D &V) /* Calculate and return the magnitude of a vector V. */ { - return sqrtl(V.x * V.x + V.y * V.y + V.z * V.z); + return sqrtl(V.x * V.x + V.y * V.y + V.z * V.z); } /* long double vecMag */ /******************************************************************************/ -Vec3D vecNormalize (const Vec3D& V) +Vec3D vecNormalize(const Vec3D &V) /* Return the normalized form of V. */ { - Vec3D C; + Vec3D C; - C.x = V.x/sqrtl(V.x * V.x + V.y * V.y + V.z * V.z); - C.y = V.y/sqrtl(V.x * V.x + V.y * V.y + V.z * V.z); - C.z = V.z/sqrtl(V.x * V.x + V.y * V.y + V.z * V.z); + C.x = V.x / sqrtl(V.x * V.x + V.y * V.y + V.z * V.z); + C.y = V.y / sqrtl(V.x * V.x + V.y * V.y + V.z * V.z); + C.z = V.z / sqrtl(V.x * V.x + V.y * V.y + V.z * V.z); - return C; + return C; -} /* Vec3D vecNormalize */ +} /* Vec3D vecNormalize */ /******************************************************************************/ -long double vecDot (const Vec3D& A, const Vec3D& B) +long double vecDot(const Vec3D &A, const Vec3D &B) /* Calculate and return the dot product of two vectors. */ { - return A.x * B.x + A.y * B.y + A.z * B.z; + return A.x * B.x + A.y * B.y + A.z * B.z; -} /* long double vecDot */ +} /* long double vecDot */ /******************************************************************************/ -long double sqrMetersToExcessD (long double area) -{ - return area * 360.0L / (4.0L * M_PI * DgGeoSphRF::earthRadiusKM() * +long double sqrMetersToExcessD(long double area) { + return area * 360.0L / + (4.0L * M_PI * DgGeoSphRF::earthRadiusKM() * DgGeoSphRF::earthRadiusKM()); } /* long double metersToExcessD */ /******************************************************************************/ -long double metersToGCDegrees (long double meters) -{ - long double earthCircum = (2.0L * M_PI * DgGeoSphRF::earthRadiusKM()); - return meters * 360.0L / earthCircum; +long double metersToGCDegrees(long double meters) { + long double earthCircum = (2.0L * M_PI * DgGeoSphRF::earthRadiusKM()); + return meters * 360.0L / earthCircum; } /* long double metersToGCDegrees */ @@ -467,73 +479,82 @@ long double maxval(long double val1, long double val2) /* return the maxmum of two variables */ - { +{ long double maxx; - if (val1>val2) maxx=val1; - else maxx=val2; + if (val1 > val2) + maxx = val1; + else + maxx = val2; return maxx; - } /* long double maxval */ +} /* long double maxval */ /******************************************************************************/ long double minval(long double val1, long double val2) /* return the minmum of two variables */ - { +{ long double minn; - if (val11.0L) v.z = 1.0L; - if (v.z<-1.0L) v.z = -1.0L; - sv.lat=asinl(v.z); - if ((sv.lat== M_PI_2) || (sv.lat==-M_PI_2)) sv.lon=0.0L; - else sv.lon=atan2l(v.y,v.x); - return sv; - } - else - { - dgcerr << "Error: in function xyzll, asin domain error.\n"; - return sv; + if (fabsl(v.z) - 1.0L < PRECISION) { + if (v.z > 1.0L) + v.z = 1.0L; + if (v.z < -1.0L) + v.z = -1.0L; + sv.lat = asinl(v.z); + if ((sv.lat == M_PI_2) || (sv.lat == -M_PI_2)) + sv.lon = 0.0L; + else + sv.lon = atan2l(v.y, v.x); + return sv; + } else { + dgcerr << "Error: in function xyzll, asin domain error.\n"; + return sv; } - } /* GeoCoord xyzll */ +} /* GeoCoord xyzll */ /******************************************************************************/ GeoCoord sphTricenpoint(GeoCoord sp[3]) @@ -541,125 +562,139 @@ GeoCoord sphTricenpoint(GeoCoord sp[3]) Calculate and return the center point of a sphere triangle Input and output unit is degree. */ - { +{ Vec3D p[3], cp, cpn; GeoCoord cv; int i; - for (i=0;i<3;i++) p[i]=llxyz(sp[i]); - cp.x=(p[0].x+p[1].x+p[2].x)/3; - cp.y=(p[0].y+p[1].y+p[2].y)/3; - cp.z=(p[0].z+p[1].z+p[2].z)/3; - cpn=vecNormalize(cp); - cv=xyzll(cpn); + for (i = 0; i < 3; i++) + p[i] = llxyz(sp[i]); + cp.x = (p[0].x + p[1].x + p[2].x) / 3; + cp.y = (p[0].y + p[1].y + p[2].y) / 3; + cp.z = (p[0].z + p[1].z + p[2].z) / 3; + cpn = vecNormalize(cp); + cv = xyzll(cpn); return cv; - } /* GeoCoord sphTricenpoint */ +} /* GeoCoord sphTricenpoint */ /******************************************************************************/ -long double chorddist(const GeoCoord& ll1, const GeoCoord& ll2) +long double chorddist(const GeoCoord &ll1, const GeoCoord &ll2) /* Calculate the chord distance between two points on the sphere Input unit is degree. Output unit regard the earth's radius as one unit. */ - { - //long double la1,lo1,la2,lo2; - Vec3D p1,p2; - p1.x=cosl(ll1.lat)*cosl(ll1.lon); - p1.y=cosl(ll1.lat)*sinl(ll1.lon); - p1.z=sinl(ll1.lat); - p2.x=cosl(ll2.lat)*cosl(ll2.lon); - p2.y=cosl(ll2.lat)*sinl(ll2.lon); - p2.z=sinl(ll2.lat); - return - sqrtl((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)+(p1.z-p2.z)*(p1.z-p2.z)); +{ + // long double la1,lo1,la2,lo2; + Vec3D p1, p2; + p1.x = cosl(ll1.lat) * cosl(ll1.lon); + p1.y = cosl(ll1.lat) * sinl(ll1.lon); + p1.z = sinl(ll1.lat); + p2.x = cosl(ll2.lat) * cosl(ll2.lon); + p2.y = cosl(ll2.lat) * sinl(ll2.lon); + p2.z = sinl(ll2.lat); + return sqrtl((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y) + + (p1.z - p2.z) * (p1.z - p2.z)); - } /* long double chorddist */ +} /* long double chorddist */ /******************************************************************************/ -long double spheredist(const GeoCoord& ll1, const GeoCoord& ll2) +long double spheredist(const GeoCoord &ll1, const GeoCoord &ll2) /* Calculate the sphere distance between two points on the sphere Input and Output unit is in radius. */ - { +{ long double sd; - sd=cosl(M_PI_2-ll1.lat)*cosl(M_PI_2-ll2.lat)+ - sinl(M_PI_2-ll1.lat)*sinl(M_PI_2-ll2.lat)*cosl(ll1.lon-ll2.lon); - if (sd>1) sd=1; - if (sd<-1) sd=-1; + sd = + cosl(M_PI_2 - ll1.lat) * cosl(M_PI_2 - ll2.lat) + + sinl(M_PI_2 - ll1.lat) * sinl(M_PI_2 - ll2.lat) * cosl(ll1.lon - ll2.lon); + if (sd > 1) + sd = 1; + if (sd < -1) + sd = -1; return acosl(sd); - }/* long double spheredist */ +} /* long double spheredist */ /******************************************************************************/ -void sphTriSolve(SphTri* tri) +void sphTriSolve(SphTri *tri) /* Input: three vertices's lat and lon in radius. Output: three edges, three angles, area. - the unit for length is kilometer, the unit for angle is degree + the unit for length is kilometer, the unit for angle is degree */ - { +{ int i; - long double l1[2],l2[2],l3[2],p,mindist; - - mindist=5; - l1[0]=tri->verts[0].lat; l1[1]=tri->verts[0].lon; - l2[0]=tri->verts[1].lat; l2[1]=tri->verts[1].lon; - l3[0]=tri->verts[2].lat; l3[1]=tri->verts[2].lon; - tri->edges[0]=acosl(cosl(M_PI_2-l2[0])*cosl(M_PI_2-l3[0])+ - sinl(M_PI_2-l2[0])*sinl(M_PI_2-l3[0])*cosl(l2[1]-l3[1])); - tri->edges[1]=acosl(cosl(M_PI_2-l1[0])*cosl(M_PI_2-l3[0])+ - sinl(M_PI_2-l1[0])*sinl(M_PI_2-l3[0])*cosl(l1[1]-l3[1])); - tri->edges[2]=acosl(cosl(M_PI_2-l2[0])*cosl(M_PI_2-l1[0])+ - sinl(M_PI_2-l2[0])*sinl(M_PI_2-l1[0])*cosl(l2[1]-l1[1])); - tri->angles[0]=acosl((cosl(tri->edges[0])-cosl(tri->edges[1])* - cosl(tri->edges[2]))/(sinl(tri->edges[1])*sinl(tri->edges[2]))); - tri->angles[1]=acosl((cosl(tri->edges[1])-cosl(tri->edges[0])* - cosl(tri->edges[2]))/(sinl(tri->edges[0])*sinl(tri->edges[2]))); - tri->angles[2]=acosl((cosl(tri->edges[2])-cosl(tri->edges[0])* - cosl(tri->edges[1]))/(sinl(tri->edges[0])*sinl(tri->edges[1]))); - for (i=0;i<3;i++) tri->edges[i]=tri->edges[i]*DgGeoSphRF::earthRadiusKM(); - if (tri->edges[0]edges[0]+tri->edges[1]+tri->edges[2])/2; - tri->area=sqrtl(p*(p-tri->edges[0])*(p-tri->edges[1])*(p-tri->edges[2])); - } - else tri->area=(tri->angles[0]+tri->angles[1]+tri->angles[2]-M_PI)* - DgGeoSphRF::earthRadiusKM()*DgGeoSphRF::earthRadiusKM(); - - } /* void sphTriSolve(SphTri* tri) */ - -/******************************************************************************/ -GeoCoord GCmidpoint(const GeoCoord& pp1, const GeoCoord& pp2) -{ - Vec3D pt1,pt2,mpt; + long double l1[2], l2[2], l3[2], p, mindist; + + mindist = 5; + l1[0] = tri->verts[0].lat; + l1[1] = tri->verts[0].lon; + l2[0] = tri->verts[1].lat; + l2[1] = tri->verts[1].lon; + l3[0] = tri->verts[2].lat; + l3[1] = tri->verts[2].lon; + tri->edges[0] = + acosl(cosl(M_PI_2 - l2[0]) * cosl(M_PI_2 - l3[0]) + + sinl(M_PI_2 - l2[0]) * sinl(M_PI_2 - l3[0]) * cosl(l2[1] - l3[1])); + tri->edges[1] = + acosl(cosl(M_PI_2 - l1[0]) * cosl(M_PI_2 - l3[0]) + + sinl(M_PI_2 - l1[0]) * sinl(M_PI_2 - l3[0]) * cosl(l1[1] - l3[1])); + tri->edges[2] = + acosl(cosl(M_PI_2 - l2[0]) * cosl(M_PI_2 - l1[0]) + + sinl(M_PI_2 - l2[0]) * sinl(M_PI_2 - l1[0]) * cosl(l2[1] - l1[1])); + tri->angles[0] = + acosl((cosl(tri->edges[0]) - cosl(tri->edges[1]) * cosl(tri->edges[2])) / + (sinl(tri->edges[1]) * sinl(tri->edges[2]))); + tri->angles[1] = + acosl((cosl(tri->edges[1]) - cosl(tri->edges[0]) * cosl(tri->edges[2])) / + (sinl(tri->edges[0]) * sinl(tri->edges[2]))); + tri->angles[2] = + acosl((cosl(tri->edges[2]) - cosl(tri->edges[0]) * cosl(tri->edges[1])) / + (sinl(tri->edges[0]) * sinl(tri->edges[1]))); + for (i = 0; i < 3; i++) + tri->edges[i] = tri->edges[i] * DgGeoSphRF::earthRadiusKM(); + if (tri->edges[0] < mindist) { + p = (tri->edges[0] + tri->edges[1] + tri->edges[2]) / 2; + tri->area = sqrtl(p * (p - tri->edges[0]) * (p - tri->edges[1]) * + (p - tri->edges[2])); + } else + tri->area = (tri->angles[0] + tri->angles[1] + tri->angles[2] - M_PI) * + DgGeoSphRF::earthRadiusKM() * DgGeoSphRF::earthRadiusKM(); + +} /* void sphTriSolve(SphTri* tri) */ + +/******************************************************************************/ +GeoCoord GCmidpoint(const GeoCoord &pp1, const GeoCoord &pp2) { + Vec3D pt1, pt2, mpt; GeoCoord mp; - pt1=llxyz(pp1); - pt2=llxyz(pp2); - mpt.x = (pt1.x+pt2.x)/2.0L; - mpt.y = (pt1.y+pt2.y)/2.0L; - mpt.z = (pt1.z+pt2.z)/2.0L; + pt1 = llxyz(pp1); + pt2 = llxyz(pp2); + mpt.x = (pt1.x + pt2.x) / 2.0L; + mpt.y = (pt1.y + pt2.y) / 2.0L; + mpt.z = (pt1.z + pt2.z) / 2.0L; mpt = vecNormalize(mpt); mp = xyzll(mpt); return mp; } /******************************************************************************/ -long double Azimuth(const GeoCoord& pt1, const GeoCoord& pt2) +long double Azimuth(const GeoCoord &pt1, const GeoCoord &pt2) /* return the azimuth of pt2 relatived to pt1 output azimuth unit is radius */ - { +{ long double az; - az=atan2l(cosl(pt2.lat)*sinl(pt2.lon-pt1.lon), cosl(pt1.lat)*sinl(pt2.lat)- - sinl(pt1.lat)*cosl(pt2.lat)*cosl(pt2.lon-pt1.lon)); + az = atan2l(cosl(pt2.lat) * sinl(pt2.lon - pt1.lon), + cosl(pt1.lat) * sinl(pt2.lat) - + sinl(pt1.lat) * cosl(pt2.lat) * cosl(pt2.lon - pt1.lon)); return az; - } /* long double Azimuth */ +} /* long double Azimuth */ /******************************************************************************/ -GeoCoord GCintersect(const GeoCoord& sv11, const GeoCoord& sv12, - const GeoCoord& sv21, const GeoCoord& sv22, int sign) +GeoCoord GCintersect(const GeoCoord &sv11, const GeoCoord &sv12, + const GeoCoord &sv21, const GeoCoord &sv22, int sign) /* return the intersect point of two great circle sign=1: two great circle segment with ends of sv11 and sv12, @@ -667,217 +702,219 @@ GeoCoord GCintersect(const GeoCoord& sv11, const GeoCoord& sv12, sign=0: two whole great circle with one pass sv11 and sv12, one pass sv21 and sv22, return the intersect point on North Hemisphere */ - { +{ GeoCoord pt; - Vec3D nn1,nn2, - p11,p12, - p21,p22, - pp,pp2; - long double a,b,maxlon,minlon; - -/* calculate the intersect point of two great circle */ - p11=llxyz(sv11); - p12=llxyz(sv12); - p21=llxyz(sv21); - p22=llxyz(sv22); - nn1.x=p11.y*p12.z-p12.y*p11.z; - nn1.y=-p11.x*p12.z+p12.x*p11.z; - nn1.z=p11.x*p12.y-p12.x*p11.y; - nn2.x=p21.y*p22.z-p22.y*p21.z; - nn2.y=-p21.x*p22.z+p22.x*p21.z; - nn2.z=p21.x*p22.y-p22.x*p21.y; - if ((nn2.z*nn1.y-nn1.z*nn2.y)!= 0.0L) { - b=(nn1.x*nn2.y-nn2.x*nn1.y)/(nn2.z*nn1.y-nn1.z*nn2.y); - a=(nn2.x*nn1.z-nn1.x*nn2.z)/(nn1.y*nn2.z-nn2.y*nn1.z); - pp.x=1/sqrtl(a*a+b*b+1); - pp.y=a*pp.x; - pp.z=b*pp.x; - } else if (((nn2.z*nn1.y-nn1.z*nn2.y)==0.0L) && - ((nn1.x*nn2.y-nn2.x*nn1.y)==0.0L) && ((nn1.x*nn2.z-nn2.x*nn1.z)==0.0L)) { + Vec3D nn1, nn2, p11, p12, p21, p22, pp, pp2; + long double a, b, maxlon, minlon; + + /* calculate the intersect point of two great circle */ + p11 = llxyz(sv11); + p12 = llxyz(sv12); + p21 = llxyz(sv21); + p22 = llxyz(sv22); + nn1.x = p11.y * p12.z - p12.y * p11.z; + nn1.y = -p11.x * p12.z + p12.x * p11.z; + nn1.z = p11.x * p12.y - p12.x * p11.y; + nn2.x = p21.y * p22.z - p22.y * p21.z; + nn2.y = -p21.x * p22.z + p22.x * p21.z; + nn2.z = p21.x * p22.y - p22.x * p21.y; + if ((nn2.z * nn1.y - nn1.z * nn2.y) != 0.0L) { + b = (nn1.x * nn2.y - nn2.x * nn1.y) / (nn2.z * nn1.y - nn1.z * nn2.y); + a = (nn2.x * nn1.z - nn1.x * nn2.z) / (nn1.y * nn2.z - nn2.y * nn1.z); + pp.x = 1 / sqrtl(a * a + b * b + 1); + pp.y = a * pp.x; + pp.z = b * pp.x; + } else if (((nn2.z * nn1.y - nn1.z * nn2.y) == 0.0L) && + ((nn1.x * nn2.y - nn2.x * nn1.y) == 0.0L) && + ((nn1.x * nn2.z - nn2.x * nn1.z) == 0.0L)) { report("Error in GCintersect: the two great circle planes are parallel.\n", DgBase::Fatal); - } else if (((nn2.z*nn1.y-nn1.z*nn2.y)==0.0L) && (nn1.z!=0.0L)) { - pp.x=0.0L; - pp.y=1.0L/sqrtl(1+nn1.y*nn1.y/nn1.z/nn1.z); - pp.z=-nn1.y/nn1.z*pp.y; - } else if (((nn2.z*nn1.y-nn1.z*nn2.y)==0.0L) && (nn2.z!=0.0L)) { - pp.x=0.0L; - pp.y=1.0L/sqrtl(1.0L+nn2.y*nn2.y/nn2.z/nn2.z); - pp.z=-nn2.y/nn2.z*pp.y; - } else if (((nn2.z*nn1.y-nn1.z*nn2.y)==0.0L) && (nn1.y!=0.0L)) { - pp.x=0.0L; - pp.z=1/sqrtl(1+nn1.z*nn1.z/nn1.y/nn1.y); - pp.y=-nn1.z/nn1.y*pp.z; - } else if (((nn2.z*nn1.y-nn1.z*nn2.y)==0.0L) && (nn2.y!=0.0L)) { - pp.x=0.0L; - pp.z=1.0L/sqrtl(1.0L+nn2.z*nn2.z/nn2.y/nn2.y); - pp.y=-nn2.z/nn2.y*pp.z; - } - - if (sign==0) { - if (pp.z<0.0L) { - pp.x=0.0L-pp.x; - pp.y=-pp.y; - pp.z=-pp.z; + } else if (((nn2.z * nn1.y - nn1.z * nn2.y) == 0.0L) && (nn1.z != 0.0L)) { + pp.x = 0.0L; + pp.y = 1.0L / sqrtl(1 + nn1.y * nn1.y / nn1.z / nn1.z); + pp.z = -nn1.y / nn1.z * pp.y; + } else if (((nn2.z * nn1.y - nn1.z * nn2.y) == 0.0L) && (nn2.z != 0.0L)) { + pp.x = 0.0L; + pp.y = 1.0L / sqrtl(1.0L + nn2.y * nn2.y / nn2.z / nn2.z); + pp.z = -nn2.y / nn2.z * pp.y; + } else if (((nn2.z * nn1.y - nn1.z * nn2.y) == 0.0L) && (nn1.y != 0.0L)) { + pp.x = 0.0L; + pp.z = 1 / sqrtl(1 + nn1.z * nn1.z / nn1.y / nn1.y); + pp.y = -nn1.z / nn1.y * pp.z; + } else if (((nn2.z * nn1.y - nn1.z * nn2.y) == 0.0L) && (nn2.y != 0.0L)) { + pp.x = 0.0L; + pp.z = 1.0L / sqrtl(1.0L + nn2.z * nn2.z / nn2.y / nn2.y); + pp.y = -nn2.z / nn2.y * pp.z; + } + + if (sign == 0) { + if (pp.z < 0.0L) { + pp.x = 0.0L - pp.x; + pp.y = -pp.y; + pp.z = -pp.z; } - pt=xyzll(pp); + pt = xyzll(pp); return pt; } else { - /* judge if the point is on the two great circle segment */ + /* judge if the point is on the two great circle segment */ - pt=xyzll(pp); - maxlon=maxval(sv11.lon,sv12.lon); - minlon=minval(sv11.lon,sv12.lon); - if ((pt.lon<=maxlon) && (pt.lon>=minlon)) - return pt; - else { - pp2.x=-pp.x; - pp2.y=-pp.y; - pp2.z=-pp.z; - pt=xyzll(pp2); - if ((pt.lon<=maxlon) && (pt.lat>=minlon)) - return pt; - else { - dgcerr << "Error of GCintersect: the point is not on great circle segment.\n"; - pt.lat=UNDEFVAL; pt.lon=UNDEFVAL; + pt = xyzll(pp); + maxlon = maxval(sv11.lon, sv12.lon); + minlon = minval(sv11.lon, sv12.lon); + if ((pt.lon <= maxlon) && (pt.lon >= minlon)) return pt; + else { + pp2.x = -pp.x; + pp2.y = -pp.y; + pp2.z = -pp.z; + pt = xyzll(pp2); + if ((pt.lon <= maxlon) && (pt.lat >= minlon)) + return pt; + else { + dgcerr << "Error of GCintersect: the point is not on great circle " + "segment.\n"; + pt.lat = UNDEFVAL; + pt.lon = UNDEFVAL; + return pt; + } } } - } } /* GeoCoord GCintersect */ /******************************************************************************/ -long double GCptlat(long double lon, const GeoCoord& sv1, const GeoCoord& sv2) +long double GCptlat(long double lon, const GeoCoord &sv1, const GeoCoord &sv2) /* return latitude of the point on great circle segment with known longtitude */ - { - long double lat,a,b,c; - Vec3D p1,p2; - - p1.x=cosl(sv1.lat)*cosl(sv1.lon); - p1.y=cosl(sv1.lat)*sinl(sv1.lon); - p1.z=sinl(sv1.lat); - p2.x=cosl(sv2.lat)*cosl(sv2.lon); - p2.y=cosl(sv2.lat)*sinl(sv2.lon); - p2.z=sinl(sv2.lat); - - a=(p1.y*p2.z-p1.z*p2.y)*cosl(lon); - b=(p1.x*p2.z-p1.z*p2.x)*sinl(lon); - c=(p1.x*p2.y-p1.y*p2.x); - if (c!=0.0L) lat=atanl((b-a)/c); +{ + long double lat, a, b, c; + Vec3D p1, p2; + + p1.x = cosl(sv1.lat) * cosl(sv1.lon); + p1.y = cosl(sv1.lat) * sinl(sv1.lon); + p1.z = sinl(sv1.lat); + p2.x = cosl(sv2.lat) * cosl(sv2.lon); + p2.y = cosl(sv2.lat) * sinl(sv2.lon); + p2.z = sinl(sv2.lat); + + a = (p1.y * p2.z - p1.z * p2.y) * cosl(lon); + b = (p1.x * p2.z - p1.z * p2.x) * sinl(lon); + c = (p1.x * p2.y - p1.y * p2.x); + if (c != 0.0L) + lat = atanl((b - a) / c); else { - lat = UNDEFVAL; - dgcerr << "Error in GCptlat: the two end points are at one longitude.\n"; + lat = UNDEFVAL; + dgcerr << "Error in GCptlat: the two end points are at one longitude.\n"; } - return(lat); - } /* long double GCptlat */ + return (lat); +} /* long double GCptlat */ /******************************************************************************/ -int ptinsphtri(const GeoCoord& pt, GeoCoord sTri[3]) +int ptinsphtri(const GeoCoord &pt, GeoCoord sTri[3]) /* decide if a point is in a spherical triangle return 1 if it is in and return 0 if not */ { - int i; - long double p0, t0; - Vec3D pp, ptri[3]; - - for (i = 0; i < 3; i++) ptri[i] = llxyz(sTri[i]); - pp = llxyz(pt); - p0 = pp.x * (ptri[1].y * ptri[2].z - ptri[2].y * ptri[1].z) - - pp.y * (ptri[1].x * ptri[2].z - ptri[2].x * ptri[1].z) + - pp.z * (ptri[1].x * ptri[2].y - ptri[2].x * ptri[1].y); + int i; + long double p0, t0; + Vec3D pp, ptri[3]; + + for (i = 0; i < 3; i++) + ptri[i] = llxyz(sTri[i]); + pp = llxyz(pt); + p0 = pp.x * (ptri[1].y * ptri[2].z - ptri[2].y * ptri[1].z) - + pp.y * (ptri[1].x * ptri[2].z - ptri[2].x * ptri[1].z) + + pp.z * (ptri[1].x * ptri[2].y - ptri[2].x * ptri[1].y); + + t0 = ptri[0].x * (ptri[1].y * ptri[2].z - ptri[2].y * ptri[1].z) - + ptri[0].y * (ptri[1].x * ptri[2].z - ptri[2].x * ptri[1].z) + + ptri[0].z * (ptri[1].x * ptri[2].y - ptri[2].x * ptri[1].y); + + if (p0 * t0 < 0.0L) { + return 0; + } else { + p0 = pp.x * (ptri[0].y * ptri[2].z - ptri[2].y * ptri[0].z) - + pp.y * (ptri[0].x * ptri[2].z - ptri[2].x * ptri[0].z) + + pp.z * (ptri[0].x * ptri[2].y - ptri[2].x * ptri[0].y); - t0 = ptri[0].x * (ptri[1].y * ptri[2].z - ptri[2].y * ptri[1].z) - - ptri[0].y * (ptri[1].x * ptri[2].z - ptri[2].x * ptri[1].z) + - ptri[0].z * (ptri[1].x * ptri[2].y - ptri[2].x * ptri[1].y); + t0 = ptri[1].x * (ptri[0].y * ptri[2].z - ptri[2].y * ptri[0].z) - + ptri[1].y * (ptri[0].x * ptri[2].z - ptri[2].x * ptri[0].z) + + ptri[1].z * (ptri[0].x * ptri[2].y - ptri[2].x * ptri[0].y); - if (p0 * t0 < 0.0L) - { + if (p0 * t0 < 0.0L) { return 0; - } - else - { - p0 = pp.x * (ptri[0].y * ptri[2].z - ptri[2].y * ptri[0].z) - - pp.y * (ptri[0].x * ptri[2].z - ptri[2].x * ptri[0].z) + - pp.z * (ptri[0].x * ptri[2].y - ptri[2].x * ptri[0].y); - - t0 = ptri[1].x * (ptri[0].y * ptri[2].z - ptri[2].y * ptri[0].z) - - ptri[1].y * (ptri[0].x * ptri[2].z - ptri[2].x * ptri[0].z) + - ptri[1].z * (ptri[0].x * ptri[2].y - ptri[2].x * ptri[0].y); - - if (p0 * t0 < 0.0L) - { - return 0; - } - else - { - p0 = pp.x * (ptri[0].y * ptri[1].z - ptri[1].y * ptri[0].z) - - pp.y * (ptri[0].x * ptri[1].z - ptri[1].x * ptri[0].z) + - pp.z * (ptri[0].x * ptri[1].y - ptri[1].x * ptri[0].y); + } else { + p0 = pp.x * (ptri[0].y * ptri[1].z - ptri[1].y * ptri[0].z) - + pp.y * (ptri[0].x * ptri[1].z - ptri[1].x * ptri[0].z) + + pp.z * (ptri[0].x * ptri[1].y - ptri[1].x * ptri[0].y); - t0 = ptri[2].x * (ptri[0].y * ptri[1].z - ptri[1].y * ptri[0].z) - - ptri[2].y * (ptri[0].x * ptri[1].z - ptri[1].x * ptri[0].z) + - ptri[2].z * (ptri[0].x * ptri[1].y - ptri[1].x * ptri[0].y); + t0 = ptri[2].x * (ptri[0].y * ptri[1].z - ptri[1].y * ptri[0].z) - + ptri[2].y * (ptri[0].x * ptri[1].z - ptri[1].x * ptri[0].z) + + ptri[2].z * (ptri[0].x * ptri[1].y - ptri[1].x * ptri[0].y); - return (p0 * t0 >= 0.0L); - } - } + return (p0 * t0 >= 0.0L); + } + } -} /* ptinsphtri */ +} /* ptinsphtri */ /******************************************************************************/ -GeoCoord GCdaz(const GeoCoord& pt, long double distance, long double az) +GeoCoord GCdaz(const GeoCoord &pt, long double distance, long double az) /* compute and return the point whose azimuth (az) and distance (distance) relative to a known point (pt) are given. input parameter distance and az are in radius, pt in degree. */ - { +{ GeoCoord pt2; long double sinlat, sinlon, coslon; - if ((fabsl(az)1.0L) sinlat=1.0L; - if (sinlat<-1.0L) sinlat=-1.0L; - pt2.lat=asinl(sinlat); - if ((pt2.lat==M_PI_2) || (pt2.lat==-M_PI_2)) pt2.lon=0.0L; - else - { - sinlon=sinl(az)*sinl(distance)/cosl(pt2.lat); - coslon=(cosl(distance)-sinl(pt.lat)*sinl(pt2.lat))/ - cosl(pt.lat)/cosl(pt2.lat); - if (sinlon>1.0L) sinlon=1.0L; - if (sinlon<-1.0L) sinlon=-1.0L; - if (coslon>1.0L) coslon=1.0L; - if (coslon<-1.0L) coslon=-1.0L; - pt2.lon=pt.lon+atan2l(sinlon,coslon); - } - if (pt2.lon>M_PI+PRECISION) pt2.lon -=2.0L*M_PI; - if (pt2.lon<-M_PI-PRECISION) pt2.lon +=2.0L*M_PI; + if ((fabsl(az) < PRECISION) || (fabsl(fabsl(az) - M_PI) < PRECISION)) { + if (fabsl(az) < PRECISION) + pt2.lat = pt.lat + distance; + else + pt2.lat = pt.lat - distance; + pt2.lon = pt.lon; + if (fabsl(pt2.lat - M_PI_2) < PRECISION) { + pt2.lat = M_PI_2; + pt2.lon = 0.0L; + } + if (fabsl(pt2.lat + M_PI_2) < PRECISION) { + pt2.lat = -M_PI; + pt2.lon = 0.0L; + } + } else { + sinlat = sinl(pt.lat) * cosl(distance) + + cosl(pt.lat) * sinl(distance) * cosl(az); + if (sinlat > 1.0L) + sinlat = 1.0L; + if (sinlat < -1.0L) + sinlat = -1.0L; + pt2.lat = asinl(sinlat); + if ((pt2.lat == M_PI_2) || (pt2.lat == -M_PI_2)) + pt2.lon = 0.0L; + else { + sinlon = sinl(az) * sinl(distance) / cosl(pt2.lat); + coslon = (cosl(distance) - sinl(pt.lat) * sinl(pt2.lat)) / cosl(pt.lat) / + cosl(pt2.lat); + if (sinlon > 1.0L) + sinlon = 1.0L; + if (sinlon < -1.0L) + sinlon = -1.0L; + if (coslon > 1.0L) + coslon = 1.0L; + if (coslon < -1.0L) + coslon = -1.0L; + pt2.lon = pt.lon + atan2l(sinlon, coslon); + } + if (pt2.lon > M_PI + PRECISION) + pt2.lon -= 2.0L * M_PI; + if (pt2.lon < -M_PI - PRECISION) + pt2.lon += 2.0L * M_PI; } return pt2; - } /* GeoCoord GCdaz */ +} /* GeoCoord GCdaz */ /******************************************************************************/ diff --git a/src/lib/dglib/lib/DgProjGnomonicRF.cpp b/src/lib/dglib/lib/DgProjGnomonicRF.cpp index 0438ee8c..50fa34d3 100644 --- a/src/lib/dglib/lib/DgProjGnomonicRF.cpp +++ b/src/lib/dglib/lib/DgProjGnomonicRF.cpp @@ -24,147 +24,145 @@ #include -#include -#include #include "proj4.h" +#include +#include #define EPS10 1.e-10L #define N_POLE 0 #define S_POLE 1 -#define EQUIT 2 -#define OBLIQ 3 +#define EQUIT 2 +#define OBLIQ 3 //////////////////////////////////////////////////////////////////////////////// -DgProjGnomonicRF::DgProjGnomonicRF(DgRFNetwork& networkIn, const string& nameIn, - const DgGeoCoord& proj0In, long double x0In, long double y0In, long double k0In, - long double to_meterIn, long double fr_meterIn) - : DgGeoProjRF (networkIn, nameIn, proj0In, x0In, y0In, k0In, - to_meterIn, fr_meterIn) -{ - if (fabs(fabs(phi0()) - M_PI_2) < EPS10) - mode_ = phi0() < 0.0L ? S_POLE : N_POLE; - else if (fabs(phi0()) < EPS10) - mode_ = EQUIT; - else { - mode_ = OBLIQ; - sinph0_ = sin(phi0()); - cosph0_ = cos(phi0()); - } +DgProjGnomonicRF::DgProjGnomonicRF(DgRFNetwork &networkIn, const string &nameIn, + const DgGeoCoord &proj0In, long double x0In, + long double y0In, long double k0In, + long double to_meterIn, + long double fr_meterIn) + : DgGeoProjRF(networkIn, nameIn, proj0In, x0In, y0In, k0In, to_meterIn, + fr_meterIn) { + if (fabs(fabs(phi0()) - M_PI_2) < EPS10) + mode_ = phi0() < 0.0L ? S_POLE : N_POLE; + else if (fabs(phi0()) < EPS10) + mode_ = EQUIT; + else { + mode_ = OBLIQ; + sinph0_ = sin(phi0()); + cosph0_ = cos(phi0()); + } } // DgProjGnomonicRF::DgProjGnomonicRF //////////////////////////////////////////////////////////////////////////////// -DgDVec2D -DgProjGnomonicRF::projForward (const DgGeoCoord& addIn, - const DgEllipsoidRF&) const -{ - //cout << "gnom projForward: " << *this << " coord: " << addIn << endl; - - DgDVec2D xy(DgDVec2D::undefDgDVec2D); - - long double coslam, cosphi, sinphi; - - sinphi = sin(addIn.lat()); - cosphi = cos(addIn.lat()); - coslam = cos(addIn.lon()); - - switch (mode_) { - case EQUIT: - xy.setY(cosphi * coslam); - break; - case OBLIQ: - xy.setY(sinph0_ * sinphi + cosph0_ * cosphi * coslam); - break; - case S_POLE: - xy.setY(- sinphi); - break; - case N_POLE: - xy.setY(sinphi); - break; - } - - if (xy.y() <= EPS10) - { - ::report(string("DgProjGnomonicRF::projForward() point out of range\n") + - string("proj0: ") + string(proj0()) + - string("\nprojecting point: ") + string(addIn), DgBase::Fatal); - } - - xy.setY(1.0L / xy.y()); - xy.setX(xy.y() * cosphi * sin(addIn.lon())); - - switch (mode_) { - case EQUIT: - xy.setY(xy.y() * sinphi); - break; - case OBLIQ: - xy.setY(xy.y() * (cosph0_ * sinphi - sinph0_ * cosphi * coslam)); - break; - case N_POLE: - coslam = - coslam; - FALLTHROUGH - case S_POLE: - xy.setY(xy.y() * cosphi * coslam); - break; - } - - //cout << " => " << xy << endl; - - return (xy); +DgDVec2D DgProjGnomonicRF::projForward(const DgGeoCoord &addIn, + const DgEllipsoidRF &) const { + // cout << "gnom projForward: " << *this << " coord: " << addIn << endl; + + DgDVec2D xy(DgDVec2D::undefDgDVec2D); + + long double coslam, cosphi, sinphi; + + sinphi = sin(addIn.lat()); + cosphi = cos(addIn.lat()); + coslam = cos(addIn.lon()); + + switch (mode_) { + case EQUIT: + xy.setY(cosphi * coslam); + break; + case OBLIQ: + xy.setY(sinph0_ * sinphi + cosph0_ * cosphi * coslam); + break; + case S_POLE: + xy.setY(-sinphi); + break; + case N_POLE: + xy.setY(sinphi); + break; + } + + if (xy.y() <= EPS10) { + ::report(string("DgProjGnomonicRF::projForward() point out of range\n") + + string("proj0: ") + string(proj0()) + + string("\nprojecting point: ") + string(addIn), + DgBase::Fatal); + } + + xy.setY(1.0L / xy.y()); + xy.setX(xy.y() * cosphi * sin(addIn.lon())); + + switch (mode_) { + case EQUIT: + xy.setY(xy.y() * sinphi); + break; + case OBLIQ: + xy.setY(xy.y() * (cosph0_ * sinphi - sinph0_ * cosphi * coslam)); + break; + case N_POLE: + coslam = -coslam; + FALLTHROUGH + case S_POLE: + xy.setY(xy.y() * cosphi * coslam); + break; + } + + // cout << " => " << xy << endl; + + return (xy); } // DgDVec2D DgProjGnomonicRF::projForward //////////////////////////////////////////////////////////////////////////////// -DgGeoCoord -DgProjGnomonicRF::projInverse (const DgDVec2D& addIn, - const DgEllipsoidRF&) const +DgGeoCoord DgProjGnomonicRF::projInverse(const DgDVec2D &addIn, + const DgEllipsoidRF &) const // // spheroid only; needs to be verified at some point // { - DgDVec2D xy(addIn); - DgGeoCoord lp; - - long double rh, cosz, sinz; - - rh = hypot(xy.x(), xy.y()); - lp.setLat(atan(rh)); - sinz = sin(lp.lat()); - cosz = sqrt(1. - sinz * sinz); - if (fabs(rh) <= EPS10) { - lp.setLat(phi0()); - lp.setLon(0.0L); - } else { - switch (mode_) { - case OBLIQ: - lp.setLat(cosz * sinph0_ + xy.y() * sinz * cosph0_ / rh); - if (fabs(lp.lat()) >= 1.0L) - lp.setLat(lp.lat() > 0.0L ? M_PI_2 : - M_PI_2); - else - lp.setLat(asin(lp.lat())); - xy.setY((cosz - sinph0_ * sin(lp.lat())) * rh); - xy.setX(xy.x() * sinz * cosph0_); - break; - case EQUIT: - lp.setLat(xy.y() * sinz / rh); - if (fabs(lp.lat()) >= 1.0L) - lp.setLat(lp.lat() > 0.0L ? M_PI_2 : - M_PI_2); - else - lp.setLat(asin(lp.lat())); - xy.setY(cosz * rh); - xy.setX(xy.x() * sinz); - break; - case S_POLE: - lp.setLat(lp.lat() - M_PI_2); - break; - case N_POLE: - lp.setLat(M_PI_2 - lp.lat()); - xy.setY(-xy.y()); - break; - } - lp.setLon(atan2(xy.x(), xy.y())); - } - return (lp); + DgDVec2D xy(addIn); + DgGeoCoord lp; + + long double rh, cosz, sinz; + + rh = hypot(xy.x(), xy.y()); + lp.setLat(atan(rh)); + sinz = sin(lp.lat()); + cosz = sqrt(1. - sinz * sinz); + if (fabs(rh) <= EPS10) { + lp.setLat(phi0()); + lp.setLon(0.0L); + } else { + switch (mode_) { + case OBLIQ: + lp.setLat(cosz * sinph0_ + xy.y() * sinz * cosph0_ / rh); + if (fabs(lp.lat()) >= 1.0L) + lp.setLat(lp.lat() > 0.0L ? M_PI_2 : -M_PI_2); + else + lp.setLat(asin(lp.lat())); + xy.setY((cosz - sinph0_ * sin(lp.lat())) * rh); + xy.setX(xy.x() * sinz * cosph0_); + break; + case EQUIT: + lp.setLat(xy.y() * sinz / rh); + if (fabs(lp.lat()) >= 1.0L) + lp.setLat(lp.lat() > 0.0L ? M_PI_2 : -M_PI_2); + else + lp.setLat(asin(lp.lat())); + xy.setY(cosz * rh); + xy.setX(xy.x() * sinz); + break; + case S_POLE: + lp.setLat(lp.lat() - M_PI_2); + break; + case N_POLE: + lp.setLat(M_PI_2 - lp.lat()); + xy.setY(-xy.y()); + break; + } + lp.setLon(atan2(xy.x(), xy.y())); + } + return (lp); } // DgGeoCoord DgProjGnomonicRF::projInverse