antipodes.cpp (2861B)
1 #include<cstdio> 2 #include<iostream> 3 #include<string> 4 #include<vector> 5 #include<utility> 6 #include<cmath> 7 8 /* 9 * How to use: 10 * take airports.dat from http://openflights.org/data.html and do: 11 * rev airports.dat | cut -d, -f5,6,7,8 | rev | tr -d '"' | 12 * tr ',' ' ' | sed // 's/ /-/' | ./antipodes 13 * (the use of rev is to avoid problems with ',' in earlier fields) 14 * results are in euclidean (pairs with euclidean distance in tore) 15 * and haversine (distance in m according to haversine formula) 16 * WARNING: output files are *not* sorted 17 */ 18 19 using namespace std; 20 21 vector<string> a; 22 vector<pair<double, double> > p; 23 vector<double> h; 24 25 #define PI 3.14159265358979 26 #define EPS 0.000000000001 27 // https://en.wikipedia.org/wiki/Earth_radius#Mean_radius 28 #define R 6371009 29 #define THRESH 19800000. 30 31 double euclidean(double la1, double lo1, double la2, double lo2) { 32 double lo2b = lo2+180.; 33 if (lo2b > 180.) 34 lo2b -= 360; 35 double la2b = -la2; 36 double dlo = lo1 - lo2b; 37 double dla = la1 - la2b; 38 dlo *= dlo; 39 dla *= dla; 40 return R * PI * (1 - sqrt(dlo + dla)/360); 41 } 42 43 // https://en.wikipedia.org/wiki/Haversine_formula#The_haversine_formula 44 double haversine(double la1, double lo1, double la2, double lo2) { 45 // degrees to radians 46 lo1 *= PI/180; 47 lo2 *= PI/180; 48 la1 *= PI/180; 49 la2 *= PI/180; 50 51 double tmp1 = sin(.5*(lo2 - lo1)); 52 tmp1 *= tmp1; 53 double tmp2 = sin(.5*(la2 - la1)); 54 tmp2 *= tmp2; 55 return R * 2. * asin(sqrt(tmp2 + cos(la1) * cos(la2) * tmp1)); 56 } 57 58 int main() { 59 while (1) { 60 char s[10]; 61 int ret = scanf("%s", s); 62 if (ret != 1) 63 break; 64 double la, lo, hh; 65 scanf("%lf%lf%lf", &la, &lo, &hh); 66 string s2(s); 67 a.push_back(s2); 68 h.push_back(hh); 69 p.push_back(make_pair(la, lo)); 70 //printf("%ld %lf %lf\n", p.size()-1, p[p.size()-1].first, p[p.size()-1].second); 71 } 72 73 FILE *feuclidean = fopen("euclidean", "w"); 74 FILE *fhaversine = fopen("haversine", "w"); 75 FILE *fellipsoid = fopen("ellipsoid_in", "w"); 76 77 for (unsigned int i = 0; i < p.size(); i++) { 78 if (!(i % 100)) 79 fprintf(stderr, "%d of %ld\n", i, p.size()); 80 for (unsigned int j = i+1; j < p.size(); j++) { 81 double de = euclidean(p[i].first, p[i].second, p[j].first, p[j].second); 82 if (de > THRESH) { 83 fprintf(feuclidean, "%d-%s %d-%s %lf\n", i, a[i].c_str(), j, a[j].c_str(), de); 84 } 85 double d = haversine(p[i].first, p[i].second, p[j].first, p[j].second); 86 if (d > THRESH) { 87 printf("%lf\n", abs(d - de)); 88 fprintf(fhaversine, "%d-%s %d-%s %lf\n", i, a[i].c_str(), j, a[j].c_str(), d); 89 fprintf(fellipsoid, "%d-%s %d-%s %lf %lf %lf %lf %lf %lf %lf\n", i, a[i].c_str(), j, a[j].c_str(), d, 90 p[i].first, p[i].second, h[i], p[j].first, p[j].second, h[j]); 91 } 92 } 93 } 94 95 fclose(fhaversine); 96 fclose(fellipsoid); 97 fclose(feuclidean); 98 99 return 0; 100 } 101