#include <stdio.h>
#include <math.h>
#include <string.h>
#define RHO 6367
#define PI 3.141592654

float convert_DD(float *, float *, float *);
float calc_SphereCoord_lat( char *, float *);
float calc_SphereCoord_long(char *, float *);
float deg2rad(float *);
float straight_D(float *, float *, float *, float *, float *, float *);

struct coordinate
{
   float fDegLat;
   float fMinLat;
   float fSecLat;
   float fDegLon;
   float fMinLon;
   float fSecLon;
   char  szDirectionLat[10];
   char  szDirectionLon[10];
};

struct st3DCoord
{
   float x;
   float y;
   float z;
};

float calc_3DCoord(float phi, float theta, struct st3DCoord *st3dCoord )
{

   st3dCoord -> x = (RHO * (sin(phi)) * (cos(theta)));
   st3dCoord -> y = (RHO * (sin(phi)) * (sin(theta)));
   st3dCoord -> z = (RHO * (cos(phi)));

   return(0);
}

void read_GPS( char * szMessage, struct coordinate *coord )
{
   printf("Please enter in lattitudinal co-ordinates for %s:\n", szMessage);
   scanf("%f%f%f%s", &coord -> fDegLat,
                     &coord -> fMinLat,
                     &coord -> fSecLat,
                     coord -> szDirectionLat );

   printf("Please enter in longitudinal co-ordinates for %s:\n", szMessage);
   scanf("%f%f%f%s", &coord -> fDegLon,
                     &coord -> fMinLon,
                     &coord -> fSecLon,
                     coord -> szDirectionLon );
   return;
}

void get3DCoordinate( struct coordinate *coord,
                      struct st3DCoord *st3dcoord )
{
   float fDecDegLat;
   float fDecDegLon;
   float phi, thet;
   float tmp;

   fDecDegLat = convert_DD( &coord -> fDegLat,
                            &coord -> fMinLat,
                            &coord -> fSecLat );

   fDecDegLon = convert_DD( &coord -> fDegLon,
                            &coord -> fMinLon,
                            &coord -> fSecLon );

   tmp = calc_SphereCoord_lat( coord -> szDirectionLat,
                               &fDecDegLat );
   phi = deg2rad( &tmp );

   tmp =  calc_SphereCoord_long( coord -> szDirectionLon,
                                          &fDecDegLon );
   thet = deg2rad( &tmp );

   calc_3DCoord( phi, thet, st3dcoord );
}

float surfaceDistance( struct coordinate *start,
                       struct coordinate *finish )
{
   struct st3DCoord st3dstart;
   struct st3DCoord st3dfinish;
   float t1, t2, t3;
   float d, sd;
   float theta;

   get3DCoordinate( start, &st3dstart );
   get3DCoordinate( finish, &st3dfinish  );

   t1 = (st3dfinish.x - st3dstart.x) * (st3dfinish.x - st3dstart.x);
   t2 = (st3dfinish.y - st3dstart.y) * (st3dfinish.y - st3dstart.y);
   t3 = (st3dfinish.z - st3dstart.z) * (st3dfinish.z - st3dstart.z);

   d = sqrt( t1 + t2 + t3 );

//   sd = 2 * (sin( d / 2 ));
   theta = acos((RHO*RHO+RHO*RHO-d*d)/(2*RHO*RHO));
   sd = theta * RHO;

   return( sd );
}

int main(void)

{
   struct coordinate start;
   struct coordinate finish;
   float  fDistance = 0;

   printf("Please enter co-ordinates in format of 123 45 30 North\n" );
   printf("where 123 = Degrees 45 = minutes 30 = seconds and North = direction.\n");

   read_GPS( "start point", &start );
   read_GPS( "finish point", &finish );

   fDistance = surfaceDistance( &start, &finish );

   printf( "Surface distance is [%f]\n", fDistance );
   printf("%c", getchar());
   return(0);
}


float convert_DD(float *degrees, float *minutes, float *seconds)
{
   float dec_deg,totalmins;

   totalmins = *minutes + (*seconds / 60);
   dec_deg = (totalmins /60) + *degrees;
   printf("The decimal degree is %f", dec_deg);

   return(dec_deg);

}

float calc_SphereCoord_lat(char *direction, float *dec_deg)
{
   int temp_1;
   float phi;
   char s5[] = "North";
   temp_1 = strcasecmp(direction,s5);
   if (temp_1 == 0)
      phi = 90 - *dec_deg;
   else
      phi = 90 + *dec_deg;

   printf("The spherical co-ordinate phi is: %f\n", phi);

   return(phi);
}

float calc_SphereCoord_long(char *direction, float *dec_deg)
{
   int temp;
   float theta;
   char s6[] = "East";
   temp = strcasecmp(direction,s6);
   if (temp == 0)
      theta = *dec_deg;
   else
      theta = (-*dec_deg);

   printf("The spherical co-ordinate theta is: %f\n", theta);

   return(theta);
}

float deg2rad(float *dec_deg)
{
   float radians;
   radians = (*dec_deg * 2 * PI) / 360;
   printf("%f as radians is: %f\n", *dec_deg, radians);
   return(radians);
}