Fortran library for Geodesics  1.49
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-2017) <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 tstinv()
171  double precision tstdat(20, 12)
172  common /tstcom/ tstdat
173  double precision lat1, lon1, azi1, lat2, lon2, azi2,
174  + s12, a12, m12, mm12, mm21, ss12
175  double precision azi1a, azi2a, s12a, a12a,
176  + m12a, mm12a, mm21a, ss12a
177  double precision a, f
178  integer r, assert, i, omask
179  include 'geodesic.inc'
180 
181 * WGS84 values
182  a = 6378137d0
183  f = 1/298.257223563d0
184  omask = 1 + 2 + 4 + 8
185  r = 0
186 
187  do i = 1,20
188  lat1 = tstdat(i, 1)
189  lon1 = tstdat(i, 2)
190  azi1 = tstdat(i, 3)
191  lat2 = tstdat(i, 4)
192  lon2 = tstdat(i, 5)
193  azi2 = tstdat(i, 6)
194  s12 = tstdat(i, 7)
195  a12 = tstdat(i, 8)
196  m12 = tstdat(i, 9)
197  mm12 = tstdat(i, 10)
198  mm21 = tstdat(i, 11)
199  ss12 = tstdat(i, 12)
200  call invers(a, f, lat1, lon1, lat2, lon2,
201  + s12a, azi1a, azi2a, omask, a12a, m12a, mm12a, mm21a, ss12a)
202  r = r + assert(azi1, azi1a, 1d-13)
203  r = r + assert(azi2, azi2a, 1d-13)
204  r = r + assert(s12, s12a, 1d-8)
205  r = r + assert(a12, a12a, 1d-13)
206  r = r + assert(m12, m12a, 1d-8)
207  r = r + assert(mm12, mm12a, 1d-15)
208  r = r + assert(mm21, mm21a, 1d-15)
209  r = r + assert(ss12, ss12a, 0.1d0)
210  end do
211 
212  tstinv = r
213  return
214  end
215 
216  integer function tstdir()
217  double precision tstdat(20, 12)
218  common /tstcom/ tstdat
219  double precision lat1, lon1, azi1, lat2, lon2, azi2,
220  + s12, a12, m12, mm12, mm21, ss12
221  double precision lat2a, lon2a, azi2a, a12a,
222  + m12a, mm12a, mm21a, ss12a
223  double precision a, f
224  integer r, assert, i, omask, flags
225  include 'geodesic.inc'
226 
227 * WGS84 values
228  a = 6378137d0
229  f = 1/298.257223563d0
230  omask = 1 + 2 + 4 + 8
231  flags = 2
232  r = 0
233 
234  do i = 1,20
235  lat1 = tstdat(i, 1)
236  lon1 = tstdat(i, 2)
237  azi1 = tstdat(i, 3)
238  lat2 = tstdat(i, 4)
239  lon2 = tstdat(i, 5)
240  azi2 = tstdat(i, 6)
241  s12 = tstdat(i, 7)
242  a12 = tstdat(i, 8)
243  m12 = tstdat(i, 9)
244  mm12 = tstdat(i, 10)
245  mm21 = tstdat(i, 11)
246  ss12 = tstdat(i, 12)
247  call direct(a, f, lat1, lon1, azi1, s12, flags,
248  + lat2a, lon2a, azi2a, omask, a12a, m12a, mm12a, mm21a, ss12a)
249  r = r + assert(lat2, lat2a, 1d-13)
250  r = r + assert(lon2, lon2a, 1d-13)
251  r = r + assert(azi2, azi2a, 1d-13)
252  r = r + assert(a12, a12a, 1d-13)
253  r = r + assert(m12, m12a, 1d-8)
254  r = r + assert(mm12, mm12a, 1d-15)
255  r = r + assert(mm21, mm21a, 1d-15)
256  r = r + assert(ss12, ss12a, 0.1d0)
257  end do
258 
259  tstdir = r
260  return
261  end
262 
263  integer function tstarc()
264  double precision tstdat(20, 12)
265  common /tstcom/ tstdat
266  double precision lat1, lon1, azi1, lat2, lon2, azi2,
267  + s12, a12, m12, mm12, mm21, ss12
268  double precision lat2a, lon2a, azi2a, s12a,
269  + m12a, mm12a, mm21a, ss12a
270  double precision a, f
271  integer r, assert, i, omask, flags
272  include 'geodesic.inc'
273 
274 * WGS84 values
275  a = 6378137d0
276  f = 1/298.257223563d0
277  omask = 1 + 2 + 4 + 8
278  flags = 1 + 2
279  r = 0
280 
281  do i = 1,20
282  lat1 = tstdat(i, 1)
283  lon1 = tstdat(i, 2)
284  azi1 = tstdat(i, 3)
285  lat2 = tstdat(i, 4)
286  lon2 = tstdat(i, 5)
287  azi2 = tstdat(i, 6)
288  s12 = tstdat(i, 7)
289  a12 = tstdat(i, 8)
290  m12 = tstdat(i, 9)
291  mm12 = tstdat(i, 10)
292  mm21 = tstdat(i, 11)
293  ss12 = tstdat(i, 12)
294  call direct(a, f, lat1, lon1, azi1, a12, flags,
295  + lat2a, lon2a, azi2a, omask, s12a, m12a, mm12a, mm21a, ss12a)
296  r = r + assert(lat2, lat2a, 1d-13)
297  r = r + assert(lon2, lon2a, 1d-13)
298  r = r + assert(azi2, azi2a, 1d-13)
299  r = r + assert(s12, s12a, 1d-8)
300  r = r + assert(m12, m12a, 1d-8)
301  r = r + assert(mm12, mm12a, 1d-15)
302  r = r + assert(mm21, mm21a, 1d-15)
303  r = r + assert(ss12, ss12a, 0.1d0)
304  end do
305 
306  tstarc = r
307  return
308  end
309 
310  integer function notnan(x)
311  double precision x
312  if (x .eq. x) then
313  notnan = 1
314  else
315  notnan = 0
316  end if
317 
318  return
319  end
320 
321  integer function tstg0()
322  double precision azi1, azi2, s12, a12, m12, mm12, mm21, ss12
323  double precision a, f
324  integer r, assert, omask
325  include 'geodesic.inc'
326 
327 * WGS84 values
328  a = 6378137d0
329  f = 1/298.257223563d0
330  omask = 0
331  r = 0
332  call invers(a, f, 40.6d0, -73.8d0, 49.01666667d0, 2.55d0,
333  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
334  r = r + assert(azi1, 53.47022d0, 0.5d-5)
335  r = r + assert(azi2, 111.59367d0, 0.5d-5)
336  r = r + assert(s12, 5853226d0, 0.5d0)
337 
338  tstg0 = r
339  return
340  end
341 
342  integer function tstg1()
343  double precision lat2, lon2, azi2, a12, m12, mm12, mm21, ss12
344  double precision a, f
345  integer r, assert, omask, flags
346  include 'geodesic.inc'
347 
348 * WGS84 values
349  a = 6378137d0
350  f = 1/298.257223563d0
351  omask = 0
352  flags = 0
353  r = 0
354  call direct(a, f, 40.63972222d0, -73.77888889d0, 53.5d0, 5850d3,
355  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
356  r = r + assert(lat2, 49.01467d0, 0.5d-5)
357  r = r + assert(lon2, 2.56106d0, 0.5d-5)
358  r = r + assert(azi2, 111.62947d0, 0.5d-5)
359 
360  tstg1 = r
361  return
362  end
363 
364  integer function tstg2()
365 * Check fix for antipodal prolate bug found 2010-09-04
366  double precision azi1, azi2, s12, a12, m12, mm12, mm21, ss12
367  double precision a, f
368  integer r, assert, omask
369  include 'geodesic.inc'
370 
371  a = 6.4d6
372  f = -1/150d0
373  omask = 0
374  r = 0
375  call invers(a, f, 0.07476d0, 0d0, -0.07476d0, 180d0,
376  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
377  r = r + assert(azi1, 90.00078d0, 0.5d-5)
378  r = r + assert(azi2, 90.00078d0, 0.5d-5)
379  r = r + assert(s12, 20106193d0, 0.5d0)
380  call invers(a, f, 0.1d0, 0d0, -0.1d0, 180d0,
381  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
382  r = r + assert(azi1, 90.00105d0, 0.5d-5)
383  r = r + assert(azi2, 90.00105d0, 0.5d-5)
384  r = r + assert(s12, 20106193d0, 0.5d0)
385 
386  tstg2 = r
387  return
388  end
389 
390  integer function tstg4()
391 * Check fix for short line bug found 2010-05-21
392  double precision azi1, azi2, s12, a12, m12, mm12, mm21, ss12
393  double precision a, f
394  integer r, assert, omask
395  include 'geodesic.inc'
396 
397 * WGS84 values
398  a = 6378137d0
399  f = 1/298.257223563d0
400  omask = 0
401  r = 0
402  call invers(a, f,
403  + 36.493349428792d0, 0d0, 36.49334942879201d0, .0000008d0,
404  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
405  r = r + assert(s12, 0.072d0, 0.5d-3)
406 
407  tstg4 = r
408  return
409  end
410 
411  integer function tstg5()
412  double precision lat2, lon2, azi2, a12, m12, mm12, mm21, ss12
413  double precision a, f
414  integer r, assert, omask, flags
415  include 'geodesic.inc'
416 
417 * WGS84 values
418  a = 6378137d0
419  f = 1/298.257223563d0
420  omask = 0
421  flags = 0
422  r = 0
423  call direct(a, f, 0.01777745589997d0, 30d0, 0d0, 10d6,
424  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
425  if (lon2 .lt. 0) then
426  r = r + assert(lon2, -150d0, 0.5d-5)
427  r = r + assert(abs(azi2), 180d0, 0.5d-5)
428  else
429  r = r + assert(lon2, 30d0, 0.5d-5)
430  r = r + assert(azi2, 0d0, 0.5d-5)
431  end if
432 
433  tstg5 = r
434  return
435  end
436 
437  integer function tstg6()
438 * Check fix for volatile sbet12a bug found 2011-06-25 (gcc 4.4d0.4d0
439 * x86 -O3). Found again on 2012-03-27 with tdm-mingw32 (g++ 4.6d0.1d0).
440  double precision azi1, azi2, s12, a12, m12, mm12, mm21, ss12
441  double precision a, f
442  integer r, assert, omask
443  include 'geodesic.inc'
444 
445 * WGS84 values
446  a = 6378137d0
447  f = 1/298.257223563d0
448  omask = 0
449  r = 0
450  call invers(a, f, 88.202499451857d0, 0d0,
451  + -88.202499451857d0, 179.981022032992859592d0,
452  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
453  r = r + assert(s12, 20003898.214d0, 0.5d-3)
454  call invers(a, f, 89.262080389218d0, 0d0,
455  + -89.262080389218d0, 179.992207982775375662d0,
456  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
457  r = r + assert(s12, 20003925.854d0, 0.5d-3)
458  call invers(a, f, 89.333123580033d0, 0d0,
459  + -89.333123580032997687d0, 179.99295812360148422d0,
460  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
461  r = r + assert(s12, 20003926.881d0, 0.5d-3)
462 
463  tstg6 = r
464  return
465  end
466 
467  integer function tstg9()
468 * Check fix for volatile x bug found 2011-06-25 (gcc 4.4d0.4d0 x86 -O3)
469  double precision azi1, azi2, s12, a12, m12, mm12, mm21, ss12
470  double precision a, f
471  integer r, assert, omask
472  include 'geodesic.inc'
473 
474 * WGS84 values
475  a = 6378137d0
476  f = 1/298.257223563d0
477  omask = 0
478  r = 0
479  call invers(a, f, 56.320923501171d0, 0d0,
480  + -56.320923501171d0, 179.664747671772880215d0,
481  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
482  r = r + assert(s12, 19993558.287d0, 0.5d-3)
483 
484  tstg9 = r
485  return
486  end
487 
488  integer function tstg10()
489 * Check fix for adjust tol1_ bug found 2011-06-25 (Visual Studio 10 rel
490 * + debug)
491  double precision azi1, azi2, s12, a12, m12, mm12, mm21, ss12
492  double precision a, f
493  integer r, assert, omask
494  include 'geodesic.inc'
495 
496 * WGS84 values
497  a = 6378137d0
498  f = 1/298.257223563d0
499  omask = 0
500  r = 0
501  call invers(a, f, 52.784459512564d0, 0d0,
502  + -52.784459512563990912d0, 179.634407464943777557d0,
503  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
504  r = r + assert(s12, 19991596.095d0, 0.5d-3)
505 
506  tstg10 = r
507  return
508  end
509 
510  integer function tstg11()
511 * Check fix for bet2 = -bet1 bug found 2011-06-25 (Visual Studio 10 rel
512 * + debug)
513  double precision azi1, azi2, s12, a12, m12, mm12, mm21, ss12
514  double precision a, f
515  integer r, assert, omask
516  include 'geodesic.inc'
517 
518 * WGS84 values
519  a = 6378137d0
520  f = 1/298.257223563d0
521  omask = 0
522  r = 0
523  call invers(a, f, 48.522876735459d0, 0d0,
524  + -48.52287673545898293d0, 179.599720456223079643d0,
525  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
526  r = r + assert(s12, 19989144.774d0, 0.5d-3)
527 
528  tstg11 = r
529  return
530  end
531 
532  integer function tstg12()
533 * Check fix for inverse geodesics on extreme prolate/oblate ellipsoids
534 * Reported 2012-08-29 Stefan Guenther <stefan.gunther@embl.de>; fixed
535 * 2012-10-07
536  double precision azi1, azi2, s12, a12, m12, mm12, mm21, ss12
537  double precision a, f
538  integer r, assert, omask
539  include 'geodesic.inc'
540 
541  a = 89.8d0
542  f = -1.83d0
543  omask = 0
544  r = 0
545  call invers(a, f, 0d0, 0d0, -10d0, 160d0,
546  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
547  r = r + assert(azi1, 120.27d0, 1d-2)
548  r = r + assert(azi2, 105.15d0, 1d-2)
549  r = r + assert(s12, 266.7d0, 1d-1)
550 
551  tstg12 = r
552  return
553  end
554 
555  integer function tstg14()
556 * Check fix for inverse ignoring lon12 = nan
557  double precision azi1, azi2, s12, a12, m12, mm12, mm21, ss12
558  double precision a, f, latfix
559  integer r, notnan, omask
560  include 'geodesic.inc'
561 
562 * WGS84 values
563  a = 6378137d0
564  f = 1/298.257223563d0
565  omask = 0
566  r = 0
567  call invers(a, f, 0d0, 0d0, 1d0, latfix(91d0),
568  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
569  r = r + notnan(azi1)
570  r = r + notnan(azi2)
571  r = r + notnan(s12)
572  tstg14 = r
573  return
574  end
575 
576  integer function tstg15()
577 * Initial implementation of Math::eatanhe was wrong for e^2 < 0. This
578 * checks that this is fixed.
579  double precision lat2, lon2, azi2, a12, m12, mm12, mm21, ss12
580  double precision a, f
581  integer r, assert, omask, flags
582  include 'geodesic.inc'
583 
584  a = 6.4d6
585  f = -1/150.0d0
586  omask = 8
587  flags = 0
588  r = 0
589  call direct(a, f, 1d0, 2d0, 3d0, 4d0,
590  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
591  r = r + assert(ss12, 23700d0, 0.5d0)
592 
593  tstg15 = r
594  return
595  end
596 
597  integer function tstg17()
598 * Check fix for LONG_UNROLL bug found on 2015-05-07
599  double precision lat2, lon2, azi2, a12, m12, mm12, mm21, ss12
600  double precision a, f
601  integer r, assert, omask, flags
602  include 'geodesic.inc'
603 
604 * WGS84 values
605  a = 6378137d0
606  f = 1/298.257223563d0
607  omask = 0
608  flags = 2
609  r = 0
610  call direct(a, f, 40d0, -75d0, -10d0, 2d7,
611  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
612  r = r + assert(lat2, -39d0, 1d0)
613  r = r + assert(lon2, -254d0, 1d0)
614  r = r + assert(azi2, -170d0, 1d0)
615  flags = 0
616  call direct(a, f, 40d0, -75d0, -10d0, 2d7,
617  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
618  r = r + assert(lat2, -39d0, 1d0)
619  r = r + assert(lon2, 105d0, 1d0)
620  r = r + assert(azi2, -170d0, 1d0)
621 
622  tstg17 = r
623  return
624  end
625 
626  integer function tstg26()
627 * Check 0/0 problem with area calculation on sphere 2015-09-08
628  double precision azi1, azi2, s12, a12, m12, mm12, mm21, ss12
629  double precision a, f
630  integer r, assert, omask
631  include 'geodesic.inc'
632 
633  a = 6.4d6
634  f = 0
635  omask = 8
636  r = 0
637  call invers(a, f, 1d0, 2d0, 3d0, 4d0,
638  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
639  r = r + assert(ss12, 49911046115.0d0, 0.5d0)
640 
641  tstg26 = r
642  return
643  end
644 
645  integer function tstg28()
646 * Check fix for LONG_UNROLL bug found on 2015-05-07
647  double precision lat2, lon2, azi2, a12, m12, mm12, mm21, ss12
648  double precision a, f
649  integer r, assert, omask, flags
650  include 'geodesic.inc'
651 
652  a = 6.4d6
653  f = 0.1d0
654  omask = 1 + 2 + 4 + 8
655  flags = 0
656  r = 0
657  call direct(a, f, 1d0, 2d0, 10d0, 5d6,
658  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
659  r = r + assert(a12, 48.55570690d0, 0.5d-8)
660 
661  tstg28 = r
662  return
663  end
664 
665  integer function tstg33()
666 * Check max(-0.0,+0.0) issues 2015-08-22 (triggered by bugs in Octave --
667 * sind(-0.0) = +0.0 -- and in some version of Visual Studio --
668 * fmod(-0.0, 360.0) = +0.0.
669  double precision azi1, azi2, s12, a12, m12, mm12, mm21, ss12
670  double precision a, f
671  integer r, assert, omask
672  include 'geodesic.inc'
673 
674 * WGS84 values
675  a = 6378137d0
676  f = 1/298.257223563d0
677  omask = 0
678  r = 0
679  call invers(a, f, 0d0, 0d0, 0d0, 179d0,
680  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
681  r = r + assert(azi1, 90.00000d0, 0.5d-5)
682  r = r + assert(azi2, 90.00000d0, 0.5d-5)
683  r = r + assert(s12, 19926189d0, 0.5d0)
684  call invers(a, f, 0d0, 0d0, 0d0, 179.5d0,
685  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
686  r = r + assert(azi1, 55.96650d0, 0.5d-5)
687  r = r + assert(azi2, 124.03350d0, 0.5d-5)
688  r = r + assert(s12, 19980862d0, 0.5d0)
689  call invers(a, f, 0d0, 0d0, 0d0, 180d0,
690  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
691  r = r + assert(azi1, 0.00000d0, 0.5d-5)
692  r = r + assert(abs(azi2), 180.00000d0, 0.5d-5)
693  r = r + assert(s12, 20003931d0, 0.5d0)
694  call invers(a, f, 0d0, 0d0, 1d0, 180d0,
695  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
696  r = r + assert(azi1, 0.00000d0, 0.5d-5)
697  r = r + assert(abs(azi2), 180.00000d0, 0.5d-5)
698  r = r + assert(s12, 19893357d0, 0.5d0)
699  a = 6.4d6
700  f = 0
701  call invers(a, f, 0d0, 0d0, 0d0, 179d0,
702  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
703  r = r + assert(azi1, 90.00000d0, 0.5d-5)
704  r = r + assert(azi2, 90.00000d0, 0.5d-5)
705  r = r + assert(s12, 19994492d0, 0.5d0)
706  call invers(a, f, 0d0, 0d0, 0d0, 180d0,
707  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
708  r = r + assert(azi1, 0.00000d0, 0.5d-5)
709  r = r + assert(abs(azi2), 180.00000d0, 0.5d-5)
710  r = r + assert(s12, 20106193d0, 0.5d0)
711  call invers(a, f, 0d0, 0d0, 1d0, 180d0,
712  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
713  r = r + assert(azi1, 0.00000d0, 0.5d-5)
714  r = r + assert(abs(azi2), 180.00000d0, 0.5d-5)
715  r = r + assert(s12, 19994492d0, 0.5d0)
716  a = 6.4d6
717  f = -1/300.0d0
718  call invers(a, f, 0d0, 0d0, 0d0, 179d0,
719  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
720  r = r + assert(azi1, 90.00000d0, 0.5d-5)
721  r = r + assert(azi2, 90.00000d0, 0.5d-5)
722  r = r + assert(s12, 19994492d0, 0.5d0)
723  call invers(a, f, 0d0, 0d0, 0d0, 180d0,
724  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
725  r = r + assert(azi1, 90.00000d0, 0.5d-5)
726  r = r + assert(azi2, 90.00000d0, 0.5d-5)
727  r = r + assert(s12, 20106193d0, 0.5d0)
728  call invers(a, f, 0d0, 0d0, 0.5d0, 180d0,
729  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
730  r = r + assert(azi1, 33.02493d0, 0.5d-5)
731  r = r + assert(azi2, 146.97364d0, 0.5d-5)
732  r = r + assert(s12, 20082617d0, 0.5d0)
733  call invers(a, f, 0d0, 0d0, 1d0, 180d0,
734  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
735  r = r + assert(azi1, 0.00000d0, 0.5d-5)
736  r = r + assert(abs(azi2), 180.00000d0, 0.5d-5)
737  r = r + assert(s12, 20027270d0, 0.5d0)
738 
739  tstg33 = r
740  return
741  end
742 
743  integer function tstg55()
744 * Check fix for nan + point on equator or pole not returning all nans in
745 * Geodesic::Inverse, found 2015-09-23.
746  double precision azi1, azi2, s12, a12, m12, mm12, mm21, ss12
747  double precision a, f
748  integer r, notnan, omask
749  include 'geodesic.inc'
750 
751 * WGS84 values
752  a = 6378137d0
753  f = 1/298.257223563d0
754  omask = 0
755  r = 0
756  call invers(a, f, 91d0, 0d0, 0d0, 90d0,
757  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
758  r = r + notnan(azi1)
759  r = r + notnan(azi2)
760  r = r + notnan(s12)
761  call invers(a, f, 91d0, 0d0, 90d0, 9d0,
762  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
763  r = r + notnan(azi1)
764  r = r + notnan(azi2)
765  r = r + notnan(s12)
766  tstg55 = r
767  return
768  end
769 
770  integer function tstg59()
771 * Check for points close with longitudes close to 180 deg apart.
772  double precision azi1, azi2, s12, a12, m12, mm12, mm21, ss12
773  double precision a, f
774  integer r, assert, omask
775  include 'geodesic.inc'
776 
777 * WGS84 values
778  a = 6378137d0
779  f = 1/298.257223563d0
780  omask = 0
781  r = 0
782  call invers(a, f, 5d0, 0.00000000000001d0, 10d0, 180d0,
783  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
784  r = r + assert(azi1, 0.000000000000035d0, 1.5d-14);
785  r = r + assert(azi2, 179.99999999999996d0, 1.5d-14);
786  r = r + assert(s12, 18345191.174332713d0, 2.5d-9);
787  tstg59 = r
788  return
789  end
790 
791  integer function tstg61()
792 * Make sure small negative azimuths are west-going
793  double precision lat2, lon2, azi2, a12, m12, mm12, mm21, ss12
794  double precision a, f
795  integer r, assert, omask, flags
796  include 'geodesic.inc'
797 
798 * WGS84 values
799  a = 6378137d0
800  f = 1/298.257223563d0
801  omask = 0
802  flags = 2
803  r = 0
804  call direct(a, f, 45d0, 0d0, -0.000000000000000003d0, 1d7,
805  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
806  r = r + assert(lat2, 45.30632d0, 0.5d-5)
807  r = r + assert(lon2, -180d0, 0.5d-5)
808  r = r + assert(abs(azi2), 180d0, 0.5d-5)
809 
810  tstg61 = r
811  return
812  end
813 
814  integer function tstg73()
815 * Check for backwards from the pole bug reported by Anon on 2016-02-13.
816 * This only affected the Java implementation. It was introduced in Java
817 * version 1.44 and fixed in 1.46-SNAPSHOT on 2016-01-17.
818  double precision lat2, lon2, azi2, a12, m12, mm12, mm21, ss12
819  double precision a, f
820  integer r, assert, omask, flags
821  include 'geodesic.inc'
822 
823 * WGS84 values
824  a = 6378137d0
825  f = 1/298.257223563d0
826  omask = 0
827  flags = 0
828  r = 0
829  call direct(a, f, 90d0, 10d0, 180d0, -1d6,
830  + flags, lat2, lon2, azi2, omask, a12, m12, mm12, mm21, ss12)
831  r = r + assert(lat2, 81.04623d0, 0.5d-5)
832  r = r + assert(lon2, -170d0, 0.5d-5)
833  r = r + assert(azi2, 0d0, 0.5d-5)
834 
835  tstg73 = r
836  return
837  end
838 
839  integer function tstg74()
840 * Check fix for inaccurate areas, bug introduced in v1.46, fixed
841 * 2015-10-16.
842  double precision azi1, azi2, s12, a12, m12, mm12, mm21, ss12
843  double precision a, f
844  integer r, assert, omask
845  include 'geodesic.inc'
846 
847 * WGS84 values
848  a = 6378137d0
849  f = 1/298.257223563d0
850  omask = 1 + 2 + 4 + 8
851  r = 0
852  call invers(a, f, 54.1589d0, 15.3872d0, 54.1591d0, 15.3877d0,
853  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
854  r = r + assert(azi1, 55.723110355d0, 5d-9);
855  r = r + assert(azi2, 55.723515675d0, 5d-9);
856  r = r + assert(s12, 39.527686385d0, 5d-9);
857  r = r + assert(a12, 0.000355495d0, 5d-9);
858  r = r + assert(m12, 39.527686385d0, 5d-9);
859  r = r + assert(mm12, 0.999999995d0, 5d-9);
860  r = r + assert(mm21, 0.999999995d0, 5d-9);
861  r = r + assert(ss12, 286698586.30197d0, 5d-4);
862  tstg74 = r
863  return
864  end
865 
866  integer function tstg76()
867 * The distance from Wellington and Salamanca (a classic failure of
868 * Vincenty
869  double precision azi1, azi2, s12, a12, m12, mm12, mm21, ss12
870  double precision a, f
871  integer r, assert, omask
872  include 'geodesic.inc'
873 
874 * WGS84 values
875  a = 6378137d0
876  f = 1/298.257223563d0
877  omask = 0
878  r = 0
879  call invers(a, f,
880  + -(41+19/60d0), 174+49/60d0, 40+58/60d0, -(5+30/60d0),
881  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
882  r = r + assert(azi1, 160.39137649664d0, 0.5d-11)
883  r = r + assert(azi2, 19.50042925176d0, 0.5d-11)
884  r = r + assert(s12, 19960543.857179d0, 0.5d-6)
885  tstg76 = r
886  return
887  end
888 
889  integer function tstg78()
890 * An example where the NGS calculator fails to converge
891  double precision azi1, azi2, s12, a12, m12, mm12, mm21, ss12
892  double precision a, f
893  integer r, assert, omask
894  include 'geodesic.inc'
895 
896 * WGS84 values
897  a = 6378137d0
898  f = 1/298.257223563d0
899  omask = 0
900  r = 0
901  call invers(a, f, 27.2d0, 0d0, -27.1d0, 179.5d0,
902  + s12, azi1, azi2, omask, a12, m12, mm12, mm21, ss12)
903  r = r + assert(azi1, 45.82468716758d0, 0.5d-11)
904  r = r + assert(azi2, 134.22776532670d0, 0.5d-11)
905  r = r + assert(s12, 19974354.765767d0, 0.5d-6)
906  tstg78 = r
907  return
908  end
909 
910  integer function tstp0()
911 * Check fix for pole-encircling bug found 2011-03-16
912  double precision lata(4), lona(4)
913  data lata / 89d0, 89d0, 89d0, 89d0 /
914  data lona / 0d0, 90d0, 180d0, 270d0 /
915  double precision latb(4), lonb(4)
916  data latb / -89d0, -89d0, -89d0, -89d0 /
917  data lonb / 0d0, 90d0, 180d0, 270d0 /
918  double precision latc(4), lonc(4)
919  data latc / 0d0, -1d0, 0d0, 1d0 /
920  data lonc / -1d0, 0d0, 1d0, 0d0 /
921  double precision latd(3), lond(3)
922  data latd / 90d0, 0d0, 0d0 /
923  data lond / 0d0, 0d0, 90d0 /
924  double precision a, f, aa, pp
925  integer r, assert
926  include 'geodesic.inc'
927 
928 * WGS84 values
929  a = 6378137d0
930  f = 1/298.257223563d0
931  r = 0
932 
933  call area(a, f, lata, lona, 4, aa, pp)
934  r = r + assert(pp, 631819.8745d0, 1d-4)
935  r = r + assert(aa, 24952305678.0d0, 1d0)
936 
937  call area(a, f, latb, lonb, 4, aa, pp)
938  r = r + assert(pp, 631819.8745d0, 1d-4)
939  r = r + assert(aa, -24952305678.0d0, 1d0)
940 
941  call area(a, f, latc, lonc, 4, aa, pp)
942  r = r + assert(pp, 627598.2731d0, 1d-4)
943  r = r + assert(aa, 24619419146.0d0, 1d0)
944 
945  call area(a, f, latd, lond, 3, aa, pp)
946  r = r + assert(pp, 30022685d0, 1d0)
947  r = r + assert(aa, 63758202715511.0d0, 1d0)
948 
949  tstp0 = r
950  return
951  end
952 
953  integer function tstp5()
954 * Check fix for Planimeter pole crossing bug found 2011-06-24
955  double precision lat(3), lon(3)
956  data lat / 89d0, 89d0, 89d0 /
957  data lon / 0.1d0, 90.1d0, -179.9d0 /
958  double precision a, f, aa, pp
959  integer r, assert
960  include 'geodesic.inc'
961 
962 * WGS84 values
963  a = 6378137d0
964  f = 1/298.257223563d0
965  r = 0
966 
967  call area(a, f, lat, lon, 3, aa, pp)
968  r = r + assert(pp, 539297d0, 1d0)
969  r = r + assert(aa, 12476152838.5d0, 1d0)
970 
971  tstp5 = r
972  return
973  end
974 
975  integer function tstp6()
976 * Check fix for pole-encircling bug found 2011-03-16
977  double precision lata(3), lona(3)
978  data lata / 9d0, 9d0, 9d0 /
979  data lona / -0.00000000000001d0, 180d0, 0d0 /
980  double precision latb(3), lonb(3)
981  data latb / 9d0, 9d0, 9d0 /
982  data lonb / 0.00000000000001d0, 0d0, 180d0 /
983  double precision latc(3), lonc(3)
984  data latc / 9d0, 9d0, 9d0 /
985  data lonc / 0.00000000000001d0, 180d0, 0d0 /
986  double precision latd(3), lond(3)
987  data latd / 9d0, 9d0, 9d0 /
988  data lond / -0.00000000000001d0, 0d0, 180d0 /
989  double precision a, f, aa, pp
990  integer r, assert
991  include 'geodesic.inc'
992 
993 * WGS84 values
994  a = 6378137d0
995  f = 1/298.257223563d0
996  r = 0
997 
998  call area(a, f, lata, lona, 3, aa, pp)
999  r = r + assert(pp, 36026861d0, 1d0)
1000  r = r + assert(aa, 0d0, 1d0)
1001 
1002  tstp6 = r
1003  return
1004  end
1005 
1006  integer function tstp12()
1007 * AA of arctic circle (not really -- adjunct to rhumb-AA test)
1008  double precision lat(2), lon(2)
1009  data lat / 66.562222222d0, 66.562222222d0 /
1010  data lon / 0d0, 180d0 /
1011  double precision a, f, aa, pp
1012  integer r, assert
1013  include 'geodesic.inc'
1014 
1015 * WGS84 values
1016  a = 6378137d0
1017  f = 1/298.257223563d0
1018  r = 0
1019 
1020  call area(a, f, lat, lon, 2, aa, pp)
1021  r = r + assert(pp, 10465729d0, 1d0)
1022  r = r + assert(aa, 0d0, 1d0)
1023 
1024  tstp12 = r
1025  return
1026  end
1027 
1028  integer function tstp13()
1029 * Check encircling pole twice
1030  double precision lat(6), lon(6)
1031  data lat / 89d0, 89d0, 89d0, 89d0, 89d0, 89d0 /
1032  data lon / -360d0, -240d0, -120d0, 0d0, 120d0, 240d0 /
1033  double precision a, f, aa, pp
1034  integer r, assert
1035  include 'geodesic.inc'
1036 
1037 * WGS84 values
1038  a = 6378137d0
1039  f = 1/298.257223563d0
1040  r = 0
1041 
1042  call area(a, f, lat, lon, 6, aa, pp)
1043  r = r + assert(pp, 1160741d0, 1d0)
1044  r = r + assert(aa, 32415230256.0d0, 1d0)
1045 
1046  tstp13 = r
1047  return
1048  end
1049 
1050  program geodtest
1051  integer n, i
1052  integer tstinv, tstdir, tstarc,
1053  + tstg0, tstg1, tstg2, tstg5, tstg6, tstg9, tstg10, tstg11,
1054  + tstg12, tstg14, tstg15, tstg17, tstg26, tstg28, tstg33,
1055  + tstg55, tstg59, tstg61, tstg73, tstg74, tstg76, tstg78,
1056  + tstp0, tstp5, tstp6, tstp12, tstp13
1057 
1058  n = 0
1059  i = tstinv()
1060  if (i .gt. 0) then
1061  n = n + 1
1062  print *, 'tstinv fail:', i
1063  end if
1064  i = tstdir()
1065  if (i .gt. 0) then
1066  n = n + 1
1067  print *, 'tstdir fail:', i
1068  end if
1069  i = tstarc()
1070  if (i .gt. 0) then
1071  n = n + 1
1072  print *, 'tstarc fail:', i
1073  end if
1074  i = tstg0()
1075  if (i .gt. 0) then
1076  n = n + 1
1077  print *, 'tstg0 fail:', i
1078  end if
1079  i = tstg1()
1080  if (i .gt. 0) then
1081  n = n + 1
1082  print *, 'tstg1 fail:', i
1083  end if
1084  i = tstg2()
1085  if (i .gt. 0) then
1086  n = n + 1
1087  print *, 'tstg2 fail:', i
1088  end if
1089  i = tstg5()
1090  if (i .gt. 0) then
1091  n = n + 1
1092  print *, 'tstg5 fail:', i
1093  end if
1094  i = tstg6()
1095  if (i .gt. 0) then
1096  n = n + 1
1097  print *, 'tstg6 fail:', i
1098  end if
1099  i = tstg9()
1100  if (i .gt. 0) then
1101  n = n + 1
1102  print *, 'tstg9 fail:', i
1103  end if
1104  i = tstg10()
1105  if (i .gt. 0) then
1106  n = n + 1
1107  print *, 'tstg10 fail:', i
1108  end if
1109  i = tstg11()
1110  if (i .gt. 0) then
1111  n = n + 1
1112  print *, 'tstg11 fail:', i
1113  end if
1114  i = tstg12()
1115  if (i .gt. 0) then
1116  n = n + 1
1117  print *, 'tstg12 fail:', i
1118  end if
1119  i = tstg14()
1120  if (i .gt. 0) then
1121  n = n + 1
1122  print *, 'tstg14 fail:', i
1123  end if
1124  i = tstg15()
1125  if (i .gt. 0) then
1126  n = n + 1
1127  print *, 'tstg15 fail:', i
1128  end if
1129  i = tstg17()
1130  if (i .gt. 0) then
1131  n = n + 1
1132  print *, 'tstg17 fail:', i
1133  end if
1134  i = tstg26()
1135  if (i .gt. 0) then
1136  n = n + 1
1137  print *, 'tstg26 fail:', i
1138  end if
1139  i = tstg28()
1140  if (i .gt. 0) then
1141  n = n + 1
1142  print *, 'tstg28 fail:', i
1143  end if
1144  i = tstg33()
1145  if (i .gt. 0) then
1146  n = n + 1
1147  print *, 'tstg33 fail:', i
1148  end if
1149  i = tstg55()
1150  if (i .gt. 0) then
1151  n = n + 1
1152  print *, 'tstg55 fail:', i
1153  end if
1154  i = tstg59()
1155  if (i .gt. 0) then
1156  n = n + 1
1157  print *, 'tstg59 fail:', i
1158  end if
1159  i = tstg61()
1160  if (i .gt. 0) then
1161  n = n + 1
1162  print *, 'tstg61 fail:', i
1163  end if
1164  i = tstg73()
1165  if (i .gt. 0) then
1166  n = n + 1
1167  print *, 'tstg73 fail:', i
1168  end if
1169  i = tstg74()
1170  if (i .gt. 0) then
1171  n = n + 1
1172  print *, 'tstg74 fail:', i
1173  end if
1174  i = tstg76()
1175  if (i .gt. 0) then
1176  n = n + 1
1177  print *, 'tstg76 fail:', i
1178  end if
1179  i = tstg78()
1180  if (i .gt. 0) then
1181  n = n + 1
1182  print *, 'tstg78 fail:', i
1183  end if
1184  i = tstp0()
1185  if (i .gt. 0) then
1186  n = n + 1
1187  print *, 'tstp0 fail:', i
1188  end if
1189  i = tstp5()
1190  if (i .gt. 0) then
1191  n = n + 1
1192  print *, 'tstp5 fail:', i
1193  end if
1194  i = tstp6()
1195  if (i .gt. 0) then
1196  n = n + 1
1197  print *, 'tstp6 fail:', i
1198  end if
1199  i = tstp12()
1200  if (i .gt. 0) then
1201  n = n + 1
1202  print *, 'tstp12 fail:', i
1203  end if
1204  i = tstp13()
1205  if (i .gt. 0) then
1206  n = n + 1
1207  print *, 'tstp13 fail:', i
1208  end if
1209 
1210  if (n .gt. 0) then
1211  stop 1
1212  end if
1213 
1214  stop
1215  end
1216 
1217 *> @endcond SKIP