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)
