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 tstp0() 1062* Check fix for pole-encircling bug found 2011-03-16 1063 double precision lata(4), lona(4) 1064 data lata / 89d0, 89d0, 89d0, 89d0 / 1065 data lona / 0d0, 90d0, 180d0, 270d0 / 1066 double precision latb(4), lonb(4) 1067 data latb / -89d0, -89d0, -89d0, -89d0 / 1068 data lonb / 0d0, 90d0, 180d0, 270d0 / 1069 double precision latc(4), lonc(4) 1070 data latc / 0d0, -1d0, 0d0, 1d0 / 1071 data lonc / -1d0, 0d0, 1d0, 0d0 / 1072 double precision latd(3), lond(3) 1073 data latd / 90d0, 0d0, 0d0 / 1074 data lond / 0d0, 0d0, 90d0 / 1075 double precision a, f, AA, PP 1076 integer r, assert 1077 include 'geodesic.inc' 1078 1079* WGS84 values 1080 a = 6378137d0 1081 f = 1/298.257223563d0 1082 r = 0 1083 1084 call area(a, f, lata, lona, 4, AA, PP) 1085 r = r + assert(PP, 631819.8745d0, 1d-4) 1086 r = r + assert(AA, 24952305678.0d0, 1d0) 1087 1088 call area(a, f, latb, lonb, 4, AA, PP) 1089 r = r + assert(PP, 631819.8745d0, 1d-4) 1090 r = r + assert(AA, -24952305678.0d0, 1d0) 1091 1092 call area(a, f, latc, lonc, 4, AA, PP) 1093 r = r + assert(PP, 627598.2731d0, 1d-4) 1094 r = r + assert(AA, 24619419146.0d0, 1d0) 1095 1096 call area(a, f, latd, lond, 3, AA, PP) 1097 r = r + assert(PP, 30022685d0, 1d0) 1098 r = r + assert(AA, 63758202715511.0d0, 1d0) 1099 1100 tstp0 = r 1101 return 1102 end 1103 1104 integer function tstp5() 1105* Check fix for Planimeter pole crossing bug found 2011-06-24 1106 double precision lat(3), lon(3) 1107 data lat / 89d0, 89d0, 89d0 / 1108 data lon / 0.1d0, 90.1d0, -179.9d0 / 1109 double precision a, f, AA, PP 1110 integer r, assert 1111 include 'geodesic.inc' 1112 1113* WGS84 values 1114 a = 6378137d0 1115 f = 1/298.257223563d0 1116 r = 0 1117 1118 call area(a, f, lat, lon, 3, AA, PP) 1119 r = r + assert(PP, 539297d0, 1d0) 1120 r = r + assert(AA, 12476152838.5d0, 1d0) 1121 1122 tstp5 = r 1123 return 1124 end 1125 1126 integer function tstp6() 1127* Check fix for pole-encircling bug found 2011-03-16 1128 double precision lata(3), lona(3) 1129 data lata / 9d0, 9d0, 9d0 / 1130 data lona / -0.00000000000001d0, 180d0, 0d0 / 1131 double precision latb(3), lonb(3) 1132 data latb / 9d0, 9d0, 9d0 / 1133 data lonb / 0.00000000000001d0, 0d0, 180d0 / 1134 double precision latc(3), lonc(3) 1135 data latc / 9d0, 9d0, 9d0 / 1136 data lonc / 0.00000000000001d0, 180d0, 0d0 / 1137 double precision latd(3), lond(3) 1138 data latd / 9d0, 9d0, 9d0 / 1139 data lond / -0.00000000000001d0, 0d0, 180d0 / 1140 double precision a, f, AA, PP 1141 integer r, assert 1142 include 'geodesic.inc' 1143 1144* WGS84 values 1145 a = 6378137d0 1146 f = 1/298.257223563d0 1147 r = 0 1148 1149 call area(a, f, lata, lona, 3, AA, PP) 1150 r = r + assert(PP, 36026861d0, 1d0) 1151 r = r + assert(AA, 0d0, 1d0) 1152 1153 tstp6 = r 1154 return 1155 end 1156 1157 integer function tstp12() 1158* AA of arctic circle (not really -- adjunct to rhumb-AA test) 1159 double precision lat(2), lon(2) 1160 data lat / 66.562222222d0, 66.562222222d0 / 1161 data lon / 0d0, 180d0 / 1162 double precision a, f, AA, PP 1163 integer r, assert 1164 include 'geodesic.inc' 1165 1166* WGS84 values 1167 a = 6378137d0 1168 f = 1/298.257223563d0 1169 r = 0 1170 1171 call area(a, f, lat, lon, 2, AA, PP) 1172 r = r + assert(PP, 10465729d0, 1d0) 1173 r = r + assert(AA, 0d0, 1d0) 1174 1175 tstp12 = r 1176 return 1177 end 1178 1179 integer function tstp13() 1180* Check encircling pole twice 1181 double precision lat(6), lon(6) 1182 data lat / 89d0, 89d0, 89d0, 89d0, 89d0, 89d0 / 1183 data lon / -360d0, -240d0, -120d0, 0d0, 120d0, 240d0 / 1184 double precision a, f, AA, PP 1185 integer r, assert 1186 include 'geodesic.inc' 1187 1188* WGS84 values 1189 a = 6378137d0 1190 f = 1/298.257223563d0 1191 r = 0 1192 1193 call area(a, f, lat, lon, 6, AA, PP) 1194 r = r + assert(PP, 1160741d0, 1d0) 1195 r = r + assert(AA, 32415230256.0d0, 1d0) 1196 1197 tstp13 = r 1198 return 1199 end 1200 1201 integer function tstp15() 1202* Coverage tests, includes Planimeter15 - Planimeter18 (combinations of 1203* reverse and sign). But flags aren't supported in the Fortran 1204* implementation. 1205 double precision lat(3), lon(3) 1206 data lat / 2d0, 1d0, 3d0 / 1207 data lon / 1d0, 2d0, 3d0 / 1208 double precision a, f, AA, PP 1209 integer r, assert 1210 include 'geodesic.inc' 1211 1212* WGS84 values 1213 a = 6378137d0 1214 f = 1/298.257223563d0 1215 r = 0 1216 1217 call area(a, f, lat, lon, 3, AA, PP) 1218 r = r + assert(AA, 18454562325.45119d0, 1d0) 1219* Interchanging lat and lon is equivalent to traversing the polygon 1220* backwards. 1221 call area(a, f, lon, lat, 3, AA, PP) 1222 r = r + assert(AA, -18454562325.45119d0, 1d0) 1223 1224 tstp15 = r 1225 return 1226 end 1227 1228 integer function tstp19() 1229* Coverage tests, includes Planimeter19 - Planimeter20 (degenerate 1230* polygons). 1231 double precision lat(1), lon(1) 1232 data lat / 1d0 / 1233 data lon / 1d0 / 1234 double precision a, f, AA, PP 1235 integer r, assert 1236 include 'geodesic.inc' 1237 1238* WGS84 values 1239 a = 6378137d0 1240 f = 1/298.257223563d0 1241 r = 0 1242 1243 call area(a, f, lat, lon, 1, AA, PP) 1244 r = r + assert(AA, 0d0, 0d0) 1245 r = r + assert(PP, 0d0, 0d0) 1246 1247 tstp19 = r 1248 return 1249 end 1250 1251 integer function tstp21() 1252* Some test to add code coverage: multiple circlings of pole (includes 1253* Planimeter21 - Planimeter28). 1254 double precision lat(12), lon(12), lonr(12) 1255 data lat / 12*45d0 / 1256 data lon / 60d0, 180d0, -60d0, 1257 + 60d0, 180d0, -60d0, 1258 + 60d0, 180d0, -60d0, 1259 + 60d0, 180d0, -60d0 / 1260 data lonr / -60d0, 180d0, 60d0, 1261 + -60d0, 180d0, 60d0, 1262 + -60d0, 180d0, 60d0, 1263 + -60d0, 180d0, 60d0 / 1264 double precision a, f, AA, PP, AA1 1265 integer r, assert 1266 include 'geodesic.inc' 1267 1268* WGS84 values 1269 a = 6378137d0 1270 f = 1/298.257223563d0 1271 r = 0 1272* Area for one circuit 1273 AA1 = 39433884866571.4277d0 1274 1275 do 10 i = 3,4 1276 call area(a, f, lat, lon, 3*i, AA, PP) 1277 r = r + assert(AA, AA1*i, 0.5d0) 1278 call area(a, f, lat, lonr, 3*i, AA, PP) 1279 r = r + assert(AA, -AA1*i, 0.5d0) 1280 10 continue 1281 1282 tstp21 = r 1283 return 1284 end 1285 1286 program geodtest 1287 integer n, i 1288 integer tstinv, tstdir, tstarc, 1289 + tstg0, tstg1, tstg2, tstg5, tstg6, tstg9, tstg10, tstg11, 1290 + tstg12, tstg14, tstg15, tstg17, tstg26, tstg28, tstg33, 1291 + tstg55, tstg59, tstg61, tstg73, tstg74, tstg76, tstg78, 1292 + tstg80, tstg84, tstg92, 1293 + tstp0, tstp5, tstp6, tstp12, tstp13, tstp15, tstp19, tstp21 1294 1295 n = 0 1296 i = tstinv() 1297 if (i .gt. 0) then 1298 n = n + 1 1299 print *, 'tstinv fail:', i 1300 end if 1301 i = tstdir() 1302 if (i .gt. 0) then 1303 n = n + 1 1304 print *, 'tstdir fail:', i 1305 end if 1306 i = tstarc() 1307 if (i .gt. 0) then 1308 n = n + 1 1309 print *, 'tstarc fail:', i 1310 end if 1311 i = tstg0() 1312 if (i .gt. 0) then 1313 n = n + 1 1314 print *, 'tstg0 fail:', i 1315 end if 1316 i = tstg1() 1317 if (i .gt. 0) then 1318 n = n + 1 1319 print *, 'tstg1 fail:', i 1320 end if 1321 i = tstg2() 1322 if (i .gt. 0) then 1323 n = n + 1 1324 print *, 'tstg2 fail:', i 1325 end if 1326 i = tstg5() 1327 if (i .gt. 0) then 1328 n = n + 1 1329 print *, 'tstg5 fail:', i 1330 end if 1331 i = tstg6() 1332 if (i .gt. 0) then 1333 n = n + 1 1334 print *, 'tstg6 fail:', i 1335 end if 1336 i = tstg9() 1337 if (i .gt. 0) then 1338 n = n + 1 1339 print *, 'tstg9 fail:', i 1340 end if 1341 i = tstg10() 1342 if (i .gt. 0) then 1343 n = n + 1 1344 print *, 'tstg10 fail:', i 1345 end if 1346 i = tstg11() 1347 if (i .gt. 0) then 1348 n = n + 1 1349 print *, 'tstg11 fail:', i 1350 end if 1351 i = tstg12() 1352 if (i .gt. 0) then 1353 n = n + 1 1354 print *, 'tstg12 fail:', i 1355 end if 1356 i = tstg14() 1357 if (i .gt. 0) then 1358 n = n + 1 1359 print *, 'tstg14 fail:', i 1360 end if 1361 i = tstg15() 1362 if (i .gt. 0) then 1363 n = n + 1 1364 print *, 'tstg15 fail:', i 1365 end if 1366 i = tstg17() 1367 if (i .gt. 0) then 1368 n = n + 1 1369 print *, 'tstg17 fail:', i 1370 end if 1371 i = tstg26() 1372 if (i .gt. 0) then 1373 n = n + 1 1374 print *, 'tstg26 fail:', i 1375 end if 1376 i = tstg28() 1377 if (i .gt. 0) then 1378 n = n + 1 1379 print *, 'tstg28 fail:', i 1380 end if 1381 i = tstg33() 1382 if (i .gt. 0) then 1383 n = n + 1 1384 print *, 'tstg33 fail:', i 1385 end if 1386 i = tstg55() 1387 if (i .gt. 0) then 1388 n = n + 1 1389 print *, 'tstg55 fail:', i 1390 end if 1391 i = tstg59() 1392 if (i .gt. 0) then 1393 n = n + 1 1394 print *, 'tstg59 fail:', i 1395 end if 1396 i = tstg61() 1397 if (i .gt. 0) then 1398 n = n + 1 1399 print *, 'tstg61 fail:', i 1400 end if 1401 i = tstg73() 1402 if (i .gt. 0) then 1403 n = n + 1 1404 print *, 'tstg73 fail:', i 1405 end if 1406 i = tstg74() 1407 if (i .gt. 0) then 1408 n = n + 1 1409 print *, 'tstg74 fail:', i 1410 end if 1411 i = tstg76() 1412 if (i .gt. 0) then 1413 n = n + 1 1414 print *, 'tstg76 fail:', i 1415 end if 1416 i = tstg78() 1417 if (i .gt. 0) then 1418 n = n + 1 1419 print *, 'tstg78 fail:', i 1420 end if 1421 i = tstg80() 1422 if (i .gt. 0) then 1423 n = n + 1 1424 print *, 'tstg80 fail:', i 1425 end if 1426 i = tstg84() 1427 if (i .gt. 0) then 1428 n = n + 1 1429 print *, 'tstg84 fail:', i 1430 end if 1431 i = tstg92() 1432 if (i .gt. 0) then 1433 n = n + 1 1434 print *, 'tstg92 fail:', i 1435 end if 1436 i = tstp0() 1437 if (i .gt. 0) then 1438 n = n + 1 1439 print *, 'tstp0 fail:', i 1440 end if 1441 i = tstp5() 1442 if (i .gt. 0) then 1443 n = n + 1 1444 print *, 'tstp5 fail:', i 1445 end if 1446 i = tstp6() 1447 if (i .gt. 0) then 1448 n = n + 1 1449 print *, 'tstp6 fail:', i 1450 end if 1451 i = tstp12() 1452 if (i .gt. 0) then 1453 n = n + 1 1454 print *, 'tstp12 fail:', i 1455 end if 1456 i = tstp13() 1457 if (i .gt. 0) then 1458 n = n + 1 1459 print *, 'tstp13 fail:', i 1460 end if 1461 i = tstp15() 1462 if (i .gt. 0) then 1463 n = n + 1 1464 print *, 'tstp15 fail:', i 1465 end if 1466 i = tstp19() 1467 if (i .gt. 0) then 1468 n = n + 1 1469 print *, 'tstp19 fail:', i 1470 end if 1471 i = tstp21() 1472 if (i .gt. 0) then 1473 n = n + 1 1474 print *, 'tstp21 fail:', i 1475 end if 1476 1477 if (n .gt. 0) then 1478 stop 1 1479 end if 1480 1481 stop 1482 end 1483 1484*> @endcond SKIP 1485