bigworld

compute diameter of flights graph
git clone https://a3nm.net/git/bigworld/
Log | Files | Refs

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