Fortran library for Geodesics  2.0
geodtest.for
Go to the documentation of this file.
1 *> @file geodtest.for
2 *! @brief Test suite for the geodesic routines in Fortran
3 *!
4 *! Run these tests by configuring with cmake and running "make test".
5 *!
6 *! Copyright (c) Charles Karney (2015-2021) <charles@karney.com> and
7 *! licensed under the MIT/X11 License. For more information, see
8 *! https://geographiclib.sourceforge.io/
9 
10 *> @cond SKIP
11 
12  block data tests
13 
14  integer j
15  double precision tstdat(20, 12)
16  common /tstcom/ tstdat
17  data (tstdat(1,j), j = 1,12) /
18  + 35.60777d0,-139.44815d0,111.098748429560326d0,
19  + -11.17491d0,-69.95921d0,129.289270889708762d0,
20  + 8935244.5604818305d0,80.50729714281974d0,6273170.2055303837d0,
21  + 0.16606318447386067d0,0.16479116945612937d0,
22  + 12841384694976.432d0 /
23  data (tstdat(2,j), j = 1,12) /
24  + 55.52454d0,106.05087d0,22.020059880982801d0,
25  + 77.03196d0,197.18234d0,109.112041110671519d0,
26  + 4105086.1713924406d0,36.892740690445894d0,
27  + 3828869.3344387607d0,
28  + 0.80076349608092607d0,0.80101006984201008d0,
29  + 61674961290615.615d0 /
30  data (tstdat(3,j), j = 1,12) /
31  + -21.97856d0,142.59065d0,-32.44456876433189d0,
32  + 41.84138d0,98.56635d0,-41.84359951440466d0,
33  + 8394328.894657671d0,75.62930491011522d0,6161154.5773110616d0,
34  + 0.24816339233950381d0,0.24930251203627892d0,
35  + -6637997720646.717d0 /
36  data (tstdat(4,j), j = 1,12) /
37  + -66.99028d0,112.2363d0,173.73491240878403d0,
38  + -12.70631d0,285.90344d0,2.512956620913668d0,
39  + 11150344.2312080241d0,100.278634181155759d0,
40  + 6289939.5670446687d0,
41  + -0.17199490274700385d0,-0.17722569526345708d0,
42  + -121287239862139.744d0 /
43  data (tstdat(5,j), j = 1,12) /
44  + -17.42761d0,173.34268d0,-159.033557661192928d0,
45  + -15.84784d0,5.93557d0,-20.787484651536988d0,
46  + 16076603.1631180673d0,144.640108810286253d0,
47  + 3732902.1583877189d0,
48  + -0.81273638700070476d0,-0.81299800519154474d0,
49  + 97825992354058.708d0 /
50  data (tstdat(6,j), j = 1,12) /
51  + 32.84994d0,48.28919d0,150.492927788121982d0,
52  + -56.28556d0,202.29132d0,48.113449399816759d0,
53  + 16727068.9438164461d0,150.565799985466607d0,
54  + 3147838.1910180939d0,
55  + -0.87334918086923126d0,-0.86505036767110637d0,
56  + -72445258525585.010d0 /
57  data (tstdat(7,j), j = 1,12) /
58  + 6.96833d0,52.74123d0,92.581585386317712d0,
59  + -7.39675d0,206.17291d0,90.721692165923907d0,
60  + 17102477.2496958388d0,154.147366239113561d0,
61  + 2772035.6169917581d0,
62  + -0.89991282520302447d0,-0.89986892177110739d0,
63  + -1311796973197.995d0 /
64  data (tstdat(8,j), j = 1,12) /
65  + -50.56724d0,-16.30485d0,-105.439679907590164d0,
66  + -33.56571d0,-94.97412d0,-47.348547835650331d0,
67  + 6455670.5118668696d0,58.083719495371259d0,
68  + 5409150.7979815838d0,
69  + 0.53053508035997263d0,0.52988722644436602d0,
70  + 41071447902810.047d0 /
71  data (tstdat(9,j), j = 1,12) /
72  + -58.93002d0,-8.90775d0,140.965397902500679d0,
73  + -8.91104d0,133.13503d0,19.255429433416599d0,
74  + 11756066.0219864627d0,105.755691241406877d0,
75  + 6151101.2270708536d0,
76  + -0.26548622269867183d0,-0.27068483874510741d0,
77  + -86143460552774.735d0 /
78  data (tstdat(10,j), j = 1,12) /
79  + -68.82867d0,-74.28391d0,93.774347763114881d0,
80  + -50.63005d0,-8.36685d0,34.65564085411343d0,
81  + 3956936.926063544d0,35.572254987389284d0,3708890.9544062657d0,
82  + 0.81443963736383502d0,0.81420859815358342d0,
83  + -41845309450093.787d0 /
84  data (tstdat(11,j), j = 1,12) /
85  + -10.62672d0,-32.0898d0,-86.426713286747751d0,
86  + 5.883d0,-134.31681d0,-80.473780971034875d0,
87  + 11470869.3864563009d0,103.387395634504061d0,
88  + 6184411.6622659713d0,
89  + -0.23138683500430237d0,-0.23155097622286792d0,
90  + 4198803992123.548d0 /
91  data (tstdat(12,j), j = 1,12) /
92  + -21.76221d0,166.90563d0,29.319421206936428d0,
93  + 48.72884d0,213.97627d0,43.508671946410168d0,
94  + 9098627.3986554915d0,81.963476716121964d0,
95  + 6299240.9166992283d0,
96  + 0.13965943368590333d0,0.14152969707656796d0,
97  + 10024709850277.476d0 /
98  data (tstdat(13,j), j = 1,12) /
99  + -19.79938d0,-174.47484d0,71.167275780171533d0,
100  + -11.99349d0,-154.35109d0,65.589099775199228d0,
101  + 2319004.8601169389d0,20.896611684802389d0,
102  + 2267960.8703918325d0,
103  + 0.93427001867125849d0,0.93424887135032789d0,
104  + -3935477535005.785d0 /
105  data (tstdat(14,j), j = 1,12) /
106  + -11.95887d0,-116.94513d0,92.712619830452549d0,
107  + 4.57352d0,7.16501d0,78.64960934409585d0,
108  + 13834722.5801401374d0,124.688684161089762d0,
109  + 5228093.177931598d0,
110  + -0.56879356755666463d0,-0.56918731952397221d0,
111  + -9919582785894.853d0 /
112  data (tstdat(15,j), j = 1,12) /
113  + -87.85331d0,85.66836d0,-65.120313040242748d0,
114  + 66.48646d0,16.09921d0,-4.888658719272296d0,
115  + 17286615.3147144645d0,155.58592449699137d0,
116  + 2635887.4729110181d0,
117  + -0.90697975771398578d0,-0.91095608883042767d0,
118  + 42667211366919.534d0 /
119  data (tstdat(16,j), j = 1,12) /
120  + 1.74708d0,128.32011d0,-101.584843631173858d0,
121  + -11.16617d0,11.87109d0,-86.325793296437476d0,
122  + 12942901.1241347408d0,116.650512484301857d0,
123  + 5682744.8413270572d0,
124  + -0.44857868222697644d0,-0.44824490340007729d0,
125  + 10763055294345.653d0 /
126  data (tstdat(17,j), j = 1,12) /
127  + -25.72959d0,-144.90758d0,-153.647468693117198d0,
128  + -57.70581d0,-269.17879d0,-48.343983158876487d0,
129  + 9413446.7452453107d0,84.664533838404295d0,
130  + 6356176.6898881281d0,
131  + 0.09492245755254703d0,0.09737058264766572d0,
132  + 74515122850712.444d0 /
133  data (tstdat(18,j), j = 1,12) /
134  + -41.22777d0,122.32875d0,14.285113402275739d0,
135  + -7.57291d0,130.37946d0,10.805303085187369d0,
136  + 3812686.035106021d0,34.34330804743883d0,3588703.8812128856d0,
137  + 0.82605222593217889d0,0.82572158200920196d0,
138  + -2456961531057.857d0 /
139  data (tstdat(19,j), j = 1,12) /
140  + 11.01307d0,138.25278d0,79.43682622782374d0,
141  + 6.62726d0,247.05981d0,103.708090215522657d0,
142  + 11911190.819018408d0,107.341669954114577d0,
143  + 6070904.722786735d0,
144  + -0.29767608923657404d0,-0.29785143390252321d0,
145  + 17121631423099.696d0 /
146  data (tstdat(20,j), j = 1,12) /
147  + -29.47124d0,95.14681d0,-163.779130441688382d0,
148  + -27.46601d0,-69.15955d0,-15.909335945554969d0,
149  + 13487015.8381145492d0,121.294026715742277d0,
150  + 5481428.9945736388d0,
151  + -0.51527225545373252d0,-0.51556587964721788d0,
152  + 104679964020340.318d0 /
153  end
154 
155  integer function assert(x, y, d)
156  double precision x, y, d
157 
158  if (abs(x - y) .le. d) then
159  assert = 0
160  else
161  assert = 1
162  print 10, x, y, d
163  10 format(1x, 'assert fails: ',
164  + g14.7, ' != ', g14.7, ' +/- ', g10.3)
165  end if
166 
167  return
168  end
169 
170  integer function chknan(x)
171  double precision x
172 
173  if (x .ne. x) then
174  chknan = 0
175  else
176  chknan = 1
177  end if
178 
179  return
180  end
181 
182  integer function tstinv()
183  double precision tstdat(20, 12)
184  common /tstcom/ tstdat
185  double precision lat1, lon1, azi1, lat2, lon2, azi2,
186  + s12, a12, m12, MM12, MM21, SS12
187  double precision azi1a, azi2a, s12a, a12a,
188  + m12a, MM12a, MM21a, SS12a
189  double precision a, f
190  integer r, assert, i, omask
191  include 'geodesic.inc'
192 
193 * WGS84 values
194  a = 6378137d0
195  f = 1/298.257223563d0
196  omask = 1 + 2 + 4 + 8
197  r = 0
198 
199  do 10 i = 1,20
200  lat1 = tstdat(i, 1)
201  lon1 = tstdat(i, 2)
202  azi1 = tstdat(i, 3)
203  lat2 = tstdat(i, 4)
204  lon2 = tstdat(i, 5)
205  azi2 = tstdat(i, 6)
206  s12 = tstdat(i, 7)
207  a12 = tstdat(i, 8)
208  m12 = tstdat(i, 9)
209  mm12 = tstdat(i, 10)
210  mm21 = tstdat(i, 11)
211  ss12 = tstdat(i, 12)
212  call invers(a, f, lat1, lon1, lat2, lon2,
213  + s12a, azi1a, azi2a, omask, a12a, m12a, mm12a, mm21a, ss12a)
214  r = r + assert(azi1, azi1a, 1d-13)
215  r = r + assert(azi2, azi2a, 1d-13)
216  r = r + assert(s12, s12a, 1d-8)
217  r = r + assert(a12, a12a, 1d-13)
218  r = r + assert(m12, m12a, 1d-8)
219  r = r + assert(mm12, mm12a, 1d-15)
220  r = r + assert(mm21, mm21a, 1d-15)
221  r = r + assert(ss12, ss12a, 0.1d0)
222  10 continue
223 
224  tstinv = r
225  return
226  end
227 
228  integer function tstdir()
229  double precision tstdat(20, 12)
230  common /tstcom/ tstdat
231  double precision lat1, lon1, azi1, lat2, lon2, azi2,
232  + s12, a12, m12, MM12, MM21, SS12
233  double precision lat2a, lon2a, azi2a, a12a,
234  + m12a, MM12a, MM21a, SS12a
235  double precision a, f
236  integer r, assert, i, omask, flags
237  include 'geodesic.inc'
238 
239 * WGS84 values
240  a = 6378137d0
241  f = 1/298.257223563d0
242  omask = 1 + 2 + 4 + 8
243  flags = 2
244  r = 0
245 
246  do 10 i = 1,20
247  lat1 = tstdat(i, 1)
248  lon1 = tstdat(i, 2)
249  azi1 = tstdat(i, 3)
250  lat2 = tstdat(i, 4)
251  lon2 = tstdat(i, 5)
252  azi2 = tstdat(i, 6)
253  s12 = tstdat(i, 7)
254  a12 = tstdat(i, 8)
255  m12 = tstdat(i, 9)
256  mm12 = tstdat(i, 10)
257  mm21 = tstdat(i, 11)
258  ss12 = tstdat(i, 12)
259  call direct(a, f, lat1, lon1, azi1, s12, flags,
260  + lat2a, lon2a, azi2a, omask, a12a, m12a, mm12a, mm21a, ss12a)
261  r = r + assert(lat2, lat2a, 1d-13)
262  r = r + assert(lon2, lon2a, 1d-13)
263  r = r + assert(azi2, azi2a, 1d-13)
264  r = r + assert(a12, a12a, 1d-13)
265  r = r + assert(m12, m12a, 1d-8)
266  r = r + assert(mm12, mm12a, 1d-15)
267  r = r + assert(mm21, mm21a, 1d-15)
268  r = r + assert(ss12, ss12a, 0.1d0)
269  10 continue
270 
271  tstdir = r
272  return
273  end
274 
275  integer function tstarc()
276  double precision tstdat(20, 12)
277  common /tstcom/ tstdat
278  double precision lat1, lon1, azi1, lat2, lon2, azi2,
279  + s12, a12, m12, MM12, MM21, SS12
280  double precision lat2a, lon2a, azi2a, s12a,
281  + m12a, MM12a, MM21a, SS12a
282  double precision a, f
283  integer r, assert, i, omask, flags
284  include 'geodesic.inc'
285 
286 * WGS84 values
287  a = 6378137d0
288  f = 1/298.257223563d0
289  omask = 1 + 2 + 4 + 8
290  flags = 1 + 2
291  r = 0
292 
293  do 10 i = 1,20
294  lat1 = tstdat(i, 1)
295  lon1 = tstdat(i, 2)
296  azi1 = tstdat(i, 3)
297  lat2 = tstdat(i, 4)
298  lon2 = tstdat(i, 5)
299  azi2 = tstdat(i, 6)
300  s12 = tstdat(i, 7)
301  a12 = tstdat(i, 8)
302  m12 = tstdat(i, 9)
303  mm12 = tstdat(i, 10)
304  mm21 = tstdat(i, 11)
305  ss12 = tstdat(i, 12)
306  call direct(a, f, lat1, lon1, azi1, a12, flags,
307  + lat2a, lon2a, azi2a, omask, s12a, m12a, mm12a, mm21a, ss12a)
308  r = r + assert(lat2, lat2a, 1d-13)
309  r = r + assert(lon2, lon2a, 1d-13)
310  r = r + assert(azi2, azi2a, 1d-13)
311  r = r + assert(s12, s12a, 1d-8)
312  r = r + assert(m12, m12a, 1d-8)
313  r = r + assert(mm12, mm12a, 1d-15)
314  r = r + assert(mm21, mm21a, 1d-15)
315  r = r + assert(ss12, ss12a, 0.1d0)
316  10 continue
317 
318  tstarc = r
319  return
320  end
321 
322  integer function tstg0()
323  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
324  double precision a, f
325  integer r, assert, omask
326  include 'geodesic.inc'
327 
328 * WGS84 values
329  a = 6378137d0
330  f = 1/298.257223563d0
331  omask = 0
332  r = 0
333  call invers(a, f, 40.6d0, -73.8d0, 49.01666667d0, 2.55d0,
334  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
335  r = r + assert(azi1, 53.47022d0, 0.5d-5)
336  r = r + assert(azi2, 111.59367d0, 0.5d-5)
337  r = r + assert(s12, 5853226d0, 0.5d0)
338 
339  tstg0 = r
340  return
341  end
342 
343  integer function tstg1()
344  double precision lat2, lon2, azi2, a12, m12, MM12, MM21, SS12
345  double precision a, f
346  integer r, assert, omask, flags
347  include 'geodesic.inc'
348 
349 * WGS84 values
350  a = 6378137d0
351  f = 1/298.257223563d0
352  omask = 0
353  flags = 0
354  r = 0
355  call direct(a, f, 40.63972222d0, -73.77888889d0, 53.5d0, 5850d3,
356  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
357  r = r + assert(lat2, 49.01467d0, 0.5d-5)
358  r = r + assert(lon2, 2.56106d0, 0.5d-5)
359  r = r + assert(azi2, 111.62947d0, 0.5d-5)
360 
361  tstg1 = r
362  return
363  end
364 
365  integer function tstg2()
366 * Check fix for antipodal prolate bug found 2010-09-04
367  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
368  double precision a, f
369  integer r, assert, omask
370  include 'geodesic.inc'
371 
372  a = 6.4d6
373  f = -1/150d0
374  omask = 0
375  r = 0
376  call invers(a, f, 0.07476d0, 0d0, -0.07476d0, 180d0,
377  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
378  r = r + assert(azi1, 90.00078d0, 0.5d-5)
379  r = r + assert(azi2, 90.00078d0, 0.5d-5)
380  r = r + assert(s12, 20106193d0, 0.5d0)
381  call invers(a, f, 0.1d0, 0d0, -0.1d0, 180d0,
382  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
383  r = r + assert(azi1, 90.00105d0, 0.5d-5)
384  r = r + assert(azi2, 90.00105d0, 0.5d-5)
385  r = r + assert(s12, 20106193d0, 0.5d0)
386 
387  tstg2 = r
388  return
389  end
390 
391  integer function tstg4()
392 * Check fix for short line bug found 2010-05-21
393  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
394  double precision a, f
395  integer r, assert, omask
396  include 'geodesic.inc'
397 
398 * WGS84 values
399  a = 6378137d0
400  f = 1/298.257223563d0
401  omask = 0
402  r = 0
403  call invers(a, f,
404  + 36.493349428792d0, 0d0, 36.49334942879201d0, .0000008d0,
405  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
406  r = r + assert(s12, 0.072d0, 0.5d-3)
407 
408  tstg4 = r
409  return
410  end
411 
412  integer function tstg5()
413  double precision lat2, lon2, azi2, a12, m12, MM12, MM21, SS12
414  double precision a, f
415  integer r, assert, omask, flags
416  include 'geodesic.inc'
417 
418 * WGS84 values
419  a = 6378137d0
420  f = 1/298.257223563d0
421  omask = 0
422  flags = 0
423  r = 0
424  call direct(a, f, 0.01777745589997d0, 30d0, 0d0, 10d6,
425  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
426  if (lon2 .lt. 0) then
427  r = r + assert(lon2, -150d0, 0.5d-5)
428  r = r + assert(abs(azi2), 180d0, 0.5d-5)
429  else
430  r = r + assert(lon2, 30d0, 0.5d-5)
431  r = r + assert(azi2, 0d0, 0.5d-5)
432  end if
433 
434  tstg5 = r
435  return
436  end
437 
438  integer function tstg6()
439 * Check fix for volatile sbet12a bug found 2011-06-25 (gcc 4.4d0.4d0
440 * x86 -O3). Found again on 2012-03-27 with tdm-mingw32 (g++ 4.6d0.1d0).
441  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
442  double precision a, f
443  integer r, assert, omask
444  include 'geodesic.inc'
445 
446 * WGS84 values
447  a = 6378137d0
448  f = 1/298.257223563d0
449  omask = 0
450  r = 0
451  call invers(a, f, 88.202499451857d0, 0d0,
452  + -88.202499451857d0, 179.981022032992859592d0,
453  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
454  r = r + assert(s12, 20003898.214d0, 0.5d-3)
455  call invers(a, f, 89.262080389218d0, 0d0,
456  + -89.262080389218d0, 179.992207982775375662d0,
457  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
458  r = r + assert(s12, 20003925.854d0, 0.5d-3)
459  call invers(a, f, 89.333123580033d0, 0d0,
460  + -89.333123580032997687d0, 179.99295812360148422d0,
461  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
462  r = r + assert(s12, 20003926.881d0, 0.5d-3)
463 
464  tstg6 = r
465  return
466  end
467 
468  integer function tstg9()
469 * Check fix for volatile x bug found 2011-06-25 (gcc 4.4d0.4d0 x86 -O3)
470  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
471  double precision a, f
472  integer r, assert, omask
473  include 'geodesic.inc'
474 
475 * WGS84 values
476  a = 6378137d0
477  f = 1/298.257223563d0
478  omask = 0
479  r = 0
480  call invers(a, f, 56.320923501171d0, 0d0,
481  + -56.320923501171d0, 179.664747671772880215d0,
482  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
483  r = r + assert(s12, 19993558.287d0, 0.5d-3)
484 
485  tstg9 = r
486  return
487  end
488 
489  integer function tstg10()
490 * Check fix for adjust tol1_ bug found 2011-06-25 (Visual Studio 10 rel
491 * + debug)
492  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
493  double precision a, f
494  integer r, assert, omask
495  include 'geodesic.inc'
496 
497 * WGS84 values
498  a = 6378137d0
499  f = 1/298.257223563d0
500  omask = 0
501  r = 0
502  call invers(a, f, 52.784459512564d0, 0d0,
503  + -52.784459512563990912d0, 179.634407464943777557d0,
504  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
505  r = r + assert(s12, 19991596.095d0, 0.5d-3)
506 
507  tstg10 = r
508  return
509  end
510 
511  integer function tstg11()
512 * Check fix for bet2 = -bet1 bug found 2011-06-25 (Visual Studio 10 rel
513 * + debug)
514  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
515  double precision a, f
516  integer r, assert, omask
517  include 'geodesic.inc'
518 
519 * WGS84 values
520  a = 6378137d0
521  f = 1/298.257223563d0
522  omask = 0
523  r = 0
524  call invers(a, f, 48.522876735459d0, 0d0,
525  + -48.52287673545898293d0, 179.599720456223079643d0,
526  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
527  r = r + assert(s12, 19989144.774d0, 0.5d-3)
528 
529  tstg11 = r
530  return
531  end
532 
533  integer function tstg12()
534 * Check fix for inverse geodesics on extreme prolate/oblate ellipsoids
535 * Reported 2012-08-29 Stefan Guenther <stefan.gunther@embl.de>; fixed
536 * 2012-10-07
537  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
538  double precision a, f
539  integer r, assert, omask
540  include 'geodesic.inc'
541 
542  a = 89.8d0
543  f = -1.83d0
544  omask = 0
545  r = 0
546  call invers(a, f, 0d0, 0d0, -10d0, 160d0,
547  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
548  r = r + assert(azi1, 120.27d0, 1d-2)
549  r = r + assert(azi2, 105.15d0, 1d-2)
550  r = r + assert(s12, 266.7d0, 1d-1)
551 
552  tstg12 = r
553  return
554  end
555 
556  integer function tstg14()
557 * Check fix for inverse ignoring lon12 = nan
558  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
559  double precision a, f, LatFix
560  integer r, chknan, omask
561  include 'geodesic.inc'
562 
563 * WGS84 values
564  a = 6378137d0
565  f = 1/298.257223563d0
566  omask = 0
567  r = 0
568  call invers(a, f, 0d0, 0d0, 1d0, latfix(91d0),
569  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
570  r = r + chknan(azi1)
571  r = r + chknan(azi2)
572  r = r + chknan(s12)
573 
574  tstg14 = r
575  return
576  end
577 
578  integer function tstg15()
579 * Initial implementation of Math::eatanhe was wrong for e^2 < 0. This
580 * checks that this is fixed.
581  double precision lat2, lon2, azi2, a12, m12, MM12, MM21, SS12
582  double precision a, f
583  integer r, assert, omask, flags
584  include 'geodesic.inc'
585 
586  a = 6.4d6
587  f = -1/150.0d0
588  omask = 8
589  flags = 0
590  r = 0
591  call direct(a, f, 1d0, 2d0, 3d0, 4d0,
592  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
593  r = r + assert(ss12, 23700d0, 0.5d0)
594 
595  tstg15 = r
596  return
597  end
598 
599  integer function tstg17()
600 * Check fix for LONG_UNROLL bug found on 2015-05-07
601  double precision lat2, lon2, azi2, a12, m12, MM12, MM21, SS12
602  double precision a, f
603  integer r, assert, omask, flags
604  include 'geodesic.inc'
605 
606 * WGS84 values
607  a = 6378137d0
608  f = 1/298.257223563d0
609  omask = 0
610  flags = 2
611  r = 0
612  call direct(a, f, 40d0, -75d0, -10d0, 2d7,
613  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
614  r = r + assert(lat2, -39d0, 1d0)
615  r = r + assert(lon2, -254d0, 1d0)
616  r = r + assert(azi2, -170d0, 1d0)
617  flags = 0
618  call direct(a, f, 40d0, -75d0, -10d0, 2d7,
619  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
620  r = r + assert(lat2, -39d0, 1d0)
621  r = r + assert(lon2, 105d0, 1d0)
622  r = r + assert(azi2, -170d0, 1d0)
623 
624  tstg17 = r
625  return
626  end
627 
628  integer function tstg26()
629 * Check 0/0 problem with area calculation on sphere 2015-09-08
630  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
631  double precision a, f
632  integer r, assert, omask
633  include 'geodesic.inc'
634 
635  a = 6.4d6
636  f = 0
637  omask = 8
638  r = 0
639  call invers(a, f, 1d0, 2d0, 3d0, 4d0,
640  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
641  r = r + assert(ss12, 49911046115.0d0, 0.5d0)
642 
643  tstg26 = r
644  return
645  end
646 
647  integer function tstg28()
648 * Check fix for LONG_UNROLL bug found on 2015-05-07
649  double precision lat2, lon2, azi2, a12, m12, MM12, MM21, SS12
650  double precision a, f
651  integer r, assert, omask, flags
652  include 'geodesic.inc'
653 
654  a = 6.4d6
655  f = 0.1d0
656  omask = 1 + 2 + 4 + 8
657  flags = 0
658  r = 0
659  call direct(a, f, 1d0, 2d0, 10d0, 5d6,
660  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
661  r = r + assert(a12, 48.55570690d0, 0.5d-8)
662 
663  tstg28 = r
664  return
665  end
666 
667  integer function tstg33()
668 * Check max(-0.0,+0.0) issues 2015-08-22 (triggered by bugs in Octave --
669 * sind(-0.0) = +0.0 -- and in some version of Visual Studio --
670 * fmod(-0.0, 360.0) = +0.0.
671  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
672  double precision a, f
673  integer r, assert, omask
674  include 'geodesic.inc'
675 
676 * WGS84 values
677  a = 6378137d0
678  f = 1/298.257223563d0
679  omask = 0
680  r = 0
681  call invers(a, f, 0d0, 0d0, 0d0, 179d0,
682  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
683  r = r + assert(azi1, 90.00000d0, 0.5d-5)
684  r = r + assert(azi2, 90.00000d0, 0.5d-5)
685  r = r + assert(s12, 19926189d0, 0.5d0)
686  call invers(a, f, 0d0, 0d0, 0d0, 179.5d0,
687  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
688  r = r + assert(azi1, 55.96650d0, 0.5d-5)
689  r = r + assert(azi2, 124.03350d0, 0.5d-5)
690  r = r + assert(s12, 19980862d0, 0.5d0)
691  call invers(a, f, 0d0, 0d0, 0d0, 180d0,
692  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
693  r = r + assert(azi1, 0.00000d0, 0.5d-5)
694  r = r + assert(abs(azi2), 180.00000d0, 0.5d-5)
695  r = r + assert(s12, 20003931d0, 0.5d0)
696  call invers(a, f, 0d0, 0d0, 1d0, 180d0,
697  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
698  r = r + assert(azi1, 0.00000d0, 0.5d-5)
699  r = r + assert(abs(azi2), 180.00000d0, 0.5d-5)
700  r = r + assert(s12, 19893357d0, 0.5d0)
701  a = 6.4d6
702  f = 0
703  call invers(a, f, 0d0, 0d0, 0d0, 179d0,
704  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
705  r = r + assert(azi1, 90.00000d0, 0.5d-5)
706  r = r + assert(azi2, 90.00000d0, 0.5d-5)
707  r = r + assert(s12, 19994492d0, 0.5d0)
708  call invers(a, f, 0d0, 0d0, 0d0, 180d0,
709  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
710  r = r + assert(azi1, 0.00000d0, 0.5d-5)
711  r = r + assert(abs(azi2), 180.00000d0, 0.5d-5)
712  r = r + assert(s12, 20106193d0, 0.5d0)
713  call invers(a, f, 0d0, 0d0, 1d0, 180d0,
714  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
715  r = r + assert(azi1, 0.00000d0, 0.5d-5)
716  r = r + assert(abs(azi2), 180.00000d0, 0.5d-5)
717  r = r + assert(s12, 19994492d0, 0.5d0)
718  a = 6.4d6
719  f = -1/300.0d0
720  call invers(a, f, 0d0, 0d0, 0d0, 179d0,
721  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
722  r = r + assert(azi1, 90.00000d0, 0.5d-5)
723  r = r + assert(azi2, 90.00000d0, 0.5d-5)
724  r = r + assert(s12, 19994492d0, 0.5d0)
725  call invers(a, f, 0d0, 0d0, 0d0, 180d0,
726  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
727  r = r + assert(azi1, 90.00000d0, 0.5d-5)
728  r = r + assert(azi2, 90.00000d0, 0.5d-5)
729  r = r + assert(s12, 20106193d0, 0.5d0)
730  call invers(a, f, 0d0, 0d0, 0.5d0, 180d0,
731  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
732  r = r + assert(azi1, 33.02493d0, 0.5d-5)
733  r = r + assert(azi2, 146.97364d0, 0.5d-5)
734  r = r + assert(s12, 20082617d0, 0.5d0)
735  call invers(a, f, 0d0, 0d0, 1d0, 180d0,
736  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
737  r = r + assert(azi1, 0.00000d0, 0.5d-5)
738  r = r + assert(abs(azi2), 180.00000d0, 0.5d-5)
739  r = r + assert(s12, 20027270d0, 0.5d0)
740 
741  tstg33 = r
742  return
743  end
744 
745  integer function tstg55()
746 * Check fix for nan + point on equator or pole not returning all nans in
747 * Geodesic::Inverse, found 2015-09-23.
748  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
749  double precision a, f
750  integer r, chknan, omask
751  include 'geodesic.inc'
752 
753 * WGS84 values
754  a = 6378137d0
755  f = 1/298.257223563d0
756  omask = 0
757  r = 0
758  call invers(a, f, 91d0, 0d0, 0d0, 90d0,
759  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
760  r = r + chknan(azi1)
761  r = r + chknan(azi2)
762  r = r + chknan(s12)
763  call invers(a, f, 91d0, 0d0, 90d0, 9d0,
764  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
765  r = r + chknan(azi1)
766  r = r + chknan(azi2)
767  r = r + chknan(s12)
768 
769  tstg55 = r
770  return
771  end
772 
773  integer function tstg59()
774 * Check for points close with longitudes close to 180 deg apart.
775  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
776  double precision a, f
777  integer r, assert, omask
778  include 'geodesic.inc'
779 
780 * WGS84 values
781  a = 6378137d0
782  f = 1/298.257223563d0
783  omask = 0
784  r = 0
785  call invers(a, f, 5d0, 0.00000000000001d0, 10d0, 180d0,
786  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
787  r = r + assert(azi1, 0.000000000000035d0, 1.5d-14)
788  r = r + assert(azi2, 179.99999999999996d0, 1.5d-14)
789  r = r + assert(s12, 18345191.174332713d0, 5d-9)
790 
791  tstg59 = r
792  return
793  end
794 
795  integer function tstg61()
796 * Make sure small negative azimuths are west-going
797  double precision lat2, lon2, azi2, a12, m12, MM12, MM21, SS12
798  double precision a, f
799  integer r, assert, omask, flags
800  include 'geodesic.inc'
801 
802 * WGS84 values
803  a = 6378137d0
804  f = 1/298.257223563d0
805  omask = 0
806  flags = 2
807  r = 0
808  call direct(a, f, 45d0, 0d0, -0.000000000000000003d0, 1d7,
809  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
810  r = r + assert(lat2, 45.30632d0, 0.5d-5)
811  r = r + assert(lon2, -180d0, 0.5d-5)
812  r = r + assert(abs(azi2), 180d0, 0.5d-5)
813 
814  tstg61 = r
815  return
816  end
817 
818  integer function tstg73()
819 * Check for backwards from the pole bug reported by Anon on 2016-02-13.
820 * This only affected the Java implementation. It was introduced in Java
821 * version 1.44 and fixed in 1.46-SNAPSHOT on 2016-01-17.
822 * Also the + sign on azi2 is a check on the normalizing of azimuths
823 * (converting -0.0 to +0.0).
824  double precision lat2, lon2, azi2, a12, m12, MM12, MM21, SS12
825  double precision a, f
826  integer r, assert, omask, flags
827  include 'geodesic.inc'
828 
829 * WGS84 values
830  a = 6378137d0
831  f = 1/298.257223563d0
832  omask = 0
833  flags = 0
834  r = 0
835  call direct(a, f, 90d0, 10d0, 180d0, -1d6,
836  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
837  r = r + assert(lat2, 81.04623d0, 0.5d-5)
838  r = r + assert(lon2, -170d0, 0.5d-5)
839  r = r + assert(azi2, 0d0, 0d0)
840  r = r + assert(sign(1d0, azi2), 1d0, 0d0)
841 
842  tstg73 = r
843  return
844  end
845 
846  integer function tstg74()
847 * Check fix for inaccurate areas, bug introduced in v1.46, fixed
848 * 2015-10-16.
849  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
850  double precision a, f
851  integer r, assert, omask
852  include 'geodesic.inc'
853 
854 * WGS84 values
855  a = 6378137d0
856  f = 1/298.257223563d0
857  omask = 1 + 2 + 4 + 8
858  r = 0
859  call invers(a, f, 54.1589d0, 15.3872d0, 54.1591d0, 15.3877d0,
860  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
861  r = r + assert(azi1, 55.723110355d0, 5d-9)
862  r = r + assert(azi2, 55.723515675d0, 5d-9)
863  r = r + assert(s12, 39.527686385d0, 5d-9)
864  r = r + assert(a12, 0.000355495d0, 5d-9)
865  r = r + assert(m12, 39.527686385d0, 5d-9)
866  r = r + assert(mm12, 0.999999995d0, 5d-9)
867  r = r + assert(mm21, 0.999999995d0, 5d-9)
868  r = r + assert(ss12, 286698586.30197d0, 5d-4)
869 
870  tstg74 = r
871  return
872  end
873 
874  integer function tstg76()
875 * The distance from Wellington and Salamanca (a classic failure of
876 * Vincenty
877  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
878  double precision a, f
879  integer r, assert, omask
880  include 'geodesic.inc'
881 
882 * WGS84 values
883  a = 6378137d0
884  f = 1/298.257223563d0
885  omask = 0
886  r = 0
887  call invers(a, f,
888  + -(41+19/60d0), 174+49/60d0, 40+58/60d0, -(5+30/60d0),
889  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
890  r = r + assert(azi1, 160.39137649664d0, 0.5d-11)
891  r = r + assert(azi2, 19.50042925176d0, 0.5d-11)
892  r = r + assert(s12, 19960543.857179d0, 0.5d-6)
893 
894  tstg76 = r
895  return
896  end
897 
898  integer function tstg78()
899 * An example where the NGS calculator fails to converge
900  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
901  double precision a, f
902  integer r, assert, omask
903  include 'geodesic.inc'
904 
905 * WGS84 values
906  a = 6378137d0
907  f = 1/298.257223563d0
908  omask = 0
909  r = 0
910  call invers(a, f, 27.2d0, 0d0, -27.1d0, 179.5d0,
911  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
912  r = r + assert(azi1, 45.82468716758d0, 0.5d-11)
913  r = r + assert(azi2, 134.22776532670d0, 0.5d-11)
914  r = r + assert(s12, 19974354.765767d0, 0.5d-6)
915 
916  tstg78 = r
917  return
918  end
919 
920  integer function tstg80()
921 * Some tests to add code coverage: computing scale in special cases + zero
922 * length geodesic (includes GeodSolve80 - GeodSolve83).
923  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
924  double precision a, f
925  integer r, assert, omask
926  include 'geodesic.inc'
927 
928 * WGS84 values
929  a = 6378137d0
930  f = 1/298.257223563d0
931  omask = 4
932  r = 0
933 
934  call invers(a, f, 0d0, 0d0, 0d0, 90d0,
935  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
936  r = r + assert(mm12, -0.00528427534d0, 0.5d-10)
937  r = r + assert(mm21, -0.00528427534d0, 0.5d-10)
938 
939  call invers(a, f, 0d0, 0d0, 1d-6, 1d-6,
940  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
941  r = r + assert(mm12, 1d0, 0.5d-10)
942  r = r + assert(mm21, 1d0, 0.5d-10)
943 
944  omask = 15
945  call invers(a, f, 20.001d0, 0d0, 20.001d0, 0d0,
946  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
947  r = r + assert(a12, 0d0, 1d-13)
948  r = r + assert(s12, 0d0, 1d-8)
949  r = r + assert(azi1, 180d0, 1d-13)
950  r = r + assert(azi2, 180d0, 1d-13)
951  r = r + assert(m12, 0d0, 1d-8)
952  r = r + assert(mm12, 1d0, 1d-15)
953  r = r + assert(mm21, 1d0, 1d-15)
954  r = r + assert(ss12, 0d0, 1d-10)
955  r = r + assert(sign(1d0, a12), 1d0, 0d0)
956  r = r + assert(sign(1d0, s12), 1d0, 0d0)
957  r = r + assert(sign(1d0, m12), 1d0, 0d0)
958 
959  call invers(a, f, 90d0, 0d0, 90d0, 180d0,
960  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
961  r = r + assert(a12, 0d0, 1d-13)
962  r = r + assert(s12, 0d0, 1d-8)
963  r = r + assert(azi1, 0d0, 1d-13)
964  r = r + assert(azi2, 180d0, 1d-13)
965  r = r + assert(m12, 0d0, 1d-8)
966  r = r + assert(mm12, 1d0, 1d-15)
967  r = r + assert(mm21, 1d0, 1d-15)
968  r = r + assert(ss12, 127516405431022d0, 0.5d0)
969 
970  tstg80 = r
971  return
972  end
973 
974  integer function tstg84()
975 * Tests for python implementation to check fix for range errors with
976 * {fmod,sin,cos}(inf) (includes GeodSolve84 - GeodSolve86).
977  double precision lat2, lon2, azi2, a12, m12, MM12, MM21, SS12
978  double precision a, f, nan, inf, LatFix
979  integer r, assert, chknan, omask, flags
980  include 'geodesic.inc'
981 
982 * WGS84 values
983  a = 6378137d0
984  f = 1/298.257223563d0
985  omask = 0
986  flags = 0
987  inf = 1d0/latfix(0d0)
988  nan = latfix(91d0)
989  r = 0
990  call direct(a, f, 0d0, 0d0, 90d0, inf,
991  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
992  r = r + chknan(lat2)
993  r = r + chknan(lon2)
994  r = r + chknan(azi2)
995  call direct(a, f, 0d0, 0d0, 90d0, nan,
996  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
997  r = r + chknan(lat2)
998  r = r + chknan(lon2)
999  r = r + chknan(azi2)
1000  call direct(a, f, 0d0, 0d0, inf, 1000d0,
1001  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
1002  r = r + chknan(lat2)
1003  r = r + chknan(lon2)
1004  r = r + chknan(azi2)
1005  call direct(a, f, 0d0, 0d0, nan, 1000d0,
1006  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
1007  r = r + chknan(lat2)
1008  r = r + chknan(lon2)
1009  r = r + chknan(azi2)
1010  call direct(a, f, 0d0, inf, 90d0, 1000d0,
1011  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
1012  r = r + assert(lat2, 0d0, 0d0)
1013  r = r + chknan(lon2)
1014  r = r + assert(azi2, 90d0, 0d0)
1015  call direct(a, f, 0d0, nan, 90d0, 1000d0,
1016  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
1017  r = r + assert(lat2, 0d0, 0d0)
1018  r = r + chknan(lon2)
1019  r = r + assert(azi2, 90d0, 0d0)
1020  call direct(a, f, inf, 0d0, 90d0, 1000d0,
1021  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
1022  r = r + chknan(lat2)
1023  r = r + chknan(lon2)
1024  r = r + chknan(azi2)
1025  call direct(a, f, nan, 0d0, 90d0, 1000d0,
1026  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
1027  r = r + chknan(lat2)
1028  r = r + chknan(lon2)
1029  r = r + chknan(azi2)
1030 
1031  tstg84 = r
1032  return
1033  end
1034 
1035  integer function tstg92()
1036 * Check fix for inaccurate hypot with python 3.[89]. Problem reported
1037 * by agdhruv https://github.com/geopy/geopy/issues/466 ; see
1038 * https://bugs.python.org/issue43088
1039  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
1040  double precision a, f
1041  integer r, assert, omask
1042  include 'geodesic.inc'
1043 
1044 * WGS84 values
1045  a = 6378137d0
1046  f = 1/298.257223563d0
1047  omask = 0
1048  r = 0
1049  call invers(a, f,
1050  + 37.757540000000006d0, -122.47018d0,
1051  + 37.75754d0, -122.470177d0,
1052  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
1053  r = r + assert(azi1, 89.99999923d0, 1d-7 )
1054  r = r + assert(azi2, 90.00000106d0, 1d-7 )
1055  r = r + assert(s12, 0.264d0, 0.5d-3)
1056 
1057  tstg92 = r
1058  return
1059  end
1060 
1061  integer function tstg94()
1062 * Check fix for lat2 = nan being treated as lat2 = 0 (bug found 2021-07-26)
1063 
1064  double precision azi1, azi2, s12, a12, m12, MM12, MM21, SS12
1065  double precision a, f, LatFix
1066  integer r, chknan, omask
1067  include 'geodesic.inc'
1068 
1069 * WGS84 values
1070  a = 6378137d0
1071  f = 1/298.257223563d0
1072  omask = 0
1073  r = 0
1074  call invers(a, f, 0d0, 0d0, latfix(91d0), 90d0,
1075  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
1076  r = r + chknan(azi1)
1077  r = r + chknan(azi2)
1078  r = r + chknan(s12)
1079 
1080  tstg94 = r
1081  return
1082  end
1083 
1084  integer function tstp0()
1085 * Check fix for pole-encircling bug found 2011-03-16
1086  double precision lata(4), lona(4)
1087  data lata / 89d0, 89d0, 89d0, 89d0 /
1088  data lona / 0d0, 90d0, 180d0, 270d0 /
1089  double precision latb(4), lonb(4)
1090  data latb / -89d0, -89d0, -89d0, -89d0 /
1091  data lonb / 0d0, 90d0, 180d0, 270d0 /
1092  double precision latc(4), lonc(4)
1093  data latc / 0d0, -1d0, 0d0, 1d0 /
1094  data lonc / -1d0, 0d0, 1d0, 0d0 /
1095  double precision latd(3), lond(3)
1096  data latd / 90d0, 0d0, 0d0 /
1097  data lond / 0d0, 0d0, 90d0 /
1098  double precision a, f, AA, PP
1099  integer r, assert
1100  include 'geodesic.inc'
1101 
1102 * WGS84 values
1103  a = 6378137d0
1104  f = 1/298.257223563d0
1105  r = 0
1106 
1107  call area(a, f, lata, lona, 4, aa, pp)
1108  r = r + assert(pp, 631819.8745d0, 1d-4)
1109  r = r + assert(aa, 24952305678.0d0, 1d0)
1110 
1111  call area(a, f, latb, lonb, 4, aa, pp)
1112  r = r + assert(pp, 631819.8745d0, 1d-4)
1113  r = r + assert(aa, -24952305678.0d0, 1d0)
1114 
1115  call area(a, f, latc, lonc, 4, aa, pp)
1116  r = r + assert(pp, 627598.2731d0, 1d-4)
1117  r = r + assert(aa, 24619419146.0d0, 1d0)
1118 
1119  call area(a, f, latd, lond, 3, aa, pp)
1120  r = r + assert(pp, 30022685d0, 1d0)
1121  r = r + assert(aa, 63758202715511.0d0, 1d0)
1122 
1123  tstp0 = r
1124  return
1125  end
1126 
1127  integer function tstp5()
1128 * Check fix for Planimeter pole crossing bug found 2011-06-24
1129  double precision lat(3), lon(3)
1130  data lat / 89d0, 89d0, 89d0 /
1131  data lon / 0.1d0, 90.1d0, -179.9d0 /
1132  double precision a, f, AA, PP
1133  integer r, assert
1134  include 'geodesic.inc'
1135 
1136 * WGS84 values
1137  a = 6378137d0
1138  f = 1/298.257223563d0
1139  r = 0
1140 
1141  call area(a, f, lat, lon, 3, aa, pp)
1142  r = r + assert(pp, 539297d0, 1d0)
1143  r = r + assert(aa, 12476152838.5d0, 1d0)
1144 
1145  tstp5 = r
1146  return
1147  end
1148 
1149  integer function tstp6()
1150 * Check fix for pole-encircling bug found 2011-03-16
1151  double precision lata(3), lona(3)
1152  data lata / 9d0, 9d0, 9d0 /
1153  data lona / -0.00000000000001d0, 180d0, 0d0 /
1154  double precision latb(3), lonb(3)
1155  data latb / 9d0, 9d0, 9d0 /
1156  data lonb / 0.00000000000001d0, 0d0, 180d0 /
1157  double precision latc(3), lonc(3)
1158  data latc / 9d0, 9d0, 9d0 /
1159  data lonc / 0.00000000000001d0, 180d0, 0d0 /
1160  double precision latd(3), lond(3)
1161  data latd / 9d0, 9d0, 9d0 /
1162  data lond / -0.00000000000001d0, 0d0, 180d0 /
1163  double precision a, f, AA, PP
1164  integer r, assert
1165  include 'geodesic.inc'
1166 
1167 * WGS84 values
1168  a = 6378137d0
1169  f = 1/298.257223563d0
1170  r = 0
1171 
1172  call area(a, f, lata, lona, 3, aa, pp)
1173  r = r + assert(pp, 36026861d0, 1d0)
1174  r = r + assert(aa, 0d0, 1d0)
1175 
1176  tstp6 = r
1177  return
1178  end
1179 
1180  integer function tstp12()
1181 * AA of arctic circle (not really -- adjunct to rhumb-AA test)
1182  double precision lat(2), lon(2)
1183  data lat / 66.562222222d0, 66.562222222d0 /
1184  data lon / 0d0, 180d0 /
1185  double precision a, f, AA, PP
1186  integer r, assert
1187  include 'geodesic.inc'
1188 
1189 * WGS84 values
1190  a = 6378137d0
1191  f = 1/298.257223563d0
1192  r = 0
1193 
1194  call area(a, f, lat, lon, 2, aa, pp)
1195  r = r + assert(pp, 10465729d0, 1d0)
1196  r = r + assert(aa, 0d0, 1d0)
1197 
1198  tstp12 = r
1199  return
1200  end
1201 
1202  integer function tstp13()
1203 * Check encircling pole twice
1204  double precision lat(6), lon(6)
1205  data lat / 89d0, 89d0, 89d0, 89d0, 89d0, 89d0 /
1206  data lon / -360d0, -240d0, -120d0, 0d0, 120d0, 240d0 /
1207  double precision a, f, AA, PP
1208  integer r, assert
1209  include 'geodesic.inc'
1210 
1211 * WGS84 values
1212  a = 6378137d0
1213  f = 1/298.257223563d0
1214  r = 0
1215 
1216  call area(a, f, lat, lon, 6, aa, pp)
1217  r = r + assert(pp, 1160741d0, 1d0)
1218  r = r + assert(aa, 32415230256.0d0, 1d0)
1219 
1220  tstp13 = r
1221  return
1222  end
1223 
1224  integer function tstp15()
1225 * Coverage tests, includes Planimeter15 - Planimeter18 (combinations of
1226 * reverse and sign). But flags aren't supported in the Fortran
1227 * implementation.
1228  double precision lat(3), lon(3)
1229  data lat / 2d0, 1d0, 3d0 /
1230  data lon / 1d0, 2d0, 3d0 /
1231  double precision a, f, AA, PP
1232  integer r, assert
1233  include 'geodesic.inc'
1234 
1235 * WGS84 values
1236  a = 6378137d0
1237  f = 1/298.257223563d0
1238  r = 0
1239 
1240  call area(a, f, lat, lon, 3, aa, pp)
1241  r = r + assert(aa, 18454562325.45119d0, 1d0)
1242 * Interchanging lat and lon is equivalent to traversing the polygon
1243 * backwards.
1244  call area(a, f, lon, lat, 3, aa, pp)
1245  r = r + assert(aa, -18454562325.45119d0, 1d0)
1246 
1247  tstp15 = r
1248  return
1249  end
1250 
1251  integer function tstp19()
1252 * Coverage tests, includes Planimeter19 - Planimeter20 (degenerate
1253 * polygons).
1254  double precision lat(1), lon(1)
1255  data lat / 1d0 /
1256  data lon / 1d0 /
1257  double precision a, f, AA, PP
1258  integer r, assert
1259  include 'geodesic.inc'
1260 
1261 * WGS84 values
1262  a = 6378137d0
1263  f = 1/298.257223563d0
1264  r = 0
1265 
1266  call area(a, f, lat, lon, 1, aa, pp)
1267  r = r + assert(aa, 0d0, 0d0)
1268  r = r + assert(pp, 0d0, 0d0)
1269 
1270  tstp19 = r
1271  return
1272  end
1273 
1274  integer function tstp21()
1275 * Some test to add code coverage: multiple circlings of pole (includes
1276 * Planimeter21 - Planimeter28).
1277  double precision lat(12), lon(12), lonr(12)
1278  data lat / 12*45d0 /
1279  data lon / 60d0, 180d0, -60d0,
1280  + 60d0, 180d0, -60d0,
1281  + 60d0, 180d0, -60d0,
1282  + 60d0, 180d0, -60d0 /
1283  data lonr / -60d0, 180d0, 60d0,
1284  + -60d0, 180d0, 60d0,
1285  + -60d0, 180d0, 60d0,
1286  + -60d0, 180d0, 60d0 /
1287  double precision a, f, AA, PP, AA1
1288  integer i, r, assert
1289  include 'geodesic.inc'
1290 
1291 * WGS84 values
1292  a = 6378137d0
1293  f = 1/298.257223563d0
1294  r = 0
1295 * Area for one circuit
1296  aa1 = 39433884866571.4277d0
1297 
1298  do 10 i = 3,4
1299  call area(a, f, lat, lon, 3*i, aa, pp)
1300  r = r + assert(aa, aa1*i, 0.5d0)
1301  call area(a, f, lat, lonr, 3*i, aa, pp)
1302  r = r + assert(aa, -aa1*i, 0.5d0)
1303  10 continue
1304 
1305  tstp21 = r
1306  return
1307  end
1308 
1309  program geodtest
1310  integer n, i
1311  integer tstinv, tstdir, tstarc,
1312  + tstg0, tstg1, tstg2, tstg5, tstg6, tstg9, tstg10, tstg11,
1313  + tstg12, tstg14, tstg15, tstg17, tstg26, tstg28, tstg33,
1314  + tstg55, tstg59, tstg61, tstg73, tstg74, tstg76, tstg78,
1315  + tstg80, tstg84, tstg92, tstg94,
1316  + tstp0, tstp5, tstp6, tstp12, tstp13, tstp15, tstp19, tstp21
1317 
1318  n = 0
1319  i = tstinv()
1320  if (i .gt. 0) then
1321  n = n + 1
1322  print *, 'tstinv fail:', i
1323  end if
1324  i = tstdir()
1325  if (i .gt. 0) then
1326  n = n + 1
1327  print *, 'tstdir fail:', i
1328  end if
1329  i = tstarc()
1330  if (i .gt. 0) then
1331  n = n + 1
1332  print *, 'tstarc fail:', i
1333  end if
1334  i = tstg0()
1335  if (i .gt. 0) then
1336  n = n + 1
1337  print *, 'tstg0 fail:', i
1338  end if
1339  i = tstg1()
1340  if (i .gt. 0) then
1341  n = n + 1
1342  print *, 'tstg1 fail:', i
1343  end if
1344  i = tstg2()
1345  if (i .gt. 0) then
1346  n = n + 1
1347  print *, 'tstg2 fail:', i
1348  end if
1349  i = tstg5()
1350  if (i .gt. 0) then
1351  n = n + 1
1352  print *, 'tstg5 fail:', i
1353  end if
1354  i = tstg6()
1355  if (i .gt. 0) then
1356  n = n + 1
1357  print *, 'tstg6 fail:', i
1358  end if
1359  i = tstg9()
1360  if (i .gt. 0) then
1361  n = n + 1
1362  print *, 'tstg9 fail:', i
1363  end if
1364  i = tstg10()
1365  if (i .gt. 0) then
1366  n = n + 1
1367  print *, 'tstg10 fail:', i
1368  end if
1369  i = tstg11()
1370  if (i .gt. 0) then
1371  n = n + 1
1372  print *, 'tstg11 fail:', i
1373  end if
1374  i = tstg12()
1375  if (i .gt. 0) then
1376  n = n + 1
1377  print *, 'tstg12 fail:', i
1378  end if
1379  i = tstg14()
1380  if (i .gt. 0) then
1381  n = n + 1
1382  print *, 'tstg14 fail:', i
1383  end if
1384  i = tstg15()
1385  if (i .gt. 0) then
1386  n = n + 1
1387  print *, 'tstg15 fail:', i
1388  end if
1389  i = tstg17()
1390  if (i .gt. 0) then
1391  n = n + 1
1392  print *, 'tstg17 fail:', i
1393  end if
1394  i = tstg26()
1395  if (i .gt. 0) then
1396  n = n + 1
1397  print *, 'tstg26 fail:', i
1398  end if
1399  i = tstg28()
1400  if (i .gt. 0) then
1401  n = n + 1
1402  print *, 'tstg28 fail:', i
1403  end if
1404  i = tstg33()
1405  if (i .gt. 0) then
1406  n = n + 1
1407  print *, 'tstg33 fail:', i
1408  end if
1409  i = tstg55()
1410  if (i .gt. 0) then
1411  n = n + 1
1412  print *, 'tstg55 fail:', i
1413  end if
1414  i = tstg59()
1415  if (i .gt. 0) then
1416  n = n + 1
1417  print *, 'tstg59 fail:', i
1418  end if
1419  i = tstg61()
1420  if (i .gt. 0) then
1421  n = n + 1
1422  print *, 'tstg61 fail:', i
1423  end if
1424  i = tstg73()
1425  if (i .gt. 0) then
1426  n = n + 1
1427  print *, 'tstg73 fail:', i
1428  end if
1429  i = tstg74()
1430  if (i .gt. 0) then
1431  n = n + 1
1432  print *, 'tstg74 fail:', i
1433  end if
1434  i = tstg76()
1435  if (i .gt. 0) then
1436  n = n + 1
1437  print *, 'tstg76 fail:', i
1438  end if
1439  i = tstg78()
1440  if (i .gt. 0) then
1441  n = n + 1
1442  print *, 'tstg78 fail:', i
1443  end if
1444  i = tstg80()
1445  if (i .gt. 0) then
1446  n = n + 1
1447  print *, 'tstg80 fail:', i
1448  end if
1449  i = tstg84()
1450  if (i .gt. 0) then
1451  n = n + 1
1452  print *, 'tstg84 fail:', i
1453  end if
1454  i = tstg92()
1455  if (i .gt. 0) then
1456  n = n + 1
1457  print *, 'tstg92 fail:', i
1458  end if
1459  i = tstg94()
1460  if (i .gt. 0) then
1461  n = n + 1
1462  print *, 'tstg94 fail:', i
1463  end if
1464  i = tstp0()
1465  if (i .gt. 0) then
1466  n = n + 1
1467  print *, 'tstp0 fail:', i
1468  end if
1469  i = tstp5()
1470  if (i .gt. 0) then
1471  n = n + 1
1472  print *, 'tstp5 fail:', i
1473  end if
1474  i = tstp6()
1475  if (i .gt. 0) then
1476  n = n + 1
1477  print *, 'tstp6 fail:', i
1478  end if
1479  i = tstp12()
1480  if (i .gt. 0) then
1481  n = n + 1
1482  print *, 'tstp12 fail:', i
1483  end if
1484  i = tstp13()
1485  if (i .gt. 0) then
1486  n = n + 1
1487  print *, 'tstp13 fail:', i
1488  end if
1489  i = tstp15()
1490  if (i .gt. 0) then
1491  n = n + 1
1492  print *, 'tstp15 fail:', i
1493  end if
1494  i = tstp19()
1495  if (i .gt. 0) then
1496  n = n + 1
1497  print *, 'tstp19 fail:', i
1498  end if
1499  i = tstp21()
1500  if (i .gt. 0) then
1501  n = n + 1
1502  print *, 'tstp21 fail:', i
1503  end if
1504 
1505  if (n .gt. 0) then
1506  stop 1
1507  end if
1508 
1509  stop
1510  end
1511 
1512 *> @endcond SKIP
direct
Definition: geodesic.inc:15
area
Definition: geodesic.inc:31
invers
Definition: geodesic.inc:23
Generated by doxygen 1.8.20