GeodesicLine.java
Go to the documentation of this file.
1 
8 package net.sf.geographiclib;
9 
93 public class GeodesicLine {
94 
95  private static final int nC1_ = Geodesic.nC1_;
96  private static final int nC1p_ = Geodesic.nC1p_;
97  private static final int nC2_ = Geodesic.nC2_;
98  private static final int nC3_ = Geodesic.nC3_;
99  private static final int nC4_ = Geodesic.nC4_;
100 
101  private double _lat1, _lon1, _azi1;
102  private double _a, _f, _b, _c2, _f1, _salp0, _calp0, _k2,
103  _salp1, _calp1, _ssig1, _csig1, _dn1, _stau1, _ctau1, _somg1, _comg1,
104  _A1m1, _A2m1, _A3c, _B11, _B21, _B31, _A4, _B41;
105  private double _a13, _s13;
106  // index zero elements of _C1a, _C1pa, _C2a, _C3a are unused
107  private double _C1a[], _C1pa[], _C2a[], _C3a[],
108  _C4a[]; // all the elements of _C4a are used
109  private int _caps;
110 
128  double lat1, double lon1, double azi1) {
129  this(g, lat1, lon1, azi1, GeodesicMask.ALL);
130  }
131 
178  double lat1, double lon1, double azi1,
179  int caps) {
180  azi1 = GeoMath.AngNormalize(azi1);
181  double salp1, calp1;
182  // Guard against underflow in salp0
183  { Pair p = GeoMath.sincosd(GeoMath.AngRound(azi1));
184  salp1 = p.first; calp1 = p.second; }
185  LineInit(g, lat1, lon1, azi1, salp1, calp1, caps);
186  }
187 
188  private void LineInit(Geodesic g,
189  double lat1, double lon1,
190  double azi1, double salp1, double calp1,
191  int caps) {
192  _a = g._a;
193  _f = g._f;
194  _b = g._b;
195  _c2 = g._c2;
196  _f1 = g._f1;
197  // Always allow latitude and azimuth and unrolling the longitude
198  _caps = caps | GeodesicMask.LATITUDE | GeodesicMask.AZIMUTH |
200 
201  _lat1 = GeoMath.LatFix(lat1);
202  _lon1 = lon1;
203  _azi1 = azi1; _salp1 = salp1; _calp1 = calp1;
204  double cbet1, sbet1;
205  { Pair p = GeoMath.sincosd(GeoMath.AngRound(_lat1));
206  sbet1 = _f1 * p.first; cbet1 = p.second; }
207  // Ensure cbet1 = +epsilon at poles
208  { Pair p = GeoMath.norm(sbet1, cbet1);
209  sbet1 = p.first; cbet1 = Math.max(Geodesic.tiny_, p.second); }
210  _dn1 = Math.sqrt(1 + g._ep2 * GeoMath.sq(sbet1));
211 
212  // Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0),
213  _salp0 = _salp1 * cbet1; // alp0 in [0, pi/2 - |bet1|]
214  // Alt: calp0 = hypot(sbet1, calp1 * cbet1). The following
215  // is slightly better (consider the case salp1 = 0).
216  _calp0 = GeoMath.hypot(_calp1, _salp1 * sbet1);
217  // Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1).
218  // sig = 0 is nearest northward crossing of equator.
219  // With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line).
220  // With bet1 = pi/2, alp1 = -pi, sig1 = pi/2
221  // With bet1 = -pi/2, alp1 = 0 , sig1 = -pi/2
222  // Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1).
223  // With alp0 in (0, pi/2], quadrants for sig and omg coincide.
224  // No atan2(0,0) ambiguity at poles since cbet1 = +epsilon.
225  // With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi.
226  _ssig1 = sbet1; _somg1 = _salp0 * sbet1;
227  _csig1 = _comg1 = sbet1 != 0 || _calp1 != 0 ? cbet1 * _calp1 : 1;
228  { Pair p = GeoMath.norm(_ssig1, _csig1);
229  _ssig1 = p.first; _csig1 = p.second; } // sig1 in (-pi, pi]
230  // GeoMath.norm(_somg1, _comg1); -- don't need to normalize!
231 
232  _k2 = GeoMath.sq(_calp0) * g._ep2;
233  double eps = _k2 / (2 * (1 + Math.sqrt(1 + _k2)) + _k2);
234 
235  if ((_caps & GeodesicMask.CAP_C1) != 0) {
236  _A1m1 = Geodesic.A1m1f(eps);
237  _C1a = new double[nC1_ + 1];
238  Geodesic.C1f(eps, _C1a);
239  _B11 = Geodesic.SinCosSeries(true, _ssig1, _csig1, _C1a);
240  double s = Math.sin(_B11), c = Math.cos(_B11);
241  // tau1 = sig1 + B11
242  _stau1 = _ssig1 * c + _csig1 * s;
243  _ctau1 = _csig1 * c - _ssig1 * s;
244  // Not necessary because C1pa reverts C1a
245  // _B11 = -SinCosSeries(true, _stau1, _ctau1, _C1pa, nC1p_);
246  }
247 
248  if ((_caps & GeodesicMask.CAP_C1p) != 0) {
249  _C1pa = new double[nC1p_ + 1];
250  Geodesic.C1pf(eps, _C1pa);
251  }
252 
253  if ((_caps & GeodesicMask.CAP_C2) != 0) {
254  _C2a = new double[nC2_ + 1];
255  _A2m1 = Geodesic.A2m1f(eps);
256  Geodesic.C2f(eps, _C2a);
257  _B21 = Geodesic.SinCosSeries(true, _ssig1, _csig1, _C2a);
258  }
259 
260  if ((_caps & GeodesicMask.CAP_C3) != 0) {
261  _C3a = new double[nC3_];
262  g.C3f(eps, _C3a);
263  _A3c = -_f * _salp0 * g.A3f(eps);
264  _B31 = Geodesic.SinCosSeries(true, _ssig1, _csig1, _C3a);
265  }
266 
267  if ((_caps & GeodesicMask.CAP_C4) != 0) {
268  _C4a = new double[nC4_];
269  g.C4f(eps, _C4a);
270  // Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0)
271  _A4 = GeoMath.sq(_a) * _calp0 * _salp0 * g._e2;
272  _B41 = Geodesic.SinCosSeries(false, _ssig1, _csig1, _C4a);
273  }
274  }
275 
277  double lat1, double lon1,
278  double azi1, double salp1, double calp1,
279  int caps, boolean arcmode, double s13_a13) {
280  LineInit(g, lat1, lon1, azi1, salp1, calp1, caps);
281  GenSetDistance(arcmode, s13_a13);
282  }
283 
292  private GeodesicLine() { _caps = 0; }
293 
311  public GeodesicData Position(double s12) {
312  return Position(false, s12, GeodesicMask.STANDARD);
313  }
334  public GeodesicData Position(double s12, int outmask) {
335  return Position(false, s12, outmask);
336  }
337 
355  public GeodesicData ArcPosition(double a12) {
356  return Position(true, a12, GeodesicMask.STANDARD);
357  }
374  public GeodesicData ArcPosition(double a12, int outmask) {
375  return Position(true, a12, outmask);
376  }
377 
423  public GeodesicData Position(boolean arcmode, double s12_a12,
424  int outmask) {
425  outmask &= _caps & GeodesicMask.OUT_MASK;
426  GeodesicData r = new GeodesicData();
427  if (!( Init() &&
428  (arcmode ||
429  (_caps & (GeodesicMask.OUT_MASK & GeodesicMask.DISTANCE_IN)) != 0)
430  ))
431  // Uninitialized or impossible distance calculation requested
432  return r;
433  r.lat1 = _lat1; r.azi1 = _azi1;
434  r.lon1 = ((outmask & GeodesicMask.LONG_UNROLL) != 0) ? _lon1 :
435  GeoMath.AngNormalize(_lon1);
436 
437  // Avoid warning about uninitialized B12.
438  double sig12, ssig12, csig12, B12 = 0, AB1 = 0;
439  if (arcmode) {
440  // Interpret s12_a12 as spherical arc length
441  r.a12 = s12_a12;
442  sig12 = Math.toRadians(s12_a12);
443  { Pair p = GeoMath.sincosd(s12_a12);
444  ssig12 = p.first; csig12 = p.second; }
445  } else {
446  // Interpret s12_a12 as distance
447  r.s12 = s12_a12;
448  double
449  tau12 = s12_a12 / (_b * (1 + _A1m1)),
450  s = Math.sin(tau12),
451  c = Math.cos(tau12);
452  // tau2 = tau1 + tau12
453  B12 = - Geodesic.SinCosSeries(true,
454  _stau1 * c + _ctau1 * s,
455  _ctau1 * c - _stau1 * s,
456  _C1pa);
457  sig12 = tau12 - (B12 - _B11);
458  ssig12 = Math.sin(sig12); csig12 = Math.cos(sig12);
459  if (Math.abs(_f) > 0.01) {
460  // Reverted distance series is inaccurate for |f| > 1/100, so correct
461  // sig12 with 1 Newton iteration. The following table shows the
462  // approximate maximum error for a = WGS_a() and various f relative to
463  // GeodesicExact.
464  // erri = the error in the inverse solution (nm)
465  // errd = the error in the direct solution (series only) (nm)
466  // errda = the error in the direct solution
467  // (series + 1 Newton) (nm)
468  //
469  // f erri errd errda
470  // -1/5 12e6 1.2e9 69e6
471  // -1/10 123e3 12e6 765e3
472  // -1/20 1110 108e3 7155
473  // -1/50 18.63 200.9 27.12
474  // -1/100 18.63 23.78 23.37
475  // -1/150 18.63 21.05 20.26
476  // 1/150 22.35 24.73 25.83
477  // 1/100 22.35 25.03 25.31
478  // 1/50 29.80 231.9 30.44
479  // 1/20 5376 146e3 10e3
480  // 1/10 829e3 22e6 1.5e6
481  // 1/5 157e6 3.8e9 280e6
482  double
483  ssig2 = _ssig1 * csig12 + _csig1 * ssig12,
484  csig2 = _csig1 * csig12 - _ssig1 * ssig12;
485  B12 = Geodesic.SinCosSeries(true, ssig2, csig2, _C1a);
486  double serr = (1 + _A1m1) * (sig12 + (B12 - _B11)) - s12_a12 / _b;
487  sig12 = sig12 - serr / Math.sqrt(1 + _k2 * GeoMath.sq(ssig2));
488  ssig12 = Math.sin(sig12); csig12 = Math.cos(sig12);
489  // Update B12 below
490  }
491  r.a12 = Math.toDegrees(sig12);
492  }
493 
494  double ssig2, csig2, sbet2, cbet2, salp2, calp2;
495  // sig2 = sig1 + sig12
496  ssig2 = _ssig1 * csig12 + _csig1 * ssig12;
497  csig2 = _csig1 * csig12 - _ssig1 * ssig12;
498  double dn2 = Math.sqrt(1 + _k2 * GeoMath.sq(ssig2));
499  if ((outmask & (GeodesicMask.DISTANCE | GeodesicMask.REDUCEDLENGTH |
500  GeodesicMask.GEODESICSCALE)) != 0) {
501  if (arcmode || Math.abs(_f) > 0.01)
502  B12 = Geodesic.SinCosSeries(true, ssig2, csig2, _C1a);
503  AB1 = (1 + _A1m1) * (B12 - _B11);
504  }
505  // sin(bet2) = cos(alp0) * sin(sig2)
506  sbet2 = _calp0 * ssig2;
507  // Alt: cbet2 = hypot(csig2, salp0 * ssig2);
508  cbet2 = GeoMath.hypot(_salp0, _calp0 * csig2);
509  if (cbet2 == 0)
510  // I.e., salp0 = 0, csig2 = 0. Break the degeneracy in this case
511  cbet2 = csig2 = Geodesic.tiny_;
512  // tan(alp0) = cos(sig2)*tan(alp2)
513  salp2 = _salp0; calp2 = _calp0 * csig2; // No need to normalize
514 
515  if ((outmask & GeodesicMask.DISTANCE) != 0 && arcmode)
516  r.s12 = _b * ((1 + _A1m1) * sig12 + AB1);
517 
518  if ((outmask & GeodesicMask.LONGITUDE) != 0) {
519  // tan(omg2) = sin(alp0) * tan(sig2)
520  double somg2 = _salp0 * ssig2, comg2 = csig2, // No need to normalize
521  E = GeoMath.copysign(1, _salp0); // east or west going?
522  // omg12 = omg2 - omg1
523  double omg12 = ((outmask & GeodesicMask.LONG_UNROLL) != 0)
524  ? E * (sig12
525  - (Math.atan2( ssig2, csig2) - Math.atan2( _ssig1, _csig1))
526  + (Math.atan2(E*somg2, comg2) - Math.atan2(E*_somg1, _comg1)))
527  : Math.atan2(somg2 * _comg1 - comg2 * _somg1,
528  comg2 * _comg1 + somg2 * _somg1);
529  double lam12 = omg12 + _A3c *
530  ( sig12 + (Geodesic.SinCosSeries(true, ssig2, csig2, _C3a)
531  - _B31));
532  double lon12 = Math.toDegrees(lam12);
533  r.lon2 = ((outmask & GeodesicMask.LONG_UNROLL) != 0) ? _lon1 + lon12 :
535  }
536 
537  if ((outmask & GeodesicMask.LATITUDE) != 0)
538  r.lat2 = GeoMath.atan2d(sbet2, _f1 * cbet2);
539 
540  if ((outmask & GeodesicMask.AZIMUTH) != 0)
541  r.azi2 = GeoMath.atan2d(salp2, calp2);
542 
543  if ((outmask &
545  double
546  B22 = Geodesic.SinCosSeries(true, ssig2, csig2, _C2a),
547  AB2 = (1 + _A2m1) * (B22 - _B21),
548  J12 = (_A1m1 - _A2m1) * sig12 + (AB1 - AB2);
549  if ((outmask & GeodesicMask.REDUCEDLENGTH) != 0)
550  // Add parens around (_csig1 * ssig2) and (_ssig1 * csig2) to ensure
551  // accurate cancellation in the case of coincident points.
552  r.m12 = _b * ((dn2 * (_csig1 * ssig2) - _dn1 * (_ssig1 * csig2))
553  - _csig1 * csig2 * J12);
554  if ((outmask & GeodesicMask.GEODESICSCALE) != 0) {
555  double t = _k2 * (ssig2 - _ssig1) * (ssig2 + _ssig1) / (_dn1 + dn2);
556  r.M12 = csig12 + (t * ssig2 - csig2 * J12) * _ssig1 / _dn1;
557  r.M21 = csig12 - (t * _ssig1 - _csig1 * J12) * ssig2 / dn2;
558  }
559  }
560 
561  if ((outmask & GeodesicMask.AREA) != 0) {
562  double
563  B42 = Geodesic.SinCosSeries(false, ssig2, csig2, _C4a);
564  double salp12, calp12;
565  if (_calp0 == 0 || _salp0 == 0) {
566  // alp12 = alp2 - alp1, used in atan2 so no need to normalize
567  salp12 = salp2 * _calp1 - calp2 * _salp1;
568  calp12 = calp2 * _calp1 + salp2 * _salp1;
569  } else {
570  // tan(alp) = tan(alp0) * sec(sig)
571  // tan(alp2-alp1) = (tan(alp2) -tan(alp1)) / (tan(alp2)*tan(alp1)+1)
572  // = calp0 * salp0 * (csig1-csig2) / (salp0^2 + calp0^2 * csig1*csig2)
573  // If csig12 > 0, write
574  // csig1 - csig2 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1)
575  // else
576  // csig1 - csig2 = csig1 * (1 - csig12) + ssig12 * ssig1
577  // No need to normalize
578  salp12 = _calp0 * _salp0 *
579  (csig12 <= 0 ? _csig1 * (1 - csig12) + ssig12 * _ssig1 :
580  ssig12 * (_csig1 * ssig12 / (1 + csig12) + _ssig1));
581  calp12 = GeoMath.sq(_salp0) + GeoMath.sq(_calp0) * _csig1 * csig2;
582  }
583  r.S12 = _c2 * Math.atan2(salp12, calp12) + _A4 * (B42 - _B41);
584  }
585 
586  return r;
587  }
588 
598  public void SetDistance(double s13) {
599  _s13 = s13;
600  GeodesicData g = Position(false, _s13, 0);
601  _a13 = g.a12;
602  }
603 
613  void SetArc(double a13) {
614  _a13 = a13;
616  _s13 = g.s12;
617  }
618 
630  public void GenSetDistance(boolean arcmode, double s13_a13) {
631  if (arcmode)
632  SetArc(s13_a13);
633  else
634  SetDistance(s13_a13);
635  }
636 
640  private boolean Init() { return _caps != 0; }
641 
645  public double Latitude()
646  { return Init() ? _lat1 : Double.NaN; }
647 
651  public double Longitude()
652  { return Init() ? _lon1 : Double.NaN; }
653 
657  public double Azimuth()
658  { return Init() ? _azi1 : Double.NaN; }
659 
664  public Pair AzimuthCosines() {
665  return new Pair(Init() ? _salp1 : Double.NaN,
666  Init() ? _calp1 : Double.NaN);
667  }
668 
673  public double EquatorialAzimuth() {
674  return Init() ?
675  GeoMath.atan2d(_salp0, _calp0) : Double.NaN;
676  }
677 
683  return new Pair(Init() ? _salp0 : Double.NaN,
684  Init() ? _calp0 : Double.NaN);
685  }
686 
691  public double EquatorialArc() {
692  return Init() ?
693  GeoMath.atan2d(_ssig1, _csig1) : Double.NaN;
694  }
695 
700  public double MajorRadius()
701  { return Init() ? _a : Double.NaN; }
702 
707  public double Flattening()
708  { return Init() ? _f : Double.NaN; }
709 
714  public int Capabilities() { return _caps; }
715 
720  public boolean Capabilities(int testcaps) {
721  testcaps &= GeodesicMask.OUT_ALL;
722  return (_caps & testcaps) == testcaps;
723  }
724 
733  public double GenDistance(boolean arcmode)
734  { return Init() ? (arcmode ? _a13 : _s13) : Double.NaN; }
735 
739  public double Distance() { return GenDistance(false); }
740 
744  public double Arc() { return GenDistance(true); }
745 
746 }
static double A1m1f(double eps)
Definition: Geodesic.java:1676
GeodesicData Position(boolean arcmode, double s12_a12, int outmask)
GeodesicData Position(double s12)
static void C1f(double eps, double c[])
Definition: Geodesic.java:1687
boolean Capabilities(int testcaps)
static double hypot(double x, double y)
Definition: GeoMath.java:49
GeodesicData ArcPosition(double a12)
static void C2f(double eps, double c[])
Definition: Geodesic.java:1754
static final double tiny_
Definition: Geodesic.java:226
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
void GenSetDistance(boolean arcmode, double s13_a13)
void C3f(double eps, double c[])
Definition: Geodesic.java:1649
void g(const string &key, int i)
Definition: testBTree.cpp:41
void C4f(double eps, double c[])
Definition: Geodesic.java:1662
static double atan2d(double y, double x)
Definition: GeoMath.java:274
static double A2m1f(double eps)
Definition: Geodesic.java:1743
static Pair norm(double sinx, double cosx)
Definition: GeoMath.java:120
GeodesicData ArcPosition(double a12, int outmask)
GeodesicLine(Geodesic g, double lat1, double lon1, double azi1, int caps)
GeodesicData Position(double s12, int outmask)
RealScalar s
static double SinCosSeries(boolean sinp, double sinx, double cosx, double c[])
Definition: Geodesic.java:1229
static double LatFix(double x)
Definition: GeoMath.java:203
static double copysign(double x, double y)
Definition: GeoMath.java:104
DiscreteKey E(5, 2)
static Pair sincosd(double x)
Definition: GeoMath.java:240
float * p
static double AngRound(double x)
Definition: GeoMath.java:168
void LineInit(Geodesic g, double lat1, double lon1, double azi1, double salp1, double calp1, int caps)
static double sq(double x)
Definition: GeoMath.java:39
static void C1pf(double eps, double c[])
Definition: Geodesic.java:1715
static double AngNormalize(double x)
Definition: GeoMath.java:191
double A3f(double eps)
Definition: Geodesic.java:1644
double GenDistance(boolean arcmode)
GeodesicLine(Geodesic g, double lat1, double lon1, double azi1, double salp1, double calp1, int caps, boolean arcmode, double s13_a13)
GeodesicLine(Geodesic g, double lat1, double lon1, double azi1)
Point2 t(10, 10)


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:18