text: /*****************************************************************************
  *                                                                            *
  *  ------------------------------- arclen.c -------------------------------  *
  *                                                                            *
  *****************************************************************************/
  
  #include <math.h>
  
  #include "geometry.h"
  
  /*****************************************************************************
  *                                                                            *
  *  -------------------------------- arclen --------------------------------  *
  *                                                                            *
  *****************************************************************************/
  
  void arclen(SPoint p1, SPoint p2, double *length) {
  
  Point              p1_rct,
                     p2_rct;
  
  double             alpha,
                     dot;
  
  /*****************************************************************************
  *                                                                            *
  *  Convert the spherical coordinates to rectilinear coordinates.             *
  *                                                                            *
  *****************************************************************************/
  
  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);
  
  /*****************************************************************************
  *                                                                            *
  *  Get the angle between the line segments from the origin to each point.    *
  *                                                                            *
  *****************************************************************************/
  
  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));
  
  /*****************************************************************************
  *                                                                            *
  *  Compute the length of the arc along the spherical surface.                *
  *                                                                            *
  *****************************************************************************/
  
  *length = alpha * p1.rho;
  
  return;
  
  }
  