GeodesicTest.java
Go to the documentation of this file.
1 package net.sf.geographiclib.test;
2 
3 import static org.junit.Assert.assertEquals;
4 import static org.junit.Assert.assertTrue;
5 import org.junit.Test;
6 import net.sf.geographiclib.*;
7 
8 public class GeodesicTest {
9 
10  private static boolean isNaN(double x) { return x != x; }
11  private static final PolygonArea polygon =
12  new PolygonArea(Geodesic.WGS84, false);
13  private static final PolygonArea polyline =
14  new PolygonArea(Geodesic.WGS84, true);
15 
16  private static PolygonResult Planimeter(double points[][]) {
17  polygon.Clear();
18  for (int i = 0; i < points.length; ++i) {
19  polygon.AddPoint(points[i][0], points[i][1]);
20  }
21  return polygon.Compute(false, true);
22  }
23 
24  private static PolygonResult PolyLength(double points[][]) {
25  polyline.Clear();
26  for (int i = 0; i < points.length; ++i) {
27  polyline.AddPoint(points[i][0], points[i][1]);
28  }
29  return polyline.Compute(false, true);
30  }
31 
32  private static final double testcases[][] = {
33  {35.60777, -139.44815, 111.098748429560326,
34  -11.17491, -69.95921, 129.289270889708762,
35  8935244.5604818305, 80.50729714281974, 6273170.2055303837,
36  0.16606318447386067, 0.16479116945612937, 12841384694976.432},
37  {55.52454, 106.05087, 22.020059880982801,
38  77.03196, 197.18234, 109.112041110671519,
39  4105086.1713924406, 36.892740690445894, 3828869.3344387607,
40  0.80076349608092607, 0.80101006984201008, 61674961290615.615},
41  {-21.97856, 142.59065, -32.44456876433189,
42  41.84138, 98.56635, -41.84359951440466,
43  8394328.894657671, 75.62930491011522, 6161154.5773110616,
44  0.24816339233950381, 0.24930251203627892, -6637997720646.717},
45  {-66.99028, 112.2363, 173.73491240878403,
46  -12.70631, 285.90344, 2.512956620913668,
47  11150344.2312080241, 100.278634181155759, 6289939.5670446687,
48  -0.17199490274700385, -0.17722569526345708, -121287239862139.744},
49  {-17.42761, 173.34268, -159.033557661192928,
50  -15.84784, 5.93557, -20.787484651536988,
51  16076603.1631180673, 144.640108810286253, 3732902.1583877189,
52  -0.81273638700070476, -0.81299800519154474, 97825992354058.708},
53  {32.84994, 48.28919, 150.492927788121982,
54  -56.28556, 202.29132, 48.113449399816759,
55  16727068.9438164461, 150.565799985466607, 3147838.1910180939,
56  -0.87334918086923126, -0.86505036767110637, -72445258525585.010},
57  {6.96833, 52.74123, 92.581585386317712,
58  -7.39675, 206.17291, 90.721692165923907,
59  17102477.2496958388, 154.147366239113561, 2772035.6169917581,
60  -0.89991282520302447, -0.89986892177110739, -1311796973197.995},
61  {-50.56724, -16.30485, -105.439679907590164,
62  -33.56571, -94.97412, -47.348547835650331,
63  6455670.5118668696, 58.083719495371259, 5409150.7979815838,
64  0.53053508035997263, 0.52988722644436602, 41071447902810.047},
65  {-58.93002, -8.90775, 140.965397902500679,
66  -8.91104, 133.13503, 19.255429433416599,
67  11756066.0219864627, 105.755691241406877, 6151101.2270708536,
68  -0.26548622269867183, -0.27068483874510741, -86143460552774.735},
69  {-68.82867, -74.28391, 93.774347763114881,
70  -50.63005, -8.36685, 34.65564085411343,
71  3956936.926063544, 35.572254987389284, 3708890.9544062657,
72  0.81443963736383502, 0.81420859815358342, -41845309450093.787},
73  {-10.62672, -32.0898, -86.426713286747751,
74  5.883, -134.31681, -80.473780971034875,
75  11470869.3864563009, 103.387395634504061, 6184411.6622659713,
76  -0.23138683500430237, -0.23155097622286792, 4198803992123.548},
77  {-21.76221, 166.90563, 29.319421206936428,
78  48.72884, 213.97627, 43.508671946410168,
79  9098627.3986554915, 81.963476716121964, 6299240.9166992283,
80  0.13965943368590333, 0.14152969707656796, 10024709850277.476},
81  {-19.79938, -174.47484, 71.167275780171533,
82  -11.99349, -154.35109, 65.589099775199228,
83  2319004.8601169389, 20.896611684802389, 2267960.8703918325,
84  0.93427001867125849, 0.93424887135032789, -3935477535005.785},
85  {-11.95887, -116.94513, 92.712619830452549,
86  4.57352, 7.16501, 78.64960934409585,
87  13834722.5801401374, 124.688684161089762, 5228093.177931598,
88  -0.56879356755666463, -0.56918731952397221, -9919582785894.853},
89  {-87.85331, 85.66836, -65.120313040242748,
90  66.48646, 16.09921, -4.888658719272296,
91  17286615.3147144645, 155.58592449699137, 2635887.4729110181,
92  -0.90697975771398578, -0.91095608883042767, 42667211366919.534},
93  {1.74708, 128.32011, -101.584843631173858,
94  -11.16617, 11.87109, -86.325793296437476,
95  12942901.1241347408, 116.650512484301857, 5682744.8413270572,
96  -0.44857868222697644, -0.44824490340007729, 10763055294345.653},
97  {-25.72959, -144.90758, -153.647468693117198,
98  -57.70581, -269.17879, -48.343983158876487,
99  9413446.7452453107, 84.664533838404295, 6356176.6898881281,
100  0.09492245755254703, 0.09737058264766572, 74515122850712.444},
101  {-41.22777, 122.32875, 14.285113402275739,
102  -7.57291, 130.37946, 10.805303085187369,
103  3812686.035106021, 34.34330804743883, 3588703.8812128856,
104  0.82605222593217889, 0.82572158200920196, -2456961531057.857},
105  {11.01307, 138.25278, 79.43682622782374,
106  6.62726, 247.05981, 103.708090215522657,
107  11911190.819018408, 107.341669954114577, 6070904.722786735,
108  -0.29767608923657404, -0.29785143390252321, 17121631423099.696},
109  {-29.47124, 95.14681, -163.779130441688382,
110  -27.46601, -69.15955, -15.909335945554969,
111  13487015.8381145492, 121.294026715742277, 5481428.9945736388,
112  -0.51527225545373252, -0.51556587964721788, 104679964020340.318}};
113 
114  @Test
115  public void InverseCheck() {
116  for (int i = 0; i < testcases.length; ++i) {
117  double
118  lat1 = testcases[i][0], lon1 = testcases[i][1], azi1 = testcases[i][2],
119  lat2 = testcases[i][3], lon2 = testcases[i][4], azi2 = testcases[i][5],
120  s12 = testcases[i][6], a12 = testcases[i][7], m12 = testcases[i][8],
121  M12 = testcases[i][9], M21 = testcases[i][10], S12 = testcases[i][11];
122  GeodesicData inv = Geodesic.WGS84.Inverse(lat1, lon1, lat2, lon2,
123  GeodesicMask.ALL |
125  assertEquals(lon2, inv.lon2, 1e-13);
126  assertEquals(azi1, inv.azi1, 1e-13);
127  assertEquals(azi2, inv.azi2, 1e-13);
128  assertEquals(s12, inv.s12, 1e-8);
129  assertEquals(a12, inv.a12, 1e-13);
130  assertEquals(m12, inv.m12, 1e-8);
131  assertEquals(M12, inv.M12, 1e-15);
132  assertEquals(M21, inv.M21, 1e-15);
133  assertEquals(S12, inv.S12, 0.1);
134  }
135  }
136 
137  @Test
138  public void DirectCheck() {
139  for (int i = 0; i < testcases.length; ++i) {
140  double
141  lat1 = testcases[i][0], lon1 = testcases[i][1], azi1 = testcases[i][2],
142  lat2 = testcases[i][3], lon2 = testcases[i][4], azi2 = testcases[i][5],
143  s12 = testcases[i][6], a12 = testcases[i][7], m12 = testcases[i][8],
144  M12 = testcases[i][9], M21 = testcases[i][10], S12 = testcases[i][11];
145  GeodesicData dir = Geodesic.WGS84.Direct(lat1, lon1, azi1, s12,
146  GeodesicMask.ALL |
148  assertEquals(lat2, dir.lat2, 1e-13);
149  assertEquals(lon2, dir.lon2, 1e-13);
150  assertEquals(azi2, dir.azi2, 1e-13);
151  assertEquals(a12, dir.a12, 1e-13);
152  assertEquals(m12, dir.m12, 1e-8);
153  assertEquals(M12, dir.M12, 1e-15);
154  assertEquals(M21, dir.M21, 1e-15);
155  assertEquals(S12, dir.S12, 0.1);
156  }
157  }
158 
159  @Test
160  public void ArcDirectCheck() {
161  for (int i = 0; i < testcases.length; ++i) {
162  double
163  lat1 = testcases[i][0], lon1 = testcases[i][1], azi1 = testcases[i][2],
164  lat2 = testcases[i][3], lon2 = testcases[i][4], azi2 = testcases[i][5],
165  s12 = testcases[i][6], a12 = testcases[i][7], m12 = testcases[i][8],
166  M12 = testcases[i][9], M21 = testcases[i][10], S12 = testcases[i][11];
167  GeodesicData dir = Geodesic.WGS84.ArcDirect(lat1, lon1, azi1, a12,
168  GeodesicMask.ALL |
170  assertEquals(lat2, dir.lat2, 1e-13);
171  assertEquals(lon2, dir.lon2, 1e-13);
172  assertEquals(azi2, dir.azi2, 1e-13);
173  assertEquals(s12, dir.s12, 1e-8);
174  assertEquals(m12, dir.m12, 1e-8);
175  assertEquals(M12, dir.M12, 1e-15);
176  assertEquals(M21, dir.M21, 1e-15);
177  assertEquals(S12, dir.S12, 0.1);
178  }
179  }
180 
181  @Test
182  public void GeodSolve0() {
183  GeodesicData inv = Geodesic.WGS84.Inverse(40.6, -73.8,
184  49.01666667, 2.55);
185  assertEquals(inv.azi1, 53.47022, 0.5e-5);
186  assertEquals(inv.azi2, 111.59367, 0.5e-5);
187  assertEquals(inv.s12, 5853226, 0.5);
188  }
189 
190  @Test
191  public void GeodSolve1() {
192  GeodesicData dir = Geodesic.WGS84.Direct(40.63972222, -73.77888889,
193  53.5, 5850e3);
194  assertEquals(dir.lat2, 49.01467, 0.5e-5);
195  assertEquals(dir.lon2, 2.56106, 0.5e-5);
196  assertEquals(dir.azi2, 111.62947, 0.5e-5);
197  }
198 
199  @Test
200  public void GeodSolve2() {
201  // Check fix for antipodal prolate bug found 2010-09-04
202  Geodesic geod = new Geodesic(6.4e6, -1/150.0);
203  GeodesicData inv = geod.Inverse(0.07476, 0, -0.07476, 180);
204  assertEquals(inv.azi1, 90.00078, 0.5e-5);
205  assertEquals(inv.azi2, 90.00078, 0.5e-5);
206  assertEquals(inv.s12, 20106193, 0.5);
207  inv = geod.Inverse(0.1, 0, -0.1, 180);
208  assertEquals(inv.azi1, 90.00105, 0.5e-5);
209  assertEquals(inv.azi2, 90.00105, 0.5e-5);
210  assertEquals(inv.s12, 20106193, 0.5);
211  }
212 
213  @Test
214  public void GeodSolve4() {
215  // Check fix for short line bug found 2010-05-21
216  GeodesicData inv = Geodesic.WGS84.Inverse(36.493349428792, 0,
217  36.49334942879201, .0000008);
218  assertEquals(inv.s12, 0.072, 0.5e-3);
219  }
220 
221  @Test
222  public void GeodSolve5() {
223  // Check fix for point2=pole bug found 2010-05-03
224  GeodesicData dir = Geodesic.WGS84.Direct(0.01777745589997, 30, 0, 10e6);
225  assertEquals(dir.lat2, 90, 0.5e-5);
226  if (dir.lon2 < 0) {
227  assertEquals(dir.lon2, -150, 0.5e-5);
228  assertEquals(Math.abs(dir.azi2), 180, 0.5e-5);
229  } else {
230  assertEquals(dir.lon2, 30, 0.5e-5);
231  assertEquals(dir.azi2, 0, 0.5e-5);
232  }
233  }
234 
235  @Test
236  public void GeodSolve6() {
237  // Check fix for volatile sbet12a bug found 2011-06-25 (gcc 4.4.4
238  // x86 -O3). Found again on 2012-03-27 with tdm-mingw32 (g++ 4.6.1).
239  GeodesicData inv =
240  Geodesic.WGS84.Inverse(88.202499451857, 0,
241  -88.202499451857, 179.981022032992859592);
242  assertEquals(inv.s12, 20003898.214, 0.5e-3);
243  inv = Geodesic.WGS84.Inverse(89.262080389218, 0,
244  -89.262080389218, 179.992207982775375662);
245  assertEquals(inv.s12, 20003925.854, 0.5e-3);
246  inv = Geodesic.WGS84.Inverse(89.333123580033, 0, -89.333123580032997687,
247  179.99295812360148422);
248  assertEquals(inv.s12, 20003926.881, 0.5e-3);
249  }
250 
251  @Test
252  public void GeodSolve9() {
253  // Check fix for volatile x bug found 2011-06-25 (gcc 4.4.4 x86 -O3)
254  GeodesicData inv =
255  Geodesic.WGS84.Inverse(56.320923501171, 0,
256  -56.320923501171, 179.664747671772880215);
257  assertEquals(inv.s12, 19993558.287, 0.5e-3);
258  }
259 
260  @Test
261  public void GeodSolve10() {
262  // Check fix for adjust tol1_ bug found 2011-06-25 (Visual Studio
263  // 10 rel + debug)
264  GeodesicData inv =
265  Geodesic.WGS84.Inverse(52.784459512564, 0,
266  -52.784459512563990912, 179.634407464943777557);
267  assertEquals(inv.s12, 19991596.095, 0.5e-3);
268  }
269 
270  @Test
271  public void GeodSolve11() {
272  // Check fix for bet2 = -bet1 bug found 2011-06-25 (Visual Studio
273  // 10 rel + debug)
274  GeodesicData inv =
275  Geodesic.WGS84.Inverse(48.522876735459, 0,
276  -48.52287673545898293, 179.599720456223079643);
277  assertEquals(inv.s12, 19989144.774, 0.5e-3);
278  }
279 
280  @Test
281  public void GeodSolve12() {
282  // Check fix for inverse geodesics on extreme prolate/oblate
283  // ellipsoids Reported 2012-08-29 Stefan Guenther
284  // <stefan.gunther@embl.de>; fixed 2012-10-07
285  Geodesic geod = new Geodesic(89.8, -1.83);
286  GeodesicData inv = geod.Inverse(0, 0, -10, 160);
287  assertEquals(inv.azi1, 120.27, 1e-2);
288  assertEquals(inv.azi2, 105.15, 1e-2);
289  assertEquals(inv.s12, 266.7, 1e-1);
290  }
291 
292  @Test
293  public void GeodSolve14() {
294  // Check fix for inverse ignoring lon12 = nan
295  GeodesicData inv = Geodesic.WGS84.Inverse(0, 0, 1, Double.NaN);
296  assertTrue(isNaN(inv.azi1));
297  assertTrue(isNaN(inv.azi2));
298  assertTrue(isNaN(inv.s12));
299  }
300 
301  @Test
302  public void GeodSolve15() {
303  // Initial implementation of Math::eatanhe was wrong for e^2 < 0. This
304  // checks that this is fixed.
305  Geodesic geod = new Geodesic(6.4e6, -1/150.0);
306  GeodesicData dir = geod.Direct(1, 2, 3, 4, GeodesicMask.AREA);
307  assertEquals(dir.S12, 23700, 0.5);
308  }
309 
310  @Test
311  public void GeodSolve17() {
312  // Check fix for LONG_UNROLL bug found on 2015-05-07
313  GeodesicData dir =
314  Geodesic.WGS84.Direct(40, -75, -10, 2e7,
316  assertEquals(dir.lat2, -39, 1);
317  assertEquals(dir.lon2, -254, 1);
318  assertEquals(dir.azi2, -170, 1);
319  GeodesicLine line = Geodesic.WGS84.Line(40, -75, -10);
321  assertEquals(dir.lat2, -39, 1);
322  assertEquals(dir.lon2, -254, 1);
323  assertEquals(dir.azi2, -170, 1);
324  dir = Geodesic.WGS84.Direct(40, -75, -10, 2e7);
325  assertEquals(dir.lat2, -39, 1);
326  assertEquals(dir.lon2, 105, 1);
327  assertEquals(dir.azi2, -170, 1);
328  dir = line.Position(2e7);
329  assertEquals(dir.lat2, -39, 1);
330  assertEquals(dir.lon2, 105, 1);
331  assertEquals(dir.azi2, -170, 1);
332  }
333 
334  @Test
335  public void GeodSolve26() {
336  // Check 0/0 problem with area calculation on sphere 2015-09-08
337  Geodesic geod = new Geodesic(6.4e6, 0);
338  GeodesicData inv = geod.Inverse(1, 2, 3, 4, GeodesicMask.AREA);
339  assertEquals(inv.S12, 49911046115.0, 0.5);
340  }
341 
342  @Test
343  public void GeodSolve28() {
344  // Check for bad placement of assignment of r.a12 with |f| > 0.01 (bug in
345  // Java implementation fixed on 2015-05-19).
346  Geodesic geod = new Geodesic(6.4e6, 0.1);
347  GeodesicData dir = geod.Direct(1, 2, 10, 5e6);
348  assertEquals(dir.a12, 48.55570690, 0.5e-8);
349  }
350 
351  @Test
352  public void GeodSolve29() {
353  // Check longitude unrolling with inverse calculation 2015-09-16
354  GeodesicData dir = Geodesic.WGS84.Inverse(0, 539, 0, 181);
355  assertEquals(dir.lon1, 179, 1e-10);
356  assertEquals(dir.lon2, -179, 1e-10);
357  assertEquals(dir.s12, 222639, 0.5);
358  dir = Geodesic.WGS84.Inverse(0, 539, 0, 181,
361  assertEquals(dir.lon1, 539, 1e-10);
362  assertEquals(dir.lon2, 541, 1e-10);
363  assertEquals(dir.s12, 222639, 0.5);
364  }
365 
366  @Test
367  public void GeodSolve33() {
368  // Check max(-0.0,+0.0) issues 2015-08-22 (triggered by bugs in Octave --
369  // sind(-0.0) = +0.0 -- and in some version of Visual Studio --
370  // fmod(-0.0, 360.0) = +0.0.
371  GeodesicData inv = Geodesic.WGS84.Inverse(0, 0, 0, 179);
372  assertEquals(inv.azi1, 90.00000, 0.5e-5);
373  assertEquals(inv.azi2, 90.00000, 0.5e-5);
374  assertEquals(inv.s12, 19926189, 0.5);
375  inv = Geodesic.WGS84.Inverse(0, 0, 0, 179.5);
376  assertEquals(inv.azi1, 55.96650, 0.5e-5);
377  assertEquals(inv.azi2, 124.03350, 0.5e-5);
378  assertEquals(inv.s12, 19980862, 0.5);
379  inv = Geodesic.WGS84.Inverse(0, 0, 0, 180);
380  assertEquals(inv.azi1, 0.00000, 0.5e-5);
381  assertEquals(Math.abs(inv.azi2), 180.00000, 0.5e-5);
382  assertEquals(inv.s12, 20003931, 0.5);
383  inv = Geodesic.WGS84.Inverse(0, 0, 1, 180);
384  assertEquals(inv.azi1, 0.00000, 0.5e-5);
385  assertEquals(Math.abs(inv.azi2), 180.00000, 0.5e-5);
386  assertEquals(inv.s12, 19893357, 0.5);
387  Geodesic geod = new Geodesic(6.4e6, 0);
388  inv = geod.Inverse(0, 0, 0, 179);
389  assertEquals(inv.azi1, 90.00000, 0.5e-5);
390  assertEquals(inv.azi2, 90.00000, 0.5e-5);
391  assertEquals(inv.s12, 19994492, 0.5);
392  inv = geod.Inverse(0, 0, 0, 180);
393  assertEquals(inv.azi1, 0.00000, 0.5e-5);
394  assertEquals(Math.abs(inv.azi2), 180.00000, 0.5e-5);
395  assertEquals(inv.s12, 20106193, 0.5);
396  inv = geod.Inverse(0, 0, 1, 180);
397  assertEquals(inv.azi1, 0.00000, 0.5e-5);
398  assertEquals(Math.abs(inv.azi2), 180.00000, 0.5e-5);
399  assertEquals(inv.s12, 19994492, 0.5);
400  geod = new Geodesic(6.4e6, -1/300.0);
401  inv = geod.Inverse(0, 0, 0, 179);
402  assertEquals(inv.azi1, 90.00000, 0.5e-5);
403  assertEquals(inv.azi2, 90.00000, 0.5e-5);
404  assertEquals(inv.s12, 19994492, 0.5);
405  inv = geod.Inverse(0, 0, 0, 180);
406  assertEquals(inv.azi1, 90.00000, 0.5e-5);
407  assertEquals(inv.azi2, 90.00000, 0.5e-5);
408  assertEquals(inv.s12, 20106193, 0.5);
409  inv = geod.Inverse(0, 0, 0.5, 180);
410  assertEquals(inv.azi1, 33.02493, 0.5e-5);
411  assertEquals(inv.azi2, 146.97364, 0.5e-5);
412  assertEquals(inv.s12, 20082617, 0.5);
413  inv = geod.Inverse(0, 0, 1, 180);
414  assertEquals(inv.azi1, 0.00000, 0.5e-5);
415  assertEquals(Math.abs(inv.azi2), 180.00000, 0.5e-5);
416  assertEquals(inv.s12, 20027270, 0.5);
417  }
418 
419  @Test
420  public void GeodSolve55() {
421  // Check fix for nan + point on equator or pole not returning all nans in
422  // Geodesic::Inverse, found 2015-09-23.
423  GeodesicData inv = Geodesic.WGS84.Inverse(Double.NaN, 0, 0, 90);
424  assertTrue(isNaN(inv.azi1));
425  assertTrue(isNaN(inv.azi2));
426  assertTrue(isNaN(inv.s12));
427  inv = Geodesic.WGS84.Inverse(Double.NaN, 0, 90, 3);
428  assertTrue(isNaN(inv.azi1));
429  assertTrue(isNaN(inv.azi2));
430  assertTrue(isNaN(inv.s12));
431  }
432 
433  @Test
434  public void GeodSolve59() {
435  // Check for points close with longitudes close to 180 deg apart.
436  GeodesicData inv = Geodesic.WGS84.Inverse(5, 0.00000000000001, 10, 180);
437  assertEquals(inv.azi1, 0.000000000000035, 1.5e-14);
438  assertEquals(inv.azi2, 179.99999999999996, 1.5e-14);
439  assertEquals(inv.s12, 18345191.174332713, 4e-9);
440  }
441 
442  @Test
443  public void GeodSolve61() {
444  // Make sure small negative azimuths are west-going
445  GeodesicData dir =
446  Geodesic.WGS84.Direct(45, 0, -0.000000000000000003, 1e7,
448  assertEquals(dir.lat2, 45.30632, 0.5e-5);
449  assertEquals(dir.lon2, -180, 0.5e-5);
450  assertEquals(Math.abs(dir.azi2), 180, 0.5e-5);
451  GeodesicLine line = Geodesic.WGS84.InverseLine(45, 0, 80,
452  -0.000000000000000003);
453  dir = line.Position(1e7, GeodesicMask.STANDARD | GeodesicMask.LONG_UNROLL);
454  assertEquals(dir.lat2, 45.30632, 0.5e-5);
455  assertEquals(dir.lon2, -180, 0.5e-5);
456  assertEquals(Math.abs(dir.azi2), 180, 0.5e-5);
457  }
458 
459  @Test
460  public void GeodSolve65() {
461  // Check for bug in east-going check in GeodesicLine (needed to check for
462  // sign of 0) and sign error in area calculation due to a bogus override
463  // of the code for alp12. Found/fixed on 2015-12-19.
464  GeodesicLine line = Geodesic.WGS84.InverseLine(30, -0.000000000000000001,
465  -31, 180);
466  GeodesicData dir =
468  assertEquals(dir.lat1, 30.00000 , 0.5e-5);
469  assertEquals(dir.lon1, -0.00000 , 0.5e-5);
470  assertEquals(Math.abs(dir.azi1), 180.00000, 0.5e-5);
471  assertEquals(dir.lat2, -60.23169 , 0.5e-5);
472  assertEquals(dir.lon2, -0.00000 , 0.5e-5);
473  assertEquals(Math.abs(dir.azi2), 180.00000, 0.5e-5);
474  assertEquals(dir.s12 , 10000000 , 0.5);
475  assertEquals(dir.a12 , 90.06544 , 0.5e-5);
476  assertEquals(dir.m12 , 6363636 , 0.5);
477  assertEquals(dir.M12 , -0.0012834, 0.5e7);
478  assertEquals(dir.M21 , 0.0013749 , 0.5e-7);
479  assertEquals(dir.S12 , 0 , 0.5);
480  dir = line.Position(2e7, GeodesicMask.ALL | GeodesicMask.LONG_UNROLL);
481  assertEquals(dir.lat1, 30.00000 , 0.5e-5);
482  assertEquals(dir.lon1, -0.00000 , 0.5e-5);
483  assertEquals(Math.abs(dir.azi1), 180.00000, 0.5e-5);
484  assertEquals(dir.lat2, -30.03547 , 0.5e-5);
485  assertEquals(dir.lon2, -180.00000, 0.5e-5);
486  assertEquals(dir.azi2, -0.00000 , 0.5e-5);
487  assertEquals(dir.s12 , 20000000 , 0.5);
488  assertEquals(dir.a12 , 179.96459 , 0.5e-5);
489  assertEquals(dir.m12 , 54342 , 0.5);
490  assertEquals(dir.M12 , -1.0045592, 0.5e7);
491  assertEquals(dir.M21 , -0.9954339, 0.5e-7);
492  assertEquals(dir.S12 , 127516405431022.0, 0.5);
493  }
494 
495  @Test
496  public void GeodSolve69() {
497  // Check for InverseLine if line is slightly west of S and that s13 is
498  // correctly set.
499  GeodesicLine line =
500  Geodesic.WGS84.InverseLine(-5, -0.000000000000002, -10, 180);
501  GeodesicData dir =
503  assertEquals(dir.lat2, 4.96445 , 0.5e-5);
504  assertEquals(dir.lon2, -180.00000, 0.5e-5);
505  assertEquals(dir.azi2, -0.00000 , 0.5e-5);
506  dir = line.Position(0.5 * line.Distance(),
508  assertEquals(dir.lat2, -87.52461 , 0.5e-5);
509  assertEquals(dir.lon2, -0.00000 , 0.5e-5);
510  assertEquals(dir.azi2, -180.00000, 0.5e-5);
511  }
512 
513  @Test
514  public void GeodSolve71() {
515  // Check that DirectLine sets s13.
516  GeodesicLine line = Geodesic.WGS84.DirectLine(1, 2, 45, 1e7);
517  GeodesicData dir =
518  line.Position(0.5 * line.Distance(),
520  assertEquals(dir.lat2, 30.92625, 0.5e-5);
521  assertEquals(dir.lon2, 37.54640, 0.5e-5);
522  assertEquals(dir.azi2, 55.43104, 0.5e-5);
523  }
524 
525  @Test
526  public void GeodSolve73() {
527  // Check for backwards from the pole bug reported by Anon on 2016-02-13.
528  // This only affected the Java implementation. It was introduced in Java
529  // version 1.44 and fixed in 1.46-SNAPSHOT on 2016-01-17.
530  GeodesicData dir = Geodesic.WGS84.Direct(90, 10, 180, -1e6);
531  assertEquals(dir.lat2, 81.04623, 0.5e-5);
532  assertEquals(dir.lon2, -170, 0.5e-5);
533  assertEquals(dir.azi2, 0, 0.5e-5);
534  }
535 
536  @Test
537  public void GeodSolve74() {
538  // Check fix for inaccurate areas, bug introduced in v1.46, fixed
539  // 2015-10-16.
540  GeodesicData inv = Geodesic.WGS84.Inverse(54.1589, 15.3872,
541  54.1591, 15.3877,
542  GeodesicMask.ALL);
543  assertEquals(inv.azi1, 55.723110355, 5e-9);
544  assertEquals(inv.azi2, 55.723515675, 5e-9);
545  assertEquals(inv.s12, 39.527686385, 5e-9);
546  assertEquals(inv.a12, 0.000355495, 5e-9);
547  assertEquals(inv.m12, 39.527686385, 5e-9);
548  assertEquals(inv.M12, 0.999999995, 5e-9);
549  assertEquals(inv.M21, 0.999999995, 5e-9);
550  assertEquals(inv.S12, 286698586.30197, 5e-4);
551  }
552 
553  @Test
554  public void GeodSolve76() {
555  // The distance from Wellington and Salamanca (a classic failure of
556  // Vincenty)
557  GeodesicData inv = Geodesic.WGS84.Inverse(-(41+19/60.0), 174+49/60.0,
558  40+58/60.0, -(5+30/60.0));
559  assertEquals(inv.azi1, 160.39137649664, 0.5e-11);
560  assertEquals(inv.azi2, 19.50042925176, 0.5e-11);
561  assertEquals(inv.s12, 19960543.857179, 0.5e-6);
562  }
563 
564  @Test
565  public void GeodSolve78() {
566  // An example where the NGS calculator fails to converge */
567  GeodesicData inv = Geodesic.WGS84.Inverse(27.2, 0.0, -27.1, 179.5);
568  assertEquals(inv.azi1, 45.82468716758, 0.5e-11);
569  assertEquals(inv.azi2, 134.22776532670, 0.5e-11);
570  assertEquals(inv.s12, 19974354.765767, 0.5e-6);
571  }
572 
573  @Test
574  public void Planimeter0() {
575  // Check fix for pole-encircling bug found 2011-03-16
576  double pa[][] = {{89, 0}, {89, 90}, {89, 180}, {89, 270}};
577  PolygonResult a = Planimeter(pa);
578  assertEquals(a.perimeter, 631819.8745, 1e-4);
579  assertEquals(a.area, 24952305678.0, 1);
580 
581  double pb[][] = {{-89, 0}, {-89, 90}, {-89, 180}, {-89, 270}};
582  a = Planimeter(pb);
583  assertEquals(a.perimeter, 631819.8745, 1e-4);
584  assertEquals(a.area, -24952305678.0, 1);
585 
586  double pc[][] = {{0, -1}, {-1, 0}, {0, 1}, {1, 0}};
587  a = Planimeter(pc);
588  assertEquals(a.perimeter, 627598.2731, 1e-4);
589  assertEquals(a.area, 24619419146.0, 1);
590 
591  double pd[][] = {{90, 0}, {0, 0}, {0, 90}};
592  a = Planimeter(pd);
593  assertEquals(a.perimeter, 30022685, 1);
594  assertEquals(a.area, 63758202715511.0, 1);
595  a = PolyLength(pd);
596  assertEquals(a.perimeter, 20020719, 1);
597  assertTrue(isNaN(a.area));
598  }
599 
600  @Test
601  public void Planimeter5() {
602  // Check fix for Planimeter pole crossing bug found 2011-06-24
603  double points[][] = {{89, 0.1}, {89, 90.1}, {89, -179.9}};
604  PolygonResult a = Planimeter(points);
605  assertEquals(a.perimeter, 539297, 1);
606  assertEquals(a.area, 12476152838.5, 1);
607  }
608 
609  @Test
610  public void Planimeter6() {
611  // Check fix for Planimeter lon12 rounding bug found 2012-12-03
612  double pa[][] = {{9, -0.00000000000001}, {9, 180}, {9, 0}};
613  PolygonResult a = Planimeter(pa);
614  assertEquals(a.perimeter, 36026861, 1);
615  assertEquals(a.area, 0, 1);
616  double pb[][] = {{9, 0.00000000000001}, {9, 0}, {9, 180}};
617  a = Planimeter(pb);
618  assertEquals(a.perimeter, 36026861, 1);
619  assertEquals(a.area, 0, 1);
620  double pc[][] = {{9, 0.00000000000001}, {9, 180}, {9, 0}};
621  a = Planimeter(pc);
622  assertEquals(a.perimeter, 36026861, 1);
623  assertEquals(a.area, 0, 1);
624  double pd[][] = {{9, -0.00000000000001}, {9, 0}, {9, 180}};
625  a = Planimeter(pd);
626  assertEquals(a.perimeter, 36026861, 1);
627  assertEquals(a.area, 0, 1);
628  }
629 
630  @Test
631  public void Planimeter12() {
632  // Area of arctic circle (not really -- adjunct to rhumb-area test)
633  double points[][] = {{66.562222222, 0}, {66.562222222, 180}};
634  PolygonResult a = Planimeter(points);
635  assertEquals(a.perimeter, 10465729, 1);
636  assertEquals(a.area, 0, 1);
637  }
638 
639  @Test
640  public void Planimeter13() {
641  // Check encircling pole twice
642  double points[][] = {{89,-360}, {89,-240}, {89,-120},
643  {89,0}, {89,120}, {89,240}};
644  PolygonResult a = Planimeter(points);
645  assertEquals(a.perimeter, 1160741, 1);
646  assertEquals(a.area, 32415230256.0, 1);
647  }
648 
649 }
GeodesicData Position(double s12)
static final Geodesic WGS84
Definition: Geodesic.java:1207
static PolygonResult Planimeter(double points[][])
static final PolygonArea polygon
int RealScalar int RealScalar int RealScalar * pc
static PolygonResult PolyLength(double points[][])
GeodesicLine InverseLine(double lat1, double lon1, double lat2, double lon2)
Definition: Geodesic.java:1077
void AddPoint(double lat, double lon)
GeodesicData Direct(double lat1, double lon1, double azi1, double s12)
Definition: Geodesic.java:313
static final PolygonArea polyline
GeodesicData Inverse(double lat1, double lon1, double lat2, double lon2)
Definition: Geodesic.java:615
GeodesicLine DirectLine(double lat1, double lon1, double azi1, double s12)
Definition: Geodesic.java:477
Array< double, 1, 3 > e(1./3., 0.5, 2.)
GeodesicLine Line(double lat1, double lon1, double azi1)
Definition: Geodesic.java:1130
GeodesicData ArcDirect(double lat1, double lon1, double azi1, double a12)
Definition: Geodesic.java:366
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


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