'===========================================================================
' Subject: Find angle between Lat and Long    Date: 09-15-02 (  :  )       
'  Author: David Williams                     Code: Qbasic,QB,PDS          
'  Origin: david.williams@ablelink.org      Packet: ALGOR.ABC
'===========================================================================
' NAVIG.BAS 
' David Williams  2001, 2002 
  
' david.williams@ablelink.org 
  
' Calculates distances and compass bearings between points, 
' given their latitudes and longitudes. Can also be used to 
' find true north, if two points are visible from each other 
' and have both been accurately located in latitude and 
' longitude, e.g. by using GPS. The compass bearing of one 
' point from the other can be calculated, so the line of sight 
' between them has a known compass bearing, allowing north to 
' be found. 
  
DECLARE FUNCTION InRad# () 
DECLARE FUNCTION Norm# (A#) 
DECLARE FUNCTION ATan# (Y#, X#) 
  
DEFDBL A-Z 
  
DIM SHARED PI 
PI = 4 * ATN(1) 
  
DIM LAT(1 TO 2), LON(1 TO 2) 
  
PRad = 6371 ' Earth radius in km. Change to change units or planets 
  
CLS 
PRINT "This program calculates the compass bearings and distance" 
PRINT "between any two points 'A' and 'B' on the earth's surface." 
PRINT "Input longitudes in degrees [and minutes] west of Greenwich" 
PRINT "and latitudes in degrees [and minutes] north of the equator." 
PRINT "Minutes can be omitted or substituted by fractional degrees." 
PRINT "Use negative numbers (of degrees) for opposite directions." 
PRINT "Decimal fractions of degrees or minutes are permitted." 
PRINT 
PRINT "Note that the bearings are TRUE, not magnetic. If you want" 
PRINT "to calculate a magnetic bearing, you must add the local" 
PRINT "magnetic variation." 
PRINT 
PRINT "Bearings are rounded to the nearest minute (1/60 degree)." 
PRINT 
  
StA: 
PRINT "Latitude of Point A"; 
LAT(1) = InRad 
PRINT "Longitude of Point A"; 
LON(1) = InRad 
StB: 
PRINT "Latitude of Point B"; 
LAT(2) = InRad 
PRINT "Longitude of Point B"; 
LON(2) = InRad 
  
PRINT 
  
FOR P% = 1 TO 2 
  Q% = 3 - P% 
  LG = LON(P%) - LON(Q%) 
  Y = SIN(LAT(Q%)) 
  T = COS(LAT(Q%)) 
  X = T * SIN(LG) 
  Z = T * COS(LG) 
  R = SQR(Y * Y + Z * Z) 
  A = ATan(Y, Z) - LAT(P%) 
  Z = R * SIN(A) 
  B% = Norm(ATan(X, Z)) * 10800 / PI 
  IF B% = 21600 THEN B% = 0 
  
  DG% = B% \ 60 
  MN% = B% MOD 60 
  H$ = RIGHT$(STR$(INT(MN% * 5 / 3 + 100.5)), 2) 
  PRINT "Bearing from "; CHR$(64 + P%); " toward "; CHR$(64 + Q%); 
  PRINT ": "; DG%; "degrees,"; MN%; "minutes. ("; 
  PRINT LTRIM$(STR$(DG%)); "."; H$; " degrees.)" 
  
NEXT 
  
Y = -R * COS(A) 
E = ATan(Y, SQR(1 - Y * Y)) 
D! = PRad * (PI / 2 + E) 
  
PRINT 
PRINT "Distance between points: "; D!; "kilometres," 
PRINT TAB(26); D! * .62137; "statute miles,"; 
PRINT TAB(26); D! * .54; "nautical miles." 
  
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 InRad 
  DO 
    LINE INPUT " (Degrees[, Minutes])? "; X$ 
    X$ = RTRIM$(LTRIM$(X$)) 
    C% = INSTR(X$, ",") 
    IF C% <> 1 AND C% < LEN(X$) THEN EXIT DO 
    BEEP 
    PRINT "Illegal entry. Try again:"; 
  LOOP 
  IF C% THEN 
    D$ = LEFT$(X$, C% - 1) 
    M = VAL(MID$(X$, C% + 1)) 
  ELSE 
    D$ = X$ 
    M = 0 
  END IF 
  IF LEFT$(D$, 1) = "-" THEN M = -ABS(M) 
  D = (VAL(D$) + M / 60) * PI / 180 
  InRad = Norm(D) 
END FUNCTION 
  
FUNCTION Norm (A) 
  P = 2 * PI 
  Norm = A - P * INT(A / P) 
END FUNCTION