Code:

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

#include "geometry.h"
#include "list.h"

#define EARTH_RADIUS 3440.065



void arclen(SPoint p1, SPoint p2, double *length) {

Point              p1_rct,
                   p2_rct;

double             alpha,
                   dot;



p1_rct.x = p1.rho * sin(p1.phi) * cos(p1.theta);
p1_rct.y = p1.rho * sin(p1.phi) * sin(p1.theta);
p1_rct.z = p1.rho * cos(p1.phi);

p2_rct.x = p2.rho * sin(p2.phi) * cos(p2.theta);
p2_rct.y = p2.rho * sin(p2.phi) * sin(p2.theta);
p2_rct.z = p2.rho * cos(p2.phi);



dot = (p1_rct.x * p2_rct.x) + (p1_rct.y * p2_rct.y) + (p1_rct.z * p2_rct.z);
alpha = acos(dot / pow(p1.rho, 2.0));



*length = alpha * p1.rho;

return;

}

void list_init(List *list, void (*destroy)(void *data)) {



list->size = 0;
list->destroy = destroy;
list->head = NULL;
list->tail = NULL;

return;

}



void list_destroy(List *list) {

void               *data;



while (list_size(list) > 0) {

   if (list_rem_next(list, NULL, (void **)&data) == 0 && list->destroy !=
      NULL) {

      list->destroy(data);

   }

}



memset(list, 0, sizeof(List));

return;

}



int list_ins_next(List *list, ListElmt *element, const void *data) {

ListElmt           *new_element;



if ((new_element = (ListElmt *)malloc(sizeof(ListElmt))) == NULL)
   return -1;



new_element->data = (void *)data;

if (element == NULL) {



   if (list_size(list) == 0)
      list->tail = new_element;

   new_element->next = list->head;
   list->head = new_element;

   }

else {



   if (element->next == NULL)
      list->tail = new_element;

   new_element->next = element->next;
   element->next = new_element;

}



list->size++;

return 0;

}


int list_rem_next(List *list, ListElmt *element, void **data) {

ListElmt           *old_element;



if (list_size(list) == 0)
   return -1;



if (element == NULL) {



   *data = list->head->data;
   old_element = list->head;
   list->head = list->head->next;

   if (list_size(list) == 1)
      list->tail = NULL;

   }

else {



   if (element->next == NULL)
      return -1;

   *data = element->next->data;
   old_element = element->next;
   element->next = element->next->next;

   if (element->next == NULL)
      list->tail = element;

}



free(old_element);



list->size--;

return 0;

}

int cvxhull(const List *P, List *polygon) {

ListElmt           *element;

Point              *min,
                   *low,
                   *p0,
                   *pi,
                   *pc;

double             z,
                   length1,
                   length2;

int                count;


min =(Point*) list_data(list_head(P));

for (element = list_head(P); element != NULL; element = list_next(element)) {

   p0 =(Point*) list_data(element);


   if (p0->y < min->y) {

      min = p0;
      low =(Point*) list_data(element);

      }

   else {


      if (p0->y == min->y && p0->x < min->x) {

         min = p0;
         low =(Point*) list_data(element);

      }

   }

}

list_init(polygon, NULL);


p0 = low;

do {



   if (list_ins_next(polygon, list_tail(polygon), p0) != 0) {

      list_destroy(polygon);
      return -1;

   }

   count = 0;

   for (element = list_head(P); element != NULL; element =
      list_next(element)) {


      if ((pi =(Point*) list_data(element)) == p0)
         continue;

      count++;


      if (count == 1) {

         pc =(Point*) list_data(element);
         continue;

      }


      if ((z = ((pi->x - p0->x) * (pc->y - p0->y)) - ((pi->y - p0->y) * (pc->x
         - p0->x))) > 0) {



         pc = pi;

         }

      else if (z == 0) {



         length1 = sqrt(pow(pi->x - p0->x, 2.0) + pow(pi->y - p0->y, 2.0));
         length2 = sqrt(pow(pc->x - p0->x, 2.0) + pow(pc->y - p0->y, 2.0));

         if (length1 > length2) {



            pc = pi;

         }

      }

   }



   p0 = pc;


} while (p0 != low);

return 0;

}


/*****************************************************************************
*                                                                            *
*  Define data for computing convex hulls.                                   *
*                                                                            *
*****************************************************************************/

#define            CVXPCT               10

