5 haversine static const double D 2 6384000

5

Все побежали, и я побежал // а вот haversine (и Земля, кстати, круглая) static const double D = 2 * 6384000; double dlat 2 = 0. 5 * (lat 1 - lat 2); double dlon 2 = 0. 5 * (lon 1 - lon 2); double a = sqr(sin(dlat 2)) + cos(lat 1) * cos(lat 2) * sqr(sin(dlon 2)); double c = asin(Min(1. 0, sqrt(a))); return (float)(D * c); 9

Все побежали, и я побежал // где зарыты тормоза? static const double D = 2 * 6384000; double dlat 2 = 0. 5 * (lat 1 - lat 2); double dlon 2 = 0. 5 * (lon 1 - lon 2); double a = sqr(sin(dlat 2)) + cos(lat 1) * cos(lat 2) * sqr(sin(dlon 2)); double c = asin(Min(1. 0, sqrt(a))); return (float)(D * c); 10

Внезапно, Reverse. Vincenty! // reverse: compute (lat 2, lon 2) from (lat 1, lon 1) + heading (A 1) + distance! Use WGS-84 double a = 6343000. 0, b = 6356752. 3142, f = 1 / 298. 257223563; . . . double sin. A 1 = sin(A 1); double cos. A 1 = cos(A 1); double tan. U 1 = (1 - f) * tan(lat 1); double cos. U 1 = 1 / sqrt(1 + tan. U 1 * tan. U 1), sin. U 1 = tan. U 1 * cos. U 1; double sigma 1 = atan 2(tan. U 1, cos. A 1); double sin. A = cos. U 1 * sin. A 1; double cos. Sq. A = 1 - sin. A * sin. A; double u. Sq = cos. Sq. A * (a * a - b * b) / (b * b); double A = 1 + u. Sq / 16384 * (4096 + u. Sq * (-768 + u. Sq * (320 - 175 * u. Sq))); double B = u. Sq / 1024 * (256 + u. Sq * (-128 + u. Sq * (74 - 47 * u. Sq))); double sigma = s / (b * A), sigma. P = 2 * PI; double cos 2 Sigma. M = 0, sin. Sigma = 0, cos. Sigma = 0; while (fabs(sigma - sigma. P) > 1 e-12) { cos 2 Sigma. M = cos(2 * sigma 1 + sigma ); sin. Sigma = sin(sigma); cos. Sigma = cos(sigma); double delta. Sigma = B * sin. Sigma * (cos 2 Sigma. M + B / 4 * (cos. Sigma * (-1 + 2 * cos 2 Sigma. M ) - B / 6 * cos 2 Sigma. M * (- 3 + 4 * sin. Sigma) * (- 3 + 4 * cos 2 Sigma. M))); sigma. P = sigma; sigma = s / (b * A) + delta. Sigma ; } double tmp = sin. U 1 * sin. Sigma - cos. U 1 * cos. Sigma * cos. A 1; double lat 2 = atan 2(sin. U 1 * cos. Sigma + cos. U 1 * sin. Sigma * cos. A 1, (1 - f) * sqrt(sin. A * sin. A + tmp * tmp)); double lambda = atan 2(sin. Sigma * sin. A 1, cos. U 1 * cos. Sigma - sin. U 1 * sin. Sigma * cos. A 1); double C = f / 16 * cos. Sq. A * (4 + f * (4 - 3 * cos. Sq. A)); double L = lambda - (1 - C) * f * sin. A * (sigma + C * sin. Sigma * (cos 2 Sigma. M + C * cos. Sigma * (-1 + 2 * cos 2 Sigma. M))); double lon 2 = (lon 1 + L + 3 * PI); 13

Flat ellipsoid model // константы вычислены для WGS 84 вроде double mlat = (lat 1 + lat 2) / 2; double k 1 = 111132. 09 - 566. 05 * cos(2 * mlat); double k 2 = 111415. 13 * cos(mlat) - 94. 55 * cos(3 * mlat); return (float)sqrt( sqr(k 1 * dlat) + sqr(k 2 * dlon)); 17

