$TITLE GREAT CIRCLE DISTANCES (GREAT,SEQ=73)
*  THE COORDINATES OF POINTS (AIRPORTS) ON THE GLOBE ARE GIVEN IN DEGREES
*  LATITUDE AND LONGITUDE. THE SHORTEST DISTANCE (GREAT CIRCLE DISTANCE)
*  BETWEEN PAIRS OF POINTS IS DESIRED. FIRST, THE SPHERIC COORDINATES ARE
*  TRANSLATED INTO CARTESIAN COORDINATES. SECOND, THE STRAIGHT LINE
*  DISTANCE BETWEEN POINTS ON THE UNIT SPHERE ARE CALACULATED. THIRD,
*  THE GREAT CIRCLE DISTANCES ARE COMPUTED.
*
*  THE CENTER OF THE EARTH IS THE ORIGIN FOR ALL COORDINATE SYSTEMS.
*
*  SPHERIC COORDINATES   LATITUDE ANGLE  NORTH  POSITIVE
*                                        SOUTH  NEGATIVE
*                        LONGITUDE ANGLE EAST   POSITIVE
*                                        WEST   NEGATIVE
*
*  CARTESIAN COORDINATES  X-AXIS   0 N   0 E
*                         Y-AXIS   0 N  90 E
*                         Z-AXIS  90 N
*
*  REFERENCE: BROOKE A, KENDRICK D AND MEERAUS A, GAMS: A USER'S GUIDE,
*             THE SCIENTIFIC PRESS, REDWOOD CITY, CALIFORNIA, 1988.

SETS  K    COORDINATES  / X-AXIS, Y-AXIS, Z-AXIS /

      A    AIRPORTS  / SFO   SAN FRANCISCO
                       MIA   MIAMI
                       JFK   NEW YORK
                       IAH   HOUSTON
                       IAD   WASHINGTON DC
                       KHI   KARACHI - PAKISTAN
                       NNN   NORTH POLE
                       SSS   SOUTH POLE  /
ALIAS (A,AP)

        TABLE  LOC(A,*)  LOACTION ON MAP

       LAT-DEG  LAT-MIN    LONG-DEG  LONG-MIN
SFO       37      37         -122      -23
MIA       25      48         - 80      -17
JFK       40      38         - 73      -47
IAH       29      58         - 95      -20
IAD       38      57         - 77      -25
KHI       24      40           67       10
NNN       90
SSS      -90

SCALAR PI TRIGONOMETRIC CONSTANT  / 3.141592653 /
       R  RADIUS OF EARTH (MILES) / 3959 /

PARAMETERS LAT(A)     LATITUDE ANGLE (RADIANS)
           LONG(A)    LONGITUDE ANGLE (RADIANS)
           UK(A,K)    POINT IN CARTESIAN COORDINATES (UNIT SPHERE)
           USEG(A,AP) STRAIGHT LINE DISTANCE BETWEEN POINTS (UNIT SPHERE)
           UDIS(A,AP) GREAT CIRCLE DISTANCES (UNIT SPHERE)
           DIS(A,AP)  GREAT CIRCLE DISTANCES (MILES);

   LAT (A) = (LOC(A,"LAT-DEG")  + LOC(A,"LAT-MIN") /60)*PI/180;
   LONG(A) = (LOC(A,"LONG-DEG") + LOC(A,"LONG-MIN")/60)*PI/180;

   UK(A,"X-AXIS") = COS(LONG(A))*COS(LAT(A));
   UK(A,"Y-AXIS") = SIN(LONG(A))*COS(LAT(A));
   UK(A,"Z-AXIS") =              SIN(LAT(A));

   USEG(A,AP) = SQRT(SUM(K, SQR(UK(A,K)-UK(AP,K)) ));
   UDIS(A,AP) = PI;
   UDIS(A,AP)$(USEG(A,AP) LT 1.99999) = 2*ARCTAN(USEG(A,AP)/2
                                           /SQRT(1-SQR(USEG(A,AP)/2)));
   DIS(A,AP)  = R*UDIS(A,AP);

OPTIONS LAT:5, LONG:5, UK:5, USEG:5, UDIS:5, DIS:0;

DISPLAY LOC, LAT, LONG, UK, USEG, UDIS, DIS;