static Point CvxTestP[CVXPCT] = {

   { 0.0,  1.0, 0.0},
   {-3.0,  5.0, 0.0},
   {-2.0, -3.0, 0.0},
   { 0.0,  0.0, 0.0},
   { 3.0,  2.0, 0.0},
   {-5.0, -2.0, 0.0},
   { 5.0,  3.0, 0.0},
   {-3.0,  3.0, 0.0},
   { 2.0,  3.0, 0.0},
   {-3.0, -1.0, 0.0}

};


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

List               P,
                   polygon;

ListElmt           *element;

Point              p1_rct,
                   p2_rct,
                   p3_rct,
                   p4_rct,
                   *point;

SPoint             p1_sph,
                   p2_sph;

double             length;

int                actual,
                   i;


/*
fprintf(stdout, "Determining whether two line segments intersect\n");

p1_rct.x = 2.0;
p1_rct.y = -1.0;
p1_rct.z = 0.0;

p2_rct.x = 4.0;
p2_rct.y = -1.0;
p2_rct.z = 0.0;

p3_rct.x = 1.0;
p3_rct.y = -2.0;
p3_rct.z = 0.0;

p4_rct.x = 3.0;
p4_rct.y = 2.0;
p4_rct.z = 0.0;

if (lint(p1_rct, p2_rct, p3_rct, p4_rct)) {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) Y", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

   }

else {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) N", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

}

fprintf(stdout, " (N=OK)\n");

p1_rct.x = -4.0;
p1_rct.y = -3.0;
p1_rct.z = 0.0;

p2_rct.x = -2.0;
p2_rct.y = -1.0;
p2_rct.z = 0.0;

p3_rct.x = -1.0;
p3_rct.y = 0.0;
p3_rct.z = 0.0;

p4_rct.x = 1.0;
p4_rct.y = 2.0;
p4_rct.z = 0.0;

if (lint(p1_rct, p2_rct, p3_rct, p4_rct)) {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) Y", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

   }

else {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) N", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

}

fprintf(stdout, " (N=OK)\n");

p1_rct.x = -4.0;
p1_rct.y = -3.0;
p1_rct.z = 0.0;

p2_rct.x = -2.0;
p2_rct.y = -1.0;
p2_rct.z = 0.0;

p3_rct.x = -4.0;
p3_rct.y = -1.0;
p3_rct.z = 0.0;

p4_rct.x = -3.0;
p4_rct.y = -2.0;
p4_rct.z = 0.0;

if (lint(p1_rct, p2_rct, p3_rct, p4_rct)) {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) Y", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

   }

else {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) N", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

}

fprintf(stdout, " (Y=OK)\n");

p1_rct.x = -4.0;
p1_rct.y = -3.0;
p1_rct.z = 0.0;

p2_rct.x = -2.0;
p2_rct.y = -1.0;
p2_rct.z = 0.0;

p3_rct.x = -3.0;
p3_rct.y = -2.0;
p3_rct.z = 0.0;

p4_rct.x = -3.0;
p4_rct.y = -2.0;
p4_rct.z = 0.0;

if (lint(p1_rct, p2_rct, p3_rct, p4_rct)) {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) Y", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

   }

else {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) N", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

}

fprintf(stdout, " (Y=OK)\n");

p1_rct.x = -3.0;
p1_rct.y = -2.0;
p1_rct.z = 0.0;

p2_rct.x = -3.0;
p2_rct.y = -2.0;
p2_rct.z = 0.0;

p3_rct.x = -3.0;
p3_rct.y = -2.0;
p3_rct.z = 0.0;

p4_rct.x = -3.0;
p4_rct.y = -2.0;
p4_rct.z = 0.0;

if (lint(p1_rct, p2_rct, p3_rct, p4_rct)) {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) Y", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

   }

else {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) N", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

}

fprintf(stdout, " (Y=OK)\n");

p1_rct.x = -3.0;
p1_rct.y = -2.0;
p1_rct.z = 0.0;

p2_rct.x = 3.0;
p2_rct.y = 4.0;
p2_rct.z = 0.0;

p3_rct.x = -1.0;
p3_rct.y = 2.0;
p3_rct.z = 0.0;

p4_rct.x = 3.0;
p4_rct.y = -1.0;
p4_rct.z = 0.0;

if (lint(p1_rct, p2_rct, p3_rct, p4_rct)) {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) Y", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

   }

else {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) N", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

}

fprintf(stdout, " (Y=OK)\n");

p1_rct.x = -3.0;
p1_rct.y = -2.0;
p1_rct.z = 0.0;

p2_rct.x = 3.0;
p2_rct.y = 4.0;
p2_rct.z = 0.0;

p3_rct.x = -1.0;
p3_rct.y = -2.0;
p3_rct.z = 0.0;

p4_rct.x = 3.0;
p4_rct.y = 1.0;
p4_rct.z = 0.0;

if (lint(p1_rct, p2_rct, p3_rct, p4_rct)) {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) Y", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

   }

else {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) N", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

}

fprintf(stdout, " (N=OK)\n");

p1_rct.x = -3.0;
p1_rct.y = -2.0;
p1_rct.z = 0.0;

p2_rct.x = 3.0;
p2_rct.y = 4.0;
p2_rct.z = 0.0;

p3_rct.x = -3.0;
p3_rct.y = -2.0;
p3_rct.z = 0.0;

p4_rct.x = -4.0;
p4_rct.y = 3.0;
p4_rct.z = 0.0;

if (lint(p1_rct, p2_rct, p3_rct, p4_rct)) {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) Y", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

   }

else {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) N", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

}

fprintf(stdout, " (Y=OK)\n");

p1_rct.x = -3.0;
p1_rct.y = -1.0;
p1_rct.z = 0.0;

p2_rct.x = 3.0;
p2_rct.y = 4.0;
p2_rct.z = 0.0;

p3_rct.x = -3.0;
p3_rct.y = -2.0;
p3_rct.z = 0.0;

p4_rct.x = -4.0;
p4_rct.y = 3.0;
p4_rct.z = 0.0;

if (lint(p1_rct, p2_rct, p3_rct, p4_rct)) {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) Y", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

   }

else {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) N", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

}

fprintf(stdout, " (N=OK)\n");

p1_rct.x = -3.0;
p1_rct.y = -2.0;
p1_rct.z = 0.0;

p2_rct.x = 3.0;
p2_rct.y = 4.0;
p2_rct.z = 0.0;

p3_rct.x = -3.0;
p3_rct.y = -1.0;
p3_rct.z = 0.0;

p4_rct.x = -4.0;
p4_rct.y = 3.0;
p4_rct.z = 0.0;

if (lint(p1_rct, p2_rct, p3_rct, p4_rct)) {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) Y", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

   }

else {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) N", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

}

fprintf(stdout, " (N=OK)\n");

p1_rct.x = -3.0;
p1_rct.y = 0.0;
p1_rct.z = 0.0;

p2_rct.x = -4.0;
p2_rct.y = 3.0;
p2_rct.z = 0.0;

p3_rct.x = -5.0;
p3_rct.y = 6.0;
p3_rct.z = 0.0;

p4_rct.x = -6.0;
p4_rct.y = 9.0;
p4_rct.z = 0.0;

if (lint(p1_rct, p2_rct, p3_rct, p4_rct)) {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) Y", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

   }

else {

   fprintf(stdout, "(%+.1lf,%+.1lf) to (%+.1lf,%+.1lf) and (%+.1lf,%+.1lf) "
      "to (%+.1lf,%+.1lf) N", p1_rct.x, p1_rct.y, p2_rct.x, p2_rct.y,
      p3_rct.x, p3_rct.y, p4_rct.x, p4_rct.y);

}
*/
fprintf(stdout, " (N=OK)\n");

