#include <math.h>
#include <stdio.h>
#include <time.h>
#include <stdlib.h>

#define NITER 10 * 1000

#define EARTH_RADIUS 3438.461

double distance(double lat1degrees, double lon1degrees, double lat2degrees, double lon2degrees)
{
  double lat1;
  double lon1;
  double lat2;
  double lon2;

  double delta_lon;
  double sin_delta_lon;
  double cos_delta_lon;
  double cos_lat2_times_sin_delta_lon;
  double sin_lat1_times_sin_lat2;
  double cos_lat1_times_sin_lat2;
  double sin_lat1_times_cos_lat2;
  double cos_lat1_times_cos_lat2;

  double factor1;
  double factor2;
  double y;
  double x;
  double res;

  lat1 = lat1degrees * M_PI / 180.0;
  lon1 = lon1degrees * M_PI / 180.0;
  lat2 = lat2degrees * M_PI / 180.0;
  lon2 = lon2degrees * M_PI / 180.0;

  delta_lon = lon2 - lon1;
  sin_delta_lon = sin(delta_lon);
  cos_delta_lon = cos(delta_lon);
  cos_lat2_times_sin_delta_lon = cos(lat2) * sin(delta_lon);
  sin_lat1_times_sin_lat2 = sin(lat1) * sin(lat2);
  cos_lat1_times_sin_lat2 = cos(lat1) * sin(lat2);
  sin_lat1_times_cos_lat2 = sin(lat1) * cos(lat2);
  cos_lat1_times_cos_lat2 = cos(lat1) * cos(lat2);

  factor1 = cos_lat2_times_sin_delta_lon;
  factor1 *= factor1;
  factor2 = cos_lat1_times_sin_lat2 - sin_lat1_times_cos_lat2 * cos_delta_lon;
  factor2 *= factor2;
  y = sqrt(factor1 + factor2);
  x = sin_lat1_times_sin_lat2 + cos_lat1_times_cos_lat2 * cos_delta_lon;
  res = EARTH_RADIUS * atan2(y, x);
    
  return res;
}

#include <signal.h>

quit(int sig)
{
	exit(1);
}

main(int argc, char * argv[]){

  float ressin;
  float rescos;
  float reseuclid;
  int iter;
  int samp;
  float arg;

  clock_t t1;
  clock_t t2;

  /* Quit after ten seconds */
  signal(SIGALRM, quit);
  alarm(10);

  int niter = NITER;
  if (argc > 1) {
    niter = atoi(argv[1]);
  }
  t1 = clock();
  for (iter = 0; iter < niter; iter++) {
    distance(36.12, -86.67, 33.94, -118.4);
  }
  t2 = clock();
  printf("That took %f seconds for %d iterations\n", ((double)(t2 - t1))/CLOCKS_PER_SEC, niter);

  printf("distance from BNA to LAX: %.4f\n", distance(36.12, -86.67, 33.94, -118.4));

  return 0;
}