• 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 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)