fprintf(stdout, "Computing a convex hull\n");
fprintf(stdout, "Points in P\n");

list_init(&P, NULL);

for (i = 0; i < CVXPCT; i++) {

   if (list_ins_next(&P, list_tail(&P), &CvxTestP[i]) != 0) {

      list_destroy(&P);
      return 1;

   }

   fprintf(stdout, "P[%03d]=(%+.1lf,%+.1lf,%+.1lf)\n", i, CvxTestP[i].x,
      CvxTestP[i].y, CvxTestP[i].z);

}

if (cvxhull(&P, &polygon) != 0) {

   list_destroy(&P);
   return 1;

}

fprintf(stdout, "Points in the convex hull\n");

i = 0;

for (element = list_head(&polygon); element != NULL; element =
   list_next(element)) {

   i++;
   point =(Point*) list_data(element);

   fprintf(stdout, "polygon[%03d]=(%+.1lf,%+.1lf,%+.1lf)\n", i, point->x,
      point->y, point->z);

}

list_destroy(&P);
list_destroy(&polygon);

fprintf(stdout, "Computing arc lengths on spherical surfaces\n");

p1_sph.phi = DEGTORAD(90.0);
p1_sph.theta = DEGTORAD(0.0);
p1_sph.rho = EARTH_RADIUS;

