Fortran library for Geodesics  1.50.1
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-2019) <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 
956  call invers(a, f, 90d0, 0d0, 90d0, 180d0,
957  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
958  r = r + assert(a12, 0d0, 1d-13)
959  r = r + assert(s12, 0d0, 1d-8)
960  r = r + assert(azi1, 0d0, 1d-13)
961  r = r + assert(azi2, 180d0, 1d-13)
962  r = r + assert(m12, 0d0, 1d-8)
963  r = r + assert(mm12, 1d0, 1d-15)
964  r = r + assert(mm21, 1d0, 1d-15)
965  r = r + assert(ss12, 127516405431022d0, 0.5d0)
966 
967  tstg80 = r
968  return
969  end
970 
971  integer function tstg84()
972 * Tests for python implementation to check fix for range errors with
973 * {fmod,sin,cos}(inf) (includes GeodSolve84 - GeodSolve86).
974  double precision lat2, lon2, azi2, a12, m12, MM12, MM21, SS12
975  double precision a, f, nan, inf, LatFix
976  integer r, assert, chknan, omask, flags
977  include 'geodesic.inc'
978 
979 * WGS84 values
980  a = 6378137d0
981  f = 1/298.257223563d0
982  omask = 0
983  flags = 0
984  inf = 1d0/latfix(0d0)
985  nan = latfix(91d0)
986  r = 0
987  call direct(a, f, 0d0, 0d0, 90d0, inf,
988  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
989  r = r + chknan(lat2)
990  r = r + chknan(lon2)
991  r = r + chknan(azi2)
992  call direct(a, f, 0d0, 0d0, 90d0, nan,
993  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
994  r = r + chknan(lat2)
995  r = r + chknan(lon2)
996  r = r + chknan(azi2)
997  call direct(a, f, 0d0, 0d0, inf, 1000d0,
998  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
999  r = r + chknan(lat2)
1000  r = r + chknan(lon2)
1001  r = r + chknan(azi2)
1002  call direct(a, f, 0d0, 0d0, nan, 1000d0,
1003  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
1004  r = r + chknan(lat2)
1005  r = r + chknan(lon2)
1006  r = r + chknan(azi2)
1007  call direct(a, f, 0d0, inf, 90d0, 1000d0,
1008  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
1009  r = r + assert(lat2, 0d0, 0d0)
1010  r = r + chknan(lon2)
1011  r = r + assert(azi2, 90d0, 0d0)
1012  call direct(a, f, 0d0, nan, 90d0, 1000d0,
1013  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
1014  r = r + assert(lat2, 0d0, 0d0)
1015  r = r + chknan(lon2)
1016  r = r + assert(azi2, 90d0, 0d0)
1017  call direct(a, f, inf, 0d0, 90d0, 1000d0,
1018  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
1019  r = r + chknan(lat2)
1020  r = r + chknan(lon2)
1021  r = r + chknan(azi2)
1022  call direct(a, f, nan, 0d0, 90d0, 1000d0,
1023  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
1024  r = r + chknan(lat2)
1025  r = r + chknan(lon2)
1026  r = r + chknan(azi2)
1027 
1028  tstg84 = r
1029  return
1030  end
1031 
1032  integer function tstp0()
1033 * Check fix for pole-encircling bug found 2011-03-16
1034  double precision lata(4), lona(4)
1035  data lata / 89d0, 89d0, 89d0, 89d0 /
1036  data lona / 0d0, 90d0, 180d0, 270d0 /
1037  double precision latb(4), lonb(4)
1038  data latb / -89d0, -89d0, -89d0, -89d0 /
1039  data lonb / 0d0, 90d0, 180d0, 270d0 /
1040  double precision latc(4), lonc(4)
1041  data latc / 0d0, -1d0, 0d0, 1d0 /
1042  data lonc / -1d0, 0d0, 1d0, 0d0 /
1043  double precision latd(3), lond(3)
1044  data latd / 90d0, 0d0, 0d0 /
1045  data lond / 0d0, 0d0, 90d0 /
1046  double precision a, f, AA, PP
1047  integer r, assert
1048  include 'geodesic.inc'
1049 
1050 * WGS84 values
1051  a = 6378137d0
1052  f = 1/298.257223563d0
1053  r = 0
1054 
1055  call area(a, f, lata, lona, 4, aa, pp)
1056  r = r + assert(pp, 631819.8745d0, 1d-4)
1057  r = r + assert(aa, 24952305678.0d0, 1d0)
1058 
1059  call area(a, f, latb, lonb, 4, aa, pp)
1060  r = r + assert(pp, 631819.8745d0, 1d-4)
1061  r = r + assert(aa, -24952305678.0d0, 1d0)
1062 
1063  call area(a, f, latc, lonc, 4, aa, pp)
1064  r = r + assert(pp, 627598.2731d0, 1d-4)
1065  r = r + assert(aa, 24619419146.0d0, 1d0)
1066 
1067  call area(a, f, latd, lond, 3, aa, pp)
1068  r = r + assert(pp, 30022685d0, 1d0)
1069  r = r + assert(aa, 63758202715511.0d0, 1d0)
1070 
1071  tstp0 = r
1072  return
1073  end
1074 
1075  integer function tstp5()
1076 * Check fix for Planimeter pole crossing bug found 2011-06-24
1077  double precision lat(3), lon(3)
1078  data lat / 89d0, 89d0, 89d0 /
1079  data lon / 0.1d0, 90.1d0, -179.9d0 /
1080  double precision a, f, AA, PP
1081  integer r, assert
1082  include 'geodesic.inc'
1083 
1084 * WGS84 values
1085  a = 6378137d0
1086  f = 1/298.257223563d0
1087  r = 0
1088 
1089  call area(a, f, lat, lon, 3, aa, pp)
1090  r = r + assert(pp, 539297d0, 1d0)
1091  r = r + assert(aa, 12476152838.5d0, 1d0)
1092 
1093  tstp5 = r
1094  return
1095  end
1096 
1097  integer function tstp6()
1098 * Check fix for pole-encircling bug found 2011-03-16
1099  double precision lata(3), lona(3)
1100  data lata / 9d0, 9d0, 9d0 /
1101  data lona / -0.00000000000001d0, 180d0, 0d0 /
1102  double precision latb(3), lonb(3)
1103  data latb / 9d0, 9d0, 9d0 /
1104  data lonb / 0.00000000000001d0, 0d0, 180d0 /
1105  double precision latc(3), lonc(3)
1106  data latc / 9d0, 9d0, 9d0 /
1107  data lonc / 0.00000000000001d0, 180d0, 0d0 /
1108  double precision latd(3), lond(3)
1109  data latd / 9d0, 9d0, 9d0 /
1110  data lond / -0.00000000000001d0, 0d0, 180d0 /
1111  double precision a, f, AA, PP
1112  integer r, assert
1113  include 'geodesic.inc'
1114 
1115 * WGS84 values
1116  a = 6378137d0
1117  f = 1/298.257223563d0
1118  r = 0
1119 
1120  call area(a, f, lata, lona, 3, aa, pp)
1121  r = r + assert(pp, 36026861d0, 1d0)
1122  r = r + assert(aa, 0d0, 1d0)
1123 
1124  tstp6 = r
1125  return
1126  end
1127 
1128  integer function tstp12()
1129 * AA of arctic circle (not really -- adjunct to rhumb-AA test)
1130  double precision lat(2), lon(2)
1131  data lat / 66.562222222d0, 66.562222222d0 /
1132  data lon / 0d0, 180d0 /
1133  double precision a, f, AA, PP
1134  integer r, assert
1135  include 'geodesic.inc'
1136 
1137 * WGS84 values
1138  a = 6378137d0
1139  f = 1/298.257223563d0
1140  r = 0
1141 
1142  call area(a, f, lat, lon, 2, aa, pp)
1143  r = r + assert(pp, 10465729d0, 1d0)
1144  r = r + assert(aa, 0d0, 1d0)
1145 
1146  tstp12 = r
1147  return
1148  end
1149 
1150  integer function tstp13()
1151 * Check encircling pole twice
1152  double precision lat(6), lon(6)
1153  data lat / 89d0, 89d0, 89d0, 89d0, 89d0, 89d0 /
1154  data lon / -360d0, -240d0, -120d0, 0d0, 120d0, 240d0 /
1155  double precision a, f, AA, PP
1156  integer r, assert
1157  include 'geodesic.inc'
1158 
1159 * WGS84 values
1160  a = 6378137d0
1161  f = 1/298.257223563d0
1162  r = 0
1163 
1164  call area(a, f, lat, lon, 6, aa, pp)
1165  r = r + assert(pp, 1160741d0, 1d0)
1166  r = r + assert(aa, 32415230256.0d0, 1d0)
1167 
1168  tstp13 = r
1169  return
1170  end
1171 
1172  integer function tstp15()
1173 * Coverage tests, includes Planimeter15 - Planimeter18 (combinations of
1174 * reverse and sign). But flags aren't supported in the Fortran
1175 * implementation.
1176  double precision lat(3), lon(3)
1177  data lat / 2d0, 1d0, 3d0 /
1178  data lon / 1d0, 2d0, 3d0 /
1179  double precision a, f, AA, PP
1180  integer r, assert
1181  include 'geodesic.inc'
1182 
1183 * WGS84 values
1184  a = 6378137d0
1185  f = 1/298.257223563d0
1186  r = 0
1187 
1188  call area(a, f, lat, lon, 3, aa, pp)
1189  r = r + assert(aa, 18454562325.45119d0, 1d0)
1190 * Interchanging lat and lon is equivalent to traversing the polygon
1191 * backwards.
1192  call area(a, f, lon, lat, 3, aa, pp)
1193  r = r + assert(aa, -18454562325.45119d0, 1d0)
1194 
1195  tstp15 = r
1196  return
1197  end
1198 
1199  integer function tstp19()
1200 * Coverage tests, includes Planimeter19 - Planimeter20 (degenerate
1201 * polygons).
1202  double precision lat(1), lon(1)
1203  data lat / 1d0 /
1204  data lon / 1d0 /
1205  double precision a, f, AA, PP
1206  integer r, assert
1207  include 'geodesic.inc'
1208 
1209 * WGS84 values
1210  a = 6378137d0
1211  f = 1/298.257223563d0
1212  r = 0
1213 
1214  call area(a, f, lat, lon, 1, aa, pp)
1215  r = r + assert(aa, 0d0, 0d0)
1216  r = r + assert(pp, 0d0, 0d0)
1217 
1218  tstp19 = r
1219  return
1220  end
1221 
1222  integer function tstp21()
1223 * Some test to add code coverage: multiple circlings of pole (includes
1224 * Planimeter21 - Planimeter28).
1225  double precision lat(12), lon(12), lonr(12)
1226  data lat / 12*45d0 /
1227  data lon / 60d0, 180d0, -60d0,
1228  + 60d0, 180d0, -60d0,
1229  + 60d0, 180d0, -60d0,
1230  + 60d0, 180d0, -60d0 /
1231  data lonr / -60d0, 180d0, 60d0,
1232  + -60d0, 180d0, 60d0,
1233  + -60d0, 180d0, 60d0,
1234  + -60d0, 180d0, 60d0 /
1235  double precision a, f, AA, PP, AA1
1236  integer r, assert
1237  include 'geodesic.inc'
1238 
1239 * WGS84 values
1240  a = 6378137d0
1241  f = 1/298.257223563d0
1242  r = 0
1243 * Area for one circuit
1244  aa1 = 39433884866571.4277d0
1245 
1246  do 10 i = 3,4
1247  call area(a, f, lat, lon, 3*i, aa, pp)
1248  r = r + assert(aa, aa1*i, 0.5d0)
1249  call area(a, f, lat, lonr, 3*i, aa, pp)
1250  r = r + assert(aa, -aa1*i, 0.5d0)
1251  10 continue
1252 
1253  tstp21 = r
1254  return
1255  end
1256 
1257  program geodtest
1258  integer n, i
1259  integer tstinv, tstdir, tstarc,
1260  + tstg0, tstg1, tstg2, tstg5, tstg6, tstg9, tstg10, tstg11,
1261  + tstg12, tstg14, tstg15, tstg17, tstg26, tstg28, tstg33,
1262  + tstg55, tstg59, tstg61, tstg73, tstg74, tstg76, tstg78,
1263  + tstg80, tstg84,
1264  + tstp0, tstp5, tstp6, tstp12, tstp13, tstp15, tstp19, tstp21
1265 
1266  n = 0
1267  i = tstinv()
1268  if (i .gt. 0) then
1269  n = n + 1
1270  print *, 'tstinv fail:', i
1271  end if
1272  i = tstdir()
1273  if (i .gt. 0) then
1274  n = n + 1
1275  print *, 'tstdir fail:', i
1276  end if
1277  i = tstarc()
1278  if (i .gt. 0) then
1279  n = n + 1
1280  print *, 'tstarc fail:', i
1281  end if
1282  i = tstg0()
1283  if (i .gt. 0) then
1284  n = n + 1
1285  print *, 'tstg0 fail:', i
1286  end if
1287  i = tstg1()
1288  if (i .gt. 0) then
1289  n = n + 1
1290  print *, 'tstg1 fail:', i
1291  end if
1292  i = tstg2()
1293  if (i .gt. 0) then
1294  n = n + 1
1295  print *, 'tstg2 fail:', i
1296  end if
1297  i = tstg5()
1298  if (i .gt. 0) then
1299  n = n + 1
1300  print *, 'tstg5 fail:', i
1301  end if
1302  i = tstg6()
1303  if (i .gt. 0) then
1304  n = n + 1
1305  print *, 'tstg6 fail:', i
1306  end if
1307  i = tstg9()
1308  if (i .gt. 0) then
1309  n = n + 1
1310  print *, 'tstg9 fail:', i
1311  end if
1312  i = tstg10()
1313  if (i .gt. 0) then
1314  n = n + 1
1315  print *, 'tstg10 fail:', i
1316  end if
1317  i = tstg11()
1318  if (i .gt. 0) then
1319  n = n + 1
1320  print *, 'tstg11 fail:', i
1321  end if
1322  i = tstg12()
1323  if (i .gt. 0) then
1324  n = n + 1
1325  print *, 'tstg12 fail:', i
1326  end if
1327  i = tstg14()
1328  if (i .gt. 0) then
1329  n = n + 1
1330  print *, 'tstg14 fail:', i
1331  end if
1332  i = tstg15()
1333  if (i .gt. 0) then
1334  n = n + 1
1335  print *, 'tstg15 fail:', i
1336  end if
1337  i = tstg17()
1338  if (i .gt. 0) then
1339  n = n + 1
1340  print *, 'tstg17 fail:', i
1341  end if
1342  i = tstg26()
1343  if (i .gt. 0) then
1344  n = n + 1
1345  print *, 'tstg26 fail:', i
1346  end if
1347  i = tstg28()
1348  if (i .gt. 0) then
1349  n = n + 1
1350  print *, 'tstg28 fail:', i
1351  end if
1352  i = tstg33()
1353  if (i .gt. 0) then
1354  n = n + 1
1355  print *, 'tstg33 fail:', i
1356  end if
1357  i = tstg55()
1358  if (i .gt. 0) then
1359  n = n + 1
1360  print *, 'tstg55 fail:', i
1361  end if
1362  i = tstg59()
1363  if (i .gt. 0) then
1364  n = n + 1
1365  print *, 'tstg59 fail:', i
1366  end if
1367  i = tstg61()
1368  if (i .gt. 0) then
1369  n = n + 1
1370  print *, 'tstg61 fail:', i
1371  end if
1372  i = tstg73()
1373  if (i .gt. 0) then
1374  n = n + 1
1375  print *, 'tstg73 fail:', i
1376  end if
1377  i = tstg74()
1378  if (i .gt. 0) then
1379  n = n + 1
1380  print *, 'tstg74 fail:', i
1381  end if
1382  i = tstg76()
1383  if (i .gt. 0) then
1384  n = n + 1
1385  print *, 'tstg76 fail:', i
1386  end if
1387  i = tstg78()
1388  if (i .gt. 0) then
1389  n = n + 1
1390  print *, 'tstg78 fail:', i
1391  end if
1392  i = tstg80()
1393  if (i .gt. 0) then
1394  n = n + 1
1395  print *, 'tstg80 fail:', i
1396  end if
1397  i = tstg84()
1398  if (i .gt. 0) then
1399  n = n + 1
1400  print *, 'tstg84 fail:', i
1401  end if
1402  i = tstp0()
1403  if (i .gt. 0) then
1404  n = n + 1
1405  print *, 'tstp0 fail:', i
1406  end if
1407  i = tstp5()
1408  if (i .gt. 0) then
1409  n = n + 1
1410  print *, 'tstp5 fail:', i
1411  end if
1412  i = tstp6()
1413  if (i .gt. 0) then
1414  n = n + 1
1415  print *, 'tstp6 fail:', i
1416  end if
1417  i = tstp12()
1418  if (i .gt. 0) then
1419  n = n + 1
1420  print *, 'tstp12 fail:', i
1421  end if
1422  i = tstp13()
1423  if (i .gt. 0) then
1424  n = n + 1
1425  print *, 'tstp13 fail:', i
1426  end if
1427  i = tstp15()
1428  if (i .gt. 0) then
1429  n = n + 1
1430  print *, 'tstp15 fail:', i
1431  end if
1432  i = tstp19()
1433  if (i .gt. 0) then
1434  n = n + 1
1435  print *, 'tstp19 fail:', i
1436  end if
1437  i = tstp21()
1438  if (i .gt. 0) then
1439  n = n + 1
1440  print *, 'tstp21 fail:', i
1441  end if
1442 
1443  if (n .gt. 0) then
1444  stop 1
1445  end if
1446 
1447  stop
1448  end
1449 
1450 *> @endcond SKIP
direct
Definition: geodesic.inc:15
area
Definition: geodesic.inc:31
invers
Definition: geodesic.inc:23