1function geographiclib_test 2%GEOGRAPHICLIB_TEST The test suite for the geographiclib package 3% 4% GEOGRAPHICLIB_TEST 5% 6% runs a variety of tests and produces no output it they are successful. 7 8 n = 0; 9 i = testrand; if i, n=n+1; fprintf('testrand fail: %d\n', i); end 10 i = GeoConvert0 ; if i, n=n+1; fprintf('GeoConvert0 fail: %d\n', i); end 11 i = GeoConvert8 ; if i, n=n+1; fprintf('GeoConvert8 fail: %d\n', i); end 12 i = GeoConvert16; if i, n=n+1; fprintf('GeoConvert16 fail: %d\n', i); end 13 i = GeoConvert17; if i, n=n+1; fprintf('GeoConvert17 fail: %d\n', i); end 14 i = GeodSolve0 ; if i, n=n+1; fprintf('GeodSolve0 fail: %d\n', i); end 15 i = GeodSolve1 ; if i, n=n+1; fprintf('GeodSolve1 fail: %d\n', i); end 16 i = GeodSolve2 ; if i, n=n+1; fprintf('GeodSolve2 fail: %d\n', i); end 17 i = GeodSolve4 ; if i, n=n+1; fprintf('GeodSolve4 fail: %d\n', i); end 18 i = GeodSolve5 ; if i, n=n+1; fprintf('GeodSolve5 fail: %d\n', i); end 19 i = GeodSolve6 ; if i, n=n+1; fprintf('GeodSolve6 fail: %d\n', i); end 20 i = GeodSolve9 ; if i, n=n+1; fprintf('GeodSolve9 fail: %d\n', i); end 21 i = GeodSolve10; if i, n=n+1; fprintf('GeodSolve10 fail: %d\n', i); end 22 i = GeodSolve11; if i, n=n+1; fprintf('GeodSolve11 fail: %d\n', i); end 23 i = GeodSolve12; if i, n=n+1; fprintf('GeodSolve12 fail: %d\n', i); end 24 i = GeodSolve14; if i, n=n+1; fprintf('GeodSolve14 fail: %d\n', i); end 25 i = GeodSolve15; if i, n=n+1; fprintf('GeodSolve15 fail: %d\n', i); end 26 i = GeodSolve17; if i, n=n+1; fprintf('GeodSolve17 fail: %d\n', i); end 27 i = GeodSolve26; if i, n=n+1; fprintf('GeodSolve26 fail: %d\n', i); end 28 i = GeodSolve28; if i, n=n+1; fprintf('GeodSolve28 fail: %d\n', i); end 29 i = GeodSolve33; if i, n=n+1; fprintf('GeodSolve33 fail: %d\n', i); end 30 i = GeodSolve55; if i, n=n+1; fprintf('GeodSolve55 fail: %d\n', i); end 31 i = GeodSolve59; if i, n=n+1; fprintf('GeodSolve59 fail: %d\n', i); end 32 i = GeodSolve61; if i, n=n+1; fprintf('GeodSolve61 fail: %d\n', i); end 33 i = GeodSolve73; if i, n=n+1; fprintf('GeodSolve73 fail: %d\n', i); end 34 i = GeodSolve74; if i, n=n+1; fprintf('GeodSolve74 fail: %d\n', i); end 35 i = GeodSolve76; if i, n=n+1; fprintf('GeodSolve76 fail: %d\n', i); end 36 i = GeodSolve78; if i, n=n+1; fprintf('GeodSolve78 fail: %d\n', i); end 37 i = GeodSolve80; if i, n=n+1; fprintf('GeodSolve80 fail: %d\n', i); end 38 i = GeodSolve84; if i, n=n+1; fprintf('GeodSolve84 fail: %d\n', i); end 39 i = GeodSolve92; if i, n=n+1; fprintf('GeodSolve92 fail: %d\n', i); end 40 i = Planimeter0 ; if i, n=n+1; fprintf('Planimeter0 fail: %d\n', i); end 41 i = Planimeter5 ; if i, n=n+1; fprintf('Planimeter5 fail: %d\n', i); end 42 i = Planimeter6 ; if i, n=n+1; fprintf('Planimeter6 fail: %d\n', i); end 43 i = Planimeter12; if i, n=n+1; fprintf('Planimeter12 fail: %d\n', i); end 44 i = Planimeter13; if i, n=n+1; fprintf('Planimeter13 fail: %d\n', i); end 45 i = Planimeter15; if i, n=n+1; fprintf('Planimeter15 fail: %d\n', i); end 46 i = Planimeter19; if i, n=n+1; fprintf('Planimeter19 fail: %d\n', i); end 47 i = Planimeter21; if i, n=n+1; fprintf('Planimeter21 fail: %d\n', i); end 48 i = TransverseMercatorProj1; 49 if i, n=n+1; fprintf('TransverseMercatorProj1 fail: %d\n', i); end 50 i = TransverseMercatorProj3; 51 if i, n=n+1; fprintf('TransverseMercatorProj3 fail: %d\n', i); end 52 i = TransverseMercatorProj5; 53 if i, n=n+1; fprintf('TransverseMercatorProj5 fail: %d\n', i); end 54 i = geodreckon0; if i, n=n+1; fprintf('geodreckon0 fail: %d\n', i); end 55 i = gedistance0; if i, n=n+1; fprintf('gedistance0 fail: %d\n', i); end 56 i = tranmerc0; if i, n=n+1; fprintf('tranmerc0 fail: %d\n', i); end 57 i = mgrs0; if i, n=n+1; fprintf('mgrs0 fail: %d\n', i); end 58 i = mgrs1; if i, n=n+1; fprintf('mgrs1 fail: %d\n', i); end 59 i = mgrs2; if i, n=n+1; fprintf('mgrs2 fail: %d\n', i); end 60 i = mgrs3; if i, n=n+1; fprintf('mgrs3 fail: %d\n', i); end 61 i = mgrs4; if i, n=n+1; fprintf('mgrs4 fail: %d\n', i); end 62 i = mgrs5; if i, n=n+1; fprintf('mgrs5 fail: %d\n', i); end 63 % check for suppression of "warning: division by zero" in octave 64 [~, ~, ~, ~, ~, ~, ~, s12] = geodreckon(-30, 0, 180, 90, 1); 65 assert(~isnan(s12)); 66 assert(n == 0); 67end 68 69function n = assertEquals(x, y, d) 70 n = abs(x - y) <= d; 71 n = sum(~n(:)); 72end 73 74function n = assertNaN(x) 75 n = isnan(x); 76 n = sum(~n(:)); 77end 78 79function n = testrand 80 n = 0; 81 testcases = [ 82 35.60777, -139.44815, 111.098748429560326, ... 83 -11.17491, -69.95921, 129.289270889708762, ... 84 8935244.5604818305, 80.50729714281974, 6273170.2055303837, ... 85 0.16606318447386067, 0.16479116945612937, 12841384694976.432; 86 55.52454, 106.05087, 22.020059880982801, ... 87 77.03196, 197.18234, 109.112041110671519, ... 88 4105086.1713924406, 36.892740690445894, 3828869.3344387607, ... 89 0.80076349608092607, 0.80101006984201008, 61674961290615.615; 90 -21.97856, 142.59065, -32.44456876433189, ... 91 41.84138, 98.56635, -41.84359951440466, ... 92 8394328.894657671, 75.62930491011522, 6161154.5773110616, ... 93 0.24816339233950381, 0.24930251203627892, -6637997720646.717; 94 -66.99028, 112.2363, 173.73491240878403, ... 95 -12.70631, 285.90344, 2.512956620913668, ... 96 11150344.2312080241, 100.278634181155759, 6289939.5670446687, ... 97 -0.17199490274700385, -0.17722569526345708, -121287239862139.744; 98 -17.42761, 173.34268, -159.033557661192928, ... 99 -15.84784, 5.93557, -20.787484651536988, ... 100 16076603.1631180673, 144.640108810286253, 3732902.1583877189, ... 101 -0.81273638700070476, -0.81299800519154474, 97825992354058.708; 102 32.84994, 48.28919, 150.492927788121982, ... 103 -56.28556, 202.29132, 48.113449399816759, ... 104 16727068.9438164461, 150.565799985466607, 3147838.1910180939, ... 105 -0.87334918086923126, -0.86505036767110637, -72445258525585.010; 106 6.96833, 52.74123, 92.581585386317712, ... 107 -7.39675, 206.17291, 90.721692165923907, ... 108 17102477.2496958388, 154.147366239113561, 2772035.6169917581, ... 109 -0.89991282520302447, -0.89986892177110739, -1311796973197.995; 110 -50.56724, -16.30485, -105.439679907590164, ... 111 -33.56571, -94.97412, -47.348547835650331, ... 112 6455670.5118668696, 58.083719495371259, 5409150.7979815838, ... 113 0.53053508035997263, 0.52988722644436602, 41071447902810.047; 114 -58.93002, -8.90775, 140.965397902500679, ... 115 -8.91104, 133.13503, 19.255429433416599, ... 116 11756066.0219864627, 105.755691241406877, 6151101.2270708536, ... 117 -0.26548622269867183, -0.27068483874510741, -86143460552774.735; 118 -68.82867, -74.28391, 93.774347763114881, ... 119 -50.63005, -8.36685, 34.65564085411343, ... 120 3956936.926063544, 35.572254987389284, 3708890.9544062657, ... 121 0.81443963736383502, 0.81420859815358342, -41845309450093.787; 122 -10.62672, -32.0898, -86.426713286747751, ... 123 5.883, -134.31681, -80.473780971034875, ... 124 11470869.3864563009, 103.387395634504061, 6184411.6622659713, ... 125 -0.23138683500430237, -0.23155097622286792, 4198803992123.548; 126 -21.76221, 166.90563, 29.319421206936428, ... 127 48.72884, 213.97627, 43.508671946410168, ... 128 9098627.3986554915, 81.963476716121964, 6299240.9166992283, ... 129 0.13965943368590333, 0.14152969707656796, 10024709850277.476; 130 -19.79938, -174.47484, 71.167275780171533, ... 131 -11.99349, -154.35109, 65.589099775199228, ... 132 2319004.8601169389, 20.896611684802389, 2267960.8703918325, ... 133 0.93427001867125849, 0.93424887135032789, -3935477535005.785; 134 -11.95887, -116.94513, 92.712619830452549, ... 135 4.57352, 7.16501, 78.64960934409585, ... 136 13834722.5801401374, 124.688684161089762, 5228093.177931598, ... 137 -0.56879356755666463, -0.56918731952397221, -9919582785894.853; 138 -87.85331, 85.66836, -65.120313040242748, ... 139 66.48646, 16.09921, -4.888658719272296, ... 140 17286615.3147144645, 155.58592449699137, 2635887.4729110181, ... 141 -0.90697975771398578, -0.91095608883042767, 42667211366919.534; 142 1.74708, 128.32011, -101.584843631173858, ... 143 -11.16617, 11.87109, -86.325793296437476, ... 144 12942901.1241347408, 116.650512484301857, 5682744.8413270572, ... 145 -0.44857868222697644, -0.44824490340007729, 10763055294345.653; 146 -25.72959, -144.90758, -153.647468693117198, ... 147 -57.70581, -269.17879, -48.343983158876487, ... 148 9413446.7452453107, 84.664533838404295, 6356176.6898881281, ... 149 0.09492245755254703, 0.09737058264766572, 74515122850712.444; 150 -41.22777, 122.32875, 14.285113402275739, ... 151 -7.57291, 130.37946, 10.805303085187369, ... 152 3812686.035106021, 34.34330804743883, 3588703.8812128856, ... 153 0.82605222593217889, 0.82572158200920196, -2456961531057.857; 154 11.01307, 138.25278, 79.43682622782374, ... 155 6.62726, 247.05981, 103.708090215522657, ... 156 11911190.819018408, 107.341669954114577, 6070904.722786735, ... 157 -0.29767608923657404, -0.29785143390252321, 17121631423099.696; 158 -29.47124, 95.14681, -163.779130441688382, ... 159 -27.46601, -69.15955, -15.909335945554969, ... 160 13487015.8381145492, 121.294026715742277, 5481428.9945736388, ... 161 -0.51527225545373252, -0.51556587964721788, 104679964020340.318]; 162 163 lat1 = testcases(:,1); lon1 = testcases(:,2); azi1 = testcases(:,3); 164 lat2 = testcases(:,4); lon2 = testcases(:,5); azi2 = testcases(:,6); 165 s12 = testcases(:,7); a12 = testcases(:,8); m12 = testcases(:,9); 166 M12 = testcases(:,10); M21 = testcases(:,11); S12 = testcases(:,12); 167 [s12a, azi1a, azi2a, S12a, m12a, M12a, M21a, a12a] = ... 168 geoddistance(lat1, lon1, lat2, lon2); 169 n = n + assertEquals(azi1, azi1a, 1e-13); 170 n = n + assertEquals(azi2, azi2a, 1e-13); 171 n = n + assertEquals(s12, s12a, 1e-8); 172 n = n + assertEquals(a12, a12a, 1e-13); 173 n = n + assertEquals(m12, m12a, 1e-8); 174 n = n + assertEquals(M12, M12a, 1e-15); 175 n = n + assertEquals(M21, M21a, 1e-15); 176 n = n + assertEquals(S12, S12a, 0.1); 177 178 [lat2a, lon2a, azi2a, S12a, m12a, M12a, M21a, a12a] = ... 179 geodreckon(lat1, lon1, s12, azi1, 2); 180 n = n + assertEquals(lat2, lat2a, 1e-13); 181 n = n + assertEquals(lon2, lon2a, 1e-13); 182 n = n + assertEquals(azi2, azi2a, 1e-13); 183 n = n + assertEquals(a12, a12a, 1e-13); 184 n = n + assertEquals(m12, m12a, 1e-8); 185 n = n + assertEquals(M12, M12a, 1e-15); 186 n = n + assertEquals(M21, M21a, 1e-15); 187 n = n + assertEquals(S12, S12a, 0.1); 188 189 [lat2a, lon2a, azi2a, S12a, m12a, M12a, M21a, s12a] = ... 190 geodreckon(lat1, lon1, a12, azi1, 1+2); 191 n = n + assertEquals(lat2, lat2a, 1e-13); 192 n = n + assertEquals(lon2, lon2a, 1e-13); 193 n = n + assertEquals(azi2, azi2a, 1e-13); 194 n = n + assertEquals(s12, s12a, 1e-8); 195 n = n + assertEquals(m12, m12a, 1e-8); 196 n = n + assertEquals(M12, M12a, 1e-15); 197 n = n + assertEquals(M21, M21a, 1e-15); 198 n = n + assertEquals(S12, S12a, 0.1); 199end 200 201function ell = ellipsoid(a, f) 202 ell = [a, flat2ecc(f)]; 203end 204 205function n = GeoConvert0 206 n = 0; 207 [x, y, zone, isnorth] = utmups_fwd(33.3, 44.4); 208 mgrs = mgrs_fwd(x, y, zone, isnorth, 2); 209 n = n + ~strcmp(mgrs, '38SMB4484'); 210end 211 212function n = GeoConvert8 213% Check fix to PolarStereographic es initialization blunder (2015-05-18) 214 n = 0; 215 [x, y, zone, isnorth] = utmups_fwd(86, 0); 216 n = n + ~(zone == 0); 217 n = n + ~(isnorth); 218 n = n + assertEquals(x, 2000000, 0.5e-6); 219 n = n + assertEquals(y, 1555731.570643, 0.5e-6); 220end 221 222function n = GeoConvert16 223% Check MGRS::Forward improved rounding fix, 2015-07-22 224 n = 0; 225 mgrs = mgrs_fwd(444140.6, 3684706.3, 38, true, 8); 226 n = n + ~strcmp(mgrs, '38SMB4414060084706300'); 227end 228 229function n = GeoConvert17 230% Check MGRS::Forward digit consistency fix, 2015-07-23 231 n = 0; 232 mgrs = mgrs_fwd(500000, 63.811, 38, true, 8); 233 n = n + ~strcmp(mgrs, '38NNF0000000000063811'); 234 mgrs = mgrs_fwd(500000, 63.811, 38, true, 9); 235 n = n + ~strcmp(mgrs, '38NNF000000000000638110'); 236end 237 238function n = GeodSolve0 239 n = 0; 240 [s12, azi1, azi2] = geoddistance(40.6, -73.8, 49.01666667, 2.55); 241 n = n + assertEquals(azi1, 53.47022, 0.5e-5); 242 n = n + assertEquals(azi2, 111.59367, 0.5e-5); 243 n = n + assertEquals(s12, 5853226, 0.5); 244end 245 246function n = GeodSolve1 247 n = 0; 248 [lat2, lon2, azi2] = geodreckon(40.63972222, -73.77888889, 5850e3, 53.5); 249 n = n + assertEquals(lat2, 49.01467, 0.5e-5); 250 n = n + assertEquals(lon2, 2.56106, 0.5e-5); 251 n = n + assertEquals(azi2, 111.62947, 0.5e-5); 252end 253 254function n = GeodSolve2 255% Check fix for antipodal prolate bug found 2010-09-04 256 n = 0; 257 ell = ellipsoid(6.4e6, -1/150.0); 258 [s12, azi1, azi2] = geoddistance(0.07476, 0, -0.07476, 180, ell); 259 n = n + assertEquals(azi1, 90.00078, 0.5e-5); 260 n = n + assertEquals(azi2, 90.00078, 0.5e-5); 261 n = n + assertEquals(s12, 20106193, 0.5); 262 [s12, azi1, azi2] = geoddistance(0.1, 0, -0.1, 180, ell); 263 n = n + assertEquals(azi1, 90.00105, 0.5e-5); 264 n = n + assertEquals(azi2, 90.00105, 0.5e-5); 265 n = n + assertEquals(s12, 20106193, 0.5); 266end 267 268function n = GeodSolve4 269% Check fix for short line bug found 2010-05-21 270% This also checks the MATLAB specific bug: 271% Ensure that Lengths in geoddistance is not invoked with zero-length 272% vectors, 2015-08-25. 273 n = 0; 274 s12 = geoddistance(36.493349428792, 0, 36.49334942879201, .0000008); 275 n = n + assertEquals(s12, 0.072, 0.5e-3); 276end 277 278function n = GeodSolve5 279% Check fix for point2=pole bug found 2010-05-03 280 n = 0; 281 [lat2, lon2, azi2] = geodreckon(0.01777745589997, 30, 10e6, 0); 282 n = n + assertEquals(lat2, 90, 0.5e-5); 283 if lon2 < 0 284 n = n + assertEquals(lon2, -150, 0.5e-5); 285 n = n + assertEquals(abs(azi2), 180, 0.5e-5); 286 else 287 n = n + assertEquals(lon2, 30, 0.5e-5); 288 n = n + assertEquals(azi2, 0, 0.5e-5); 289 end 290end 291 292function n = GeodSolve6 293% Check fix for volatile sbet12a bug found 2011-06-25 (gcc 4.4.4 294% x86 -O3). Found again on 2012-03-27 with tdm-mingw32 (g++ 4.6.1). 295 n = 0; 296 s12 = geoddistance(88.202499451857, 0, ... 297 -88.202499451857, 179.981022032992859592); 298 n = n + assertEquals(s12, 20003898.214, 0.5e-3); 299 s12 = geoddistance(89.262080389218, 0, ... 300 -89.262080389218, 179.992207982775375662); 301 n = n + assertEquals(s12, 20003925.854, 0.5e-3); 302 s12 = geoddistance(89.333123580033, 0, ... 303 -89.333123580032997687, 179.99295812360148422); 304 n = n + assertEquals(s12, 20003926.881, 0.5e-3); 305end 306 307function n = GeodSolve9 308% Check fix for volatile x bug found 2011-06-25 (gcc 4.4.4 x86 -O3) 309 n = 0; 310 s12 = geoddistance(56.320923501171, 0, ... 311 -56.320923501171, 179.664747671772880215); 312 n = n + assertEquals(s12, 19993558.287, 0.5e-3); 313end 314 315function n = GeodSolve10 316% Check fix for adjust tol1_ bug found 2011-06-25 (Visual Studio 317% 10 rel + debug) 318 n = 0; 319 s12 = geoddistance(52.784459512564, 0, ... 320 -52.784459512563990912, 179.634407464943777557); 321 n = n + assertEquals(s12, 19991596.095, 0.5e-3); 322end 323 324function n = GeodSolve11 325% Check fix for bet2 = -bet1 bug found 2011-06-25 (Visual Studio 326% 10 rel + debug) 327 n = 0; 328 s12 = geoddistance(48.522876735459, 0, ... 329 -48.52287673545898293, 179.599720456223079643); 330 n = n + assertEquals(s12, 19989144.774, 0.5e-3); 331end 332 333function n = GeodSolve12 334% Check fix for inverse geodesics on extreme prolate/oblate 335% ellipsoids Reported 2012-08-29 Stefan Guenther 336% <stefan.gunther@embl.de>; fixed 2012-10-07 337 n = 0; 338 ell = ellipsoid(89.8, -1.83); 339 [s12, azi1, azi2] = geoddistance(0, 0, -10, 160, ell); 340 n = n + assertEquals(azi1, 120.27, 1e-2); 341 n = n + assertEquals(azi2, 105.15, 1e-2); 342 n = n + assertEquals(s12, 266.7, 1e-1); 343end 344 345function n = GeodSolve14 346% Check fix for inverse ignoring lon12 = nan 347 n = 0; 348 [s12, azi1, azi2] = geoddistance(0, 0, 1, NaN); 349 n = n + assertNaN(azi1); 350 n = n + assertNaN(azi2); 351 n = n + assertNaN(s12); 352end 353 354function n = GeodSolve15 355% Initial implementation of Math::eatanhe was wrong for e^2 < 0. This 356% checks that this is fixed. 357 n = 0; 358 ell = ellipsoid(6.4e6, -1/150.0); 359 [~, ~, ~, S12] = geodreckon(1, 2, 4, 3, ell); 360 n = n + assertEquals(S12, 23700, 0.5); 361end 362 363function n = GeodSolve17 364% Check fix for LONG_UNROLL bug found on 2015-05-07 365 n = 0; 366 [lat2, lon2, azi2] = geodreckon(40, -75, 2e7, -10, 2); 367 n = n + assertEquals(lat2, -39, 1); 368 n = n + assertEquals(lon2, -254, 1); 369 n = n + assertEquals(azi2, -170, 1); 370 [lat2, lon2, azi2] = geodreckon(40, -75, 2e7, -10); 371 n = n + assertEquals(lat2, -39, 1); 372 n = n + assertEquals(lon2, 105, 1); 373 n = n + assertEquals(azi2, -170, 1); 374end 375 376function n = GeodSolve26 377% Check 0/0 problem with area calculation on sphere 2015-09-08 378 n = 0; 379 ell = ellipsoid(6.4e6, 0); 380 [~, ~, ~, S12] = geoddistance(1, 2, 3, 4, ell); 381 n = n + assertEquals(S12, 49911046115.0, 0.5); 382end 383 384function n = GeodSolve28 385% Check for bad placement of assignment of r.a12 with |f| > 0.01 (bug in 386% Java implementation fixed on 2015-05-19). 387 n = 0; 388 ell = ellipsoid(6.4e6, 0.1); 389 [~, ~, ~, ~, ~, ~, ~, a12] = geodreckon(1, 2, 5e6, 10, ell); 390 n = n + assertEquals(a12, 48.55570690, 0.5e-8); 391end 392 393function n = GeodSolve33 394% Check max(-0.0,+0.0) issues 2015-08-22 (triggered by bugs in Octave -- 395% sind(-0.0) = +0.0 -- and in some version of Visual Studio -- 396% fmod(-0.0, 360.0) = +0.0. 397 n = 0; 398 [s12, azi1, azi2] = geoddistance(0, 0, 0, 179); 399 n = n + assertEquals(azi1, 90.00000, 0.5e-5); 400 n = n + assertEquals(azi2, 90.00000, 0.5e-5); 401 n = n + assertEquals(s12, 19926189, 0.5); 402 [s12, azi1, azi2] = geoddistance(0, 0, 0, 179.5); 403 n = n + assertEquals(azi1, 55.96650, 0.5e-5); 404 n = n + assertEquals(azi2, 124.03350, 0.5e-5); 405 n = n + assertEquals(s12, 19980862, 0.5); 406 [s12, azi1, azi2] = geoddistance(0, 0, 0, 180); 407 n = n + assertEquals(azi1, 0.00000, 0.5e-5); 408 n = n + assertEquals(abs(azi2), 180.00000, 0.5e-5); 409 n = n + assertEquals(s12, 20003931, 0.5); 410 [s12, azi1, azi2] = geoddistance(0, 0, 1, 180); 411 n = n + assertEquals(azi1, 0.00000, 0.5e-5); 412 n = n + assertEquals(abs(azi2), 180.00000, 0.5e-5); 413 n = n + assertEquals(s12, 19893357, 0.5); 414 ell = ellipsoid(6.4e6, 0); 415 [s12, azi1, azi2] = geoddistance(0, 0, 0, 179, ell); 416 n = n + assertEquals(azi1, 90.00000, 0.5e-5); 417 n = n + assertEquals(azi2, 90.00000, 0.5e-5); 418 n = n + assertEquals(s12, 19994492, 0.5); 419 [s12, azi1, azi2] = geoddistance(0, 0, 0, 180, ell); 420 n = n + assertEquals(azi1, 0.00000, 0.5e-5); 421 n = n + assertEquals(abs(azi2), 180.00000, 0.5e-5); 422 n = n + assertEquals(s12, 20106193, 0.5); 423 [s12, azi1, azi2] = geoddistance(0, 0, 1, 180, ell); 424 n = n + assertEquals(azi1, 0.00000, 0.5e-5); 425 n = n + assertEquals(abs(azi2), 180.00000, 0.5e-5); 426 n = n + assertEquals(s12, 19994492, 0.5); 427 ell = ellipsoid(6.4e6, -1/300.0); 428 [s12, azi1, azi2] = geoddistance(0, 0, 0, 179, ell); 429 n = n + assertEquals(azi1, 90.00000, 0.5e-5); 430 n = n + assertEquals(azi2, 90.00000, 0.5e-5); 431 n = n + assertEquals(s12, 19994492, 0.5); 432 [s12, azi1, azi2] = geoddistance(0, 0, 0, 180, ell); 433 n = n + assertEquals(azi1, 90.00000, 0.5e-5); 434 n = n + assertEquals(azi2, 90.00000, 0.5e-5); 435 n = n + assertEquals(s12, 20106193, 0.5); 436 [s12, azi1, azi2] = geoddistance(0, 0, 0.5, 180, ell); 437 n = n + assertEquals(azi1, 33.02493, 0.5e-5); 438 n = n + assertEquals(azi2, 146.97364, 0.5e-5); 439 n = n + assertEquals(s12, 20082617, 0.5); 440 [s12, azi1, azi2] = geoddistance(0, 0, 1, 180, ell); 441 n = n + assertEquals(azi1, 0.00000, 0.5e-5); 442 n = n + assertEquals(abs(azi2), 180.00000, 0.5e-5); 443 n = n + assertEquals(s12, 20027270, 0.5); 444 % Check also octave-specific versions of this problem. 445 % In 1.44 this returned [-2.0004e+07, -2.0004e+07, 0.0000e+00, 0.0000e+00] 446 s12 = geoddistance(0,0,0,[179.5, 179.6, 180, 180]); 447 n = n + assertEquals(s12, [19980862, 19989165, 20003931, 20003931], 0.5); 448end 449 450function n = GeodSolve55 451% Check fix for nan + point on equator or pole not returning all nans in 452% Geodesic::Inverse, found 2015-09-23. 453 n = 0; 454 [s12, azi1, azi2] = geoddistance(NaN, 0, 0, 90); 455 n = n + assertNaN(azi1); 456 n = n + assertNaN(azi2); 457 n = n + assertNaN(s12); 458 [s12, azi1, azi2] = geoddistance(NaN, 0, 90, 9); 459 n = n + assertNaN(azi1); 460 n = n + assertNaN(azi2); 461 n = n + assertNaN(s12); 462end 463 464function n = GeodSolve59 465% Check for points close with longitudes close to 180 deg apart. 466 n = 0; 467 [s12, azi1, azi2] = geoddistance(5, 0.00000000000001, 10, 180); 468 n = n + assertEquals(azi1, 0.000000000000035, 1.5e-14); 469 n = n + assertEquals(azi2, 179.99999999999996, 1.5e-14); 470 n = n + assertEquals(s12, 18345191.174332713, 5e-9); 471end 472 473function n = GeodSolve61 474% Make sure small negative azimuths are west-going 475 n = 0; 476 [lat2, lon2, azi2] = geodreckon(45, 0, 1e7, -0.000000000000000003, 2); 477 n = n + assertEquals(lat2, 45.30632, 0.5e-5); 478 n = n + assertEquals(lon2, -180, 0.5e-5); 479 n = n + assertEquals(abs(azi2), 180, 0.5e-5); 480end 481 482function n = GeodSolve73 483% Check for backwards from the pole bug reported by Anon on 2016-02-13. 484% This only affected the Java implementation. It was introduced in Java 485% version 1.44 and fixed in 1.46-SNAPSHOT on 2016-01-17. 486% Also the + sign on azi2 is a check on the normalizing of azimuths 487% (converting -0.0 to +0.0). 488 n = 0; 489 [lat2, lon2, azi2] = geodreckon(90, 10, -1e6, 180); 490 n = n + assertEquals(lat2, 81.04623, 0.5e-5); 491 n = n + assertEquals(lon2, -170, 0.5e-5); 492 n = n + assertEquals(azi2, 0, 0); 493 n = n + assertEquals(copysignx(1, azi2), 1, 0); 494end 495 496function n = GeodSolve74 497% Check fix for inaccurate areas, bug introduced in v1.46, fixed 498% 2015-10-16. 499 n = 0; 500 [s12, azi1, azi2, S12, m12, M12, M21, a12] = ... 501 geoddistance(54.1589, 15.3872, 54.1591, 15.3877); 502 n = n + assertEquals(azi1, 55.723110355, 5e-9); 503 n = n + assertEquals(azi2, 55.723515675, 5e-9); 504 n = n + assertEquals(s12, 39.527686385, 5e-9); 505 n = n + assertEquals(a12, 0.000355495, 5e-9); 506 n = n + assertEquals(m12, 39.527686385, 5e-9); 507 n = n + assertEquals(M12, 0.999999995, 5e-9); 508 n = n + assertEquals(M21, 0.999999995, 5e-9); 509 n = n + assertEquals(S12, 286698586.30197, 5e-4); 510end 511 512function n = GeodSolve76 513% The distance from Wellington and Salamanca (a classic failure of 514% Vincenty) 515 n = 0; 516 [s12, azi1, azi2] = ... 517 geoddistance(-(41+19/60), 174+49/60, 40+58/60, -(5+30/60)); 518 n = n + assertEquals(azi1, 160.39137649664, 0.5e-11); 519 n = n + assertEquals(azi2, 19.50042925176, 0.5e-11); 520 n = n + assertEquals(s12, 19960543.857179, 0.5e-6); 521end 522 523function n = GeodSolve78 524% An example where the NGS calculator fails to converge 525 n = 0; 526 [s12, azi1, azi2] = geoddistance(27.2, 0.0, -27.1, 179.5); 527 n = n + assertEquals(azi1, 45.82468716758, 0.5e-11); 528 n = n + assertEquals(azi2, 134.22776532670, 0.5e-11); 529 n = n + assertEquals(s12, 19974354.765767, 0.5e-6); 530end 531 532function n = GeodSolve80 533% Some tests to add code coverage: computing scale in special cases + zero 534% length geodesic (includes GeodSolve80 - GeodSolve83). 535 n = 0; 536 [~, ~, ~, ~, ~, M12, M21] = geoddistance(0, 0, 0, 90); 537 n = n + assertEquals(M12, -0.00528427534, 0.5e-10); 538 n = n + assertEquals(M21, -0.00528427534, 0.5e-10); 539 540 [~, ~, ~, ~, ~, M12, M21] = geoddistance(0, 0, 1e-6, 1e-6); 541 n = n + assertEquals(M12, 1, 0.5e-10); 542 n = n + assertEquals(M21, 1, 0.5e-10); 543 544 [s12, azi1, azi2, S12, m12, M12, M21, a12] = ... 545 geoddistance(20.001, 0, 20.001, 0); 546 n = n + assertEquals(a12, 0, 1e-13); 547 n = n + assertEquals(s12, 0, 1e-8); 548 n = n + assertEquals(azi1, 180, 1e-13); 549 n = n + assertEquals(azi2, 180, 1e-13); 550 n = n + assertEquals(m12, 0, 1e-8); 551 n = n + assertEquals(M12, 1, 1e-15); 552 n = n + assertEquals(M21, 1, 1e-15); 553 n = n + assertEquals(S12, 0, 1e-10); 554 n = n + assertEquals(copysignx(1, a12), 1, 0); 555 n = n + assertEquals(copysignx(1, s12), 1, 0); 556 n = n + assertEquals(copysignx(1, m12), 1, 0); 557 558 [s12, azi1, azi2, S12, m12, M12, M21, a12] = ... 559 geoddistance(90, 0, 90, 180); 560 n = n + assertEquals(a12, 0, 1e-13); 561 n = n + assertEquals(s12, 0, 1e-8); 562 n = n + assertEquals(azi1, 0, 1e-13); 563 n = n + assertEquals(azi2, 180, 1e-13); 564 n = n + assertEquals(m12, 0, 1e-8); 565 n = n + assertEquals(M12, 1, 1e-15); 566 n = n + assertEquals(M21, 1, 1e-15); 567 n = n + assertEquals(S12, 127516405431022.0, 0.5); 568end 569 570function n = GeodSolve84 571% Tests for python implementation to check fix for range errors with 572% {fmod,sin,cos}(inf) (includes GeodSolve84 - GeodSolve86). 573 n = 0; 574 [lat2, lon2, azi2] = geodreckon(0, 0, inf, 90); 575 n = n + assertNaN(lat2); 576 n = n + assertNaN(lon2); 577 n = n + assertNaN(azi2); 578 [lat2, lon2, azi2] = geodreckon(0, 0, nan, 90); 579 n = n + assertNaN(lat2); 580 n = n + assertNaN(lon2); 581 n = n + assertNaN(azi2); 582 [lat2, lon2, azi2] = geodreckon(0, 0, 1000, inf); 583 n = n + assertNaN(lat2); 584 n = n + assertNaN(lon2); 585 n = n + assertNaN(azi2); 586 [lat2, lon2, azi2] = geodreckon(0, 0, 1000, nan); 587 n = n + assertNaN(lat2); 588 n = n + assertNaN(lon2); 589 n = n + assertNaN(azi2); 590 [lat2, lon2, azi2] = geodreckon(0, inf, 1000, 90); 591 n = n + assertEquals(lat2, 0, 0); 592 n = n + assertNaN(lon2); 593 n = n + assertEquals(azi2, 90, 0); 594 [lat2, lon2, azi2] = geodreckon(0, nan, 1000, 90); 595 n = n + assertEquals(lat2, 0, 0); 596 n = n + assertNaN(lon2); 597 n = n + assertEquals(azi2, 90, 0); 598 [lat2, lon2, azi2] = geodreckon(inf, 0, 1000, 90); 599 n = n + assertNaN(lat2); 600 n = n + assertNaN(lon2); 601 n = n + assertNaN(azi2); 602 [lat2, lon2, azi2] = geodreckon(nan, 0, 1000, 90); 603 n = n + assertNaN(lat2); 604 n = n + assertNaN(lon2); 605 n = n + assertNaN(azi2); 606end 607 608function n = GeodSolve92 609% Check fix for inaccurate hypot with python 3.[89]. Problem reported 610% by agdhruv https://github.com/geopy/geopy/issues/466 ; see 611% https://bugs.python.org/issue43088 612 n = 0; 613 [s12, azi1, azi2] = geoddistance(37.757540000000006, -122.47018, ... 614 37.75754, -122.470177); 615 n = n + assertEquals(azi1, 89.99999923, 1e-7 ); 616 n = n + assertEquals(azi2, 90.00000106, 1e-7 ); 617 n = n + assertEquals(s12, 0.264, 0.5e-3); 618end 619 620function n = Planimeter0 621% Check fix for pole-encircling bug found 2011-03-16 622 n = 0; 623 pa = [89, 0; 89, 90; 89, 180; 89, 270]; 624 pb = [-89, 0; -89, 90; -89, 180; -89, 270]; 625 pc = [0, -1; -1, 0; 0, 1; 1, 0]; 626 pd = [90, 0; 0, 0; 0, 90]; 627 628 [area, perimeter] = geodarea(pa(:,1), pa(:,2)); 629 n = n + assertEquals(perimeter, 631819.8745, 1e-4); 630 n = n + assertEquals(area, 24952305678.0, 1); 631 632 [area, perimeter] = geodarea(pb(:,1), pb(:,2)); 633 n = n + assertEquals(perimeter, 631819.8745, 1e-4); 634 n = n + assertEquals(area, -24952305678.0, 1); 635 636 [area, perimeter] = geodarea(pc(:,1), pc(:,2)); 637 n = n + assertEquals(perimeter, 627598.2731, 1e-4); 638 n = n + assertEquals(area, 24619419146.0, 1); 639 640 [area, perimeter] = geodarea(pd(:,1), pd(:,2)); 641 n = n + assertEquals(perimeter, 30022685, 1); 642 n = n + assertEquals(area, 63758202715511.0, 1); 643 644end 645 646function n = Planimeter5 647% Check fix for Planimeter pole crossing bug found 2011-06-24 648 n = 0; 649 points = [89, 0.1; 89, 90.1; 89, -179.9]; 650 [area, perimeter] = geodarea(points(:,1), points(:,2)); 651 n = n + assertEquals(perimeter, 539297, 1); 652 n = n + assertEquals(area, 12476152838.5, 1); 653end 654 655function n = Planimeter6 656% Check fix for Planimeter lon12 rounding bug found 2012-12-03 657 n = 0; 658 pa = [9, -0.00000000000001; 9, 180; 9, 0]; 659 pb = [9, 0.00000000000001; 9, 0; 9, 180]; 660 pc = [9, 0.00000000000001; 9, 180; 9, 0]; 661 pd = [9, -0.00000000000001; 9, 0; 9, 180]; 662 663 [area, perimeter] = geodarea(pa(:,1), pa(:,2)); 664 n = n + assertEquals(perimeter, 36026861, 1); 665 n = n + assertEquals(area, 0, 1); 666 [area, perimeter] = geodarea(pb(:,1), pb(:,2)); 667 n = n + assertEquals(perimeter, 36026861, 1); 668 n = n + assertEquals(area, 0, 1); 669 [area, perimeter] = geodarea(pc(:,1), pc(:,2)); 670 n = n + assertEquals(perimeter, 36026861, 1); 671 n = n + assertEquals(area, 0, 1); 672 [area, perimeter] = geodarea(pd(:,1), pd(:,2)); 673 n = n + assertEquals(perimeter, 36026861, 1); 674 n = n + assertEquals(area, 0, 1); 675end 676 677function n = Planimeter12 678% Area of arctic circle (not really -- adjunct to rhumb-area test) 679 n = 0; 680 points = [66.562222222, 0; 66.562222222, 180]; 681 [area, perimeter] = geodarea(points(:,1), points(:,2)); 682 n = n + assertEquals(perimeter, 10465729, 1); 683 n = n + assertEquals(area, 0, 1); 684end 685 686function n = Planimeter13 687% Check encircling pole twice 688 n = 0; 689 points = [89,-360; 89,-240; 89,-120; 89,0; 89,120; 89,240]; 690 [area, perimeter] = geodarea(points(:,1), points(:,2)); 691 n = n + assertEquals(perimeter, 1160741, 1); 692 n = n + assertEquals(area, 32415230256.0, 1); 693end 694 695function n = Planimeter15 696% Coverage tests, includes Planimeter15 - Planimeter18 (combinations of 697% reverse and sign). But flags aren't supported in the MATLAB 698% implementation. 699 n = 0; 700 points = [2,1; 1,2; 3,3]; 701 area = geodarea(points(:,1), points(:,2)); 702 n = n + assertEquals(area, 18454562325.45119, 1); 703 % Interchanging lat and lon is equivalent to traversing the polygon 704 % backwards. 705 area = geodarea(points(:,2), points(:,1)); 706 n = n + assertEquals(area, -18454562325.45119, 1); 707end 708 709function n = Planimeter19 710% Coverage tests, includes Planimeter19 - Planimeter20 (degenerate 711% polygons). 712 n = 0; 713 points = [1,1]; 714 [area, perimeter] = geodarea(points(:,1), points(:,2)); 715 n = n + assertEquals(area, 0, 0); 716 n = n + assertEquals(perimeter, 0, 0); 717end 718 719function n = Planimeter21 720% Some test to add code coverage: multiple circlings of pole (includes 721% Planimeter21 - Planimeter28). 722 n = 0; 723 points = [45 60;45 180;45 -60;... 724 45 60;45 180;45 -60;... 725 45 60;45 180;45 -60;... 726 45 60;45 180;45 -60;... 727 ]; 728 r = 39433884866571.4277; % Area for one circuit 729 for i = 3 : 4 730 area = geodarea(points(1:3*i,1), points(1:3*i,2)); 731 % if i ~= 4 732 n = n + assertEquals(area, i*r, 0.5); 733 % end 734 area = geodarea(points(3*i:-1:1,1), points(3*i:-1:1,2)); 735 % if i ~= 4 736 n = n + assertEquals(area, -i*r, 0.5); 737 % end 738 end 739end 740 741function n = TransverseMercatorProj1 742% Test fix to bad meridian convergence at pole with 743% TransverseMercatorExact found 2013-06-26 744 n = 0; 745 [x, y, gam, k] = tranmerc_fwd(0, 0, 90, 75); 746 n = n + assertEquals(x, 0, 0.5e-6); 747 n = n + assertEquals(y, 10001965.72935, 0.5e-4); 748 n = n + assertEquals(gam, 75, 0.5e-12); 749 n = n + assertEquals(k, 1, 0.5e-12); 750end 751 752function n = TransverseMercatorProj3 753% Test fix to bad scale at pole with TransverseMercatorExact 754% found 2013-06-30 (quarter meridian = 10001965.7293127228128889202m) 755 n = 0; 756 [lat, lon, gam, k] = tranmerc_inv(0, 0, 0, 10001965.7293127228); 757 n = n + assertEquals(lat, 90, 1e-11); 758 if abs(lon) < 90 759 n = n + assertEquals(lon, 0, 0.5e-12); 760 n = n + assertEquals(gam, 0, 0.5e-12); 761 else 762 n = n + assertEquals(abs(lon), 180, 0.5e-12); 763 n = n + assertEquals(abs(gam), 180, 0.5e-12); 764 end 765 n = n + assertEquals(k, 1.0, 0.5e-12); 766end 767 768function n = TransverseMercatorProj5 769% Generic tests for transverse Mercator added 2017-04-15 to check use of 770% complex arithmetic to do Clenshaw sum. 771 n = 0; 772 k0 = 0.9996; 773 ell = ellipsoid(6.4e6, 1/150); 774 [x, y, gam, k] = tranmerc_fwd(0, 0, 20, 30, ell); 775 n = n + assertEquals(x * k0, 3266035.453860, 0.5e-6); 776 n = n + assertEquals(y * k0, 2518371.552676, 0.5e-6); 777 n = n + assertEquals(gam, 11.207356502141, 0.5e-12); 778 n = n + assertEquals(k * k0, 1.134138960741, 0.5e-12); 779 [lat, lon, gam, k] = tranmerc_inv(0, 0, 3.3e6 / k0, 2.5e6 / k0, ell); 780 n = n + assertEquals(lat, 19.80370996793, 0.5e-11); 781 n = n + assertEquals(lon, 30.24919702282, 0.5e-11); 782 n = n + assertEquals(gam, 11.214378172893, 0.5e-12); 783 n = n + assertEquals(k * k0, 1.137025775759, 0.5e-12); 784end 785 786function n = geodreckon0 787% Check mixed array size bugs 788 n = 0; 789 % lat1 is an array, azi1 is a scalar: 2015-08-10 790 [~, ~, ~, S12] = geodreckon([10 20], 0, 0, 0); 791 if length(S12) ~= 2, n = n+1; end 792 % scalar args except s12 is empty: 2017-03-26 793 [~, ~, ~, S12] = geodreckon(10, 0, [], 0); 794 if ~isempty(S12), n = n+1; end 795 % atan2dx needs to accommodate scalar + array arguments: 2017-03-27 796 lat2 = geodreckon(3, 4, [1, 2], 90); 797 if length(lat2) ~= 2, n = n+1; end 798end 799 800function n = gedistance0 801% gedistance(0, 0, 0, 100) wrongly return nan; 2015-09-23 802 n = 0; 803 s12 = gedistance(0, 0, 0, 100); 804 n = n + assertEquals(s12, 11131949, 0.5); 805end 806 807function n = tranmerc0 808% In 1.44, tranmerc_{fwd,inv} didn't work with array arguments. 809 n = 0; 810 % This used to result in an error 811 [x, y, gam, k] = tranmerc_fwd(0, 0, [90,90;85,85], [10,20;10,20]); 812 k0 = 0.9996; 813 n = n + assertEquals(x, [0, 0; 96820.412637, 190740.935334]/k0, 0.5e-6); 814 n = n + assertEquals(y, [9997964.943021, 9997964.943021; ... 815 9448171.516284, 9473242.646190]/k0, 0.5e-6); 816 n = n + assertEquals(gam, [10, 20; ... 817 9.962710901776, 19.929896900550], 0.5e-12); 818 n = n + assertEquals(k, [0.9996, 0.9996; ... 819 0.999714504947, 1.000044424775]/k0, 0.5e-12); 820 % This used to result in NaNs 821 [x, y, gam, k] = tranmerc_fwd(0, 0, [90,90;85,85]', [10,20;10,20]'); 822 k0 = 0.9996; 823 n = n + assertEquals(x, [0, 0; 96820.412637, 190740.935334]'/k0, 0.5e-6); 824 n = n + assertEquals(y, [9997964.943021, 9997964.943021; ... 825 9448171.516284, 9473242.646190]'/k0, 0.5e-6); 826 n = n + assertEquals(gam, [10, 20; ... 827 9.962710901776, 19.929896900550]', 0.5e-12); 828 n = n + assertEquals(k, [0.9996, 0.9996; ... 829 0.999714504947, 1.000044424775]'/k0, 0.5e-12); 830end 831 832function n = mgrs0 833% In 1.43, mgrs_inv didn't detect illegal letter combinations. 834 n = 0; 835 % This used to result in finite x, y 836 [x, y, zone, northp] = mgrs_inv('38RMB'); 837 n = n + assertNaN(x) + assertNaN(y); 838 n = n + assertEquals(zone, -4, 0); 839 n = n + assertEquals(northp, false, 0); 840end 841 842function n = mgrs1 843% In 1.44, mgrs_fwd gives the wrong results with prec = 10 or 11 in octave 844 n = 0; 845 % This used to result in '38SMB-1539607552-1539607552' 846 mgrs = mgrs_fwd(450000, 3650000, 38, 1, 11); 847 n = n + assertEquals(mgrs{1}, '38SMB5000000000050000000000', 0); 848end 849 850function n = mgrs2 851% In 1.43, mgrs_inv doesn't decode prec 11 string correctly 852 n = 0; 853 % This used to result in x = y = NaN 854 [x, y, zone, northp] = mgrs_inv('38SMB5000000000050000000000'); 855 n = n + assertEquals(x, 450000.0000005, 0.5e-6); 856 n = n + assertEquals(y, 3650000.0000005, 0.5e-6); 857 n = n + assertEquals(zone, 38, 0); 858 n = n + assertEquals(northp, true, 0); 859end 860 861function n = mgrs3 862% GeoConvert16: Check MGRS::Forward improved rounding fix, 2015-07-22 863 n = 0; 864 mgrs = mgrs_fwd(444140.6, 3684706.3, 38, 1, 8); 865 n = n + assertEquals(mgrs{1}, '38SMB4414060084706300', 0); 866end 867 868function n = mgrs4 869% GeoConvert17: Check MGRS::Forward digit consistency fix, 2015-07-23 870 n = 0; 871 mgrs = mgrs_fwd(500000, 63.811, 38, 1, 8); 872 n = n + assertEquals(mgrs{1}, '38NNF0000000000063811', 0); 873 mgrs = mgrs_fwd(500000, 63.811, 38, 1, 9); 874 n = n + assertEquals(mgrs{1}, '38NNF000000000000638110', 0); 875end 876 877function n = mgrs5 878% GeoConvert19: Check prec = -6 for UPS (to check fix to Matlab mgrs_fwd, 879% 2018-03-19) 880 n = 0; 881 mgrs = mgrs_fwd(2746000, 1515000, 0, 0, [-1,0,1]); 882 n = n + assertEquals(length(mgrs{1}), 1, 0); 883 n = n + assertEquals(mgrs{1}, 'B', 0); 884 n = n + assertEquals(length(mgrs{2}), 3, 0); 885 n = n + assertEquals(mgrs{2}, 'BKH', 0); 886 n = n + assertEquals(length(mgrs{3}), 5, 0); 887 n = n + assertEquals(mgrs{3}, 'BKH41', 0); 888end 889