Geodesic.java
Go to the documentation of this file.
1 
8 package net.sf.geographiclib;
9 
202 public class Geodesic {
203 
207  protected static final int GEOGRAPHICLIB_GEODESIC_ORDER = 6;
208 
209  protected static final int nA1_ = GEOGRAPHICLIB_GEODESIC_ORDER;
210  protected static final int nC1_ = GEOGRAPHICLIB_GEODESIC_ORDER;
211  protected static final int nC1p_ = GEOGRAPHICLIB_GEODESIC_ORDER;
212  protected static final int nA2_ = GEOGRAPHICLIB_GEODESIC_ORDER;
213  protected static final int nC2_ = GEOGRAPHICLIB_GEODESIC_ORDER;
214  protected static final int nA3_ = GEOGRAPHICLIB_GEODESIC_ORDER;
215  protected static final int nA3x_ = nA3_;
216  protected static final int nC3_ = GEOGRAPHICLIB_GEODESIC_ORDER;
217  protected static final int nC3x_ = (nC3_ * (nC3_ - 1)) / 2;
218  protected static final int nC4_ = GEOGRAPHICLIB_GEODESIC_ORDER;
219  protected static final int nC4x_ = (nC4_ * (nC4_ + 1)) / 2;
220  private static final int maxit1_ = 20;
221  private static final int maxit2_ = maxit1_ + GeoMath.digits + 10;
222 
223  // Underflow guard. We require
224  // tiny_ * epsilon() > 0
225  // tiny_ + epsilon() == epsilon()
226  protected static final double tiny_ = Math.sqrt(GeoMath.min);
227  private static final double tol0_ = GeoMath.epsilon;
228  // Increase multiplier in defn of tol1_ from 100 to 200 to fix inverse case
229  // 52.784459512564 0 -52.784459512563990912 179.634407464943777557
230  // which otherwise failed for Visual Studio 10 (Release and Debug)
231  private static final double tol1_ = 200 * tol0_;
232  private static final double tol2_ = Math.sqrt(tol0_);
233  // Check on bisection interval
234  private static final double tolb_ = tol0_ * tol2_;
235  private static final double xthresh_ = 1000 * tol2_;
236 
237  protected double _a, _f, _f1, _e2, _ep2, _b, _c2;
238  private double _n, _etol2;
239  private double _A3x[], _C3x[], _C4x[];
240 
250  public Geodesic(double a, double f) {
251  _a = a;
252  _f = f;
253  _f1 = 1 - _f;
254  _e2 = _f * (2 - _f);
255  _ep2 = _e2 / GeoMath.sq(_f1); // e2 / (1 - e2)
256  _n = _f / ( 2 - _f);
257  _b = _a * _f1;
258  _c2 = (GeoMath.sq(_a) + GeoMath.sq(_b) *
259  (_e2 == 0 ? 1 :
260  (_e2 > 0 ? GeoMath.atanh(Math.sqrt(_e2)) :
261  Math.atan(Math.sqrt(-_e2))) /
262  Math.sqrt(Math.abs(_e2))))/2; // authalic radius squared
263  // The sig12 threshold for "really short". Using the auxiliary sphere
264  // solution with dnm computed at (bet1 + bet2) / 2, the relative error in
265  // the azimuth consistency check is sig12^2 * abs(f) * min(1, 1-f/2) / 2.
266  // (Error measured for 1/100 < b/a < 100 and abs(f) >= 1/1000. For a
267  // given f and sig12, the max error occurs for lines near the pole. If
268  // the old rule for computing dnm = (dn1 + dn2)/2 is used, then the error
269  // increases by a factor of 2.) Setting this equal to epsilon gives
270  // sig12 = etol2. Here 0.1 is a safety factor (error decreased by 100)
271  // and max(0.001, abs(f)) stops etol2 getting too large in the nearly
272  // spherical case.
273  _etol2 = 0.1 * tol2_ /
274  Math.sqrt( Math.max(0.001, Math.abs(_f)) *
275  Math.min(1.0, 1 - _f/2) / 2 );
276  if (!(GeoMath.isfinite(_a) && _a > 0))
277  throw new GeographicErr("Equatorial radius is not positive");
278  if (!(GeoMath.isfinite(_b) && _b > 0))
279  throw new GeographicErr("Polar semi-axis is not positive");
280  _A3x = new double[nA3x_];
281  _C3x = new double[nC3x_];
282  _C4x = new double[nC4x_];
283 
284  A3coeff();
285  C3coeff();
286  C4coeff();
287  }
288 
313  public GeodesicData Direct(double lat1, double lon1,
314  double azi1, double s12) {
315  return Direct(lat1, lon1, azi1, false, s12, GeodesicMask.STANDARD);
316  }
337  public GeodesicData Direct(double lat1, double lon1,
338  double azi1, double s12, int outmask) {
339  return Direct(lat1, lon1, azi1, false, s12, outmask);
340  }
341 
366  public GeodesicData ArcDirect(double lat1, double lon1,
367  double azi1, double a12) {
368  return Direct(lat1, lon1, azi1, true, a12, GeodesicMask.STANDARD);
369  }
370 
391  public GeodesicData ArcDirect(double lat1, double lon1,
392  double azi1, double a12, int outmask) {
393  return Direct(lat1, lon1, azi1, true, a12, outmask);
394  }
395 
452  public GeodesicData Direct(double lat1, double lon1, double azi1,
453  boolean arcmode, double s12_a12, int outmask) {
454  // Automatically supply DISTANCE_IN if necessary
455  if (!arcmode) outmask |= GeodesicMask.DISTANCE_IN;
456  return new GeodesicLine(this, lat1, lon1, azi1, outmask)
457  . // Note the dot!
458  Position(arcmode, s12_a12, outmask);
459  }
460 
477  public GeodesicLine DirectLine(double lat1, double lon1, double azi1,
478  double s12) {
479  return DirectLine(lat1, lon1, azi1, s12, GeodesicMask.ALL);
480  }
481 
502  public GeodesicLine DirectLine(double lat1, double lon1, double azi1,
503  double s12, int caps) {
504  return GenDirectLine(lat1, lon1, azi1, false, s12, caps);
505  }
506 
523  public GeodesicLine ArcDirectLine(double lat1, double lon1, double azi1,
524  double a12) {
525  return ArcDirectLine(lat1, lon1, azi1, a12, GeodesicMask.ALL);
526  }
527 
549  public GeodesicLine ArcDirectLine(double lat1, double lon1, double azi1,
550  double a12, int caps) {
551  return GenDirectLine(lat1, lon1, azi1, true, a12, caps);
552  }
553 
577  public GeodesicLine GenDirectLine(double lat1, double lon1, double azi1,
578  boolean arcmode, double s12_a12, int caps)
579  {
580  azi1 = GeoMath.AngNormalize(azi1);
581  double salp1, calp1;
582  // Guard against underflow in salp0. Also -0 is converted to +0.
583  { Pair p = GeoMath.sincosd(GeoMath.AngRound(azi1));
584  salp1 = p.first; calp1 = p.second; }
585  // Automatically supply DISTANCE_IN if necessary
586  if (!arcmode) caps |= GeodesicMask.DISTANCE_IN;
587  return new GeodesicLine(this, lat1, lon1, azi1, salp1, calp1,
588  caps, arcmode, s12_a12);
589  }
590 
615  public GeodesicData Inverse(double lat1, double lon1,
616  double lat2, double lon2) {
617  return Inverse(lat1, lon1, lat2, lon2, GeodesicMask.STANDARD);
618  }
619 
620  private class InverseData {
621  private GeodesicData g;
622  private double salp1, calp1, salp2, calp2;
623  private InverseData() {
624  g = new GeodesicData();
625  salp1 = calp1 = salp2 = calp2 = Double.NaN;
626  }
627  }
628 
629  private InverseData InverseInt(double lat1, double lon1,
630  double lat2, double lon2, int outmask) {
632  GeodesicData r = result.g;
633  // Compute longitude difference (AngDiff does this carefully). Result is
634  // in [-180, 180] but -180 is only for west-going geodesics. 180 is for
635  // east-going and meridional geodesics.
636  r.lat1 = lat1 = GeoMath.LatFix(lat1); r.lat2 = lat2 = GeoMath.LatFix(lat2);
637  // If really close to the equator, treat as on equator.
638  lat1 = GeoMath.AngRound(lat1);
639  lat2 = GeoMath.AngRound(lat2);
640  double lon12, lon12s;
641  {
642  Pair p = GeoMath.AngDiff(lon1, lon2);
643  lon12 = p.first; lon12s = p.second;
644  }
645  if ((outmask & GeodesicMask.LONG_UNROLL) != 0) {
646  r.lon1 = lon1; r.lon2 = (lon1 + lon12) + lon12s;
647  } else {
648  r.lon1 = GeoMath.AngNormalize(lon1); r.lon2 = GeoMath.AngNormalize(lon2);
649  }
650  // Make longitude difference positive.
651  int lonsign = lon12 >= 0 ? 1 : -1;
652  // If very close to being on the same half-meridian, then make it so.
653  lon12 = lonsign * GeoMath.AngRound(lon12);
654  lon12s = GeoMath.AngRound((180 - lon12) - lonsign * lon12s);
655  double
656  lam12 = Math.toRadians(lon12), slam12, clam12;
657  { Pair p = GeoMath.sincosd(lon12 > 90 ? lon12s : lon12);
658  slam12 = p.first; clam12 = (lon12 > 90 ? -1 : 1) * p.second; }
659 
660  // Swap points so that point with higher (abs) latitude is point 1
661  // If one latitude is a nan, then it becomes lat1.
662  int swapp = Math.abs(lat1) < Math.abs(lat2) ? -1 : 1;
663  if (swapp < 0) {
664  lonsign *= -1;
665  { double t = lat1; lat1 = lat2; lat2 = t; }
666  }
667  // Make lat1 <= 0
668  int latsign = lat1 < 0 ? 1 : -1;
669  lat1 *= latsign;
670  lat2 *= latsign;
671  // Now we have
672  //
673  // 0 <= lon12 <= 180
674  // -90 <= lat1 <= 0
675  // lat1 <= lat2 <= -lat1
676  //
677  // longsign, swapp, latsign register the transformation to bring the
678  // coordinates to this canonical form. In all cases, 1 means no change was
679  // made. We make these transformations so that there are few cases to
680  // check, e.g., on verifying quadrants in atan2. In addition, this
681  // enforces some symmetries in the results returned.
682 
683  double sbet1, cbet1, sbet2, cbet2, s12x, m12x;
684  s12x = m12x = Double.NaN;
685 
686  { Pair p = GeoMath.sincosd(lat1);
687  sbet1 = _f1 * p.first; cbet1 = p.second; }
688  // Ensure cbet1 = +epsilon at poles; doing the fix on beta means that sig12
689  // will be <= 2*tiny for two points at the same pole.
690  { Pair p = GeoMath.norm(sbet1, cbet1); sbet1 = p.first; cbet1 = p.second; }
691  cbet1 = Math.max(tiny_, cbet1);
692 
693  { Pair p = GeoMath.sincosd(lat2);
694  sbet2 = _f1 * p.first; cbet2 = p.second; }
695  // Ensure cbet2 = +epsilon at poles
696  { Pair p = GeoMath.norm(sbet2, cbet2); sbet2 = p.first; cbet2 = p.second; }
697  cbet2 = Math.max(tiny_, cbet2);
698 
699  // If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure of the
700  // |bet1| - |bet2|. Alternatively (cbet1 >= -sbet1), abs(sbet2) + sbet1 is
701  // a better measure. This logic is used in assigning calp2 in Lambda12.
702  // Sometimes these quantities vanish and in that case we force bet2 = +/-
703  // bet1 exactly. An example where is is necessary is the inverse problem
704  // 48.522876735459 0 -48.52287673545898293 179.599720456223079643
705  // which failed with Visual Studio 10 (Release and Debug)
706 
707  if (cbet1 < -sbet1) {
708  if (cbet2 == cbet1)
709  sbet2 = sbet2 < 0 ? sbet1 : -sbet1;
710  } else {
711  if (Math.abs(sbet2) == -sbet1)
712  cbet2 = cbet1;
713  }
714 
715  double
716  dn1 = Math.sqrt(1 + _ep2 * GeoMath.sq(sbet1)),
717  dn2 = Math.sqrt(1 + _ep2 * GeoMath.sq(sbet2));
718 
719  double a12, sig12, calp1, salp1, calp2, salp2;
720  a12 = sig12 = calp1 = salp1 = calp2 = salp2 = Double.NaN;
721  // index zero elements of these arrays are unused
722  double C1a[] = new double[nC1_ + 1];
723  double C2a[] = new double[nC2_ + 1];
724  double C3a[] = new double[nC3_];
725 
726  boolean meridian = lat1 == -90 || slam12 == 0;
727 
728  if (meridian) {
729 
730  // Endpoints are on a single full meridian, so the geodesic might lie on
731  // a meridian.
732 
733  calp1 = clam12; salp1 = slam12; // Head to the target longitude
734  calp2 = 1; salp2 = 0; // At the target we're heading north
735 
736  double
737  // tan(bet) = tan(sig) * cos(alp)
738  ssig1 = sbet1, csig1 = calp1 * cbet1,
739  ssig2 = sbet2, csig2 = calp2 * cbet2;
740 
741  // sig12 = sig2 - sig1
742  sig12 = Math.atan2(Math.max(0.0, csig1 * ssig2 - ssig1 * csig2),
743  csig1 * csig2 + ssig1 * ssig2);
744  {
745  LengthsV v =
746  Lengths(_n, sig12, ssig1, csig1, dn1,
747  ssig2, csig2, dn2, cbet1, cbet2,
749  C1a, C2a);
750  s12x = v.s12b; m12x = v.m12b;
751  if ((outmask & GeodesicMask.GEODESICSCALE) != 0) {
752  r.M12 = v.M12; r.M21 = v.M21;
753  }
754  }
755  // Add the check for sig12 since zero length geodesics might yield m12 <
756  // 0. Test case was
757  //
758  // echo 20.001 0 20.001 0 | GeodSolve -i
759  //
760  // In fact, we will have sig12 > pi/2 for meridional geodesic which is
761  // not a shortest path.
762  if (sig12 < 1 || m12x >= 0) {
763  // Need at least 2, to handle 90 0 90 180
764  if (sig12 < 3 * tiny_)
765  sig12 = m12x = s12x = 0;
766  m12x *= _b;
767  s12x *= _b;
768  a12 = Math.toDegrees(sig12);
769  } else
770  // m12 < 0, i.e., prolate and too close to anti-podal
771  meridian = false;
772  }
773 
774  double omg12 = Double.NaN, somg12 = 2, comg12 = Double.NaN;
775  if (!meridian &&
776  sbet1 == 0 && // and sbet2 == 0
777  // Mimic the way Lambda12 works with calp1 = 0
778  (_f <= 0 || lon12s >= _f * 180)) {
779 
780  // Geodesic runs along equator
781  calp1 = calp2 = 0; salp1 = salp2 = 1;
782  s12x = _a * lam12;
783  sig12 = omg12 = lam12 / _f1;
784  m12x = _b * Math.sin(sig12);
785  if ((outmask & GeodesicMask.GEODESICSCALE) != 0)
786  r.M12 = r.M21 = Math.cos(sig12);
787  a12 = lon12 / _f1;
788 
789  } else if (!meridian) {
790 
791  // Now point1 and point2 belong within a hemisphere bounded by a
792  // meridian and geodesic is neither meridional or equatorial.
793 
794  // Figure a starting point for Newton's method
795  double dnm;
796  {
797  InverseStartV v =
798  InverseStart(sbet1, cbet1, dn1, sbet2, cbet2, dn2,
799  lam12, slam12, clam12,
800  C1a, C2a);
801  sig12 = v.sig12;
802  salp1 = v.salp1; calp1 = v.calp1;
803  salp2 = v.salp2; calp2 = v.calp2;
804  dnm = v.dnm;
805  }
806 
807  if (sig12 >= 0) {
808  // Short lines (InverseStart sets salp2, calp2, dnm)
809  s12x = sig12 * _b * dnm;
810  m12x = GeoMath.sq(dnm) * _b * Math.sin(sig12 / dnm);
811  if ((outmask & GeodesicMask.GEODESICSCALE) != 0)
812  r.M12 = r.M21 = Math.cos(sig12 / dnm);
813  a12 = Math.toDegrees(sig12);
814  omg12 = lam12 / (_f1 * dnm);
815  } else {
816 
817  // Newton's method. This is a straightforward solution of f(alp1) =
818  // lambda12(alp1) - lam12 = 0 with one wrinkle. f(alp) has exactly one
819  // root in the interval (0, pi) and its derivative is positive at the
820  // root. Thus f(alp) is positive for alp > alp1 and negative for alp <
821  // alp1. During the course of the iteration, a range (alp1a, alp1b) is
822  // maintained which brackets the root and with each evaluation of
823  // f(alp) the range is shrunk, if possible. Newton's method is
824  // restarted whenever the derivative of f is negative (because the new
825  // value of alp1 is then further from the solution) or if the new
826  // estimate of alp1 lies outside (0,pi); in this case, the new starting
827  // guess is taken to be (alp1a + alp1b) / 2.
828  double ssig1, csig1, ssig2, csig2, eps, domg12;
829  ssig1 = csig1 = ssig2 = csig2 = eps = domg12 = Double.NaN;
830  int numit = 0;
831  // Bracketing range
832  double salp1a = tiny_, calp1a = 1, salp1b = tiny_, calp1b = -1;
833  for (boolean tripn = false, tripb = false; numit < maxit2_; ++numit) {
834  // the WGS84 test set: mean = 1.47, sd = 1.25, max = 16
835  // WGS84 and random input: mean = 2.85, sd = 0.60
836  double v, dv;
837  {
838  Lambda12V w =
839  Lambda12(sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
840  slam12, clam12, numit < maxit1_, C1a, C2a, C3a);
841  v = w.lam12;
842  salp2 = w.salp2; calp2 = w.calp2;
843  sig12 = w.sig12;
844  ssig1 = w.ssig1; csig1 = w.csig1;
845  ssig2 = w.ssig2; csig2 = w.csig2;
846  eps = w.eps; domg12 = w.domg12;
847  dv = w.dlam12;
848  }
849  // 2 * tol0 is approximately 1 ulp for a number in [0, pi].
850  // Reversed test to allow escape with NaNs
851  if (tripb || !(Math.abs(v) >= (tripn ? 8 : 1) * tol0_)) break;
852  // Update bracketing values
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) {
858  double
859  dalp1 = -v/dv;
860  double
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;
865  salp1 = nsalp1;
866  { Pair p = GeoMath.norm(salp1, calp1);
867  salp1 = p.first; calp1 = p.second; }
868  // In some regimes we don't get quadratic convergence because
869  // slope -> 0. So use convergence conditions based on epsilon
870  // instead of sqrt(epsilon).
871  tripn = Math.abs(v) <= 16 * tol0_;
872  continue;
873  }
874  }
875  // Either dv was not positive or updated value was outside legal
876  // range. Use the midpoint of the bracket as the next estimate.
877  // This mechanism is not needed for the WGS84 ellipsoid, but it does
878  // catch problems with more eccentric ellipsoids. Its efficacy is
879  // such for the WGS84 test set with the starting guess set to alp1 =
880  // 90deg:
881  // the WGS84 test set: mean = 5.21, sd = 3.93, max = 24
882  // WGS84 and random input: mean = 4.74, sd = 0.99
883  salp1 = (salp1a + salp1b)/2;
884  calp1 = (calp1a + calp1b)/2;
885  { Pair p = GeoMath.norm(salp1, calp1);
886  salp1 = p.first; calp1 = p.second; }
887  tripn = false;
888  tripb = (Math.abs(salp1a - salp1) + (calp1a - calp1) < tolb_ ||
889  Math.abs(salp1 - salp1b) + (calp1 - calp1b) < tolb_);
890  }
891  {
892  // Ensure that the reduced length and geodesic scale are computed in
893  // a "canonical" way, with the I2 integral.
894  int lengthmask = outmask |
895  ((outmask &
898  LengthsV v =
899  Lengths(eps, sig12,
900  ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2,
901  lengthmask, C1a, C2a);
902  s12x = v.s12b; m12x = v.m12b;
903  if ((outmask & GeodesicMask.GEODESICSCALE) != 0) {
904  r.M12 = v.M12; r.M21 = v.M21;
905  }
906  }
907  m12x *= _b;
908  s12x *= _b;
909  a12 = Math.toDegrees(sig12);
910  if ((outmask & GeodesicMask.AREA) != 0) {
911  // omg12 = lam12 - domg12
912  double sdomg12 = Math.sin(domg12), cdomg12 = Math.cos(domg12);
913  somg12 = slam12 * cdomg12 - clam12 * sdomg12;
914  comg12 = clam12 * cdomg12 + slam12 * sdomg12;
915  }
916  }
917  }
918 
919  if ((outmask & GeodesicMask.DISTANCE) != 0)
920  r.s12 = 0 + s12x; // Convert -0 to 0
921 
922  if ((outmask & GeodesicMask.REDUCEDLENGTH) != 0)
923  r.m12 = 0 + m12x; // Convert -0 to 0
924 
925  if ((outmask & GeodesicMask.AREA) != 0) {
926  double
927  // From Lambda12: sin(alp1) * cos(bet1) = sin(alp0)
928  salp0 = salp1 * cbet1,
929  calp0 = GeoMath.hypot(calp1, salp1 * sbet1); // calp0 > 0
930  double alp12;
931  if (calp0 != 0 && salp0 != 0) {
932  double
933  // From Lambda12: tan(bet) = tan(sig) * cos(alp)
934  ssig1 = sbet1, csig1 = calp1 * cbet1,
935  ssig2 = sbet2, csig2 = calp2 * cbet2,
936  k2 = GeoMath.sq(calp0) * _ep2,
937  eps = k2 / (2 * (1 + Math.sqrt(1 + k2)) + k2),
938  // Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0).
939  A4 = GeoMath.sq(_a) * calp0 * salp0 * _e2;
940  { Pair p = GeoMath.norm(ssig1, csig1);
941  ssig1 = p.first; csig1 = p.second; }
942  { Pair p = GeoMath.norm(ssig2, csig2);
943  ssig2 = p.first; csig2 = p.second; }
944  double C4a[] = new double[nC4_];
945  C4f(eps, C4a);
946  double
947  B41 = SinCosSeries(false, ssig1, csig1, C4a),
948  B42 = SinCosSeries(false, ssig2, csig2, C4a);
949  r.S12 = A4 * (B42 - B41);
950  } else
951  // Avoid problems with indeterminate sig1, sig2 on equator
952  r.S12 = 0;
953 
954  if (!meridian && somg12 > 1) {
955  somg12 = Math.sin(omg12); comg12 = Math.cos(omg12);
956  }
957 
958  if (!meridian &&
959  comg12 > -0.7071 && // Long difference not too big
960  sbet2 - sbet1 < 1.75) { // Lat difference not too big
961  // Use tan(Gamma/2) = tan(omg12/2)
962  // * (tan(bet1/2)+tan(bet2/2))/(1+tan(bet1/2)*tan(bet2/2))
963  // with tan(x/2) = sin(x)/(1+cos(x))
964  double
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 ) );
968  } else {
969  // alp12 = alp2 - alp1, used in atan2 so no need to normalize
970  double
971  salp12 = salp2 * calp1 - calp2 * salp1,
972  calp12 = calp2 * calp1 + salp2 * salp1;
973  // The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz
974  // salp12 = -0 and alp12 = -180. However this depends on the sign
975  // being attached to 0 correctly. The following ensures the correct
976  // behavior.
977  if (salp12 == 0 && calp12 < 0) {
978  salp12 = tiny_ * calp1;
979  calp12 = -1;
980  }
981  alp12 = Math.atan2(salp12, calp12);
982  }
983  r.S12 += _c2 * alp12;
984  r.S12 *= swapp * lonsign * latsign;
985  // Convert -0 to 0
986  r.S12 += 0;
987  }
988 
989  // Convert calp, salp to azimuth accounting for lonsign, swapp, latsign.
990  if (swapp < 0) {
991  { double t = salp1; salp1 = salp2; salp2 = t; }
992  { double t = calp1; calp1 = calp2; calp2 = t; }
993  if ((outmask & GeodesicMask.GEODESICSCALE) != 0)
994  { double t = r.M12; r.M12 = r.M21; r.M21 = t; }
995  }
996 
997  salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
998  salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
999 
1000  // Returned value in [0, 180]
1001  r.a12 = a12;
1002  result.salp1 = salp1; result.calp1 = calp1;
1003  result.salp2 = salp2; result.calp2 = calp2;
1004  return result;
1005  }
1006 
1049  public GeodesicData Inverse(double lat1, double lon1,
1050  double lat2, double lon2, int outmask) {
1051  outmask &= GeodesicMask.OUT_MASK;
1052  InverseData result = InverseInt(lat1, lon1, lat2, lon2, outmask);
1053  GeodesicData r = result.g;
1054  if ((outmask & GeodesicMask.AZIMUTH) != 0) {
1055  r.azi1 = GeoMath.atan2d(result.salp1, result.calp1);
1056  r.azi2 = GeoMath.atan2d(result.salp2, result.calp2);
1057  }
1058  return r;
1059  }
1060 
1077  public GeodesicLine InverseLine(double lat1, double lon1,
1078  double lat2, double lon2) {
1079  return InverseLine(lat1, lon1, lat2, lon2, GeodesicMask.ALL);
1080  }
1081 
1102  public GeodesicLine InverseLine(double lat1, double lon1,
1103  double lat2, double lon2, int caps) {
1104  InverseData result = InverseInt(lat1, lon1, lat2, lon2, 0);
1105  double salp1 = result.salp1, calp1 = result.calp1,
1106  azi1 = GeoMath.atan2d(salp1, calp1), a12 = result.g.a12;
1107  // Ensure that a12 can be converted to a distance
1108  if ((caps & (GeodesicMask.OUT_MASK & GeodesicMask.DISTANCE_IN)) != 0)
1109  caps |= GeodesicMask.DISTANCE;
1110  return new GeodesicLine(this, lat1, lon1, azi1, salp1, calp1, caps,
1111  true, a12);
1112  }
1113 
1130  public GeodesicLine Line(double lat1, double lon1, double azi1) {
1131  return Line(lat1, lon1, azi1, GeodesicMask.ALL);
1132  }
1180  public GeodesicLine Line(double lat1, double lon1, double azi1, int caps) {
1181  return new GeodesicLine(this, lat1, lon1, azi1, caps);
1182  }
1183 
1188  public double MajorRadius() { return _a; }
1189 
1194  public double Flattening() { return _f; }
1195 
1201  public double EllipsoidArea() { return 4 * Math.PI * _c2; }
1202 
1207  public static final Geodesic WGS84 =
1209 
1210  // This is a reformulation of the geodesic problem. The notation is as
1211  // follows:
1212  // - at a general point (no suffix or 1 or 2 as suffix)
1213  // - phi = latitude
1214  // - beta = latitude on auxiliary sphere
1215  // - omega = longitude on auxiliary sphere
1216  // - lambda = longitude
1217  // - alpha = azimuth of great circle
1218  // - sigma = arc length along great circle
1219  // - s = distance
1220  // - tau = scaled distance (= sigma at multiples of pi/2)
1221  // - at northwards equator crossing
1222  // - beta = phi = 0
1223  // - omega = lambda = 0
1224  // - alpha = alpha0
1225  // - sigma = s = 0
1226  // - a 12 suffix means a difference, e.g., s12 = s2 - s1.
1227  // - s and c prefixes mean sin and cos
1228 
1229  protected static double SinCosSeries(boolean sinp,
1230  double sinx, double cosx,
1231  double c[]) {
1232  // Evaluate
1233  // y = sinp ? sum(c[i] * sin( 2*i * x), i, 1, n) :
1234  // sum(c[i] * cos((2*i+1) * x), i, 0, n-1)
1235  // using Clenshaw summation. N.B. c[0] is unused for sin series
1236  // Approx operation count = (n + 5) mult and (2 * n + 2) add
1237  int
1238  k = c.length, // Point to one beyond last element
1239  n = k - (sinp ? 1 : 0);
1240  double
1241  ar = 2 * (cosx - sinx) * (cosx + sinx), // 2 * cos(2 * x)
1242  y0 = (n & 1) != 0 ? c[--k] : 0, y1 = 0; // accumulators for sum
1243  // Now n is even
1244  n /= 2;
1245  while (n-- != 0) {
1246  // Unroll loop x 2, so accumulators return to their original role
1247  y1 = ar * y0 - y1 + c[--k];
1248  y0 = ar * y1 - y0 + c[--k];
1249  }
1250  return sinp
1251  ? 2 * sinx * cosx * y0 // sin(2 * x) * y0
1252  : cosx * (y0 - y1); // cos(x) * (y0 - y1)
1253  }
1254 
1255  private class LengthsV {
1256  private double s12b, m12b, m0, M12, M21;
1257  private LengthsV() {
1258  s12b = m12b = m0 = M12 = M21 = Double.NaN;
1259  }
1260  }
1261 
1262  private LengthsV Lengths(double eps, double sig12,
1263  double ssig1, double csig1, double dn1,
1264  double ssig2, double csig2, double dn2,
1265  double cbet1, double cbet2,
1266  int outmask,
1267  // Scratch areas of the right size
1268  double C1a[], double C2a[]) {
1269  // Return m12b = (reduced length)/_b; also calculate s12b = distance/_b,
1270  // and m0 = coefficient of secular term in expression for reduced length.
1271  outmask &= GeodesicMask.OUT_MASK;
1272  LengthsV v = new LengthsV(); // To hold s12b, m12b, m0, M12, M21;
1273 
1274  double m0x = 0, J12 = 0, A1 = 0, A2 = 0;
1275  if ((outmask & (GeodesicMask.DISTANCE | GeodesicMask.REDUCEDLENGTH |
1276  GeodesicMask.GEODESICSCALE)) != 0) {
1277  A1 = A1m1f(eps);
1278  C1f(eps, C1a);
1279  if ((outmask & (GeodesicMask.REDUCEDLENGTH |
1280  GeodesicMask.GEODESICSCALE)) != 0) {
1281  A2 = A2m1f(eps);
1282  C2f(eps, C2a);
1283  m0x = A1 - A2;
1284  A2 = 1 + A2;
1285  }
1286  A1 = 1 + A1;
1287  }
1288  if ((outmask & GeodesicMask.DISTANCE) != 0) {
1289  double B1 = SinCosSeries(true, ssig2, csig2, C1a) -
1290  SinCosSeries(true, ssig1, csig1, C1a);
1291  // Missing a factor of _b
1292  v.s12b = A1 * (sig12 + B1);
1293  if ((outmask & (GeodesicMask.REDUCEDLENGTH |
1294  GeodesicMask.GEODESICSCALE)) != 0) {
1295  double B2 = SinCosSeries(true, ssig2, csig2, C2a) -
1296  SinCosSeries(true, ssig1, csig1, C2a);
1297  J12 = m0x * sig12 + (A1 * B1 - A2 * B2);
1298  }
1299  } else if ((outmask & (GeodesicMask.REDUCEDLENGTH |
1300  GeodesicMask.GEODESICSCALE)) != 0) {
1301  // Assume here that nC1_ >= nC2_
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) -
1305  SinCosSeries(true, ssig1, csig1, C2a));
1306  }
1307  if ((outmask & GeodesicMask.REDUCEDLENGTH) != 0) {
1308  v.m0 = m0x;
1309  // Missing a factor of _b.
1310  // Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure
1311  // accurate cancellation in the case of coincident points.
1312  v.m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
1313  csig1 * csig2 * J12;
1314  }
1315  if ((outmask & GeodesicMask.GEODESICSCALE) != 0) {
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;
1320  }
1321  return v;
1322  }
1323 
1324  private static double Astroid(double x, double y) {
1325  // Solve k^4+2*k^3-(x^2+y^2-1)*k^2-2*y^2*k-y^2 = 0 for positive root k.
1326  // This solution is adapted from Geocentric::Reverse.
1327  double k;
1328  double
1329  p = GeoMath.sq(x),
1330  q = GeoMath.sq(y),
1331  r = (p + q - 1) / 6;
1332  if ( !(q == 0 && r <= 0) ) {
1333  double
1334  // Avoid possible division by zero when r = 0 by multiplying equations
1335  // for s and t by r^3 and r, resp.
1336  S = p * q / 4, // S = r^3 * s
1337  r2 = GeoMath.sq(r),
1338  r3 = r * r2,
1339  // The discriminant of the quadratic equation for T3. This is zero on
1340  // the evolute curve p^(1/3)+q^(1/3) = 1
1341  disc = S * (S + 2 * r3);
1342  double u = r;
1343  if (disc >= 0) {
1344  double T3 = S + r3;
1345  // Pick the sign on the sqrt to maximize abs(T3). This minimizes loss
1346  // of precision due to cancellation. The result is unchanged because
1347  // of the way the T is used in definition of u.
1348  T3 += T3 < 0 ? -Math.sqrt(disc) : Math.sqrt(disc); // T3 = (r * t)^3
1349  // N.B. cbrt always returns the double root. cbrt(-8) = -2.
1350  double T = GeoMath.cbrt(T3); // T = r * t
1351  // T can be zero; but then r2 / T -> 0.
1352  u += T + (T != 0 ? r2 / T : 0);
1353  } else {
1354  // T is complex, but the way u is defined the result is double.
1355  double ang = Math.atan2(Math.sqrt(-disc), -(S + r3));
1356  // There are three possible cube roots. We choose the root which
1357  // avoids cancellation. Note that disc < 0 implies that r < 0.
1358  u += 2 * r * Math.cos(ang / 3);
1359  }
1360  double
1361  v = Math.sqrt(GeoMath.sq(u) + q), // guaranteed positive
1362  // Avoid loss of accuracy when u < 0.
1363  uv = u < 0 ? q / (v - u) : u + v, // u+v, guaranteed positive
1364  w = (uv - q) / (2 * v); // positive?
1365  // Rearrange expression for k to avoid loss of accuracy due to
1366  // subtraction. Division by 0 not possible because uv > 0, w >= 0.
1367  k = uv / (Math.sqrt(uv + GeoMath.sq(w)) + w); // guaranteed positive
1368  } else { // q == 0 && r <= 0
1369  // y = 0 with |x| <= 1. Handle this case directly.
1370  // for y small, positive root is k = abs(y)/sqrt(1-x^2)
1371  k = 0;
1372  }
1373  return k;
1374  }
1375 
1376  private class InverseStartV {
1377  private double sig12, salp1, calp1,
1378  // Only updated if return val >= 0
1379  salp2, calp2,
1380  // Only updated for short lines
1381  dnm;
1382  private InverseStartV() {
1383  sig12 = salp1 = calp1 = salp2 = calp2 = dnm = Double.NaN;
1384  }
1385  }
1386 
1387  private InverseStartV InverseStart(double sbet1, double cbet1, double dn1,
1388  double sbet2, double cbet2, double dn2,
1389  double lam12,
1390  double slam12, double clam12,
1391  // Scratch areas of the right size
1392  double C1a[], double C2a[]) {
1393  // Return a starting point for Newton's method in salp1 and calp1 (function
1394  // value is -1). If Newton's method doesn't need to be used, return also
1395  // salp2 and calp2 and function value is sig12.
1396 
1397  // To hold sig12, salp1, calp1, salp2, calp2, dnm.
1398  InverseStartV w = new InverseStartV();
1399  w.sig12 = -1; // Return value
1400  double
1401  // bet12 = bet2 - bet1 in [0, pi); bet12a = bet2 + bet1 in (-pi, 0]
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;
1408  if (shortline) {
1409  double sbetm2 = GeoMath.sq(sbet1 + sbet2);
1410  // sin((bet1+bet2)/2)^2
1411  // = (sbet1 + sbet2)^2 / ((sbet1 + sbet2)^2 + (cbet1 + cbet2)^2)
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);
1416  } else {
1417  somg12 = slam12; comg12 = clam12;
1418  }
1419 
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);
1424 
1425  double
1426  ssig12 = GeoMath.hypot(w.salp1, w.calp1),
1427  csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
1428 
1429  if (shortline && ssig12 < _etol2) {
1430  // really short lines
1431  w.salp2 = cbet1 * somg12;
1432  w.calp2 = sbet12 - cbet1 * sbet2 *
1433  (comg12 >= 0 ? GeoMath.sq(somg12) / (1 + comg12) : 1 - comg12);
1434  { Pair p = GeoMath.norm(w.salp2, w.calp2);
1435  w.salp2 = p.first; w.calp2 = p.second; }
1436  // Set return value
1437  w.sig12 = Math.atan2(ssig12, csig12);
1438  } else if (Math.abs(_n) > 0.1 || // Skip astroid calc if too eccentric
1439  csig12 >= 0 ||
1440  ssig12 >= 6 * Math.abs(_n) * Math.PI * GeoMath.sq(cbet1)) {
1441  // Nothing to do, zeroth order spherical approximation is OK
1442  } else {
1443  // Scale lam12 and bet2 to x, y coordinate system where antipodal point
1444  // is at origin and singular point is at y = 0, x = -1.
1445  double y, lamscale, betscale;
1446  // In C++ volatile declaration needed to fix inverse case
1447  // 56.320923501171 0 -56.320923501171 179.664747671772880215
1448  // which otherwise fails with g++ 4.4.4 x86 -O3
1449  double x;
1450  double lam12x = Math.atan2(-slam12, -clam12); // lam12 - pi
1451  if (_f >= 0) { // In fact f == 0 does not get here
1452  // x = dlong, y = dlat
1453  {
1454  double
1455  k2 = GeoMath.sq(sbet1) * _ep2,
1456  eps = k2 / (2 * (1 + Math.sqrt(1 + k2)) + k2);
1457  lamscale = _f * cbet1 * A3f(eps) * Math.PI;
1458  }
1459  betscale = lamscale * cbet1;
1460 
1461  x = lam12x / lamscale;
1462  y = sbet12a / betscale;
1463  } else { // _f < 0
1464  // x = dlat, y = dlong
1465  double
1466  cbet12a = cbet2 * cbet1 - sbet2 * sbet1,
1467  bet12a = Math.atan2(sbet12a, cbet12a);
1468  double m12b, m0;
1469  // In the case of lon12 = 180, this repeats a calculation made in
1470  // Inverse.
1471  LengthsV v =
1472  Lengths(_n, Math.PI + bet12a,
1473  sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
1474  cbet1, cbet2, GeodesicMask.REDUCEDLENGTH, C1a, C2a);
1475  m12b = v.m12b; m0 = v.m0;
1476 
1477  x = -1 + m12b / (cbet1 * cbet2 * m0 * Math.PI);
1478  betscale = x < -0.01 ? sbet12a / x :
1479  -_f * GeoMath.sq(cbet1) * Math.PI;
1480  lamscale = betscale / cbet1;
1481  y = lam12x / lamscale;
1482  }
1483 
1484  if (y > -tol1_ && x > -1 - xthresh_) {
1485  // strip near cut
1486  if (_f >= 0) {
1487  w.salp1 = Math.min(1.0, -x);
1488  w.calp1 = - Math.sqrt(1 - GeoMath.sq(w.salp1));
1489  } else {
1490  w.calp1 = Math.max(x > -tol1_ ? 0.0 : -1.0, x);
1491  w.salp1 = Math.sqrt(1 - GeoMath.sq(w.calp1));
1492  }
1493  } else {
1494  // Estimate alp1, by solving the astroid problem.
1495  //
1496  // Could estimate alpha1 = theta + pi/2, directly, i.e.,
1497  // calp1 = y/k; salp1 = -x/(1+k); for _f >= 0
1498  // calp1 = x/(1+k); salp1 = -y/k; for _f < 0 (need to check)
1499  //
1500  // However, it's better to estimate omg12 from astroid and use
1501  // spherical formula to compute alp1. This reduces the mean number of
1502  // Newton iterations for astroid cases from 2.24 (min 0, max 6) to 2.12
1503  // (min 0 max 5). The changes in the number of iterations are as
1504  // follows:
1505  //
1506  // change percent
1507  // 1 5
1508  // 0 78
1509  // -1 16
1510  // -2 0.6
1511  // -3 0.04
1512  // -4 0.002
1513  //
1514  // The histogram of iterations is (m = number of iterations estimating
1515  // alp1 directly, n = number of iterations estimating via omg12, total
1516  // number of trials = 148605):
1517  //
1518  // iter m n
1519  // 0 148 186
1520  // 1 13046 13845
1521  // 2 93315 102225
1522  // 3 36189 32341
1523  // 4 5396 7
1524  // 5 455 1
1525  // 6 56 0
1526  //
1527  // Because omg12 is near pi, estimate work with omg12a = pi - omg12
1528  double k = Astroid(x, y);
1529  double
1530  omg12a = lamscale * ( _f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k );
1531  somg12 = Math.sin(omg12a); comg12 = -Math.cos(omg12a);
1532  // Update spherical estimate of alp1 using omg12 instead of lam12
1533  w.salp1 = cbet2 * somg12;
1534  w.calp1 = sbet12a - cbet2 * sbet1 * GeoMath.sq(somg12) / (1 - comg12);
1535  }
1536  }
1537  // Sanity check on starting guess. Backwards check allows NaN through.
1538  if (!(w.salp1 <= 0))
1539  { Pair p = GeoMath.norm(w.salp1, w.calp1);
1540  w.salp1 = p.first; w.calp1 = p.second; }
1541  else {
1542  w.salp1 = 1; w.calp1 = 0;
1543  }
1544  return w;
1545  }
1546 
1547  private class Lambda12V {
1548  private double lam12, salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
1549  eps, domg12, dlam12;
1550  private Lambda12V() {
1551  lam12 = salp2 = calp2 = sig12 = ssig1 = csig1 = ssig2 = csig2
1552  = eps = domg12 = dlam12 = Double.NaN;
1553  }
1554  }
1555 
1556  private Lambda12V Lambda12(double sbet1, double cbet1, double dn1,
1557  double sbet2, double cbet2, double dn2,
1558  double salp1, double calp1,
1559  double slam120, double clam120,
1560  boolean diffp,
1561  // Scratch areas of the right size
1562  double C1a[], double C2a[], double C3a[]) {
1563  // Object to hold lam12, salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
1564  // eps, domg12, dlam12;
1565 
1566  Lambda12V w = new Lambda12V();
1567 
1568  if (sbet1 == 0 && calp1 == 0)
1569  // Break degeneracy of equatorial line. This case has already been
1570  // handled.
1571  calp1 = -tiny_;
1572 
1573  double
1574  // sin(alp1) * cos(bet1) = sin(alp0)
1575  salp0 = salp1 * cbet1,
1576  calp0 = GeoMath.hypot(calp1, salp1 * sbet1); // calp0 > 0
1577 
1578  double somg1, comg1, somg2, comg2, somg12, comg12;
1579  // tan(bet1) = tan(sig1) * cos(alp1)
1580  // tan(omg1) = sin(alp0) * tan(sig1) = tan(omg1)=tan(alp1)*sin(bet1)
1581  w.ssig1 = sbet1; somg1 = salp0 * sbet1;
1582  w.csig1 = comg1 = calp1 * cbet1;
1583  { Pair p = GeoMath.norm(w.ssig1, w.csig1);
1584  w.ssig1 = p.first; w.csig1 = p.second; }
1585  // GeoMath.norm(somg1, comg1); -- don't need to normalize!
1586 
1587  // Enforce symmetries in the case abs(bet2) = -bet1. Need to be careful
1588  // about this case, since this can yield singularities in the Newton
1589  // iteration.
1590  // sin(alp2) * cos(bet2) = sin(alp0)
1591  w.salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1;
1592  // calp2 = sqrt(1 - sq(salp2))
1593  // = sqrt(sq(calp0) - sq(sbet2)) / cbet2
1594  // and subst for calp0 and rearrange to give (choose positive sqrt
1595  // to give alp2 in [0, pi/2]).
1596  w.calp2 = cbet2 != cbet1 || Math.abs(sbet2) != -sbet1 ?
1597  Math.sqrt(GeoMath.sq(calp1 * cbet1) +
1598  (cbet1 < -sbet1 ?
1599  (cbet2 - cbet1) * (cbet1 + cbet2) :
1600  (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :
1601  Math.abs(calp1);
1602  // tan(bet2) = tan(sig2) * cos(alp2)
1603  // tan(omg2) = sin(alp0) * tan(sig2).
1604  w.ssig2 = sbet2; somg2 = salp0 * sbet2;
1605  w.csig2 = comg2 = w.calp2 * cbet2;
1606  { Pair p = GeoMath.norm(w.ssig2, w.csig2);
1607  w.ssig2 = p.first; w.csig2 = p.second; }
1608  // GeoMath.norm(somg2, comg2); -- don't need to normalize!
1609 
1610  // sig12 = sig2 - sig1, limit to [0, pi]
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);
1613 
1614  // omg12 = omg2 - omg1, limit to [0, pi]
1615  somg12 = Math.max(0.0, comg1 * somg2 - somg1 * comg2);
1616  comg12 = comg1 * comg2 + somg1 * somg2;
1617  // eta = omg12 - lam120
1618  double eta = Math.atan2(somg12 * clam120 - comg12 * slam120,
1619  comg12 * clam120 + somg12 * slam120);
1620  double B312;
1621  double k2 = GeoMath.sq(calp0) * _ep2;
1622  w.eps = k2 / (2 * (1 + Math.sqrt(1 + k2)) + k2);
1623  C3f(w.eps, C3a);
1624  B312 = (SinCosSeries(true, w.ssig2, w.csig2, C3a) -
1625  SinCosSeries(true, w.ssig1, w.csig1, C3a));
1626  w.domg12 = -_f * A3f(w.eps) * salp0 * (w.sig12 + B312);
1627  w.lam12 = eta + w.domg12;
1628 
1629  if (diffp) {
1630  if (w.calp2 == 0)
1631  w.dlam12 = - 2 * _f1 * dn1 / sbet1;
1632  else {
1633  LengthsV v =
1634  Lengths(w.eps, w.sig12, w.ssig1, w.csig1, dn1, w.ssig2, w.csig2, dn2,
1635  cbet1, cbet2, GeodesicMask.REDUCEDLENGTH, C1a, C2a);
1636  w.dlam12 = v.m12b;
1637  w.dlam12 *= _f1 / (w.calp2 * cbet2);
1638  }
1639  }
1640 
1641  return w;
1642  }
1643 
1644  protected double A3f(double eps) {
1645  // Evaluate A3
1646  return GeoMath.polyval(nA3_ - 1, _A3x, 0, eps);
1647  }
1648 
1649  protected void C3f(double eps, double c[]) {
1650  // Evaluate C3 coeffs
1651  // Elements c[1] thru c[nC3_ - 1] are set
1652  double mult = 1;
1653  int o = 0;
1654  for (int l = 1; l < nC3_; ++l) { // l is index of C3[l]
1655  int m = nC3_ - l - 1; // order of polynomial in eps
1656  mult *= eps;
1657  c[l] = mult * GeoMath.polyval(m, _C3x, o, eps);
1658  o += m + 1;
1659  }
1660  }
1661 
1662  protected void C4f(double eps, double c[]) {
1663  // Evaluate C4 coeffs
1664  // Elements c[0] thru c[nC4_ - 1] are set
1665  double mult = 1;
1666  int o = 0;
1667  for (int l = 0; l < nC4_; ++l) { // l is index of C4[l]
1668  int m = nC4_ - l - 1; // order of polynomial in eps
1669  c[l] = mult * GeoMath.polyval(m, _C4x, o, eps);
1670  o += m + 1;
1671  mult *= eps;
1672  }
1673  }
1674 
1675  // The scale factor A1-1 = mean value of (d/dsigma)I1 - 1
1676  protected static double A1m1f(double eps) {
1677  final double coeff[] = {
1678  // (1-eps)*A1-1, polynomial in eps2 of order 3
1679  1, 4, 64, 0, 256,
1680  };
1681  int m = nA1_/2;
1682  double t = GeoMath.polyval(m, coeff, 0, GeoMath.sq(eps)) / coeff[m + 1];
1683  return (t + eps) / (1 - eps);
1684  }
1685 
1686  // The coefficients C1[l] in the Fourier expansion of B1
1687  protected static void C1f(double eps, double c[]) {
1688  final double coeff[] = {
1689  // C1[1]/eps^1, polynomial in eps2 of order 2
1690  -1, 6, -16, 32,
1691  // C1[2]/eps^2, polynomial in eps2 of order 2
1692  -9, 64, -128, 2048,
1693  // C1[3]/eps^3, polynomial in eps2 of order 1
1694  9, -16, 768,
1695  // C1[4]/eps^4, polynomial in eps2 of order 1
1696  3, -5, 512,
1697  // C1[5]/eps^5, polynomial in eps2 of order 0
1698  -7, 1280,
1699  // C1[6]/eps^6, polynomial in eps2 of order 0
1700  -7, 2048,
1701  };
1702  double
1703  eps2 = GeoMath.sq(eps),
1704  d = eps;
1705  int o = 0;
1706  for (int l = 1; l <= nC1_; ++l) { // l is index of C1p[l]
1707  int m = (nC1_ - l) / 2; // order of polynomial in eps^2
1708  c[l] = d * GeoMath.polyval(m, coeff, o, eps2) / coeff[o + m + 1];
1709  o += m + 2;
1710  d *= eps;
1711  }
1712  }
1713 
1714  // The coefficients C1p[l] in the Fourier expansion of B1p
1715  protected static void C1pf(double eps, double c[]) {
1716  final double coeff[] = {
1717  // C1p[1]/eps^1, polynomial in eps2 of order 2
1718  205, -432, 768, 1536,
1719  // C1p[2]/eps^2, polynomial in eps2 of order 2
1720  4005, -4736, 3840, 12288,
1721  // C1p[3]/eps^3, polynomial in eps2 of order 1
1722  -225, 116, 384,
1723  // C1p[4]/eps^4, polynomial in eps2 of order 1
1724  -7173, 2695, 7680,
1725  // C1p[5]/eps^5, polynomial in eps2 of order 0
1726  3467, 7680,
1727  // C1p[6]/eps^6, polynomial in eps2 of order 0
1728  38081, 61440,
1729  };
1730  double
1731  eps2 = GeoMath.sq(eps),
1732  d = eps;
1733  int o = 0;
1734  for (int l = 1; l <= nC1p_; ++l) { // l is index of C1p[l]
1735  int m = (nC1p_ - l) / 2; // order of polynomial in eps^2
1736  c[l] = d * GeoMath.polyval(m, coeff, o, eps2) / coeff[o + m + 1];
1737  o += m + 2;
1738  d *= eps;
1739  }
1740  }
1741 
1742  // The scale factor A2-1 = mean value of (d/dsigma)I2 - 1
1743  protected static double A2m1f(double eps) {
1744  final double coeff[] = {
1745  // (eps+1)*A2-1, polynomial in eps2 of order 3
1746  -11, -28, -192, 0, 256,
1747  };
1748  int m = nA2_/2;
1749  double t = GeoMath.polyval(m, coeff, 0, GeoMath.sq(eps)) / coeff[m + 1];
1750  return (t - eps) / (1 + eps);
1751  }
1752 
1753  // The coefficients C2[l] in the Fourier expansion of B2
1754  protected static void C2f(double eps, double c[]) {
1755  final double coeff[] = {
1756  // C2[1]/eps^1, polynomial in eps2 of order 2
1757  1, 2, 16, 32,
1758  // C2[2]/eps^2, polynomial in eps2 of order 2
1759  35, 64, 384, 2048,
1760  // C2[3]/eps^3, polynomial in eps2 of order 1
1761  15, 80, 768,
1762  // C2[4]/eps^4, polynomial in eps2 of order 1
1763  7, 35, 512,
1764  // C2[5]/eps^5, polynomial in eps2 of order 0
1765  63, 1280,
1766  // C2[6]/eps^6, polynomial in eps2 of order 0
1767  77, 2048,
1768  };
1769  double
1770  eps2 = GeoMath.sq(eps),
1771  d = eps;
1772  int o = 0;
1773  for (int l = 1; l <= nC2_; ++l) { // l is index of C2[l]
1774  int m = (nC2_ - l) / 2; // order of polynomial in eps^2
1775  c[l] = d * GeoMath.polyval(m, coeff, o, eps2) / coeff[o + m + 1];
1776  o += m + 2;
1777  d *= eps;
1778  }
1779  }
1780 
1781  // The scale factor A3 = mean value of (d/dsigma)I3
1782  protected void A3coeff() {
1783  final double coeff[] = {
1784  // A3, coeff of eps^5, polynomial in n of order 0
1785  -3, 128,
1786  // A3, coeff of eps^4, polynomial in n of order 1
1787  -2, -3, 64,
1788  // A3, coeff of eps^3, polynomial in n of order 2
1789  -1, -3, -1, 16,
1790  // A3, coeff of eps^2, polynomial in n of order 2
1791  3, -1, -2, 8,
1792  // A3, coeff of eps^1, polynomial in n of order 1
1793  1, -1, 2,
1794  // A3, coeff of eps^0, polynomial in n of order 0
1795  1, 1,
1796  };
1797  int o = 0, k = 0;
1798  for (int j = nA3_ - 1; j >= 0; --j) { // coeff of eps^j
1799  int m = Math.min(nA3_ - j - 1, j); // order of polynomial in n
1800  _A3x[k++] = GeoMath.polyval(m, coeff, o, _n) / coeff[o + m + 1];
1801  o += m + 2;
1802  }
1803  }
1804 
1805  // The coefficients C3[l] in the Fourier expansion of B3
1806  protected void C3coeff() {
1807  final double coeff[] = {
1808  // C3[1], coeff of eps^5, polynomial in n of order 0
1809  3, 128,
1810  // C3[1], coeff of eps^4, polynomial in n of order 1
1811  2, 5, 128,
1812  // C3[1], coeff of eps^3, polynomial in n of order 2
1813  -1, 3, 3, 64,
1814  // C3[1], coeff of eps^2, polynomial in n of order 2
1815  -1, 0, 1, 8,
1816  // C3[1], coeff of eps^1, polynomial in n of order 1
1817  -1, 1, 4,
1818  // C3[2], coeff of eps^5, polynomial in n of order 0
1819  5, 256,
1820  // C3[2], coeff of eps^4, polynomial in n of order 1
1821  1, 3, 128,
1822  // C3[2], coeff of eps^3, polynomial in n of order 2
1823  -3, -2, 3, 64,
1824  // C3[2], coeff of eps^2, polynomial in n of order 2
1825  1, -3, 2, 32,
1826  // C3[3], coeff of eps^5, polynomial in n of order 0
1827  7, 512,
1828  // C3[3], coeff of eps^4, polynomial in n of order 1
1829  -10, 9, 384,
1830  // C3[3], coeff of eps^3, polynomial in n of order 2
1831  5, -9, 5, 192,
1832  // C3[4], coeff of eps^5, polynomial in n of order 0
1833  7, 512,
1834  // C3[4], coeff of eps^4, polynomial in n of order 1
1835  -14, 7, 512,
1836  // C3[5], coeff of eps^5, polynomial in n of order 0
1837  21, 2560,
1838  };
1839  int o = 0, k = 0;
1840  for (int l = 1; l < nC3_; ++l) { // l is index of C3[l]
1841  for (int j = nC3_ - 1; j >= l; --j) { // coeff of eps^j
1842  int m = Math.min(nC3_ - j - 1, j); // order of polynomial in n
1843  _C3x[k++] = GeoMath.polyval(m, coeff, o, _n) / coeff[o + m + 1];
1844  o += m + 2;
1845  }
1846  }
1847  }
1848 
1849  protected void C4coeff() {
1850  final double coeff[] = {
1851  // C4[0], coeff of eps^5, polynomial in n of order 0
1852  97, 15015,
1853  // C4[0], coeff of eps^4, polynomial in n of order 1
1854  1088, 156, 45045,
1855  // C4[0], coeff of eps^3, polynomial in n of order 2
1856  -224, -4784, 1573, 45045,
1857  // C4[0], coeff of eps^2, polynomial in n of order 3
1858  -10656, 14144, -4576, -858, 45045,
1859  // C4[0], coeff of eps^1, polynomial in n of order 4
1860  64, 624, -4576, 6864, -3003, 15015,
1861  // C4[0], coeff of eps^0, polynomial in n of order 5
1862  100, 208, 572, 3432, -12012, 30030, 45045,
1863  // C4[1], coeff of eps^5, polynomial in n of order 0
1864  1, 9009,
1865  // C4[1], coeff of eps^4, polynomial in n of order 1
1866  -2944, 468, 135135,
1867  // C4[1], coeff of eps^3, polynomial in n of order 2
1868  5792, 1040, -1287, 135135,
1869  // C4[1], coeff of eps^2, polynomial in n of order 3
1870  5952, -11648, 9152, -2574, 135135,
1871  // C4[1], coeff of eps^1, polynomial in n of order 4
1872  -64, -624, 4576, -6864, 3003, 135135,
1873  // C4[2], coeff of eps^5, polynomial in n of order 0
1874  8, 10725,
1875  // C4[2], coeff of eps^4, polynomial in n of order 1
1876  1856, -936, 225225,
1877  // C4[2], coeff of eps^3, polynomial in n of order 2
1878  -8448, 4992, -1144, 225225,
1879  // C4[2], coeff of eps^2, polynomial in n of order 3
1880  -1440, 4160, -4576, 1716, 225225,
1881  // C4[3], coeff of eps^5, polynomial in n of order 0
1882  -136, 63063,
1883  // C4[3], coeff of eps^4, polynomial in n of order 1
1884  1024, -208, 105105,
1885  // C4[3], coeff of eps^3, polynomial in n of order 2
1886  3584, -3328, 1144, 315315,
1887  // C4[4], coeff of eps^5, polynomial in n of order 0
1888  -128, 135135,
1889  // C4[4], coeff of eps^4, polynomial in n of order 1
1890  -2560, 832, 405405,
1891  // C4[5], coeff of eps^5, polynomial in n of order 0
1892  128, 99099,
1893  };
1894  int o = 0, k = 0;
1895  for (int l = 0; l < nC4_; ++l) { // l is index of C4[l]
1896  for (int j = nC4_ - 1; j >= l; --j) { // coeff of eps^j
1897  int m = nC4_ - j - 1; // order of polynomial in n
1898  _C4x[k++] = GeoMath.polyval(m, coeff, o, _n) / coeff[o + m + 1];
1899  o += m + 2;
1900  }
1901  }
1902  }
1903 }
net.sf.geographiclib.GeodesicMask.OUT_MASK
static final int OUT_MASK
Definition: GeodesicMask.java:33
net.sf.geographiclib.GeodesicData.a12
double a12
Definition: GeodesicData.java:54
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
net.sf.geographiclib.GeodesicData.M12
double M12
Definition: GeodesicData.java:62
net.sf.geographiclib.GeoMath.AngRound
static double AngRound(double x)
Definition: GeoMath.java:168
net.sf.geographiclib.Geodesic.InverseStartV.sig12
double sig12
Definition: Geodesic.java:1377
net.sf.geographiclib.Geodesic.C3f
void C3f(double eps, double c[])
Definition: Geodesic.java:1649
net.sf.geographiclib.GeodesicMask.GEODESICSCALE
static final int GEODESICSCALE
Definition: GeodesicMask.java:75
net.sf.geographiclib.Geodesic.Direct
GeodesicData Direct(double lat1, double lon1, double azi1, double s12)
Definition: Geodesic.java:313
net.sf.geographiclib.Geodesic.Line
GeodesicLine Line(double lat1, double lon1, double azi1)
Definition: Geodesic.java:1130
net.sf.geographiclib.Geodesic.nA2_
static final int nA2_
Definition: Geodesic.java:212
net.sf.geographiclib.Geodesic.InverseData
Definition: Geodesic.java:620
net.sf.geographiclib.GeodesicMask.AREA
static final int AREA
Definition: GeodesicMask.java:79
net.sf.geographiclib.GeodesicMask.DISTANCE_IN
static final int DISTANCE_IN
Definition: GeodesicMask.java:67
net.sf.geographiclib.GeodesicLine
Definition: GeodesicLine.java:93
net.sf.geographiclib.Geodesic._A3x
double _A3x[]
Definition: Geodesic.java:239
net.sf.geographiclib.Geodesic.DirectLine
GeodesicLine DirectLine(double lat1, double lon1, double azi1, double s12)
Definition: Geodesic.java:477
d
static const double d[K][N]
Definition: igam.h:11
net.sf.geographiclib.Geodesic.nC3x_
static final int nC3x_
Definition: Geodesic.java:217
net.sf.geographiclib.GeodesicData.azi1
double azi1
Definition: GeodesicData.java:33
net.sf.geographiclib.GeoMath.LatFix
static double LatFix(double x)
Definition: GeoMath.java:203
net.sf.geographiclib.Geodesic.Astroid
static double Astroid(double x, double y)
Definition: Geodesic.java:1324
net.sf.geographiclib.Geodesic.Direct
GeodesicData Direct(double lat1, double lon1, double azi1, double s12, int outmask)
Definition: Geodesic.java:337
net.sf.geographiclib.Geodesic.nA3_
static final int nA3_
Definition: Geodesic.java:214
net.sf.geographiclib.Geodesic.Lambda12
Lambda12V Lambda12(double sbet1, double cbet1, double dn1, double sbet2, double cbet2, double dn2, double salp1, double calp1, double slam120, double clam120, boolean diffp, double C1a[], double C2a[], double C3a[])
Definition: Geodesic.java:1556
r2
static const double r2
Definition: testSmartRangeFactor.cpp:32
net.sf.geographiclib.Geodesic.C1pf
static void C1pf(double eps, double c[])
Definition: Geodesic.java:1715
net.sf.geographiclib.Geodesic.LengthsV.LengthsV
LengthsV()
Definition: Geodesic.java:1257
net.sf.geographiclib.Geodesic.nC1_
static final int nC1_
Definition: Geodesic.java:210
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
net.sf.geographiclib.Geodesic
Definition: Geodesic.java:202
net.sf.geographiclib.Geodesic.LengthsV.s12b
double s12b
Definition: Geodesic.java:1256
net.sf.geographiclib.Geodesic.Inverse
GeodesicData Inverse(double lat1, double lon1, double lat2, double lon2)
Definition: Geodesic.java:615
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
net.sf.geographiclib.GeoMath.digits
static final int digits
Definition: GeoMath.java:21
net.sf.geographiclib.GeodesicMask.NONE
static final int NONE
Definition: GeodesicMask.java:38
net.sf.geographiclib.Geodesic.xthresh_
static final double xthresh_
Definition: Geodesic.java:235
net.sf.geographiclib.Geodesic.tol2_
static final double tol2_
Definition: Geodesic.java:232
T3
static const Pose3 T3(Rot3::Rodrigues(-90, 0, 0), Point3(1, 2, 3))
net.sf.geographiclib.Geodesic.tol0_
static final double tol0_
Definition: Geodesic.java:227
net.sf.geographiclib.GeodesicData.S12
double S12
Definition: GeodesicData.java:70
T
Eigen::Triplet< double > T
Definition: Tutorial_sparse_example.cpp:6
m0
static const DiscreteKey m0(M(0), 2)
net.sf.geographiclib.Geodesic.ArcDirectLine
GeodesicLine ArcDirectLine(double lat1, double lon1, double azi1, double a12)
Definition: Geodesic.java:523
net.sf.geographiclib.Geodesic.Lambda12V
Definition: Geodesic.java:1547
net.sf.geographiclib.GeodesicMask.ALL
static final int ALL
Definition: GeodesicMask.java:84
net.sf.geographiclib.Geodesic.C1f
static void C1f(double eps, double c[])
Definition: Geodesic.java:1687
net.sf.geographiclib.Geodesic.ArcDirect
GeodesicData ArcDirect(double lat1, double lon1, double azi1, double a12)
Definition: Geodesic.java:366
net.sf.geographiclib.GeoMath.atanh
static double atanh(double x)
Definition: GeoMath.java:90
net.sf.geographiclib.GeoMath.isfinite
static boolean isfinite(double x)
Definition: GeoMath.java:304
net.sf.geographiclib.Geodesic._a
double _a
Definition: Geodesic.java:237
result
Values result
Definition: OdometryOptimize.cpp:8
net.sf.geographiclib.Constants.WGS84_a
static final double WGS84_a
Definition: Constants.java:19
net.sf.geographiclib.GeoMath.AngDiff
static Pair AngDiff(double x, double y)
Definition: GeoMath.java:221
net.sf.geographiclib.Geodesic.A3f
double A3f(double eps)
Definition: Geodesic.java:1644
net.sf.geographiclib.GeoMath.sincosd
static Pair sincosd(double x)
Definition: GeoMath.java:240
y0
double y0(double x)
Definition: j0.c:220
net.sf.geographiclib.Geodesic.ArcDirect
GeodesicData ArcDirect(double lat1, double lon1, double azi1, double a12, int outmask)
Definition: Geodesic.java:391
net.sf.geographiclib.Geodesic.nA3x_
static final int nA3x_
Definition: Geodesic.java:215
net.sf.geographiclib.Geodesic.A3coeff
void A3coeff()
Definition: Geodesic.java:1782
net.sf.geographiclib.Geodesic.Flattening
double Flattening()
Definition: Geodesic.java:1194
n
int n
Definition: BiCGSTAB_simple.cpp:1
net.sf.geographiclib.Geodesic.nA1_
static final int nA1_
Definition: Geodesic.java:209
net.sf.geographiclib.GeoMath.hypot
static double hypot(double x, double y)
Definition: GeoMath.java:49
net.sf.geographiclib.Geodesic.Lengths
LengthsV Lengths(double eps, double sig12, double ssig1, double csig1, double dn1, double ssig2, double csig2, double dn2, double cbet1, double cbet2, int outmask, double C1a[], double C2a[])
Definition: Geodesic.java:1262
net.sf.geographiclib.Geodesic.C4f
void C4f(double eps, double c[])
Definition: Geodesic.java:1662
net.sf.geographiclib.Geodesic.tol1_
static final double tol1_
Definition: Geodesic.java:231
y1
double y1(double x)
Definition: j1.c:199
net.sf.geographiclib.Geodesic.WGS84
static final Geodesic WGS84
Definition: Geodesic.java:1207
net.sf.geographiclib.Geodesic.nC3_
static final int nC3_
Definition: Geodesic.java:216
net.sf.geographiclib.Constants.WGS84_f
static final double WGS84_f
Definition: Constants.java:23
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
net.sf.geographiclib.GeoMath.AngNormalize
static double AngNormalize(double x)
Definition: GeoMath.java:191
net.sf.geographiclib.GeoMath.atan2d
static double atan2d(double y, double x)
Definition: GeoMath.java:274
net.sf.geographiclib.Geodesic.Lambda12V.Lambda12V
Lambda12V()
Definition: Geodesic.java:1550
net.sf.geographiclib.Geodesic.InverseInt
InverseData InverseInt(double lat1, double lon1, double lat2, double lon2, int outmask)
Definition: Geodesic.java:629
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:1984
net.sf.geographiclib.Geodesic.C2f
static void C2f(double eps, double c[])
Definition: Geodesic.java:1754
net.sf.geographiclib.Geodesic.LengthsV
Definition: Geodesic.java:1255
net.sf.geographiclib.Geodesic.C3coeff
void C3coeff()
Definition: Geodesic.java:1806
net.sf.geographiclib.Geodesic.nC4_
static final int nC4_
Definition: Geodesic.java:218
net.sf.geographiclib.Geodesic.A1m1f
static double A1m1f(double eps)
Definition: Geodesic.java:1676
l
static const Line3 l(Rot3(), 1, 1)
net.sf.geographiclib.GeodesicData.lon1
double lon1
Definition: GeodesicData.java:29
net.sf.geographiclib.Geodesic.InverseLine
GeodesicLine InverseLine(double lat1, double lon1, double lat2, double lon2, int caps)
Definition: Geodesic.java:1102
A2
static const double A2[]
Definition: expn.h:7
net.sf.geographiclib.Constants
Definition: Constants.java:15
net.sf.geographiclib.GeoMath.polyval
static double polyval(int N, double p[], int s, double x)
Definition: GeoMath.java:162
net.sf.geographiclib.GeoMath.epsilon
static final double epsilon
Definition: GeoMath.java:26
net.sf.geographiclib.Geodesic.EllipsoidArea
double EllipsoidArea()
Definition: Geodesic.java:1201
net.sf.geographiclib.Geodesic.Direct
GeodesicData Direct(double lat1, double lon1, double azi1, boolean arcmode, double s12_a12, int outmask)
Definition: Geodesic.java:452
net.sf.geographiclib.GeodesicMask.LONG_UNROLL
static final int LONG_UNROLL
Definition: GeodesicMask.java:88
r3
static const double r3
Definition: testSmartRangeFactor.cpp:32
net.sf.geographiclib.GeodesicData.M21
double M21
Definition: GeodesicData.java:66
net.sf.geographiclib.Geodesic.nC2_
static final int nC2_
Definition: Geodesic.java:213
net.sf.geographiclib.Geodesic.A2m1f
static double A2m1f(double eps)
Definition: Geodesic.java:1743
net.sf.geographiclib.GeodesicMask
Definition: GeodesicMask.java:23
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
net.sf.geographiclib.Pair
Definition: Pair.java:15
Eigen::Triplet< double >
net.sf.geographiclib.Geodesic.DirectLine
GeodesicLine DirectLine(double lat1, double lon1, double azi1, double s12, int caps)
Definition: Geodesic.java:502
net.sf.geographiclib.Geodesic.InverseData.g
GeodesicData g
Definition: Geodesic.java:621
net.sf.geographiclib.GeodesicData
Definition: GeodesicData.java:21
net.sf.geographiclib.Geodesic.Inverse
GeodesicData Inverse(double lat1, double lon1, double lat2, double lon2, int outmask)
Definition: Geodesic.java:1049
net.sf.geographiclib.GeoMath.norm
static Pair norm(double sinx, double cosx)
Definition: GeoMath.java:120
net.sf.geographiclib.GeodesicData.lon2
double lon2
Definition: GeodesicData.java:41
y
Scalar * y
Definition: level1_cplx_impl.h:124
tree::f
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Definition: testExpression.cpp:218
net.sf.geographiclib.GeodesicData.lat1
double lat1
Definition: GeodesicData.java:25
net.sf.geographiclib.Geodesic.nC4x_
static final int nC4x_
Definition: Geodesic.java:219
net.sf.geographiclib.GeodesicMask.AZIMUTH
static final int AZIMUTH
Definition: GeodesicMask.java:53
net.sf.geographiclib.GeodesicMask.STANDARD
static final int STANDARD
Definition: GeodesicMask.java:61
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
net.sf.geographiclib.GeoMath
Definition: GeoMath.java:16
net.sf.geographiclib.Geodesic.InverseData.salp1
double salp1
Definition: Geodesic.java:622
net.sf.geographiclib.GeodesicData.lat2
double lat2
Definition: GeodesicData.java:37
net.sf.geographiclib.Geodesic.InverseData.InverseData
InverseData()
Definition: Geodesic.java:623
net.sf.geographiclib.GeodesicData.azi2
double azi2
Definition: GeodesicData.java:45
net.sf.geographiclib.GeographicErr
Definition: GeographicErr.java:16
net.sf.geographiclib.Geodesic.InverseStartV
Definition: Geodesic.java:1376
net.sf.geographiclib.Geodesic.MajorRadius
double MajorRadius()
Definition: Geodesic.java:1188
net.sf.geographiclib.Geodesic._n
double _n
Definition: Geodesic.java:238
net.sf.geographiclib.Geodesic.tolb_
static final double tolb_
Definition: Geodesic.java:234
net.sf.geographiclib.GeodesicMask.REDUCEDLENGTH
static final int REDUCEDLENGTH
Definition: GeodesicMask.java:71
net.sf.geographiclib.Geodesic.SinCosSeries
static double SinCosSeries(boolean sinp, double sinx, double cosx, double c[])
Definition: Geodesic.java:1229
p
float * p
Definition: Tutorial_Map_using.cpp:9
net.sf.geographiclib.GeoMath.sq
static double sq(double x)
Definition: GeoMath.java:39
A1
static const double A1[]
Definition: expn.h:6
net.sf.geographiclib.Geodesic.ArcDirectLine
GeodesicLine ArcDirectLine(double lat1, double lon1, double azi1, double a12, int caps)
Definition: Geodesic.java:549
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
net.sf.geographiclib.Geodesic.GEOGRAPHICLIB_GEODESIC_ORDER
static final int GEOGRAPHICLIB_GEODESIC_ORDER
Definition: Geodesic.java:207
net.sf.geographiclib.Geodesic.maxit2_
static final int maxit2_
Definition: Geodesic.java:221
net.sf.geographiclib.Geodesic.Geodesic
Geodesic(double a, double f)
Definition: Geodesic.java:250
net.sf.geographiclib.Geodesic.InverseStart
InverseStartV InverseStart(double sbet1, double cbet1, double dn1, double sbet2, double cbet2, double dn2, double lam12, double slam12, double clam12, double C1a[], double C2a[])
Definition: Geodesic.java:1387
net.sf.geographiclib.GeoMath.cbrt
static double cbrt(double x)
Definition: GeoMath.java:115
net.sf.geographiclib.Geodesic.InverseLine
GeodesicLine InverseLine(double lat1, double lon1, double lat2, double lon2)
Definition: Geodesic.java:1077
net.sf.geographiclib.Geodesic.Lambda12V.lam12
double lam12
Definition: Geodesic.java:1548
net.sf.geographiclib.GeodesicData.s12
double s12
Definition: GeodesicData.java:49
net.sf.geographiclib.GeodesicMask.DISTANCE
static final int DISTANCE
Definition: GeodesicMask.java:57
net.sf.geographiclib.Geodesic.tiny_
static final double tiny_
Definition: Geodesic.java:226
A4
static const double A4[]
Definition: expn.h:9
align_3::t
Point2 t(10, 10)
net.sf.geographiclib.GeodesicData.m12
double m12
Definition: GeodesicData.java:58
net.sf.geographiclib.Geodesic.maxit1_
static final int maxit1_
Definition: Geodesic.java:220
net.sf.geographiclib.GeoMath.min
static final double min
Definition: GeoMath.java:31
net.sf.geographiclib.Geodesic.InverseStartV.InverseStartV
InverseStartV()
Definition: Geodesic.java:1382
net.sf.geographiclib.Geodesic.Line
GeodesicLine Line(double lat1, double lon1, double azi1, int caps)
Definition: Geodesic.java:1180
net.sf.geographiclib.Geodesic.GenDirectLine
GeodesicLine GenDirectLine(double lat1, double lon1, double azi1, boolean arcmode, double s12_a12, int caps)
Definition: Geodesic.java:577
net.sf.geographiclib.Geodesic.nC1p_
static final int nC1p_
Definition: Geodesic.java:211
S
DiscreteKey S(1, 2)
net.sf.geographiclib.Geodesic.C4coeff
void C4coeff()
Definition: Geodesic.java:1849


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:02:22