// chapter4_4.c - great circle distance, modified to use a 3D Point structure
//
#include <stdio.h>
#include <math.h>

#define PI (4*atan(1))

#define rho 3960 // radius of the earth

// a point in three dimensions
//
typedef struct { double x, y, z; } Point;

// compute x,y,z coordinates from latitude and longitude
//
Point xyz( double lat1, double long1)
{
   Point p;
   double phi = (90 - lat1)*(PI/180.0);
   double theta = (360 - long1)*(PI/180.0);
   p.x = rho*sin(phi)*cos(theta);
   p.y = rho*sin(phi)*sin(theta);
   p.z = rho*cos(phi);
   return p;
}

// inner product between two points
//
double prod( Point p, Point q)
{
  return p.x*q.x + p.y*q.y + p.z*q.z;
}

// great circle distance
//
double gc_distance(double lat1,double long1, double lat2,double long2)
{
   double  gamma, dot, dist1, dist2;
   Point p = xyz( lat1, long1);
   Point q = xyz( lat2, long2);
   dot = prod( p, q);
   dist1 = sqrt(prod(p,p)); 
   dist2 = sqrt(prod(q,q));
   gamma = acos(dot/(dist1*dist2));
   return gamma*rho;
}

int main(void)
{
  // New York, London, distance should be 3466 miles
  //
   double lat1 = 40.75, long1 = 74, lat2 = 51.5, long2 = 0;
   
   printf("Great Circle Distance: %.0f miles \n", gc_distance(lat1,long1,lat2,long2));

   return 0;
}