PROGRAM PLANETS;
{

                               Version 1.0

     This Program computes information relating to the position, distance,
  magnitude, etc for the major planets on a specified date and time.
  Refer to PRACTICAL ASTRONOMY WITH YOUR CALCULATOR by Peter Duffett-Smith
  for the calculation methods.

     This program is placed in the public domain "as is".

     Comments may be sent to the author:

             Larry Puhl
             6 Plum Court
             Sleepy Hollow, Ill. 60118
}

CONST
   W = 12; {Width of real printouts}
   D = 6;  {Decinal places in real printouts}
   Number_of_Planets = 9;
   CenterX = 300;
   CenterY = 92;
   Screen_Aspect = 2.7;
   Scale_for_inner_planets = 60;
   Scale_for_outer_planets = 2;

   Display_Degrees_in_decimal : boolean = false;
   Daylight_Savings_Correction_Enabled : boolean = true;
   Display_Hour_Angle_in_Degrees : boolean = false;
   Display_RA_in_Degrees : boolean = false;
   Zone_correction : integer = 6;
   Longitude : real = 88.3;
   Latitude : real = 43.1;


TYPE
   orbit =
      RECORD
         name : STRING[10];
         Epoch : integer;
         Period : real;
         Long_at_Epoch : real;
         Long_of_Per : real;
         Ecc : real;
         Axis : real;
         Inc : real;
         Long_of_Ascend_Node : real;
         Size : real;
         Brightness : real;
      END;
   MSDOS_REGISTERS =
      RECORD
         AX,BX,CX,DX,BP,SI,DI,DS,ES,FLAGS: Integer;
      END;

VAR
   Bios_Display_Mode : Integer Absolute $0000 : $0465;
   MSDOS_time : MSDOS_REGISTERS;
   ch : char;
   zero : real; {used to print 0}
   planets : ARRAY [1 .. 10] OF orbit;
   num : integer;
   day,month,year,day_of_year : integer;
   JD : real;


   planet_long_of_ascend_node : real;
   planet_mean_anomaly : real;
   planet_heliocentric_long : real;
   planet_true_anomaly : real;
   planet_radius_vector : real;
   planet_helio_ecliptic_lat : real;

   earth_long_of_ascend_node : real;
   earth_mean_anomaly : real;
   earth_heliocentric_long : real;
   earth_true_anomaly : real;
   earth_radius_vector : real;

   GMT : real; {Greenwich mean time}
   GST : real; {Greenwich Siderial time}
   LCT : real; {Local Civil time}
   LST : real; {Local Siderial time}

   Hour : integer;
   Minute : integer;
   T,B : real;

   projected_heliocentric_long : real;
   projected_radius_vector : real;

   Geocentric_ecliptic_Long : real;
   Geocentric_latitude : real;

   RA : real;{equatorial coordinates}
   DEC : real;

   Phase : real;
   Distance_from_Earth : real;
   Diameter : real;
   Magnitude : real;

   Hour_Angle : real;
   Hour_Angle_in_Degrees : real;
   Azimuth : real;
   Altitude : real;

   LST_Rise : real;
   LST_Set  : real;
   GST_Rise : real;
   GST_Set  : real;
   GMT_Rise : real;
   GMT_Set  : real;
   LCT_Rise : real;
   LCT_Set  : real;
   Azimuth_Set : real;
   Azimuth_Rise : real;

   y : real; {intermediate variable}
   x : real; {intermediate variable}
   z : real; {intermediate variable}
   zmax : real; {intermediate variable}
   time : real;
   Daylight_savings_start : real;
   Daylight_savings_end  : real;


