Subj : distance - direction pgm
To   : All
From : David Williams
Date : Tue Jan 23 2001 02:00 am

In case anyone is interested, here's a QBasic program I wrote recently
that calculates the distance between any two points on the earth's
surface, and also the compass bearing of one as seen (along the shortest
great-circle path) from the other. This compass bearing is, of course,
the direction in which an antenna should be aimed in order to
communicate with someone at the remote point.

I hope somebody finds it useful....

                                dow

---------------------------------------

' Navig.Bas  David Williams  2001

DECLARE SUB InRad (X#)
DECLARE FUNCTION Norm# (A#)
DECLARE FUNCTION ATan# (Y#, X#)
DECLARE FUNCTION ATang# (Y#, X#)

DEFDBL A-Z

COMMON SHARED PI
PI = 4 * ATN(1)

PRad = 6371 ' Earth radius in km. Change to change units or planets

CLS
PRINT "This program calculates the compass bearing and distance of"
PRINT "a point 'B' on the earth's surface, as seen from a point 'A'."
PRINT "Input longitudes in degrees and minutes west of Greenwich"
PRINT "and latitudes in degrees and minutes north of the equator."
PRINT "Use negative numbers (of degrees) for opposite directions."
PRINT

StA:
PRINT "Latitude of Point A";
CALL InRad(LTA)
PRINT "Longitude of Point A";
CALL InRad(LGA)
StB:
PRINT "Latitude of Point B";
CALL InRad(LTB)
PRINT "Longitude of Point B";
CALL InRad(LGB)

LG = LGA - LGB
Y = SIN(LTB)
T = COS(LTB)
X = T * SIN(LG)
Z = T * COS(LG)
R = SQR(Y * Y + Z * Z)
A = ATan(Y, Z) - LTA
Y = -R * COS(A)
Z = R * SIN(A)
B = CLNG(ATang(X, Z) * 10800 / PI)
IF B = 21600 THEN B = 0
DG = B \ 60
MN = B MOD 60
E = ATan(Y, SQR(1 - Y * Y))
D = PRad * (PI / 2 + E)

PRINT
PRINT "Distance of B from A: "; CSNG(D); "kilometres."
PRINT "Bearing of B from A: "; DG; "degrees,"; MN; "minutes."

PRINT
PRINT "1. Keep A. Change B"
PRINT "2. Change A and B"
PRINT "3. Quit"
PRINT
PRINT "Which (by number) ? ";
DO
  K$ = INKEY$
LOOP UNTIL K$ >= "1" AND K$ <= "3"
PRINT K$
PRINT
ON VAL(K$) GOTO StB, StA

END

FUNCTION ATan (Y, X)
  IF X = 0 THEN
     T = SGN(Y) * PI / 2
  ELSE
     T = ATN(Y / X)
     IF X < 0 THEN T = T + PI
  END IF
  ATan = T
END FUNCTION

FUNCTION ATang (Y, X)
  ATang = Norm(ATan(Y, X))
END FUNCTION

SUB InRad (X)
  PRINT " (Degrees, Minutes including comma) ";
  INPUT D$, M
  IF LEFT$(D$, 1) = "-" THEN M = -ABS(M)
  D = (VAL(D$) + M / 60) * PI / 180
  X = Norm(D)
END SUB

FUNCTION Norm (A)
  X = A
  P = 2 * PI
  DO WHILE X < 0
     X = X + P
  LOOP
  DO WHILE X >= P
    X = X - P
  LOOP
  Norm = X
END FUNCTION

------------------------------------------
 
* Origin: FidoNet: CAP/CANADA Support BBS : 416 287-0234 (1:250/710)

.