• #### distance - direction pgm

From David Williams@1:250/710 to All on Tue Jan 23 02:00:00 2001
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 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";
PRINT "Longitude of Point A";
StB:
PRINT "Latitude of Point B";
PRINT "Longitude of Point B";

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

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)