{       Trig Functions     }


   FUNCTION DEGREE (X:real) : real;
      BEGIN
         DEGREE := X*57.295779513
      END; {DEGREE}

   FUNCTION RADIAN (X:real) : real;
      BEGIN
         RADIAN := X/57.295779513
      END; {RADIAN}

   FUNCTION SIND (X:real) : real;
      BEGIN
         SIND := SIN(X/57.295779513)
      END; {SIND}

   FUNCTION COSD (X:real) : real;
      BEGIN
         COSD := COS(X/57.295779513)
      END; {COSD}

   FUNCTION TAN (X:real) : real;
      BEGIN
         TAN := SIN(X)/COS(X)
      END; {TAN}

   FUNCTION TAND (X:real) : real;
      BEGIN
         TAND := TAN(X/57.295779513)
      END; {TAND}

   FUNCTION ARCTAND (X:real) : real;
      BEGIN
         ARCTAND := 57.295779513*ARCTAN(X)
      END; {ARCTAND}

   FUNCTION ARCTANYX(Y,X : real) : real;
      VAR
         z,zmax : real;
      BEGIN
         z :=ArcTanD(y/x);
         IF (y>0) AND (x>0) THEN zmax := 90;
         IF (y>0) AND (x<0) THEN zmax := 180;
         IF (y<0) AND (x<0) THEN zmax := 270;
         IF (y<0) AND (x>0) THEN zmax := 360;
         WHILE z>zmax DO z := z -180;
         WHILE z<(zmax-90) DO z := Z + 180;
         ARCTANYX := z;
      END; {ARCTANYX}

   FUNCTION ARCSIN (X:real) : real;
      BEGIN
         IF X * X >= 1 THEN ARCSIN := 1.570796327 ELSE
         ARCSIN := ARCTAN(X/SQRT(1-X*X))
      END; {ARCSIN}

   FUNCTION ARCSIND (X:real) : real;
      BEGIN
         ARCSIND := 57.295779513 * ARCSIN(X)
      END; {ARCSIND}

   FUNCTION ARCCOS (X:real) : real;
      BEGIN
         IF X * X >= 1 THEN ARCCOS := 1 ELSE
         ARCCOS := 1.570796327-ARCTAN(X/SQRT(1-X*X))
      END; {ARCCOS}

   FUNCTION ARCCOSD (X:real) : real;
      BEGIN
         ARCCOSD := 57.295779513*ARCCOS(X)
      END;


   PROCEDURE Frame(UpperLeftX, UpperLeftY, LowerRightX, LowerRightY: Integer);
   VAR
      i: Integer;
   BEGIN
      GotoXY(UpperLeftX, UpperLeftY);  Write('Ú');
      FOR i:=UpperLeftX+1 TO LowerRightX-1 DO Write('Ä');
      Write('¿');
      FOR i:=UpperLeftY+1 TO LowerRightY-1 DO
      BEGIN
         GotoXY(UpperLeftX , i);  Write('³');
         GotoXY(LowerRightX, i);  Write('³');
      END;
            GotoXY(UpperLeftX, LowerRightY);
      Write('À');
      FOR i:=UpperLeftX+1 TO LowerRightX-1 DO Write('Ä');
      Write('Ù');
   END  { Frame };


