8 package net.sf.geographiclib;
237 protected double _a, _f, _f1, _e2, _ep2, _b, _c2;
238 private double _n, _etol2;
239 private double _A3x[], _C3x[], _C4x[];
261 Math.atan(Math.sqrt(-_e2))) /
262 Math.sqrt(Math.abs(_e2))))/2;
273 _etol2 = 0.1 *
tol2_ /
274 Math.sqrt( Math.max(0.001, Math.abs(_f)) *
275 Math.min(1.0, 1 - _f/2) / 2 );
277 throw new GeographicErr(
"Equatorial radius is not positive");
281 _C3x =
new double[
nC3x_];
282 _C4x =
new double[
nC4x_];
314 double azi1,
double s12) {
338 double azi1,
double s12,
int outmask) {
339 return Direct(lat1, lon1, azi1,
false, s12, outmask);
367 double azi1,
double a12) {
392 double azi1,
double a12,
int outmask) {
393 return Direct(lat1, lon1, azi1,
true, a12, outmask);
453 boolean arcmode,
double s12_a12,
int outmask) {
456 return new GeodesicLine(
this, lat1, lon1, azi1, outmask)
458 Position(arcmode, s12_a12, outmask);
503 double s12,
int caps) {
550 double a12,
int caps) {
578 boolean arcmode,
double s12_a12,
int caps)
584 salp1 =
p.first; calp1 =
p.second; }
587 return new GeodesicLine(
this, lat1, lon1, azi1, salp1, calp1,
588 caps, arcmode, s12_a12);
616 double lat2,
double lon2) {
622 private double salp1, calp1, salp2, calp2;
625 salp1 = calp1 = salp2 = calp2 = Double.NaN;
630 double lat2,
double lon2,
int outmask) {
640 double lon12, lon12s;
643 lon12 =
p.first; lon12s =
p.second;
646 r.
lon1 = lon1; r.
lon2 = (lon1 + lon12) + lon12s;
651 int lonsign = lon12 >= 0 ? 1 : -1;
656 lam12 = Math.toRadians(lon12), slam12, clam12;
658 slam12 =
p.first; clam12 = (lon12 > 90 ? -1 : 1) *
p.second; }
662 int swapp = Math.abs(lat1) < Math.abs(lat2) ? -1 : 1;
665 {
double t = lat1; lat1 = lat2; lat2 =
t; }
668 int latsign = lat1 < 0 ? 1 : -1;
683 double sbet1, cbet1, sbet2, cbet2, s12x, m12x;
684 s12x = m12x = Double.NaN;
687 sbet1 = _f1 *
p.first; cbet1 =
p.second; }
691 cbet1 = Math.max(
tiny_, cbet1);
694 sbet2 = _f1 *
p.first; cbet2 =
p.second; }
697 cbet2 = Math.max(
tiny_, cbet2);
707 if (cbet1 < -sbet1) {
709 sbet2 = sbet2 < 0 ? sbet1 : -sbet1;
711 if (Math.abs(sbet2) == -sbet1)
716 dn1 = Math.sqrt(1 + _ep2 *
GeoMath.
sq(sbet1)),
717 dn2 = Math.sqrt(1 + _ep2 *
GeoMath.
sq(sbet2));
719 double a12, sig12, calp1, salp1, calp2, salp2;
720 a12 = sig12 = calp1 = salp1 = calp2 = salp2 = Double.NaN;
722 double C1a[] =
new double[
nC1_ + 1];
723 double C2a[] =
new double[
nC2_ + 1];
724 double C3a[] =
new double[
nC3_];
726 boolean meridian = lat1 == -90 || slam12 == 0;
733 calp1 = clam12; salp1 = slam12;
734 calp2 = 1; salp2 = 0;
738 ssig1 = sbet1, csig1 = calp1 * cbet1,
739 ssig2 = sbet2, csig2 = calp2 * cbet2;
742 sig12 = Math.atan2(Math.max(0.0, csig1 * ssig2 - ssig1 * csig2),
743 csig1 * csig2 + ssig1 * ssig2);
747 ssig2, csig2, dn2, cbet1, cbet2,
750 s12x =
v.s12b; m12x =
v.m12b;
762 if (sig12 < 1 || m12x >= 0) {
764 if (sig12 < 3 *
tiny_)
765 sig12 = m12x = s12x = 0;
768 a12 = Math.toDegrees(sig12);
774 double omg12 = Double.NaN, somg12 = 2, comg12 = Double.NaN;
778 (_f <= 0 || lon12s >= _f * 180)) {
781 calp1 = calp2 = 0; salp1 = salp2 = 1;
783 sig12 = omg12 = lam12 / _f1;
784 m12x = _b * Math.sin(sig12);
786 r.
M12 = r.
M21 = Math.cos(sig12);
789 }
else if (!meridian) {
799 lam12, slam12, clam12,
802 salp1 =
v.salp1; calp1 =
v.calp1;
803 salp2 =
v.salp2; calp2 =
v.calp2;
809 s12x = sig12 * _b * dnm;
810 m12x =
GeoMath.
sq(dnm) * _b * Math.sin(sig12 / dnm);
812 r.
M12 = r.
M21 = Math.cos(sig12 / dnm);
813 a12 = Math.toDegrees(sig12);
814 omg12 = lam12 / (_f1 * dnm);
828 double ssig1, csig1, ssig2, csig2, eps, domg12;
829 ssig1 = csig1 = ssig2 = csig2 = eps = domg12 = Double.NaN;
832 double salp1a =
tiny_, calp1a = 1, salp1b =
tiny_, calp1b = -1;
833 for (
boolean tripn =
false, tripb =
false; numit <
maxit2_; ++numit) {
839 Lambda12(sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
840 slam12, clam12, numit <
maxit1_, C1a, C2a, C3a);
842 salp2 =
w.salp2; calp2 =
w.calp2;
844 ssig1 =
w.ssig1; csig1 =
w.csig1;
845 ssig2 =
w.ssig2; csig2 =
w.csig2;
846 eps =
w.eps; domg12 =
w.domg12;
851 if (tripb || !(Math.abs(
v) >= (tripn ? 8 : 1) *
tol0_))
break;
853 if (
v > 0 && (numit >
maxit1_ || calp1/salp1 > calp1b/salp1b))
854 { salp1b = salp1; calp1b = calp1; }
855 else if (
v < 0 && (numit >
maxit1_ || calp1/salp1 < calp1a/salp1a))
856 { salp1a = salp1; calp1a = calp1; }
857 if (numit < maxit1_ && dv > 0) {
861 sdalp1 = Math.sin(dalp1), cdalp1 = Math.cos(dalp1),
862 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1;
863 if (nsalp1 > 0 && Math.abs(dalp1) < Math.PI) {
864 calp1 = calp1 * cdalp1 - salp1 * sdalp1;
867 salp1 =
p.first; calp1 =
p.second; }
871 tripn = Math.abs(
v) <= 16 *
tol0_;
883 salp1 = (salp1a + salp1b)/2;
884 calp1 = (calp1a + calp1b)/2;
886 salp1 =
p.first; calp1 =
p.second; }
888 tripb = (Math.abs(salp1a - salp1) + (calp1a - calp1) <
tolb_ ||
889 Math.abs(salp1 - salp1b) + (calp1 - calp1b) <
tolb_);
894 int lengthmask = outmask |
900 ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2,
901 lengthmask, C1a, C2a);
902 s12x =
v.s12b; m12x =
v.m12b;
909 a12 = Math.toDegrees(sig12);
912 double sdomg12 = Math.sin(domg12), cdomg12 = Math.cos(domg12);
913 somg12 = slam12 * cdomg12 - clam12 * sdomg12;
914 comg12 = clam12 * cdomg12 + slam12 * sdomg12;
928 salp0 = salp1 * cbet1,
931 if (calp0 != 0 && salp0 != 0) {
934 ssig1 = sbet1, csig1 = calp1 * cbet1,
935 ssig2 = sbet2, csig2 = calp2 * cbet2,
937 eps = k2 / (2 * (1 + Math.sqrt(1 + k2)) + k2),
941 ssig1 =
p.first; csig1 =
p.second; }
943 ssig2 =
p.first; csig2 =
p.second; }
944 double C4a[] =
new double[
nC4_];
949 r.
S12 =
A4 * (B42 - B41);
954 if (!meridian && somg12 > 1) {
955 somg12 = Math.sin(omg12); comg12 = Math.cos(omg12);
960 sbet2 - sbet1 < 1.75) {
965 domg12 = 1 + comg12, dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
966 alp12 = 2 * Math.atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
967 domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );
971 salp12 = salp2 * calp1 - calp2 * salp1,
972 calp12 = calp2 * calp1 + salp2 * salp1;
977 if (salp12 == 0 && calp12 < 0) {
978 salp12 =
tiny_ * calp1;
981 alp12 = Math.atan2(salp12, calp12);
983 r.
S12 += _c2 * alp12;
984 r.
S12 *= swapp * lonsign * latsign;
991 {
double t = salp1; salp1 = salp2; salp2 =
t; }
992 {
double t = calp1; calp1 = calp2; calp2 =
t; }
997 salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
998 salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
1050 double lat2,
double lon2,
int outmask) {
1078 double lat2,
double lon2) {
1103 double lat2,
double lon2,
int caps) {
1110 return new GeodesicLine(
this, lat1, lon1, azi1, salp1, calp1, caps,
1230 double sinx,
double cosx,
1239 n = k - (sinp ? 1 : 0);
1241 ar = 2 * (cosx - sinx) * (cosx + sinx),
1242 y0 = (
n & 1) != 0 ?
c[--k] : 0,
y1 = 0;
1251 ? 2 * sinx * cosx *
y0
1256 private double s12b, m12b, m0, M12, M21;
1258 s12b = m12b = m0 = M12 = M21 = Double.NaN;
1263 double ssig1,
double csig1,
double dn1,
1264 double ssig2,
double csig2,
double dn2,
1265 double cbet1,
double cbet2,
1268 double C1a[],
double C2a[]) {
1274 double m0x = 0, J12 = 0,
A1 = 0,
A2 = 0;
1292 v.s12b =
A1 * (sig12 + B1);
1297 J12 = m0x * sig12 + (
A1 * B1 -
A2 * B2);
1302 for (
int l = 1;
l <=
nC2_; ++
l)
1303 C2a[
l] =
A1 * C1a[
l] -
A2 * C2a[
l];
1304 J12 = m0x * sig12 + (
SinCosSeries(
true, ssig2, csig2, C2a) -
1312 v.m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
1313 csig1 * csig2 * J12;
1316 double csig12 = csig1 * csig2 + ssig1 * ssig2;
1317 double t = _ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2);
1318 v.M12 = csig12 + (
t * ssig2 - csig2 * J12) * ssig1 / dn1;
1319 v.M21 = csig12 - (
t * ssig1 - csig1 * J12) * ssig2 / dn2;
1331 r = (
p +
q - 1) / 6;
1332 if ( !(
q == 0 && r <= 0) ) {
1341 disc =
S * (
S + 2 *
r3);
1348 T3 +=
T3 < 0 ? -Math.sqrt(disc) : Math.sqrt(disc);
1352 u +=
T + (
T != 0 ?
r2 /
T : 0);
1355 double ang = Math.atan2(Math.sqrt(-disc), -(
S +
r3));
1358 u += 2 * r * Math.cos(ang / 3);
1363 uv = u < 0 ?
q / (
v - u) : u +
v,
1364 w = (uv -
q) / (2 *
v);
1383 sig12 = salp1 = calp1 = salp2 = calp2 = dnm = Double.NaN;
1388 double sbet2,
double cbet2,
double dn2,
1390 double slam12,
double clam12,
1392 double C1a[],
double C2a[]) {
1402 sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
1403 cbet12 = cbet2 * cbet1 + sbet2 * sbet1;
1404 double sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
1405 boolean shortline = cbet12 >= 0 && sbet12 < 0.5 &&
1406 cbet2 * lam12 < 0.5;
1407 double somg12, comg12;
1409 double sbetm2 =
GeoMath.
sq(sbet1 + sbet2);
1412 sbetm2 /= sbetm2 +
GeoMath.
sq(cbet1 + cbet2);
1413 w.dnm = Math.sqrt(1 + _ep2 * sbetm2);
1414 double omg12 = lam12 / (_f1 *
w.dnm);
1415 somg12 = Math.sin(omg12); comg12 = Math.cos(omg12);
1417 somg12 = slam12; comg12 = clam12;
1420 w.salp1 = cbet2 * somg12;
1421 w.calp1 = comg12 >= 0 ?
1422 sbet12 + cbet2 * sbet1 *
GeoMath.
sq(somg12) / (1 + comg12) :
1423 sbet12a - cbet2 * sbet1 *
GeoMath.
sq(somg12) / (1 - comg12);
1427 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
1429 if (shortline && ssig12 < _etol2) {
1431 w.salp2 = cbet1 * somg12;
1432 w.calp2 = sbet12 - cbet1 * sbet2 *
1433 (comg12 >= 0 ?
GeoMath.
sq(somg12) / (1 + comg12) : 1 - comg12);
1435 w.salp2 =
p.first;
w.calp2 =
p.second; }
1437 w.sig12 = Math.atan2(ssig12, csig12);
1438 }
else if (Math.abs(
_n) > 0.1 ||
1440 ssig12 >= 6 * Math.abs(
_n) * Math.PI *
GeoMath.
sq(cbet1)) {
1445 double y, lamscale, betscale;
1450 double lam12x = Math.atan2(-slam12, -clam12);
1456 eps = k2 / (2 * (1 + Math.sqrt(1 + k2)) + k2);
1457 lamscale = _f * cbet1 *
A3f(eps) * Math.PI;
1459 betscale = lamscale * cbet1;
1461 x = lam12x / lamscale;
1462 y = sbet12a / betscale;
1466 cbet12a = cbet2 * cbet1 - sbet2 * sbet1,
1467 bet12a = Math.atan2(sbet12a, cbet12a);
1473 sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
1475 m12b =
v.m12b;
m0 =
v.m0;
1477 x = -1 + m12b / (cbet1 * cbet2 *
m0 * Math.PI);
1478 betscale =
x < -0.01 ? sbet12a /
x :
1480 lamscale = betscale / cbet1;
1481 y = lam12x / lamscale;
1487 w.salp1 = Math.min(1.0, -
x);
1490 w.calp1 = Math.max(
x > -
tol1_ ? 0.0 : -1.0,
x);
1530 omg12a = lamscale * ( _f >= 0 ? -
x * k/(1 + k) : -
y * (1 + k)/k );
1531 somg12 = Math.sin(omg12a); comg12 = -Math.cos(omg12a);
1533 w.salp1 = cbet2 * somg12;
1534 w.calp1 = sbet12a - cbet2 * sbet1 *
GeoMath.
sq(somg12) / (1 - comg12);
1538 if (!(
w.salp1 <= 0))
1540 w.salp1 =
p.first;
w.calp1 =
p.second; }
1542 w.salp1 = 1;
w.calp1 = 0;
1548 private double lam12, salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
1549 eps, domg12, dlam12;
1551 lam12 = salp2 = calp2 = sig12 = ssig1 = csig1 = ssig2 = csig2
1552 = eps = domg12 = dlam12 = Double.NaN;
1557 double sbet2,
double cbet2,
double dn2,
1558 double salp1,
double calp1,
1559 double slam120,
double clam120,
1562 double C1a[],
double C2a[],
double C3a[]) {
1568 if (sbet1 == 0 && calp1 == 0)
1575 salp0 = salp1 * cbet1,
1578 double somg1, comg1, somg2, comg2, somg12, comg12;
1581 w.ssig1 = sbet1; somg1 = salp0 * sbet1;
1582 w.csig1 = comg1 = calp1 * cbet1;
1584 w.ssig1 =
p.first;
w.csig1 =
p.second; }
1591 w.salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1;
1596 w.calp2 = cbet2 != cbet1 || Math.abs(sbet2) != -sbet1 ?
1599 (cbet2 - cbet1) * (cbet1 + cbet2) :
1600 (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :
1604 w.ssig2 = sbet2; somg2 = salp0 * sbet2;
1605 w.csig2 = comg2 =
w.calp2 * cbet2;
1607 w.ssig2 =
p.first;
w.csig2 =
p.second; }
1611 w.sig12 = Math.atan2(Math.max(0.0,
w.csig1 *
w.ssig2 -
w.ssig1 *
w.csig2),
1612 w.csig1 *
w.csig2 +
w.ssig1 *
w.ssig2);
1615 somg12 = Math.max(0.0, comg1 * somg2 - somg1 * comg2);
1616 comg12 = comg1 * comg2 + somg1 * somg2;
1618 double eta = Math.atan2(somg12 * clam120 - comg12 * slam120,
1619 comg12 * clam120 + somg12 * slam120);
1622 w.eps = k2 / (2 * (1 + Math.sqrt(1 + k2)) + k2);
1626 w.domg12 = -_f *
A3f(
w.eps) * salp0 * (
w.sig12 + B312);
1627 w.lam12 = eta +
w.domg12;
1631 w.dlam12 = - 2 * _f1 * dn1 / sbet1;
1634 Lengths(
w.eps,
w.sig12,
w.ssig1,
w.csig1, dn1,
w.ssig2,
w.csig2, dn2,
1637 w.dlam12 *= _f1 / (
w.calp2 * cbet2);
1644 protected double A3f(
double eps) {
1649 protected void C3f(
double eps,
double c[]) {
1654 for (
int l = 1;
l <
nC3_; ++
l) {
1662 protected void C4f(
double eps,
double c[]) {
1667 for (
int l = 0;
l <
nC4_; ++
l) {
1676 protected static double A1m1f(
double eps) {
1677 final double coeff[] = {
1683 return (
t + eps) / (1 - eps);
1687 protected static void C1f(
double eps,
double c[]) {
1688 final double coeff[] = {
1706 for (
int l = 1;
l <=
nC1_; ++
l) {
1715 protected static void C1pf(
double eps,
double c[]) {
1716 final double coeff[] = {
1718 205, -432, 768, 1536,
1720 4005, -4736, 3840, 12288,
1743 protected static double A2m1f(
double eps) {
1744 final double coeff[] = {
1746 -11, -28, -192, 0, 256,
1750 return (
t - eps) / (1 + eps);
1754 protected static void C2f(
double eps,
double c[]) {
1755 final double coeff[] = {
1773 for (
int l = 1;
l <=
nC2_; ++
l) {
1783 final double coeff[] = {
1798 for (
int j =
nA3_ - 1;
j >= 0; --
j) {
1799 int m = Math.min(
nA3_ -
j - 1,
j);
1807 final double coeff[] = {
1840 for (
int l = 1;
l <
nC3_; ++
l) {
1841 for (
int j =
nC3_ - 1;
j >=
l; --
j) {
1842 int m = Math.min(
nC3_ -
j - 1,
j);
1850 final double coeff[] = {
1856 -224, -4784, 1573, 45045,
1858 -10656, 14144, -4576, -858, 45045,
1860 64, 624, -4576, 6864, -3003, 15015,
1862 100, 208, 572, 3432, -12012, 30030, 45045,
1868 5792, 1040, -1287, 135135,
1870 5952, -11648, 9152, -2574, 135135,
1872 -64, -624, 4576, -6864, 3003, 135135,
1878 -8448, 4992, -1144, 225225,
1880 -1440, 4160, -4576, 1716, 225225,
1886 3584, -3328, 1144, 315315,
1895 for (
int l = 0;
l <
nC4_; ++
l) {
1896 for (
int j =
nC4_ - 1;
j >=
l; --
j) {