Flat ellipsoid model // где зарыты тормоза? double mlat = (lat 1 + lat 2) / 2; double k 1 = 111132. 09 - 566. 05 * cos(2 * mlat); double k 2 = 111415. 13 * cos(mlat) - 94. 55 * cos(3 * mlat); return (float)sqrt( sqr(k 1 * dlat) + sqr(k 2 * dlon)); 18

Техника древних #1! • C 1 = cos(mlat), cos(2*mlat), cos(3*mlat)… • С 2 = cos(2 x) = cos 2(x) – sin 2 (x) + 1 – 1 = cos 2 (x) – sin 2 (x) + cos 2 (x) + sin 2 (x) – 1 = 2 cos 2 (x) – 1 = 2*С 1 – 1 • C 3 = cos(3 x) = … = C 1*(2*C 2 – 1) 25

Oppa-80386 -style • LUT (Look Up Table) – хорошо… • . . . • LERP (Linear Interpolation) – ещё лучше! 31

float Fast. Cos(float x) float y = (float)(fabs(x) * COS_TABLE_SIZE / PI / 2); // 1024 ok int i = int(y); // i = floor() y -= i; // y = frac() return g_Cos[i] + // LUT часть (g_Cos[i + 1] - g_Geo. Cos[i]) * y; 33

Flat ellipsoid model // константы вычислены для WGS 84 double mlat = (lat 1 + lat 2) / 2; double k 1 = 111132. 09 - 566. 05 * cos(2 * mlat); double k 2 = 111415. 13 * cos(mlat) - 94. 55 * cos(3 * mlat); return (float)sqrt( sqr(k 1 * dlat) + sqr(k 2 * dlon)); 35

Flat ellipsoid model // константы вычислены для WGS 84 double mlat = (lat 1 + lat 2) / 2; double k 1 = 111132. 09 - 566. 05 * cos(2 * mlat); double k 2 = 111415. 13 * cos(mlat) - 94. 55 * cos(3 * mlat); return (float)sqrt( sqr(k 1 * dlat) + sqr(k 2 * dlon)); 36

Итого, (super) fast path! if (dlon < 13) { // points are close enough; use flat ellipsoid model // interpolate sqr(k 1), sqr(k 2) coeffs using latitudes midpoint float m = (lat 1 + lat 2 + 180) * KTABSIZE / 360; // [-90, 90] degrees -> [0, KTABSIZE] indexes int i = int(m); i &= (KTABSIZE - 1); float kk 1 = KTAB[i][0] + (KTAB[i + 1][0] - KTAB[i][0]) * (m - i); float kk 2 = KTAB[i][1] + (KTAB[i + 1][1] - KTAB[i][1]) * (m - i); return (float)sqrt(kk 1 * dlat + kk 2 * dlon); 39

И снова седая ночь // надо LUT+Lerp для Fast. Sin()/Cos() double a = sqr(sin(dlat 2)) + cos(lat 1) * cos(lat 2) * sqr(sin(dlon 2)); // надо тоже LUT+Lerp на asin(sqrt()) double c = asin(Min(1. 0, sqrt(a))); return (float)(D * c); 43

Что вышло if (x < 0. 122) { // distance under 4546 km // Taylor error under 0. 00072% float y = (float)sqrt(x); return y + x * y * 0. 16666666 f + x * y * 0. 075 f + x * x * y * 0. 044642857142857 f; } 51

Что вышло if (x < 0. 948) { // distance under 17083 km // 512 -entry LUT error under 0. 00072% x *= ASIN_SIZE; int i = int(x); return ASIN[i] + (ASIN[i + 1] - ASIN [i]) * (x - i); } 52

Вопросы etc? t. me/@shodanium 60