{      Astronomy Functions     }


   FUNCTION Minutes (time :real) : integer;
       BEGIN
            Minutes := ABS(TRUNC(60*(FRAC(time))))
       END;

   FUNCTION Seconds (time :real) : integer;
       BEGIN
            Seconds := ABS(TRUNC(60*(FRAC(60*FRAC(time)))))
       END;

   PROCEDURE deg_min_sec (VAR angle:real);
   BEGIN
      IF NOT Display_Degrees_in_Decimal
         THEN
            BEGIN
               IF angle < 0 THEN write('-') ELSE write(' ');
               write(abs(trunc(angle)):3,'ø ',minutes(angle):2,''' ',
                     seconds(angle):2,'"');
            END
         ELSE
            write(angle:8:4,'ø ');
       END;

   PROCEDURE hours_min_sec (VAR angle:real);
   BEGIN
      IF NOT Display_Degrees_in_Decimal
         THEN
            BEGIN
               IF angle < 0 THEN write('-') ELSE write(' ');
               write(abs(trunc(angle/15)):3,'h ',minutes(angle/15):2,'m ',
                     seconds(angle/15):2,'s');
            END
         ELSE
            write(angle/15:8:4,'h');
       END;

   FUNCTION Daylight_Savings : boolean;
     BEGIN
       Daylight_Savings := (((Month + Day/100) < Daylight_savings_end)  AND
                           ((Month + Day/100) > Daylight_savings_start) AND
                           Daylight_savings_Correction_Enabled);
     END;

   FUNCTION Anomaly (M,Ecc : real) : real;
       VAR
          d : real;
          E : real;
       BEGIN
          E := M;
          d := E - Ecc * SIN(E) - M;
          WHILE abs(d) > 0.000001 DO
             BEGIN
                E := E - d/(1-Ecc * COS(E));
                d := E - Ecc * SIN(E) - M
             END;
          Anomaly := 2 * ArcTan(sqrt((1+Ecc)/(1-Ecc)) * Tan(E/2))
        END;{Anomaly}

   FUNCTION JulianDay (month,day,year :integer) : real;
       VAR
          A,B,C,D : real;

       BEGIN
         IF (month=1) OR (month=2)
            THEN
               BEGIN
                  year := year-1;
                  month := month+12
               END;
         IF (year>1582) OR (year=1582) AND (month>10) OR (year=1582) AND (month=10) AND (day>15)
            THEN
               BEGIN
                  A := INT(year/100);
                  B := 2-A+INT(A/4)
               END
            ELSE
               BEGIN
                  B := 0
               END;

        C := INT(365.25*year);
        D := INT(30.6001*(month+1));
        JulianDay := B+C+D+day+1720994.5;
      END;{JulianDay}

FUNCTION ECL_TO_RA (L,B : real) : real;
   VAR
     x : real;
     y : real;
     z : real;
     zmax : real;
   BEGIN
     x := COSD(L);
     y := SIND(L) * 0.91746406 - TAND(B) * 0.397818676;
     z := ArcTanYX(y,x);
     ECL_TO_RA := Z;
   END;{ECL_TO_RA}

FUNCTION ECL_TO_DEC (L,B : real) : real;
   BEGIN
     ECL_TO_DEC := ARCSIND(SIND(B) * 0.91746406 + COSD(B) * SIND(L)
                   * 0.397818676);
   END;

FUNCTION LST_TO_LCT (LST,Long : real;
                     Year, Day_of_year : integer;
                     Zone_Corr : real) : real;
   VAR
     B,T,GST,GMT : real;
   BEGIN
        GST := LST + Long/15;
        IF GST > 24 THEN GST := GST - 24;
        IF GST < 0 THEN GST := GST + 24;
        T := (JulianDay(1,0,year) - 2415020.0)/36525.0;
        B := 24 - 6.6460656 - (2400.051262 * T) - (0.00002581 * T * T) +
             (24 * (year - 1900));
        T := GST -  Day_of_Year  * 0.0657098 + B;
        IF T > 24 THEN T := T - 24;
        IF T < 0 THEN T := T + 24;
        GMT := T * 0.99727;
        LST_TO_LCT := GMT - Zone_Corr;
    END;

PROCEDURE make_degrees_in_range(VAR n : real);
BEGIN
   WHILE n > 360 DO
      n := n - 360;
   WHILE n < 0 DO
      n := n + 360;
END;

PROCEDURE Make_Hours_in_Range(VAR n : real);
   BEGIN
      IF n > 24 THEN n := n - 24;
      IF n <  0 THEN n := n + 24;
   END;

PROCEDURE time_window;
   BEGIN
   WINDOW(1,1,80,25);
   frame(56,1,80,10);
   WINDOW(58,2,78,10);
   GotoXY(1,1);
   writeln('    ',month,'-',day,'-',year);
   JD := Julianday(month,day,year);
   daylight_savings_start :=
          4 + (30 - Round(7 * Frac((Julianday(4,30,year) + 1.5)/7))) * 0.01;
   daylight_savings_end :=
          10 + (31 - Round(7 * Frac((Julianday(10,31,year) + 1.5)/7))) * 0.01;
   writeln(' JD =',JD:10:0);
   write('Long = ');
   deg_min_sec(Longitude);
   writeln;
   write('Lat  = ');
   deg_min_sec(Latitude);
   writeln;
   LCT := Hour + Minute/60 + 1/120;
   writeln('    LCT = ',Hour:2,'h ',Minute:2,'m');
   GMT := Zone_correction + LCT;
   IF  Daylight_Savings THEN
      GMT := GMT - 1;
   Make_Hours_in_Range(GMT);
   writeln ('    GMT = ',trunc(GMT):2,'h ',minutes(GMT):2,'m');
   T := (JulianDay(1,0,year) - 2415020.0)/36525.0;
   B := 24 - 6.6460656 - (2400.051262 * T) - (0.00002581 * T * T) +
        (24 * (year - 1900));
   GST := 0.0657098 * Day_of_year - B + GMT * 1.002738;
   Make_Hours_in_Range(GST);
   writeln ('    GST = ',trunc(GST):2,'h ',minutes(GST):2,'m');
   LST := GST - Longitude/15;
   Make_Hours_in_Range(LST);
   writeln ('    LST = ',trunc(LST):2,'h ',minutes(LST):2,'m');
   END;{time_window}

PROCEDURE Locate_Position_of_Planet_in_Its_Own_Orbital_Plane(n:integer);
   BEGIN
    WITH planets[n] DO
     BEGIN
       planet_long_of_ascend_node := (360/365.2422) * (Julianday(month,day,year)
                                  - Julianday(1,0,Epoch)) / period;
       Make_Degrees_in_Range(planet_long_of_ascend_node);
       planet_mean_anomaly := planet_long_of_ascend_node + long_at_epoch
                              - long_of_per;
       planet_true_anomaly := DEGREE(ANOMALY(RADIAN(planet_mean_anomaly),Ecc));
       planet_heliocentric_long := planet_true_anomaly + long_of_per;
       Make_Degrees_in_Range(planet_heliocentric_long);
       planet_radius_vector := axis*(1-Ecc*Ecc)/(1+Ecc*COSD(planet_true_anomaly));
       planet_helio_ecliptic_lat := ARCSIND(SIND(planet_heliocentric_long -
                                    long_of_ascend_node) * SIND(Inc));
    END;

   END;


PROCEDURE Plot_planet(n,scale : integer);
FUNCTION compute_radius_vector(n,angle : integer): real ;
BEGIN
   compute_radius_vector := planets[n].axis *
                            (1 - planets[n].ecc * planets[n].ecc)
                            /(1 + planets[n].ecc * COSD(angle
                            - planets[n].long_of_per));
END;
VAR
delX,delY,Xold,Xnew,Yold,Ynew,A : integer;
BEGIN
   Locate_Position_of_Planet_in_Its_Own_Orbital_Plane(n);
   Ynew := Round(scale * planet_radius_vector * SIND(planet_heliocentric_long))
           + CenterY;
   Xnew := Round(scale * Screen_Aspect * planet_radius_vector *
           COSD(planet_heliocentric_long)) + CenterX;
   FOR delX := -3 TO 3 DO
      FOR delY := -1 TO 1 DO
          plot(Xnew + delX,Ynew + delY,1);
   GoToXY((Xnew DIV 8)+2,(Ynew DIV 8)+1);
   write(planets[n].name);

   Xnew := CenterX + Round(compute_radius_vector(n,0) * scale * Screen_Aspect);
   Ynew := CenterY;
   FOR A := 1 TO 30 DO
      BEGIN
         Xold := Xnew;
         Yold := Ynew;
         Xnew := Round(compute_radius_vector(n,A*12) * scale * Screen_Aspect
                 * COSD(A * 12)) + CenterX;
         Ynew := Round(compute_radius_vector(n,A*12) * scale * SIND(A * 12))
                 + CenterY;
         Draw(Xold,Yold,Xnew,Ynew,1);
      END;
END;

{    Orbital Data for Planets   }


PROCEDURE Orbital_Data;
BEGIN
WITH planets[1] DO
   BEGIN
      name :='Mercury   ';
      Epoch := 1980;
      Period := 0.24085;
      Long_at_Epoch := 231.2973;
      Long_of_Per := 77.1442128;
      Ecc := 0.2056306;
      Axis := 0.3870986;
      Inc := 7.0043579;
      Long_of_Ascend_Node := 48.0941733;
      Size := 6.74;
      Brightness := 0.000001918;
   END;
WITH planets[2] DO
   BEGIN
      name :='Venus     ';
      Epoch := 1980;
      Period := 0.61521;
      Long_at_Epoch := 355.73352;
      Long_of_Per := 131.2895792;
      Ecc := 0.0067826;
      Axis := 0.7233316;
      Inc := 3.3944359;
      Long_of_Ascend_Node := 76.4997524;
      Size := 16.92;
      Brightness := 0.00001721;
   END;
WITH planets[3] DO
   BEGIN
      name :='Earth     ';
      Epoch := 1980;
      Period := 1.00004;
      Long_at_Epoch := 98.833540;
      Long_of_Per := 102.596403;
      Ecc := 0.016718;
      Axis := 1.0;
      Inc := 0;
      Long_of_Ascend_Node := 0;
      Size := 17;
      Brightness := 0;
   END;
WITH planets[4] DO
   BEGIN
      name :='Mars      ';
      Epoch := 1980;
      Period := 1.88089;
      Long_at_Epoch := 126.30783;
      Long_of_Per := 335.6908166;
      Ecc := 0.0933865;
      Axis := 1.5236883;
      Inc := 1.8498011;
      Long_of_Ascend_Node := 49.4032001;
      Size := 9.36;
      Brightness := 0.00000454;
   END;
WITH planets[5] DO
   BEGIN
      name :='Jupiter   ';
      Epoch := 1980;
      Period := 11.86224;
      Long_at_Epoch := 146.966365;
      Long_of_Per := 14.0095493;
      Ecc := 0.0484658;
      Axis := 5.202561;
      Inc := 1.3041819;
      Long_of_Ascend_Node := 100.2520175;
      Size := 196.74;
      Brightness := 0.0001994;
   END;
WITH planets[6] DO
   BEGIN
      name :='Saturn    ';
      Epoch := 1980;
      Period := 29.4571;
      Long_at_Epoch := 165.322242;
      Long_of_Per := 92.6653974;
      Ecc := 0.0556155;
      Axis := 9.554747;
      Inc := 2.4893741;
      Long_of_Ascend_Node := 113.4888341;
      Size := 165.60;
      Brightness := 0.0001740;
   END;

WITH planets[7] DO
   BEGIN
      name :='Uranus    ';
      Epoch := 1980;
      Period := 84.01247;
      Long_at_Epoch := 228.0708551;
      Long_of_Per := 172.7363288;
      Ecc := 0.0463232;
      Axis := 19.21814;
      Inc := 0.7729895;
      Long_of_Ascend_Node := 73.8768642;
      Size := 65.80;
      Brightness := 0.00007768;
   END;
WITH planets[8] DO
   BEGIN
      name :='Neptune   ';
      Epoch := 1980;
      Period := 164.79558;
      Long_at_Epoch := 260.3578998;
      Long_of_Per := 47.8672148;
      Ecc := 0.0090021;
      Axis := 30.10957;
      Inc := 1.7716017;
      Long_of_Ascend_Node := 131.5606494;
      Size := 62.20;
      Brightness := 0.00007597;
   END;

WITH planets[9] DO
   BEGIN
      name :='Pluto     ';
      Epoch := 1980;
      Period := 250.9;
      Long_at_Epoch := 209.439;
      Long_of_Per := 222.972;
      Ecc := 0.25387;
      Axis := 39.78459;
      Inc := 17.137;
      Long_of_Ascend_Node := 109.941;
      Size := 8.2;
      Brightness := 0.000004073;
   END;
END; {Orbital_Data}

PROCEDURE Show_Menu;
BEGIN
   WINDOW(1,1,80,25);
   CLRSCR;
   Time_Window;
   WINDOW(1,1,80,25);
   GotoXY(1,1);

   FOR num := 1 TO Number_of_Planets DO
      BEGIN
         WITH planets[num] DO
            writeln(num,'  ',name);
      END;
   writeln;
   writeln('D  Change Date/Time');
   writeln('I  Plot Inner Planets');
   writeln('L  Change Long/Lat');
   writeln('M  Menu');
   writeln('O  Plot Outer Planets');
   writeln('S  Set_up');
   writeln('Q  Quit');
   writeln;
   writeln('Enter command:');
END; {Show_menu}

PROCEDURE Set_Up;
BEGIN
   writeln('Display Degrees in Decimal Degrees ( Y/N )? ');
   Read(KBD,Ch);
   IF (Ch = 'Y') OR (Ch = 'y')
      THEN Display_Degrees_in_Decimal := true;
   IF (Ch = 'N') OR (Ch = 'n')
      THEN Display_Degrees_in_Decimal := false;
   writeln('Correct for daylight savings ( Y/N )?');
   Read(KBD,Ch);
   IF (Ch = 'Y') OR (Ch = 'y')
      THEN Daylight_Savings_Correction_Enabled := true;
   IF (Ch = 'N') OR (Ch = 'n')
      THEN Daylight_Savings_Correction_Enabled := false;
   writeln('Display Hour Angle in Degrees ( Y/N )?');
   Read(KBD,Ch);
   IF (Ch = 'Y') OR (Ch = 'y')
      THEN Display_Hour_Angle_in_Degrees := true;
   IF (Ch = 'N') OR (Ch = 'n')
      THEN Display_Hour_Angle_in_Degrees := false;
   writeln('Display RA in Degrees ( Y/N )?');
   Read(KBD,Ch);
   IF (Ch = 'Y') OR (Ch = 'y')
      THEN Display_RA_in_Degrees := true;
   IF (Ch = 'N') OR (Ch = 'n')
      THEN Display_RA_in_Degrees := false;
END;


PROCEDURE Plot_Inner_Planets;
BEGIN
   HiRes;
   HiResColor(1);
   Window(1,1,80,25);
   GotoXY(62,1);
   writeln('    ',month,'-',day,'-',year);
   GotoXY(62,2);
   writeln(' JD =',JD:10:0);
   plot(CenterX,CenterY,1);
   FOR num := 1 TO 4 DO
      plot_planet(num,scale_for_inner_planets);
END;


PROCEDURE Plot_Outer_Planets;
BEGIN
   HiRes;
   HiResColor(1);
   Window(1,1,80,25);
   GotoXY(62,1);
   writeln('    ',month,'-',day,'-',year);
   GotoXY(62,2);
   writeln(' JD =',JD:10:0);
   plot(CenterX,CenterY,1);
   plot_planet(3,scale_for_outer_planets);
   FOR num := 5 TO 9 DO
      plot_planet(num,scale_for_outer_planets);
END;


PROCEDURE Change_Date_Time;
     BEGIN
       time_window;
       WINDOW(1,1,80,25);
       writeln;
       writeln ('Enter month  day  year');
       readln (month,day,year);
       writeln ('Enter hour minute');
       readln (hour,minute);
       day_of_year := TRUNC(Julianday(month,day,year)-Julianday(1,0,year));
       time_window;
     END;



PROCEDURE Change_Long_Lat;
    BEGIN
      time_window;
      WINDOW(1,1,80,25);
      writeln;
      writeln ('Enter Longitude  Latitude');
      readln (longitude,Latitude);
      Zone_correction := trunc(longitude/15) + 1;
      writeln ('Enter zone correction (Default = ',zone_correction,' ):');
      readln (zone_correction);
      time_window;
    END;


PROCEDURE Locate_Planet;
LABEL
End_of_Locate_Planet;
BEGIN
   Time_window;
   WINDOW(1,1,80,25);
   frame(1,1,54,9);
   WINDOW(3,2,53,10);
   GotoXY(1,1);
   Locate_Position_of_Planet_in_Its_Own_Orbital_Plane(num);
   WITH planets[num] DO
     BEGIN
       writeln('                          ',name);
       write('long of ascend node  = ');
       deg_min_sec(planet_long_of_ascend_node);
       writeln;
       write('mean anomaly         = ');
       deg_min_sec(planet_mean_anomaly);
       writeln;
       write('true anomaly         = ');
       deg_min_sec(planet_true_anomaly);
       writeln;
       write('heliocentric long    = ');
       deg_min_sec(planet_heliocentric_long);
       writeln;
       writeln('radius vector        = ',planet_radius_vector:8:3,' AU');
       write('helio ecliptic lat.  =');
       deg_min_sec(planet_helio_ecliptic_lat);
       writeln;
     END;

{  Locate Position of Earth in Its Own Orbital Plane   }

   IF num = 3 THEN GoTO End_of_Locate_Planet;
   WINDOW(40,2,53,10);
   GotoXY(1,1);
   writeln('     Earth');
   WITH planets[3] DO
     BEGIN

       earth_long_of_ascend_node := (360/365.2422) * (Julianday(month,day,year)
                                  - Julianday(1,0,Epoch))/period;
       Make_Degrees_in_Range(earth_long_of_ascend_node);
       deg_min_sec(earth_long_of_ascend_node);
       writeln;

       earth_mean_anomaly := earth_long_of_ascend_node + long_at_epoch
                              - long_of_per;
       deg_min_sec(earth_mean_anomaly);
       writeln;

       earth_true_anomaly := DEGREE(ANOMALY(RADIAN(earth_mean_anomaly),Ecc));
       deg_min_sec(earth_true_anomaly);
       writeln;

       earth_heliocentric_long := earth_true_anomaly + long_of_per;
       Make_Degrees_in_Range(earth_heliocentric_long);
       deg_min_sec(earth_heliocentric_long);
       writeln;

       earth_radius_vector := axis*(1-Ecc*Ecc)/(1+Ecc*COSD(earth_true_anomaly));
       writeln(earth_radius_vector:8:3,' AU');
       zero := 0;
       deg_min_sec(zero);
     END;


{  Project Position of Planet onto plane of ecliptic   }


   WINDOW(1,1,80,25);
   frame(1,10,54,13);
   WINDOW( 12,11,79,23);
   GotoXY(1,1);
   WITH planets[num] DO
     BEGIN
        y := SIND(planet_heliocentric_long - long_of_ascend_node) * COSD(Inc);
        x := COSD(planet_heliocentric_long - long_of_ascend_node);
        z :=ArcTanYX(y,x);
        projected_heliocentric_long := Z + long_of_ascend_node;
        Make_Degrees_in_Range(projected_heliocentric_long);
        write('Projected Longitude  = ');
        deg_min_sec(projected_heliocentric_long);
        writeln;
        Projected_Radius_vector := planet_radius_vector *
                                   COSD(planet_helio_ecliptic_lat);
        writeln('Projected Radius     = ',Projected_Radius_Vector:8:3,' AU');
      END;

{    Calculate Ecliptical Coordinates   }


   WINDOW(1,1,80,25);
   frame(1,18,25,24);
   WINDOW(3,19,24,23);
   GotoXY(1,1);
   IF (Planet_Radius_Vector < Earth_Radius_Vector)
   THEN
      Geocentric_Ecliptic_Long := Earth_Heliocentric_long + 180 +
           ArcTanD(Projected_Radius_Vector * SIND(Earth_Heliocentric_long -
           Projected_heliocentric_long)/(Earth_Radius_Vector -
           Projected_Radius_Vector * COSD(Earth_Heliocentric_long -
           Projected_heliocentric_long)))
   ELSE
      Geocentric_Ecliptic_Long := Projected_Heliocentric_Long +
           ArcTanD(Earth_Radius_Vector * SIND(Projected_heliocentric_long -
           Earth_Heliocentric_Long) / (Projected_Radius_Vector -
           Earth_Radius_Vector * COSD(Projected_heliocentric_long -
           Earth_Heliocentric_Long)));
    Make_Degrees_in_Range(Geocentric_Ecliptic_Long);
   writeln('      ECLIPTIC');
   writeln;
   write(' Long = ');
   deg_min_sec(Geocentric_Ecliptic_Long);
   writeln;
   Geocentric_latitude := ArcTanD(Projected_Radius_Vector *
                          TAND(planet_helio_ecliptic_lat) *
                          SIND(Geocentric_Ecliptic_Long -
                          Projected_Heliocentric_Long) / (Earth_Radius_Vector *
                          SIND(Projected_Heliocentric_Long -
                          Earth_Heliocentric_Long)));
   write(' Lat  = ');
   deg_min_sec(Geocentric_latitude);
   writeln;


{     Calculate Equatorial Coordinates    }

   WINDOW(1,1,80,25);
   frame(27,18,52,24);
   WINDOW(29,19,51,24);
   GotoXY(1,1);
   writeln('     EQUATORIAL');
   writeln;
   RA := ECL_TO_RA(Geocentric_Ecliptic_Long,Geocentric_Latitude);
   write(' RA  = ');
   IF Display_RA_in_Degrees
   THEN
      deg_min_sec(RA)
   ELSE
      hours_min_sec(RA);
   writeln;
   DEC := ECL_TO_DEC(Geocentric_Ecliptic_Long,Geocentric_Latitude);
   write(' DEC = ');
   deg_min_sec(DEC);
   writeln;
   Hour_Angle := LST - RA/15;
   IF Hour_Angle < 0 THEN Hour_Angle := Hour_Angle + 24;
   Hour_Angle_in_Degrees := Hour_Angle * 15;
   write(' HA  = ');
   IF Display_Hour_Angle_in_Degrees
   THEN
      deg_min_sec(Hour_Angle_in_Degrees)
   ELSE
      hours_min_sec(Hour_Angle_in_Degrees);

{    Calculate Horizontal Coordinates   }

   WINDOW(1,1,80,25);
   frame(54,18,80,24);
   WINDOW(56,19,78,23);
   GotoXY(1,1);
   writeln('      HORIZONTAL');
   writeln;
   Altitude := ArcSinD(SIND(Dec) * SIND(Latitude) +
               COSD(Dec) * COSD(Latitude) * COSD(Hour_Angle * 15));
   write(' Alt  = ');
   deg_min_sec(Altitude);
   writeln;

   Azimuth := ArcCosD((SIND(Dec) - SIND(Latitude) * SIND(Altitude))/
              (COSD(Latitude) * COSD(Altitude)));
   IF SIND(Hour_Angle) > 0 THEN
   Azimuth := 360 - Azimuth;
   write(' Azim = ');
   deg_min_sec(Azimuth);


{   Calculate Time and Position of Rise and Set   }


   WINDOW(1,1,80,25);
   frame(1,14,54,17);
   WINDOW(4,15,52,17);
   GotoXY(1,1);
   Azimuth_Rise := SIND(DEC)/COSD(Latitude);
   IF (Azimuth_Rise < -1) OR (Azimuth_Rise > 1) THEN
      IF Altitude < 0
         THEN
            writeln('never rises above horizon')
         ELSE
            writeln('never sets below horizon')
   ELSE
      BEGIN
        Azimuth_Rise := ARCCOSD(Azimuth_Rise);
        Azimuth_Set := 360 - Azimuth_Rise;
        LST_Rise := 24 + RA/15 - (ARCCOSD(-TAND(Latitude) * TAND(DEC)))/15;
        IF LST_Rise > 24 THEN LST_Rise := LST_Rise - 24;
        LST_SET := RA/15 + (ARCCOSD(-TAND(Latitude) * TAND(DEC)))/15;
        IF LST_SET > 24 THEN LST_SET := LST_SET - 24;

        LCT_SET := LST_TO_LCT(LST_SET,Longitude,Year,Day_of_Year,
                    Zone_Correction);
        IF Daylight_Savings THEN
           LCT_SET := LCT_SET + 1;
        IF LCT_SET < 0 THEN LCT_SET := LCT_SET + 24;

        LCT_RISE := LST_TO_LCT(LST_RISE,Longitude,Year,Day_of_Year,
                    Zone_Correction);
        IF Daylight_Savings THEN
           LCT_RISE := LCT_RISE + 1;
        IF LCT_RISE < 0 THEN LCT_RISE := LCT_RISE + 24;

        write('Rises at ',trunc(LCT_Rise):2,':',minutes(LCT_Rise):2,'  LCT');
        write('      Azimuth  ');
        deg_min_sec(Azimuth_Rise);
        writeln;
        write('Sets at  ',trunc(LCT_Set):2,':',minutes(LCT_Set):2,'  LCT');
        write('      Azimuth  ');
        deg_min_sec(Azimuth_Set);
      END;


{    Calculate Phase,Distance, Diameter, Magnitude   }


   Window(1,1,80,25);
   frame(56,11,80,17);
   Window(57,12,79,16);
   GotoXY(1,1);

   Phase := (1+COSD(Geocentric_ecliptic_long - planet_heliocentric_long))/2;
   writeln('Phase     = ',100 * Phase:6:2,'%');
   Distance_from_Earth := sqrt(Earth_Radius_Vector * Earth_Radius_Vector
                          + Planet_Radius_Vector * Planet_Radius_Vector
                          - 2 * Earth_Radius_Vector * Planet_Radius_Vector *
                          COSD(planet_heliocentric_long -
                          earth_heliocentric_long));
   writeln('Distance  = ',Distance_from_Earth:6:2,' AU');

   WITH Planets[num] DO
   Diameter := size / Distance_from_Earth;
   writeln('Diameter  = ',Diameter:6:2,'"');

   WITH Planets[num] DO
   Magnitude := 2.17147 * ln(Distance_from_Earth * Planet_Radius_Vector /
                Brightness * sqrt(Phase)) - 26.7;
   writeln('Magnitude = ',Magnitude:6:2);
End_of_Locate_Planet:
END; {Locate_Planet}

BEGIN {Planets}

{      Initalize Default Parameters   }

   CLRSCR;
   MSDOS_time.AX := $2C00;{time call to MSDOS}
   MsDos(MSDOS_time);
   Hour :=Hi(MSDOS_time.CX);
   Minute := Lo(MSDOS_time.CX);
   MSDOS_time.AX := $2A00; {Date call to MSDOS}
   MsDos(MSDOS_time);
   Year := MSDOS_time.CX;
   Month := Hi(MSDOS_time.DX);
   Day := Lo(MSDOS_time.DX);
   Orbital_Data;
   day_of_year := TRUNC(Julianday(month,day,year)-Julianday(1,0,year));

{     Main Program     }


Show_Menu;
REPEAT
   Read(KBD, Ch);
   num := ord(Ch) - ord('0');
   TextMode;
   Bios_Display_Mode := $2D;
   ClrScr;
   CASE Ch OF
      'I','i' :  Plot_Inner_Planets;
      'O','o' :  Plot_Outer_Planets;
      'D','d' :  Change_Date_Time;
      'L','l' :  Change_Long_Lat;
      'M','m' :  Show_Menu;
      'S','s' :  Set_Up;
      '1','2','3','4','5','6','7','8','9'  : Locate_Planet;
   END; {case}
   Window(1,1,80,25);
   GotoXY(1,25);
   write('(Q)uit (M)enu (I)nner (O)uter (D)ate (L)ong   (1-9) Planet');
UNTIL (Ch = 'Q') OR (Ch = 'q');
   TextMode;
   Bios_Display_Mode := $2D;
   ClrScr;
END.ööööööööööööööööööööööööööööööö