p2_sph.phi = DEGTORAD(90.0);
p2_sph.theta = DEGTORAD(60.0);
p2_sph.rho = EARTH_RADIUS;

arclen(p1_sph, p2_sph, &length);
actual = 3606;

fprintf(stdout, "Simple: phi=%+07.2lf, theta=%+07.2lf, rho=%.3lf\n",
   RADTODEG(p1_sph.phi), RADTODEG(p1_sph.theta), p1_sph.rho);

fprintf(stdout, "Simple: phi=%+07.2lf, theta=%+07.2lf, rho=%.3lf\n",
   RADTODEG(p2_sph.phi), RADTODEG(p2_sph.theta), p2_sph.rho);

fprintf(stdout, "length=%d, actual=%d, error=%6.4lf\n", (int)length,
   actual, fabs(1.0 - (floor(length) / actual)));

/* SFO (San Francisco) */
p1_sph.phi = DEGTORAD(90 - 37.62);
p1_sph.theta = DEGTORAD(-122.38);
p1_sph.rho = EARTH_RADIUS;

/* LAX (Los Angeles) */
p2_sph.phi = DEGTORAD(90.0 - 33.94);
p2_sph.theta = DEGTORAD(-118.41);
p2_sph.rho = EARTH_RADIUS;

arclen(p1_sph, p2_sph, &length);
actual = 293;

fprintf(stdout, "SFO: phi=%+07.2lf, theta=%+07.2lf, rho=%.3lf\n",
   RADTODEG(p1_sph.phi), RADTODEG(p1_sph.theta), p1_sph.rho);

fprintf(stdout, "LAX: phi=%+07.2lf, theta=%+07.2lf, rho=%.3lf\n",
   RADTODEG(p2_sph.phi), RADTODEG(p2_sph.theta), p2_sph.rho);

fprintf(stdout, "length=%d, actual=%d, error=%6.4lf\n", (int)length,
   actual, fabs(1.0 - (floor(length) / actual)));

/* SFO (San Francisco) */
p1_sph.phi = DEGTORAD(90 - 37.62);
p1_sph.theta = DEGTORAD(-122.38);
p1_sph.rho = EARTH_RADIUS;

/* HKG (Hong Kong) */
p2_sph.phi = DEGTORAD(90.0 - 22.31);
p2_sph.theta = DEGTORAD(113.92);
p2_sph.rho = EARTH_RADIUS;

arclen(p1_sph, p2_sph, &length);
actual = 6159;

fprintf(stdout, "SFO: phi=%+07.2lf, theta=%+07.2lf, rho=%.3lf\n",
   RADTODEG(p1_sph.phi), RADTODEG(p1_sph.theta), p1_sph.rho);

fprintf(stdout, "HKG: phi=%+07.2lf, theta=%+07.2lf, rho=%.3lf\n",
   RADTODEG(p2_sph.phi), RADTODEG(p2_sph.theta), p2_sph.rho);

fprintf(stdout, "length=%d, actual=%d, error=%6.4lf\n", (int)length,
   actual, fabs(1.0 - (floor(length) / actual)));

/* CDG (Paris) */
p1_sph.phi = DEGTORAD(90 - 49.01);
p1_sph.theta = DEGTORAD(2.55);
p1_sph.rho = EARTH_RADIUS;

/* PER (Perth) */
p2_sph.phi = DEGTORAD(90.0 + 31.94);
p2_sph.theta = DEGTORAD(115.97);
p2_sph.rho = EARTH_RADIUS;

arclen(p1_sph, p2_sph, &length);
actual = 7733;

fprintf(stdout, "CDG: phi=%+07.2lf, theta=%+07.2lf, rho=%.3lf\n",
   RADTODEG(p1_sph.phi), RADTODEG(p1_sph.theta), p1_sph.rho);

fprintf(stdout, "PER: phi=%+07.2lf, theta=%+07.2lf, rho=%.3lf\n",
   RADTODEG(p2_sph.phi), RADTODEG(p2_sph.theta), p2_sph.rho);

fprintf(stdout, "length=%d, actual=%d, error=%6.4lf\n", (int)length,
   actual, fabs(1.0 - (floor(length) / actual)));

return 0;

}