PROGRAM Bearing_Distance;
                   {Determines bearing and distance based on}
                   {algorithm compensating for earth's shape}
                   {Copyright 1993 Michael R. Owen, W9IP & }
                   {               Paul Wade, N1BWT }
                   {Released to the public domain Feb 23, 1993}

                   {NOTE: NORTH latitude, WEST longitude apre positive;   }
                   {      South & East should be input as negative numbers}

    {$N+,E+}  {Use coprocessor if it's present, otherwise emulate w/software}

    Uses CRT;

    Type Float = Extended;   {for maximum precision}

    TYPE        grid        = String[6];

    CONST
                MyLat : Float = 0;
                MyLon : Float = 0;
                KeepGoing : Boolean = True;
                FirstTime : Boolean = True;
                MyGrid   : Grid = 'None';
      Var
            error                       : BOOLEAN ;
            FirstLat,
            FirstLon,
            OtherLat,
            OtherLon                    : Float;
            FirstGrid,
            OtherGrid                   : Grid ;
            Choice                      : Char ;

{------------------------------------------------------------}
Function Tan(x:Float): Float;
     Begin
          Tan := sin(x)/cos(x);
     end;
{------------------------------------------------------------}
Function ArcCos(x:Float): Float;
     Begin
         ArcCos := (Pi/2) - ArcTan(x/sqrt(1-x*x));
     End;
{---------------------------------------------------------}
Function ATan2 (Y, X : Float) : Float;  {returns ArcTan in the rnge 0.. 2Pi}
     var Arg : Float;
     Label 10, 11, 12, 13, 14, 15, 16, 17;
     Begin
         If X <> 0 then Arg := ArcTan(Y/X);
         If X < 0 then GOTO 10
                  else If X = 0 then GOTO 14
                                else GOTO 11;
   10 :  Atan2 := Pi + Arg;
         Exit;
   11 :  If Y < 0 then GOTO 12
                  else GOTO 13;
   12 :  Atan2 := 2.0 * Pi + Arg;
         Exit;
   13 :  Atan2 := Arg;
         Exit;
   14 :  If Y < 0 then GOTO 15
                  else if Y = 0 then GOTO 16
                                else GOTO 17;
   15 :  Atan2 := 1.5 * Pi;
         Exit;
   16 :  Atan2 := 0.0;
         Exit;
   17 :  Atan2 := 0.5 * Pi;
         Exit;
   end;  {Atan2}
{------------------------------------------------------------}
    Procedure Clear_Buff;   {Clears keyboard buffer}
       Begin
          MemW[$0:$041A] := MemW[$0:$041C];
       end;  {Clear_Buff}
{---------------------------------------------------------------}
    Procedure Wait;     {Suspends program until key pressed}
    var Dump : Char;
       Begin
            Writeln;
            Writeln('Press any key to continue');
            Repeat Until KeyPressed;
            Dump := ReadKey;  {so KeyPressed will be false upon return}
       End;  {wait}
{---------------------------------------------------------------}
    Function DMS (InVal : Float) : String;  {converts to degrees, minutes, seconds}
    Var Deg, Min, Sec : Integer;
        DegStr, MinStr, SecStr : String;

        Begin
            Deg := Trunc(InVal);
            Min := Trunc(Frac(InVal) * 60);
            Sec := Round(((Frac(InVal) * 60) - Min) * 60);
            Str(Deg:2,DegStr);
            Str(Min:2,MinStr);
            Str(Sec:2,SecStr);
            DMS := DegStr + 'ø ' + MinStr + ''' ' + SecStr + '"'
        End;    {DMS}
{---------------------------------------------------------------}
    Function NS_EW (InVal : Float; LatOrLon : String) : String;  {Adds North, South...}
        Begin
            If LatOrLon = 'Lat' then
                If InVal > 0 then NS_EW := ' North'
                             else NS_EW := ' South'
            else
                If InVal > 0 then NS_EW := ' West'
                             else NS_EW := ' East';
        end; {NS_EW}
{---------------------------------------------------------------}
    PROCEDURE Check_Grid (VAR g : grid ; VAR error : boolean ) ;
              {verifies that the grid square is legitimate}
        VAR  i  : INTEGER ;

        BEGIN
            If Length(g) = 4 then g := g + 'LL';  {choose middle if only 4-character}
            error := FALSE ;
            i := 0;

            REPEAT
              Inc(i) ;
              if ((i = 3) OR (i = 4))
                then if ((ord(g[i])>=ord('0')) AND (ord(g[i])<=ord('9')))
                     then error := FALSE  {2nd two characters are numbers}
                     else error := TRUE

                else if ((ord(g[i])>=ord('A')) AND (ord(g[i])<=ord('Z')))
                     then error := FALSE  {first and last 2 characters are letters}
                     else if ((ord(g[i])>=ord('a')) AND (ord(g[i])<=ord('z')))
                        then BEGIN
                            g[i] := chr(ord(g[i]) - ord('a') + ord('A')) ;
                            error := FALSE
                        END { else if }
                     else error := TRUE
            UNTIL ( i=6 ) OR ( error ) ;

        END ; { Check_Grid }
{------------------------------------------------------------}
    PROCEDURE Grid_center (gridsq : grid ; VAR lon,lat : Float ) ;
                          {finds the lat/lon of center of the sub-square}
        VAR
            lonmin, londeg, latmin, latdeg : Float ;

        BEGIN
            lonmin := (5 * (ord(gridsq[5])-ord('A'))) + 2.5 ;   { center }
            londeg := 180 -( 20 * (ord(gridsq[1])-ord('A')))    { tens of deg }
                        - ( 2 * (ord(gridsq[3])-ord('0'))) ;     { two deg }

            lon := londeg - (lonmin/60) ;

            latdeg := -90 + ( 10 * (ord(gridsq[2])-ord('A')))   { tens of deg }
                    + (ord(gridsq[4])-ord('0'));                 { degrees }
            latmin := 2.5 * (ord(gridsq[6])-ord('A'))           { minutes }
                       + 1.25 ;                                   { for center }

            lat := latdeg + (latmin/60) ;

        END ; { Grid_center }
{------------------------------------------------------------}
    PROCEDURE Input_Grid (VAR aGrid : grid ; VAR lon, lat : Float; Var OK : Boolean);
                         {reads grid square input from keyboard}
        VAR error       : BOOLEAN ;
            Ans         : String;

        BEGIN
            error := TRUE ;
            while error do
            BEGIN
                Write(' 4 or 6 digit grid square: [Q to quit]  ') ;
                Readln(Ans);
                OK := NOT (Ans[1] in ['q', 'Q']);
                If OK then begin
                  aGrid := ans;
                  Check_Grid(aGrid, error) ;
                  Grid_Center(aGrid,lon,lat) ;
                  if error then Writeln('   *** Error in grid format ***') ;
                end else error := False; {lets you out}
            END;  { while }

        END ;  { Input_Grid }
{------------------------------------------------------------}
   Procedure LatLon_to_Grid(LocLat, LocLon : Float; Var GS : Grid);
                           {converts lat & long to grid square}
        var c : Integer;
            g4, R, M, L4 : Float;
            M1, M2, M3, M4, M5, M6 : Char;
        BEGIN
             G4 := -LocLon+180;
             C := Trunc (G4/20);
             M1 := CHR(C+65);
             R := ABS(LocLon/20);
             R := INT(((R-INT(R))*20+0.001));
             C := Trunc(R/2);
             IF LocLon >= 0 THEN C := ABS(C-9);
             M3 := CHR(C+48);
             M := ABS(LocLon*60);
             M := ((M/120)-INT(M/120))*120;
             M := INT(M+0.001);
             C := Trunc(M/5);
             IF LocLon >=0 THEN C := ABS(C-23);
             M5 := CHR(C+65);
             L4 := LocLat+90;
             C := Trunc(L4/10);
             M2 := CHR(C+65);
             R := ABS(LocLat/10);
             R := INT(((R-INT(R))*10+0.001));
             C := Trunc(R);
             IF LocLat <0 THEN C := ABS(C-9);
             M4 := CHR(C+48);
             M := ABS(LocLat*60);
             M := ((M/60)-INT(M/60))*60;
             C := Trunc(M/2.5);
             IF LocLat <0 THEN C:=ABS(C-23);
             M6 := CHR(C+65);
             GS := M1 + M2 + M3 + M4 + M5 + M6;
end;  {Lat_Lon_to_Grid}
{------------------------------------------------------------}
    Function In_Range(Value, LowVal, HighVal : Float) : Boolean;
                 {checks that the value is between Lowval and Highval}
       BEGIN
            In_Range := (Value >= LowVal) AND (Value <= HighVal);
       END;  {In_Range}
{------------------------------------------------------------}
    PROCEDURE Input_LatLon(Title: String;
                       Var Lat, Lon : Float; var OK : Boolean);
              {reads lat & lon from keyboard; Decimal values are OK}

        Var LatDeg, LatMin, LatSec,
            LonDeg, LonMin, LonSec   : Float;

        BEGIN
            ClrScr;
            HighVideo;
            GotoXY(40-Length(title) DIV 2,10);
            Writeln(Title);
            NormVideo;
            GoToXY(5,25);
            Write('North latitude is positive, South is negative.  [enter zeros to quit]');
            GotoXY(1,12);
            Repeat
              Write('Latitude degrees: ');
              ReadLn(LatDeg);
            Until In_Range(LatDeg, -90, 90);
            Repeat
              Write('Latitude minutes: ');
              ReadLn(LatMin);
            Until In_Range(LatMin, -50.9999, 59.9999);
             If (LatDeg < 0) and (LatMin > 0) then LatMin := -LatMin;
            Repeat
              Write('Latitude seconds: ');
              ReadLn(LatSec);
            Until In_Range(LatSec, -59.9999, 59.9999);
             If (LatDeg < 0) and (LatSec > 0) then LatSec := -LatSec;
            Lat := LatDeg + (LatMin/60) + (LatSec/3600);
            OK := Abs(Lat) > 0.001;  {check for zero, which lets you out}
            WriteLn;
           If OK then Begin
            GotoXY(1,25); ClrEol;
            GoToXY(15,25);
            Write('West longitude is positive, East is negative');
            GotoXY(1,16);
            Repeat
              Write('Longitude degrees: ');
              ReadLn(LonDeg);
            Until In_Range(LonDeg, -180, 179.9999);
            Repeat
              Write('Longitude minutes: ');
              ReadLn(LonMin);
            Until In_Range(LonMin, -59.9999, 59.9999);
             If (LonDeg < 0) and (LonMin > 0) then LonMin := -LonMin;
            Repeat
              Write('Longitude seconds: ');
              ReadLn(LonSec);
            Until In_Range(LonSec, -59.9999, 59.9999);
            If (LonDeg < 0) and (LonSec > 0) then LonSec := -LonSec;
            Lon := LonDeg + (LonMin/60) + (LonSec/3600);
            OK := Abs(Lon) > 0.001;
           end;
           GoToXY(1,25); ClrEol;
        END;  {Input_LatLon}
{------------------------------------------------------------}
     PROCEDURE SetUp;  {Gets and stores user's location for future use}

       Var SetUpLat, SetUpLon : Float;
           SetUpGrid : Grid;
           Ans                : Char;
           DataFile           : Text;

        BEGIN
                   ClrScr;
                   GoToXY(1,10);
                   Writeln('YOUR current location (in decimal form): ');
                   Writeln('Latitude:   ',Abs(MyLat):8:3,'ø ', NS_EW(MyLat,'Lat'));
                   Writeln('Longitude:  ',Abs(MyLon):8:3,'ø ', NS_EW(MyLon,'Lon'));
                   Writeln('Grid Square : ',MyGrid);
                   Writeln;
                   WriteLn('Do you want to change these values? ');
                   Repeat Until KeyPressed;
                   Ans := ReadKey;
                   If Ans in ['Y', 'y'] then begin
                     Input_LatLon('Your location',SetUpLat, SetUpLon, KeepGoing);
                     WriteLn;
                     LatLon_to_Grid(SetUpLat, SetUpLon, SetUpGrid);
                     Write('Grid Square: ');
                     Highvideo; Writeln(SetUpGrid); NormVideo;
                     Writeln;
                     If KeepGoing then begin
                       Writeln('Do you want to save these values? [Y/N] ');
                       Repeat Until KeyPressed;
                       Ans := ReadKey;
                       If Ans in ['Y', 'y'] then begin   {save data}
                          Assign(DataFile,'BD.DAT');
                          ReWrite(DataFile);
                          Write(DataFile, SetUpLat:8:3,' ',SetUpLon:8:3);
                          Close(DataFile);
                          MyLat := SetUpLat;
                          MyLon := SetUpLon;
                          MyGrid := SetUpGrid;
                          Writeln('Saved');
                     end;  {Otherwise don't change anything}
                     end;
                    end;   {if user didn't want to change anything}
        END; {SetUp}
{------------------------------------------------------------}
     Procedure Calc_GeoDist(Eplat, Eplon, Stlat, Stlon : Float;
                        var Az, Baz, Dist, Deg : Float);

                 {Taken directly from:       }
                 {Thomas, P.D., 1970, Spheroidal Geodesics, reference systems,}
                 {    & local geometry, U.S. Naval Oceanographic Office SP-138,}
                 {    165 pp.}

                 {assumes North Latitude and East Longitude are positive}

                 {EpLat, EpLon = MyLat, MyLon}
                 {Stlat, Stlon = HisLat, HisLon}
                 {Az, BAz = direct & reverse azimuith}
                 {Dist = Dist (km); Deg = central angle, discarded}

Const AL = 6378206.4;   {Clarke 1866 ellipsoid}
      BL = 6356583.8;
      D2R = pi/180.0;   {degrees to radians conversion factor}
      Pi2 = 2.0 * Pi;
Label 1, 2, 3, 4, 5, 6, 7, 8, 9;
var BOA, F,
    P1R, P2R,
    L1R, L2R,
    DLR,
    T1R, T2R,
    TM, DTM,
    STM, CTM,
    SDTM,CDTM,
    KL, KK,
    SDLMR, L,
    CD, DL, SD,
    T, U, V, D, X, E, Y, A,
    FF64, TDLPM,
    HAPBR, HAMBR,
    A1M2, A2M1        : Float;

begin
        BOA := BL/AL;
        F := 1.0 - BOA;
        P1R := Eplat * D2R;
        P2R := Stlat * D2R;
        L1R := Eplon * D2R;
        L2R := StLon * D2R;
        DLR := L2R - L1R;
        T1R := ArcTan(BOA * Tan(P1R));
        T2R := ArcTan(BOA * Tan(P2R));
        TM := (T1R + T2R) / 2.0;
        DTM := (T2R - T1R) / 2.0;
        STM := Sin(TM);
        CTM := Cos(TM);
        SDTM := Sin(DTM);
        CDTM := Cos(DTM);
        KL := STM * CDTM;
        KK := SDTM * CTM;
        SDLMR := Sin(DLR/2.0);
        L := SDTM * SDTM + SDLMR * SDLMR * (CDTM * CDTM - STM * STM);
        CD := 1.0 - 2.0 * L;
        DL := ArcCos(CD);
        SD := Sin(DL);
        T := DL/SD;
        U := 2.0 * KL * KL / (1.0 - L);
        V := 2.0 * KK * KK / L;
        D := 4.0 * T * T;
        X := U + V;
        E := -2.0 * CD;
        Y := U - V;
        A := -D * E;
        FF64 := F * F / 64.0;
        Dist := AL*SD*(T -(F/4.0)*(T*X-Y)+FF64*(X*(A+(T-(A+E)/2.0)*X)+Y*(-2.0*D+E*Y)+D*X*Y))/1000.0;
        Deg := Dist/111.195;
        TDLPM := Tan((DLR+(-((E*(4.0-X)+2.0*Y)*((F/2.0)*T+FF64*(32.0*T+(A-20.0*T)*X-2.0*(D+2.0)*Y))/4.0)*Tan(DLR)))/2.0);
        HAPBR := ATan2(SDTM,(CTM*TDLPM));
        HAMBR := Atan2(CDTM,(STM*TDLPM));
        A1M2 := Pi2 + HAMBR - HAPBR;
        A2M1 := Pi2 - HAMBR - HAPBR;
   1 :  If (A1M2 >= 0.0) AND (A1M2 < Pi2) then GOTO 5
                                          else GOTO 2;
   2 :  If A1M2 >= Pi2 then GOTO 3
                       else GOTO 4;
   3 :  A1M2 := A1M2 - Pi2;
        GOTO 1;
   4 :  A1M2 := A1M2 + Pi2;
        GOTO 1;
   5 :  If (A2M1 >= 0.0) AND (A2M1 < Pi2) then GOTO 9
                                          else GOTO 6;
   6 :  If A2M1 >= Pi2 then GOTO 7
                       else GOTO 8;
   7 :  A2M1 := A2M1 - Pi2;
        GOTO 5;
   8 :  A2M1 := A2M1 + Pi2;
        GOTO 5;
   9 : Az := A1M2 / D2R;
       BAZ := A2M1 / D2R;
end; {Calc_GeoDist}
{------------------------------------------------------------------}
    PROCEDURE Read_QTH_File;  {read user's QTH data file created by SetUp}

      Var DataFile : Text;
          IOR      : Integer;
          Ans      : Char;

        BEGIN
            Assign(DataFile,'BD.DAT');
            {$I-}                     {turn off I/O checking}
            Reset(DataFile);          {see if the file is there}
            {$I+}                     {I/O checking back on}
            IOR := IOResult;
            If IOR <> 0 then begin    {if it's not there...}
                Writeln;
                HighVideo;
                Writeln('You need to set up your current station location.');
                Writeln('Press Y to run Setup or any Other key to quit');
                NormVideo;
                Repeat Until KeyPressed;
                Ans := ReadKey;
                If Ans in ['Y','y'] then SetUp
                                    else begin
                                         ClrScr;
                                         Halt;  {quit program}
                                         end;
            end else begin     {read it}
                Readln(DataFile,MyLat,MyLon);
                LatLon_to_Grid(MyLat, MyLon, MyGrid);
                FirstTime := False;  {so it won't be read again}
                Close(DataFile);
           end;
        END;   {Read_QTH_File}
{------------------------------------------------------------}
     Function Format(GridIn: Grid) : Grid;  {makes last two lettters of grid lowercase}
       Begin
          GridIn[5] := Chr(Ord(GridIn[5])+32);
          GridIn[6] := Chr(Ord(GridIn[6])+32);
          Format := GridIn;
       end;   {Format}
{------------------------------------------------------------}
    PROCEDURE Write_Results(Grid1, Grid2           : Grid;
                            Lat1, Lat2, Lon1, Lon2 : Float);
                            {output of calculations}

        VAR   PathDistmi, PathDistKm, Bearing,
              ReverseBearing, Deg  : Float ;

        BEGIN
        ClrScr;
        GotoXY(1,10);
        Calc_GeoDist(Lat1, -Lon1, Lat2, -Lon2, Bearing, ReverseBearing, PathDistKm, Deg);
        PathDistMi := PathDistKm / 1.609344;
            Writeln(' From grid ',Format(Grid1),' to grid ',Format(Grid2) ) ;
            Writeln;
            Write ('       Distance = ');
            Highvideo;
            Writeln(PathDistMi:9:3,' mi, ', PathDistKm:9:3,' km');
            NormVideo;

            Write ('   True bearing = ');
              HighVideo;  Writeln(Bearing:6:2,'ø' ) ; NormVideo;
            IF Lon1 > Lon2      {Other station is west}
                THEN Bearing :=  360.0 - Bearing ;
            Write ('Reverse bearing = ');
            HighVideo; Writeln(ReverseBearing:6:2,'ø' ) ;
            NormVideo;
            wait;
            Clear_Buff;               {clear keyboard buffer}
        END;   {Write_Results}
{------------------------------------------------------------}
{------------------------------------------------------------}
{------------------------------------------------------------}
    BEGIN { main BD rogram }
Repeat
        ClrScr;
        Highvideo;
        GotoXY(20,1); Write('BD.Pas : Bearing and Distance Calculator');
        NormVideo;
        GotoXY(10,3); Write('Calculates bearing & distance, compensating for earth''s shape');
        GoToXY(2,6); Write('REF: Thomas, P.D., 1970, Spheroidal Geodesics, reference systems, & local');
        GoToXY(2,7); Write('                  geometry, U.S. Naval Oceanographic Office #SP-138, 165 pp.');
        GotoXY(10,9); Writeln('Copyright 1993 Michael R. Owen, W9IP & Paul Wade, N1BWT');
        GotoXY(27,10); Writeln('Northern Lights Software');
        Writeln; Writeln;
        If FirstTime then Read_QTH_File;  {check stored QTH}
        Writeln('Press <1> for calculations from your QTH to another grid square');
        Writeln('      <2> for calculations from your QTH to another latitude/longitude');
        Writeln('      <3> for calculations between any two grid squares');
        Writeln('      <4> for calculations between any two latitude/longitude points');
        Writeln('      <5> to convert latitude/longitude to grid square');
        Writeln('      <6> to convert grid square to latitude/longitude');
        Writeln('      <7> set your stored QTH');
        Writeln('     <Esc> to quit the program');

        Writeln; Writeln('Your choice: ');
        Repeat until keypressed;
        Choice := ReadKey;
        Case Choice of
           '1' : Repeat
                 ClrScr;
                 GotoXY(1,10);
                 Input_Grid(OtherGrid,OtherLon,OtherLat, KeepGoing);
                 If KeepGoing then Write_Results(MyGrid, OtherGrid, MyLat, OtherLat, MyLon, OtherLon);
                 Until Not KeepGoing;
           '2' : Repeat
                 Input_LatLon('Other location',OtherLat, OtherLon, KeepGoing);
                 If KeepGoing then begin
                     LatLon_to_Grid(OtherLat, OtherLon, OtherGrid);
                     Write_Results(MyGrid, OtherGrid, MyLat, OtherLat, MyLon, OtherLon);
                 end;
                 Until Not KeepGoing;
           '3' : Repeat
                   ClrScr;
                   GotoXY(1,10);
                   Write('   Enter first') ;
                   Input_Grid (FirstGrid,FirstLon,FirstLat, KeepGoing);
                   Writeln;
                   If KeepGoing then begin
                     Write('   Enter second') ;
                     Input_Grid (OtherGrid,OtherLon,OtherLat, KeepGoing);
                   end;
                 If KeepGoing then Write_Results(FirstGrid, OtherGrid, FirstLat, OtherLat, FirstLon, OtherLon);
                 Until Not KeepGoing;
           '4' : Repeat
                   Input_LatLon('Location #1',FirstLat, FirstLon, KeepGoing);
                   LatLon_to_Grid(FirstLat, FirstLon, FirstGrid);
                   If KeepGoing then Input_LatLon('Location #2',OtherLat, OtherLon, KeepGoing);
                   LatLon_to_Grid(OtherLat, OtherLon, OtherGrid);
                   If KeepGoing then Write_Results(FirstGrid, OtherGrid, FirstLat, OtherLat, FirstLon, OtherLon);
                   Until not KeepGoing;
           '5' : Repeat
                   Input_LatLon('Location?',FirstLat, FirstLon, KeepGoing);
                   If KeepGoing then begin
                     LatLon_to_Grid(FirstLat, FirstLon, FirstGrid);
                     Writeln;
                     HighVideo;
                     Writeln('Grid square: ',Format(FirstGrid));
                     LowVideo;
                     Wait;
                   end;
                 Until Not KeepGoing;
           '6' : Repeat
                   ClrScr;
                   GoToXY(1,10);
                   Input_Grid (OtherGrid,OtherLon,OtherLat, KeepGoing);
                   If KeepGoing then begin
                      Grid_center(OtherGrid,OtherLon,OtherLat);
                      WriteLn;
                      HighVideo;
                      Writeln('Grid     :   ',Format(OtherGrid));
                      Writeln('Latitide : ',Abs(OtherLat):8:3,'ø  =  ', DMS(OtherLat), NS_EW(OtherLat,'Lat'));
                      Write('Longitude: ',Abs(OtherLon):8:3,'ø  =  ', DMS(OtherLon),NS_EW(OtherLon,'Lon'));
                      LowVideo;
                      Wait;
                   end;
                 Until Not KeepGoing;
           '7' : SetUp;
           #27 : Begin         {Esc character quits program}
                  ClrScr;
                  Halt;        {Terminate program}
                 end;
         end;  {Case}

    Until 1 = 2;  {i.e. forever; the only way out is through <Esc> above}

    END.