35 # pragma warning (disable: 4701 4127) 43 : maxit2_(maxit1_ +
Math::digits() + 10)
47 , tiny_(sqrt(numeric_limits<real>::min()))
48 , tol0_(numeric_limits<real>::epsilon())
54 , tolb_(tol0_ * tol2_)
55 , xthresh_(1000 * tol2_)
60 , _ep2(_e2 /
Math::sq(_f1))
71 (_f > 0 ?
Math::asinh(sqrt(_ep2)) : atan(sqrt(-_e2))) /
83 , _etol2(0.1 * tol2_ /
84 sqrt( max(real(0.001), abs(_f)) * min(real(1), 1 - _f/2) / 2 ))
99 Math::real GeodesicExact::CosSeries(real sinx, real cosx,
100 const real c[],
int n) {
107 ar = 2 * (cosx - sinx) * (cosx + sinx),
108 y0 = n & 1 ? *--c : 0, y1 = 0;
113 y1 = ar * y0 - y1 + *--c;
114 y0 = ar * y1 - y0 + *--c;
116 return cosx * (y0 - y1);
120 unsigned caps)
const {
125 bool arcmode, real s12_a12,
127 real& lat2, real& lon2, real& azi2,
128 real& s12, real& m12,
129 real& M12, real& M21,
135 GenPosition(arcmode, s12_a12, outmask,
136 lat2, lon2, azi2, s12, m12, M12, M21, S12);
141 bool arcmode, real s12_a12,
142 unsigned caps)
const {
150 caps, arcmode, s12_a12);
155 unsigned caps)
const {
161 unsigned caps)
const {
165 Math::real GeodesicExact::GenInverse(real lat1, real lon1,
166 real lat2, real lon2,
167 unsigned outmask, real& s12,
168 real& salp1, real& calp1,
169 real& salp2, real& calp2,
170 real& m12, real& M12, real& M21,
177 int lonsign = lon12 >= 0 ? 1 : -1;
195 int swapp = abs(lat1) < abs(lat2) ? -1 : 1;
201 int latsign = lat1 < 0 ? 1 : -1;
216 real sbet1, cbet1, sbet2, cbet2, s12x, m12x;
224 Math::norm(sbet1, cbet1); cbet1 = max(tiny_, cbet1);
228 Math::norm(sbet2, cbet2); cbet2 = max(tiny_, cbet2);
238 if (cbet1 < -sbet1) {
240 sbet2 = sbet2 < 0 ? sbet1 : -sbet1;
242 if (abs(sbet2) == -sbet1)
247 dn1 = (_f >= 0 ? sqrt(1 + _ep2 *
Math::sq(sbet1)) :
248 sqrt(1 - _e2 *
Math::sq(cbet1)) / _f1),
249 dn2 = (_f >= 0 ? sqrt(1 + _ep2 *
Math::sq(sbet2)) :
250 sqrt(1 - _e2 *
Math::sq(cbet2)) / _f1);
254 bool meridian = lat1 == -90 || slam12 == 0;
261 calp1 = clam12; salp1 = slam12;
262 calp2 = 1; salp2 = 0;
266 ssig1 = sbet1, csig1 = calp1 * cbet1,
267 ssig2 = sbet2, csig2 = calp2 * cbet2;
270 sig12 = atan2(max(real(0), csig1 * ssig2 - ssig1 * csig2),
271 csig1 * csig2 + ssig1 * ssig2);
274 Lengths(E, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
276 s12x, m12x, dummy, M12, M21);
285 if (sig12 < 1 || m12x >= 0) {
287 if (sig12 < 3 * tiny_)
288 sig12 = m12x = s12x = 0;
298 real omg12 = 0, somg12 = 2, comg12 = 0;
301 (_f <= 0 || lon12s >= _f * 180)) {
304 calp1 = calp2 = 0; salp1 = salp2 = 1;
306 sig12 = omg12 = lam12 / _f1;
307 m12x = _b * sin(sig12);
309 M12 = M21 = cos(sig12);
312 }
else if (!meridian) {
319 sig12 = InverseStart(E, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
320 lam12, slam12, clam12,
321 salp1, calp1, salp2, calp2, dnm);
325 s12x = sig12 * _b * dnm;
326 m12x =
Math::sq(dnm) * _b * sin(sig12 / dnm);
328 M12 = M21 = cos(sig12 / dnm);
330 omg12 = lam12 / (_f1 * dnm);
346 real ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, domg12 = 0;
349 real salp1a = tiny_, calp1a = 1, salp1b = tiny_, calp1b = -1;
350 for (
bool tripn =
false, tripb =
false;
375 real v = Lambda12(sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
377 salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
378 E, domg12, numit < maxit1_, dv);
380 if (tripb || !(abs(v) >= (tripn ? 8 : 1) * tol0_))
break;
382 if (v > 0 && (numit > maxit1_ || calp1/salp1 > calp1b/salp1b))
383 { salp1b = salp1; calp1b = calp1; }
384 else if (v < 0 && (numit > maxit1_ || calp1/salp1 < calp1a/salp1a))
385 { salp1a = salp1; calp1a = calp1; }
386 if (numit < maxit1_ && dv > 0) {
390 sdalp1 = sin(dalp1), cdalp1 = cos(dalp1),
391 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1;
392 if (nsalp1 > 0 && abs(dalp1) <
Math::pi()) {
393 calp1 = calp1 * cdalp1 - salp1 * sdalp1;
399 tripn = abs(v) <= 16 * tol0_;
411 salp1 = (salp1a + salp1b)/2;
412 calp1 = (calp1a + calp1b)/2;
415 tripb = (abs(salp1a - salp1) + (calp1a - calp1) < tolb_ ||
416 abs(salp1 - salp1b) + (calp1 - calp1b) < tolb_);
420 Lengths(E, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
421 cbet1, cbet2, outmask, s12x, m12x, dummy, M12, M21);
426 if (outmask &
AREA) {
428 real sdomg12 = sin(domg12), cdomg12 = cos(domg12);
429 somg12 = slam12 * cdomg12 - clam12 * sdomg12;
430 comg12 = clam12 * cdomg12 + slam12 * sdomg12;
441 if (outmask &
AREA) {
444 salp0 = salp1 * cbet1,
447 if (calp0 != 0 && salp0 != 0) {
450 ssig1 = sbet1, csig1 = calp1 * cbet1,
451 ssig2 = sbet2, csig2 = calp2 * cbet2,
453 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2),
455 A4 =
Math::sq(_a) * calp0 * salp0 * _e2;
461 B41 = CosSeries(ssig1, csig1, C4a, nC4_),
462 B42 = CosSeries(ssig2, csig2, C4a, nC4_);
463 S12 = A4 * (B42 - B41);
470 somg12 = sin(omg12); comg12 = cos(omg12);
476 comg12 > -real(0.7071) &&
477 sbet2 - sbet1 < real(1.75)) {
481 real domg12 = 1 + comg12, dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
482 alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
483 domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );
487 salp12 = salp2 * calp1 - calp2 * salp1,
488 calp12 = calp2 * calp1 + salp2 * salp1;
493 if (salp12 == 0 && calp12 < 0) {
494 salp12 = tiny_ * calp1;
497 alp12 = atan2(salp12, calp12);
500 S12 *= swapp * lonsign * latsign;
513 salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
514 salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
521 real lat2, real lon2,
523 real& s12, real& azi1, real& azi2,
524 real& m12, real& M12, real& M21,
528 real salp1, calp1, salp2, calp2,
529 a12 = GenInverse(lat1, lon1, lat2, lon2,
530 outmask, s12, salp1, calp1, salp2, calp2,
540 real lat2, real lon2,
541 unsigned caps)
const {
542 real t, salp1, calp1, salp2, calp2,
543 a12 = GenInverse(lat1, lon1, lat2, lon2,
545 0u, t, salp1, calp1, salp2, calp2,
556 real ssig1, real csig1, real dn1,
557 real ssig2, real csig2, real dn2,
558 real cbet1, real cbet2,
unsigned outmask,
559 real& s12b, real& m12b, real& m0,
560 real& M12, real& M21)
const {
575 (sig12 + (E.
deltaE(ssig2, csig2, dn2) - E.
deltaE(ssig1, csig1, dn1)));
580 (sig12 + (E.
deltaD(ssig2, csig2, dn2) - E.
deltaD(ssig1, csig1, dn1)));
586 m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
590 real csig12 = csig1 * csig2 + ssig1 * ssig2;
591 real t = _ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2);
592 M12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1;
593 M21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2;
598 Math::real GeodesicExact::Astroid(real x, real y) {
606 if ( !(q == 0 && r <= 0) ) {
615 disc = S * (S + 2 * r3);
622 T3 += T3 < 0 ? -sqrt(disc) : sqrt(disc);
626 u += T + (T ? r2 / T : 0);
629 real ang = atan2(sqrt(-disc), -(S + r3));
632 u += 2 * r * cos(ang / 3);
637 uv = u < 0 ? q / (v - u) : u + v,
638 w = (uv - q) / (2 * v);
641 k = uv / (sqrt(uv +
Math::sq(w)) + w);
651 real sbet1, real cbet1, real dn1,
652 real sbet2, real cbet2, real dn2,
653 real lam12, real slam12, real clam12,
654 real& salp1, real& calp1,
656 real& salp2, real& calp2,
666 sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
667 cbet12 = cbet2 * cbet1 + sbet2 * sbet1;
668 #if defined(__GNUC__) && __GNUC__ == 4 && \ 669 (__GNUC_MINOR__ < 6 || defined(__MINGW32__)) 683 real sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
685 bool shortline = cbet12 >= 0 && sbet12 < real(0.5) &&
686 cbet2 * lam12 < real(0.5);
689 real sbetm2 =
Math::sq(sbet1 + sbet2);
692 sbetm2 /= sbetm2 +
Math::sq(cbet1 + cbet2);
693 dnm = sqrt(1 + _ep2 * sbetm2);
694 real omg12 = lam12 / (_f1 * dnm);
695 somg12 = sin(omg12); comg12 = cos(omg12);
697 somg12 = slam12; comg12 = clam12;
700 salp1 = cbet2 * somg12;
701 calp1 = comg12 >= 0 ?
702 sbet12 + cbet2 * sbet1 *
Math::sq(somg12) / (1 + comg12) :
703 sbet12a - cbet2 * sbet1 *
Math::sq(somg12) / (1 - comg12);
707 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
709 if (shortline && ssig12 < _etol2) {
711 salp2 = cbet1 * somg12;
712 calp2 = sbet12 - cbet1 * sbet2 *
713 (comg12 >= 0 ?
Math::sq(somg12) / (1 + comg12) : 1 - comg12);
716 sig12 = atan2(ssig12, csig12);
717 }
else if (abs(_n) > real(0.1) ||
724 real y, lamscale, betscale;
729 real lam12x = atan2(-slam12, -clam12);
734 E.
Reset(-k2, -_ep2, 1 + k2, 1 + _ep2);
735 lamscale = _e2/_f1 * cbet1 * 2 * E.
H();
737 betscale = lamscale * cbet1;
739 x = lam12x / lamscale;
740 y = sbet12a / betscale;
744 cbet12a = cbet2 * cbet1 - sbet2 * sbet1,
745 bet12a = atan2(sbet12a, cbet12a);
746 real m12b, m0, dummy;
750 sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
752 x = -1 + m12b / (cbet1 * cbet2 * m0 *
Math::pi());
753 betscale = x < -real(0.01) ? sbet12a / x :
755 lamscale = betscale / cbet1;
756 y = lam12x / lamscale;
759 if (y > -tol1_ && x > -1 - xthresh_) {
763 salp1 = min(real(1), -real(x)); calp1 = - sqrt(1 -
Math::sq(salp1));
765 calp1 = max(real(x > -tol1_ ? 0 : -1), real(x));
803 real k = Astroid(x, y);
805 omg12a = lamscale * ( _f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k );
806 somg12 = sin(omg12a); comg12 = -cos(omg12a);
808 salp1 = cbet2 * somg12;
809 calp1 = sbet12a - cbet2 * sbet1 *
Math::sq(somg12) / (1 - comg12);
816 salp1 = 1; calp1 = 0;
821 Math::real GeodesicExact::Lambda12(real sbet1, real cbet1, real dn1,
822 real sbet2, real cbet2, real dn2,
823 real salp1, real calp1,
824 real slam120, real clam120,
825 real& salp2, real& calp2,
827 real& ssig1, real& csig1,
828 real& ssig2, real& csig2,
831 bool diffp, real& dlam12)
const 834 if (sbet1 == 0 && calp1 == 0)
841 salp0 = salp1 * cbet1,
844 real somg1, comg1, somg2, comg2, somg12, comg12, cchi1, cchi2, lam12;
847 ssig1 = sbet1; somg1 = salp0 * sbet1;
848 csig1 = comg1 = calp1 * cbet1;
850 cchi1 = _f1 * dn1 * comg1;
859 salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1;
864 calp2 = cbet2 != cbet1 || abs(sbet2) != -sbet1 ?
867 (cbet2 - cbet1) * (cbet1 + cbet2) :
868 (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :
872 ssig2 = sbet2; somg2 = salp0 * sbet2;
873 csig2 = comg2 = calp2 * cbet2;
875 cchi2 = _f1 * dn2 * comg2;
881 sig12 = atan2(max(real(0), csig1 * ssig2 - ssig1 * csig2),
882 csig1 * csig2 + ssig1 * ssig2);
885 somg12 = max(real(0), comg1 * somg2 - somg1 * comg2);
886 comg12 = comg1 * comg2 + somg1 * somg2;
888 E.
Reset(-k2, -_ep2, 1 + k2, 1 + _ep2);
891 schi12 = max(real(0), cchi1 * somg2 - somg1 * cchi2),
892 cchi12 = cchi1 * cchi2 + somg1 * somg2;
894 real eta = atan2(schi12 * clam120 - cchi12 * slam120,
895 cchi12 * clam120 + schi12 * slam120);
896 real deta12 = -_e2/_f1 * salp0 * E.
H() / (
Math::pi() / 2) *
897 (sig12 + (E.
deltaH(ssig2, csig2, dn2) - E.
deltaH(ssig1, csig1, dn1)));
898 lam12 = eta + deta12;
900 domg12 = deta12 + atan2(schi12 * comg12 - cchi12 * somg12,
901 cchi12 * comg12 + schi12 * somg12);
904 dlam12 = - 2 * _f1 * dn1 / sbet1;
907 Lengths(E, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
909 dummy, dlam12, dummy, dummy, dummy);
910 dlam12 *= _f1 / (calp2 * cbet2);
917 void GeodesicExact::C4f(real eps, real c[])
const {
922 for (
int l = 0; l < nC4_; ++l) {
923 int m = nC4_ - l - 1;
static T AngNormalize(T x)
void Reset(real k2=0, real alpha2=0)
static bool isfinite(T x)
GeodesicLineExact InverseLine(real lat1, real lon1, real lat2, real lon2, unsigned caps=ALL) const
Math::real deltaE(real sn, real cn, real dn) const
Mathematical functions needed by GeographicLib.
static void sincosd(T x, T &sinx, T &cosx)
static T AngDiff(T x, T y, T &e)
Elliptic integrals and functions.
GeodesicLineExact GenDirectLine(real lat1, real lon1, real azi1, bool arcmode, real s12_a12, unsigned caps=ALL) const
static void norm(T &x, T &y)
#define GEOGRAPHICLIB_VOLATILE
GeodesicExact(real a, real f)
Header for GeographicLib::GeodesicLineExact class.
static T atan2d(T y, T x)
static T polyval(int N, const T p[], T x)
Namespace for GeographicLib.
Math::real deltaH(real sn, real cn, real dn) const
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
Exact geodesic calculations.
Header for GeographicLib::GeodesicExact class.
Exception handling for GeographicLib.
GeodesicLineExact DirectLine(real lat1, real lon1, real azi1, real s12, unsigned caps=ALL) const
friend class GeodesicLineExact
Math::real deltaD(real sn, real cn, real dn) const
GeodesicLineExact ArcDirectLine(real lat1, real lon1, real azi1, real a12, unsigned caps=ALL) const
Math::real GenDirect(real lat1, real lon1, real azi1, bool arcmode, real s12_a12, unsigned outmask, real &lat2, real &lon2, real &azi2, real &s12, real &m12, real &M12, real &M21, real &S12) const
GeodesicLineExact Line(real lat1, real lon1, real azi1, unsigned caps=ALL) const
#define GEOGRAPHICLIB_PANIC
static const GeodesicExact & WGS84()