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())
55 , tolb_(tol0_ * tol2_)
56 , xthresh_(1000 * tol2_)
58 , _f(f <= 1 ? f : 1/f)
61 , _ep2(_e2 /
Math::sq(_f1))
66 (_e2 > 0 ?
Math::atanh(sqrt(_e2)) : atan(sqrt(-_e2))) /
78 , _etol2(0.1 * tol2_ /
79 sqrt( max(real(0.001), abs(_f)) * min(real(1), 1 - _f/2) / 2 ))
95 const real c[],
int n) {
102 ar = 2 * (cosx - sinx) * (cosx + sinx),
103 y0 = n & 1 ? *--c : 0, y1 = 0;
108 y1 = ar * y0 - y1 + *--c;
109 y0 = ar * y1 - y0 + *--c;
111 return cosx * (y0 - y1);
115 unsigned caps)
const {
120 bool arcmode, real s12_a12,
122 real& lat2, real& lon2, real& azi2,
123 real& s12, real& m12,
124 real& M12, real& M21,
130 GenPosition(arcmode, s12_a12, outmask,
131 lat2, lon2, azi2, s12, m12, M12, M21, S12);
135 real lat2, real lon2,
137 real& s12, real& azi1, real& azi2,
138 real& m12, real& M12, real& M21,
147 lon12 = AngRound(lon12);
149 int lonsign = lon12 >= 0 ? 1 : -1;
152 lat1 = AngRound(lat1);
153 lat2 = AngRound(lat2);
155 int swapp = abs(lat1) >= abs(lat2) ? 1 : -1;
161 int latsign = lat1 < 0 ? 1 : -1;
176 real phi, sbet1, cbet1, sbet2, cbet2, s12x, m12x;
183 sbet1 = _f1 * sin(phi);
184 cbet1 = lat1 == -90 ? tiny_ : cos(phi);
185 SinCosNorm(sbet1, cbet1);
189 sbet2 = _f1 * sin(phi);
190 cbet2 = abs(lat2) == 90 ? tiny_ : cos(phi);
191 SinCosNorm(sbet2, cbet2);
201 if (cbet1 < -sbet1) {
203 sbet2 = sbet2 < 0 ? sbet1 : -sbet1;
205 if (abs(sbet2) == -sbet1)
210 dn1 = (_f >= 0 ? sqrt(1 + _ep2 *
Math::sq(sbet1)) :
211 sqrt(1 - _e2 *
Math::sq(cbet1)) / _f1),
212 dn2 = (_f >= 0 ? sqrt(1 + _ep2 *
Math::sq(sbet2)) :
213 sqrt(1 - _e2 *
Math::sq(cbet2)) / _f1);
217 slam12 = abs(lon12) == 180 ? 0 : sin(lam12),
221 real a12, sig12, calp1, salp1, calp2 = 0, salp2 = 0;
223 bool meridian = lat1 == -90 || slam12 == 0;
230 calp1 = clam12; salp1 = slam12;
231 calp2 = 1; salp2 = 0;
235 ssig1 = sbet1, csig1 = calp1 * cbet1,
236 ssig2 = sbet2, csig2 = calp2 * cbet2;
239 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, real(0)),
240 csig1 * csig2 + ssig1 * ssig2);
243 Lengths(E, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
244 cbet1, cbet2, s12x, m12x, dummy,
254 if (sig12 < 1 || m12x >= 0) {
270 calp1 = calp2 = 0; salp1 = salp2 = 1;
272 sig12 = omg12 = lam12 / _f1;
273 m12x = _b * sin(sig12);
275 M12 = M21 = cos(sig12);
278 }
else if (!meridian) {
285 sig12 = InverseStart(E, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
287 salp1, calp1, salp2, calp2, dnm);
291 s12x = sig12 * _b * dnm;
292 m12x =
Math::sq(dnm) * _b * sin(sig12 / dnm);
294 M12 = M21 = cos(sig12 / dnm);
296 omg12 = lam12 / (_f1 * dnm);
312 real ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0;
315 real salp1a = tiny_, calp1a = 1, salp1b = tiny_, calp1b = -1;
316 for (
bool tripn =
false, tripb =
false;
341 real v = Lambda12(sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
342 salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
343 E, omg12, numit < maxit1_, dv) - lam12;
346 if (tripb || !(abs(v) >= (tripn ? 8 : 2) * tol0_))
break;
348 if (v > 0 && (numit > maxit1_ || calp1/salp1 > calp1b/salp1b))
349 { salp1b = salp1; calp1b = calp1; }
350 else if (v < 0 && (numit > maxit1_ || calp1/salp1 < calp1a/salp1a))
351 { salp1a = salp1; calp1a = calp1; }
352 if (numit < maxit1_ && dv > 0) {
356 sdalp1 = sin(dalp1), cdalp1 = cos(dalp1),
357 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1;
358 if (nsalp1 > 0 && abs(dalp1) <
Math::pi()) {
359 calp1 = calp1 * cdalp1 - salp1 * sdalp1;
361 SinCosNorm(salp1, calp1);
365 tripn = abs(v) <= 16 * tol0_;
377 salp1 = (salp1a + salp1b)/2;
378 calp1 = (calp1a + calp1b)/2;
379 SinCosNorm(salp1, calp1);
381 tripb = (abs(salp1a - salp1) + (calp1a - calp1) < tolb_ ||
382 abs(salp1 - salp1b) + (calp1 - calp1b) < tolb_);
386 Lengths(E, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
387 cbet1, cbet2, s12x, m12x, dummy,
402 if (outmask &
AREA) {
405 salp0 = salp1 * cbet1,
408 if (calp0 != 0 && salp0 != 0) {
411 ssig1 = sbet1, csig1 = calp1 * cbet1,
412 ssig2 = sbet2, csig2 = calp2 * cbet2,
414 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2),
416 A4 =
Math::sq(_a) * calp0 * salp0 * _e2;
417 SinCosNorm(ssig1, csig1);
418 SinCosNorm(ssig2, csig2);
422 B41 = CosSeries(ssig1, csig1, C4a, nC4_),
423 B42 = CosSeries(ssig2, csig2, C4a, nC4_);
424 S12 = A4 * (B42 - B41);
431 sbet2 - sbet1 < real(1.75)) {
436 somg12 = sin(omg12), domg12 = 1 + cos(omg12),
437 dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
438 alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
439 domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );
443 salp12 = salp2 * calp1 - calp2 * salp1,
444 calp12 = calp2 * calp1 + salp2 * salp1;
449 if (salp12 == 0 && calp12 < 0) {
450 salp12 = tiny_ * calp1;
453 alp12 = atan2(salp12, calp12);
456 S12 *= swapp * lonsign * latsign;
469 salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
470 salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
488 bool scalep,
real& M12,
real& M21)
const {
496 (sig12 + E.
deltaD(ssig2, csig2, dn2) - E.
deltaD(ssig1, csig1, dn1));
500 m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) - csig1 * csig2 * J12;
503 (sig12 + E.
deltaE(ssig2, csig2, dn2) - E.
deltaE(ssig1, csig1, dn1));
505 real csig12 = csig1 * csig2 + ssig1 * ssig2;
506 real t = _ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2);
507 M12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1;
508 M21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2;
520 if ( !(q == 0 && r <= 0) ) {
529 disc = S * (S + 2 * r3);
536 T3 += T3 < 0 ? -sqrt(disc) : sqrt(disc);
540 u += T + (T ? r2 / T : 0);
543 real ang = atan2(sqrt(-disc), -(S + r3));
546 u += 2 * r * cos(ang / 3);
551 uv = u < 0 ? q / (v - u) : u + v,
552 w = (uv - q) / (2 * v);
555 k = uv / (sqrt(uv +
Math::sq(w)) + w);
564 Math::real GeodesicExact::InverseStart(EllipticFunction& E,
580 sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
581 cbet12 = cbet2 * cbet1 + sbet2 * sbet1;
582 #if defined(__GNUC__) && __GNUC__ == 4 && \
583 (__GNUC_MINOR__ < 6 || defined(__MINGW32__))
597 real sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
599 bool shortline = cbet12 >= 0 && sbet12 <
real(0.5) &&
600 cbet2 * lam12 <
real(0.5);
606 sbetm2 /= sbetm2 +
Math::sq(cbet1 + cbet2);
607 dnm = sqrt(1 + _ep2 * sbetm2);
610 real somg12 = sin(omg12), comg12 = cos(omg12);
612 salp1 = cbet2 * somg12;
613 calp1 = comg12 >= 0 ?
614 sbet12 + cbet2 * sbet1 *
Math::sq(somg12) / (1 + comg12) :
615 sbet12a - cbet2 * sbet1 *
Math::sq(somg12) / (1 - comg12);
619 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
621 if (shortline && ssig12 < _etol2) {
623 salp2 = cbet1 * somg12;
624 calp2 = sbet12 - cbet1 * sbet2 *
625 (comg12 >= 0 ?
Math::sq(somg12) / (1 + comg12) : 1 - comg12);
626 SinCosNorm(salp2, calp2);
628 sig12 = atan2(ssig12, csig12);
629 }
else if (abs(_n) >
real(0.1) ||
636 real y, lamscale, betscale;
645 E.Reset(-k2, -_ep2, 1 + k2, 1 + _ep2);
646 lamscale = _e2/_f1 * cbet1 * 2 * E.H();
648 betscale = lamscale * cbet1;
650 x = (lam12 -
Math::pi()) / lamscale;
651 y = sbet12a / betscale;
655 cbet12a = cbet2 * cbet1 - sbet2 * sbet1,
656 bet12a = atan2(sbet12a, cbet12a);
657 real m12b, m0, dummy;
661 sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
662 cbet1, cbet2, dummy, m12b, m0,
false,
664 x = -1 + m12b / (cbet1 * cbet2 * m0 *
Math::pi());
665 betscale = x < -
real(0.01) ? sbet12a / x :
667 lamscale = betscale / cbet1;
668 y = (lam12 -
Math::pi()) / lamscale;
671 if (y > -tol1_ && x > -1 - xthresh_) {
677 calp1 = max(
real(x > -tol1_ ? 0 : -1),
real(x));
715 real k = Astroid(x, y);
717 omg12a = lamscale * ( _f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k );
718 somg12 = sin(omg12a); comg12 = -cos(omg12a);
720 salp1 = cbet2 * somg12;
721 calp1 = sbet12a - cbet2 * sbet1 *
Math::sq(somg12) / (1 - comg12);
726 SinCosNorm(salp1, calp1);
728 salp1 = 1; calp1 = 0;
742 bool diffp,
real& dlam12)
const
745 if (sbet1 == 0 && calp1 == 0)
752 salp0 = salp1 * cbet1,
755 real somg1, comg1, somg2, comg2, cchi1, cchi2, lam12;
758 ssig1 = sbet1; somg1 = salp0 * sbet1;
759 csig1 = comg1 = calp1 * cbet1;
761 cchi1 = _f1 * dn1 * comg1;
762 SinCosNorm(ssig1, csig1);
770 salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1;
775 calp2 = cbet2 != cbet1 || abs(sbet2) != -sbet1 ?
778 (cbet2 - cbet1) * (cbet1 + cbet2) :
779 (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :
783 ssig2 = sbet2; somg2 = salp0 * sbet2;
784 csig2 = comg2 = calp2 * cbet2;
786 cchi2 = _f1 * dn2 * comg2;
787 SinCosNorm(ssig2, csig2);
792 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2,
real(0)),
793 csig1 * csig2 + ssig1 * ssig2);
796 omg12 = atan2(max(comg1 * somg2 - somg1 * comg2,
real(0)),
797 comg1 * comg2 + somg1 * somg2);
799 E.Reset(-k2, -_ep2, 1 + k2, 1 + _ep2);
800 real chi12 = atan2(max(cchi1 * somg2 - somg1 * cchi2,
real(0)),
801 cchi1 * cchi2 + somg1 * somg2);
803 _e2/_f1 * salp0 * E.H() / (
Math::pi() / 2) *
804 (sig12 + E.deltaH(ssig2, csig2, dn2) - E.deltaH(ssig1, csig1, dn1) );
808 dlam12 = - 2 * _f1 * dn1 / sbet1;
811 Lengths(E, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
812 cbet1, cbet2, dummy, dlam12, dummy,
813 false, dummy, dummy);
814 dlam12 *= _f1 / (calp2 * cbet2);
821 void GeodesicExact::C4f(
real eps,
real c[])
const {
824 for (
int j = nC4x_, k = nC4_; k; ) {
826 for (
int i = nC4_ - k + 1; i; --i)
827 t = eps * t + _C4x[--j];
832 for (
int k = 1; k < nC4_; ) {
843 void GeodesicExact::C4coeff() {
844 const real* cc = rawC4coeff();
846 for (
int m = 0, k = 0, h = 0; m < nC4_; ++m) {
848 for (
int j = m; j < nC4_; ++j) {
851 for (
int l = nC4_ - j; l--;)
852 t = _n * t + cc[h++];
853 _C4x[k++] = t/cc[h++];
static T AngNormalize(T x)
GeographicLib::Math::real real
static bool isfinite(T x)
Mathematical functions needed by GeographicLib.
Elliptic integrals and functions.
#define GEOGRAPHICLIB_VOLATILE
Math::real GenInverse(real lat1, real lon1, real lat2, real lon2, unsigned outmask, real &s12, real &azi1, real &azi2, real &m12, real &M12, real &M21, real &S12) const
GeodesicExact(real a, real f)
GeodesicLineExact Line(real lat1, real lon1, real azi1, unsigned caps=ALL) const
Header for GeographicLib::GeodesicLineExact class.
Namespace for GeographicLib.
static T AngDiff(T x, T y)
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
Exact geodesic calculations.
Math::real deltaE(real sn, real cn, real dn) const
Header for GeographicLib::GeodesicExact class.
Exception handling for GeographicLib.
friend class GeodesicLineExact
Math::real deltaD(real sn, real cn, real dn) const
#define GEOGRAPHICLIB_PANIC
static const GeodesicExact & WGS84()