seed7/prg/planets.sd7

(********************************************************************)
(*                                                                  *)
(*  planets.sd7   Display information about the planets             *)
(*  Copyright (C) 2006, 2007  Thomas Mertes                         *)
(*                                                                  *)
(*  This program is free software; you can redistribute it and/or   *)
(*  modify it under the terms of the GNU General Public License as  *)
(*  published by the Free Software Foundation; either version 2 of  *)
(*  the License, or (at your option) any later version.             *)
(*                                                                  *)
(*  This program is distributed in the hope that it will be useful, *)
(*  but WITHOUT ANY WARRANTY; without even the implied warranty of  *)
(*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the   *)
(*  GNU General Public License for more details.                    *)
(*                                                                  *)
(*  You should have received a copy of the GNU General Public       *)
(*  License along with this program; if not, write to the           *)
(*  Free Software Foundation, Inc., 51 Franklin Street,             *)
(*  Fifth Floor, Boston, MA  02110-1301, USA.                       *)
(*                                                                  *)
(********************************************************************)


$ include "seed7_05.s7i";
  include "time.s7i";
  include "duration.s7i";
  include "float.s7i";
  include "keybd.s7i";
  include "echo.s7i";
  include "line.s7i";
  include "draw.s7i";
  include "graph_file.s7i";
  include "window.s7i";
  include "stars.s7i";
  include "math.s7i";

var text: scr is STD_NULL;

const float:   julianDayOfEpoch is 2451545.0; # 2000-01-01 12:00:00
const integer: numberOfPlanets is 9;
const integer: centerX is 440;
const integer: centerY is 345;
const integer: panoramaHorizont is 700;
const integer: panoramaXMin is 62;
const float:   sizePanorama is 900.0;
const integer: windowWidth is 1024;
const integer: windowHeight is 768;
const float:   screenAspect is 1.0;
const integer: scaleForInnerPlanets is 210;
const integer: scaleForOuterPlanets is 9;
const float:   radianToDegrees is 57.295779513082320876798154814114;
const float:   degreesToRadian is 0.017453292519943295769236907684883;

var boolean: displayDegreesInDecimal is FALSE;
var boolean: displayHourAngleInDegrees is FALSE;
var boolean: displayRAinDegrees is FALSE;


const type: orbitType is new struct
    var string: name is "";
    var float: orbitalPeriod is 0.0;                   # days
    var float: meanLongitudeAtEpoch is 0.0;            # deg
    var float: longitudeOfPerihelion is 0.0;           # deg
    var float: longitudeOfPerihelionPerCentury is 0.0; # deg/century
    var float: eccentricity is 0.0;                    # rad
    var float: eccentricityPerCentury is 0.0;          # rad/century
    var float: axis is 0.0;                            # AU
    var float: axisPerCentury is 0.0;                  # AU/century
    var float: inclination is 0.0;
    var float: longitudeOfTheAscendingNode is 0.0;
    var float: argumentOfThePerihelion is 0.0;
    var float: size is 0.0;
    var float: brightness is 0.0;
  end struct;

var array orbitType: orbitDescr is numberOfPlanets times orbitType.value;

var time: aTime is time.value;
var integer: currTimeChange is 4;

const type: heliocentricPlanetType is new struct
    var float: axis is 0.0;
    var float: eccentricity is 0.0;
    var float: meanLongitude is 0.0;
    var float: longitudeOfPerihelion is 0.0;
    var float: meanAnomaly is 0.0;
    var float: heliocentricLongitude is 0.0;
    var float: trueAnomaly is 0.0;
    var float: radiusVector is 0.0;
    var float: heliocentricEclipticLatitude is 0.0;
    var float: projectedHeliocentricLongitude is 0.0;
    var float: projectedRadius is 0.0;
  end struct;

const type: eclipticalCoordinateType is new struct
    var float: geocentricEclipticLongitude is 0.0;
    var float: geocentricLatitude is 0.0;
  end struct;

const type: equatorialCoordinateType is new struct
    var float: rightAscension is 0.0;
    var float: declination is 0.0;
  end struct;

const type: horizontalCoordinateType is new struct
    var float: azimuth is 0.0;
    var float: altitude is 0.0;
  end struct;

const type: earthCoordinateType is new struct
    var float: longitude is -16.42;
    var float: latitude  is  48.25;
  end struct;

const type: geocentricPlanetType is new struct
    var eclipticalCoordinateType: eclipticalCoordinate is eclipticalCoordinateType.value;
    var equatorialCoordinateType: equatorialCoordinate is equatorialCoordinateType.value;
    var float: hourAngle is 0.0;
    var horizontalCoordinateType: horizontalCoordinate is horizontalCoordinateType.value;
  end struct;

var earthCoordinateType: posAtEarth is earthCoordinateType.value;
var float: localSiderealTime is 0.0;


(* Trigonometric functions for degrees *)


const func float: toDegrees (in float: x) is
  return x * radianToDegrees;


const func float: toRadians (in float: x) is
  return x * degreesToRadian;


const func float: sinD (in float: x) is
  return sin(x * degreesToRadian);


const func float: cosD (in float: x) is
  return cos(x * degreesToRadian);


const func float: tanD (in float: x) is
  return tan(x * degreesToRadian);


const func float: asinD (in float: x) is
  return asin(x) * radianToDegrees;


const func float: acosD (in float: x) is
  return acos(x) * radianToDegrees;


const func float: atanD (in float: x) is
  return atan(x) * radianToDegrees;


const func float: atan2D (in float: y, in float: x) is
  return atan2(y, x) * radianToDegrees;


const func float: frac (in float: x) is
  return x - flt(trunc(x));


(* Astronomy functions *)


const func float: julianDay (in time: aTime) is
  return flt(julianDayNumber(aTime)) +
      ((flt(aTime.second) / 60.0 + flt(aTime.minute)) / 60.0 + flt(aTime.hour - 12)) / 24.0;


const func integer: minutes (in float: time) is
  return abs(trunc(60.0 * (frac(time))));


const func integer: seconds (in float: time) is
  return abs(trunc(60.0 * (frac(60.0 * frac(time)))));


const func string: angleAsDegrees (in float: angle) is func
  result
    var string: degreesMinutesSeconds is "";
  begin
    if displayDegreesInDecimal then
      degreesMinutesSeconds &:= angle digits 4 lpad 9;
      degreesMinutesSeconds &:= "° ";
    else
      if angle < 0.0 then
        degreesMinutesSeconds := "-";
      else
        degreesMinutesSeconds := " ";
      end if;
      degreesMinutesSeconds &:= abs(trunc(angle)) lpad 3;
      degreesMinutesSeconds &:= "° ";
      degreesMinutesSeconds &:= minutes(angle) lpad 2;
      degreesMinutesSeconds &:= "' ";
      degreesMinutesSeconds &:= seconds(angle) lpad 2;
      degreesMinutesSeconds &:= "\"";
    end if;
  end func;


const func string: angleAsHours (in float: angle) is func
  result
    var string: hoursMinutesSeconds is "";
  begin
    if displayDegreesInDecimal then
      hoursMinutesSeconds &:= angle / 15.0 digits 4 lpad 8;
      hoursMinutesSeconds &:= "h";
    else
      if angle < 0.0 then
        hoursMinutesSeconds := "-";
      else
        hoursMinutesSeconds := " ";
      end if;
      hoursMinutesSeconds &:= abs(trunc(angle / 15.0)) lpad 3;
      hoursMinutesSeconds &:= "h ";
      hoursMinutesSeconds &:= minutes(angle / 15.0) lpad 2;
      hoursMinutesSeconds &:= "m ";
      hoursMinutesSeconds &:= seconds(angle / 15.0) lpad 2;
      hoursMinutesSeconds &:= "s";
    end if;
  end func;


const func float: anomaly (in float: m, in float: eccentricity) is func
  result
    var float: anomaly is 0.0;
  local
    var float: d is 0.0;
    var float: e is 0.0;
  begin
    e := m;
    d := -eccentricity * sin(m);
    while abs(d) > 0.000001 do
      e := e - d / (1.0 - eccentricity * cos(e));
      d := e - eccentricity * sin(e) - m;
    end while;
    anomaly := 2.0 * atan(sqrt((1.0 + eccentricity) / (1.0 - eccentricity)) * tan(e / 2.0))
  end func; (* anomaly *)


const func float: anomaly2 (in float: m, in float: eccentricity) is func
  result
    var float: anomaly is 0.0;
  local
    var float: e is 0.0;
    var float: deltaM is 0.0;
    var float: deltaE is 0.0;
  begin
    e := m - eccentricity * sin(m);
    repeat
      deltaM := m - e - eccentricity * sin(e);
      deltaE := deltaM / (1.0 - eccentricity * cos(e));
      e +:= deltaE;
    until abs(deltaE) < 0.000001;
    anomaly := e;
  end func; (* anomaly2 *)


const func float: eclipticalCoordinatesToRightAscension (in float: l, in float: b) is func
  result
    var float: rightAscension is 0.0;
  local
    var float: x is 0.0;
    var float: y is 0.0;
  begin
    x := cosD(l);
    y := sinD(l) * 0.91746406 - tanD(b) * 0.397818676;
    rightAscension := atan2D(y, x) / 15.0;
  end func; (* eclipticalCoordinatesToRightAscension *)


const func float: eclipticalCoordinatesToDeclination (in float: l, in float: b) is
  return asinD(sinD(b) * 0.91746406 + cosD(b) * sinD(l) * 0.397818676);


const proc: degreesToRange (inout float: degrees) is func
  begin
    while degrees > 360.0 do
      degrees -:= 360.0;
    end while;
    while degrees < 0.0 do
      degrees +:= 360.0;
    end while;
  end func;


const proc: hoursToRange (inout float: hours) is func
  begin
    while hours >= 24.0 do
      hours -:= 24.0;
    end while;
    while hours <  0.0 do
      hours +:= 24.0;
    end while;
  end func;


const func float: localSiderealToCivilTime (in float: localSiderealTime,
    in float: longitude, in time: aTime) is func
  result
    var float: localCivilTime is 0.0
  local
    var float: t is 0.0;
    var float: b is 0.0;
    var float: greenwichSiderealTime is 0.0;
    var float: greenwichMeanTime is 0.0;
  begin
    greenwichSiderealTime := localSiderealTime + longitude / 15.0;
    if greenwichSiderealTime > 24.0 then
      greenwichSiderealTime -:= 24.0;
    end if;
    if greenwichSiderealTime < 0.0 then
      greenwichSiderealTime +:= 24.0;
    end if;
    t := (julianDay(truncToYear(aTime) - 1 . DAYS) - 2415020.0) / 36525.0;
    b := 24.0 - 6.6460656 - (2400.051262 * t) - (0.00002581 * t * t) +
        flt(24 * (aTime.year - 1900));
    t := greenwichSiderealTime - flt(dayOfYear(aTime)) * 0.0657098 + b;
    hoursToRange(t);
    greenwichMeanTime := t * 0.99727;
    localCivilTime := greenwichMeanTime + flt(aTime.timeZone) / 60.0;
  end func;


const proc: showLocalTime (inout text: win) is func
  local
    var string: timeStri is "";
    const integer: column is 1;
  begin
    timeStri := aTime lpad 19;
    setPos(win, 1, column);
    writeln(win, timeStri);
    color(win, light_cyan);
    case currTimeChange of
      when {1}:
        setPos(win, 1, column);
        write(win, timeStri[.. 4]);
      when {2}:
        setPos(win, 1, column + 5);
        write(win, timeStri[6 .. 7]);
      when {3}:
        setPos(win, 1, column + 8);
        write(win, timeStri[9 .. 10]);
      when {4}:
        setPos(win, 1, column + 11);
        write(win, timeStri[12 .. 13]);
      when {5}:
        setPos(win, 1, column + 14);
        write(win, timeStri[15 .. 16]);
      when {6}:
        setPos(win, 1, column + 17);
        write(win, timeStri[18 .. 19]);
      when {7}:
        setPos(win, 1, column + 20);
        write(win, timeStri[21 .. 26]);
    end case;
    color(win, white);
    writeln(win);
  end func;


const func float: calculateGreenwichMeanSiderealTime (in time: gmTime) is func
  result
    var float: greenwichMeanSiderealTime is 0.0;
  local
    var float: d is 0.0;
    var float: d0 is 0.0;
    var float: h is 0.0;
    var integer: t is 0;
  begin
    d := julianDay(gmTime) - julianDayOfEpoch;
    d0 := flt(julianDayNumber(gmTime)) + 0.5 - julianDayOfEpoch;
    h := (flt(gmTime.second) / 60.0 + flt(gmTime.minute)) / 60.0 + flt(gmTime.hour);
    t := round(d) div 36525;
    greenwichMeanSiderealTime := 6.697374558 + 0.06570982441908 * d0 +
                                 1.00273790935 * h + 0.000026 * flt(t * t);
    hoursToRange(greenwichMeanSiderealTime);
  end func;


const proc: showTime is func
  local
    var text: win is STD_NULL;
    var integer: julianDayNumber is 0;
    var float: greenwichMeanSiderealTime is 0.0;
    var time: gmTime is time.value;
  begin
    win := openWindow(scr, 2, 58, 8, 32);
    showLocalTime(win);
    julianDayNumber := julianDayNumber(aTime);
    writeln(win, "julianDayNum =" <& julianDayNumber lpad 10);
    writeln(win, "Longitude    = " <& angleAsDegrees(posAtEarth.longitude));
    writeln(win, "Latitude     = " <& angleAsDegrees(posAtEarth.latitude));
    write(win, "localCivilTime        = ");
    writeln(win, strTime(aTime));
    gmTime := toUTC(aTime);
    write(win, "greenwichMeanTime     = ");
    writeln(win, strTime(gmTime));
    greenwichMeanSiderealTime := calculateGreenwichMeanSiderealTime(gmTime);
    write(win, "greenwichSiderealTime = ");
    write(win, trunc(greenwichMeanSiderealTime) lpad0 2);
    write(win, ":");
    write(win, minutes(greenwichMeanSiderealTime) lpad0 2);
    write(win, ":");
    writeln(win, seconds(greenwichMeanSiderealTime) lpad0 2);
    localSiderealTime := greenwichMeanSiderealTime - posAtEarth.longitude / 15.0;
    hoursToRange(localSiderealTime);
    write(win, "localSiderealTime     = ");
    write(win, trunc(localSiderealTime) lpad0 2);
    write(win, ":");
    write(win, minutes(localSiderealTime) lpad0 2);
    write(win, ":");
    write(win, seconds(localSiderealTime) lpad0 2);
  end func; (* showTime *)


const func heliocentricPlanetType: locatePositionOfPlanetInItsOrbitalPlane (
    in orbitType: planetDescr) is func
  result
    var heliocentricPlanetType: planet is heliocentricPlanetType.value;
  begin
    planet.axis := planetDescr.axis +
        planetDescr.axisPerCentury * (julianDay(aTime) - julianDayOfEpoch) / 36525.0;
    planet.eccentricity := planetDescr.eccentricity +
        planetDescr.eccentricityPerCentury * (julianDay(aTime) - julianDayOfEpoch) / 36525.0;
    planet.meanLongitude := planetDescr.meanLongitudeAtEpoch +
        360.0 / (planetDescr.orbitalPeriod / orbitDescr[3].orbitalPeriod  * 365.2425) *
        (julianDay(aTime) - julianDayOfEpoch);
    degreesToRange(planet.meanLongitude);
    planet.longitudeOfPerihelion := planetDescr.longitudeOfPerihelion +
        planetDescr.longitudeOfPerihelionPerCentury * (julianDay(aTime) - julianDayOfEpoch) / 36525.0;
    planet.meanAnomaly := planet.meanLongitude - planet.longitudeOfPerihelion;
    planet.trueAnomaly := toDegrees(anomaly(toRadians(planet.meanAnomaly), planet.eccentricity));
    planet.heliocentricLongitude := planet.trueAnomaly + planet.longitudeOfPerihelion;
    degreesToRange(planet.heliocentricLongitude);
    planet.radiusVector := planetDescr.axis * (1.0 - planet.eccentricity * planet.eccentricity) /
        (1.0 + planet.eccentricity * cosD(planet.trueAnomaly));
    planet.heliocentricEclipticLatitude := asinD(sinD(planet.heliocentricLongitude -
        planetDescr.longitudeOfTheAscendingNode) * sinD(planetDescr.inclination));
  end func;


const proc: plotPlanet (in integer: planetNumber, in integer: scale) is func
  local
    const func float: computeRadiusVector (in heliocentricPlanetType: planet, in integer: angle) is
      return planet.axis *
             (1.0 - planet.eccentricity * planet.eccentricity) /
             (1.0 + planet.eccentricity * cosD(flt(angle) - planet.longitudeOfPerihelion));

    var heliocentricPlanetType: planet is heliocentricPlanetType.value;
    var integer: oldX is 0;
    var integer: oldY is 0;
    var integer: newX is 0;
    var integer: newY is 0;
    var integer: angle is 0;
  begin
    planet := locatePositionOfPlanetInItsOrbitalPlane(orbitDescr[planetNumber]);
    newY := round(flt(scale) * planet.radiusVector * sinD(planet.heliocentricLongitude)) +
            centerY;
    newX := round(flt(scale) * planet.radiusVector * cosD(planet.heliocentricLongitude) *
            screenAspect) + centerX;
    fcircle(newX, windowHeight - newY, 3, light_cyan);
    setPosXY(scr, newX + 4, windowHeight - newY);
    write(scr, orbitDescr[planetNumber].name);

    newX := centerX + round(computeRadiusVector(planet, 0) * flt(scale) * screenAspect);
    newY := windowHeight - centerY;
    for angle range 1 to 60 do
      oldX := newX;
      oldY := newY;
      newX := round(computeRadiusVector(planet, angle * 6) * flt(scale) * screenAspect
              * cosD(flt(angle * 6))) + centerX;
      newY := windowHeight -
              (round(computeRadiusVector(planet, angle * 6) * flt(scale) * sinD(flt(angle * 6)))
              + centerY);
      lineTo(oldX, oldY, newX, newY, light_cyan);
    end for;
  end func; (* plotPlanet *)


const proc: showMenu is func
  local
    var integer: num is 0;
  begin
    clear(black);
    showTime;
    setPos(scr, 2, 1);
    writeln(scr, " P L A N E T S   Display information about the planets");
    setPos(scr, 4, 1);
    writeln(scr, " Copyright (C) 2006, 2007  Thomas Mertes");
    setPos(scr, 6, 1);
    writeln(scr, " This program is free software under the");
    writeln(scr, " terms of the GNU General Public License");
    setPos(scr, 9, 1);
    writeln(scr, " Planets is written in the Seed7 programming language");
    writeln(scr, " Homepage:    http://seed7.sourceforge.net");
    writeln(scr);
    writeln(scr);
    writeln(scr, " Possible commands are:");
    writeln(scr);
    writeln(scr, " M  Menu");
    writeln(scr, " I  Plot Inner Planets");
    writeln(scr, " O  Plot Outer Planets");
    writeln(scr, " D  Change Date/Time");
    writeln(scr, " L  Change Long/Lat");
    writeln(scr, " S  Set_up");
    writeln(scr, " Q  Quit");
    for num range 1 to numberOfPlanets do
      setPos(scr, 14 + num, 25);
      writeln(scr, " " <& num <& "  " <& orbitDescr[num].name);
    end for;
    setPos(scr, 15, 42);
    writeln(scr, "Horizontal cursor keys:");
    setPos(scr, 16, 42);
    writeln(scr, "  Select time step");
    setPos(scr, 18, 42);
    writeln(scr, "Vertical cursor keys:");
    setPos(scr, 19, 42);
    write(scr,   "  Add or subtract a time step");
    setPos(scr, 25, 1);
    writeln(scr, " Enter command:");
  end func; (* showMenu *)


const proc: setup is func
  local
    var char: ch is ' ';
  begin
    clear(curr_win);
    setPos(scr, 1, 1);
    writeln(scr, "Display Degrees in Decimal Degrees (Y/N)? ");
    ch := getc(KEYBOARD);
    if upper(ch) = 'Y' then
      displayDegreesInDecimal := TRUE;
    elsif upper(ch) = 'N' then
      displayDegreesInDecimal := FALSE;
    end if;
    writeln(scr, "Display Hour Angle in Degrees (Y/N)? ");
    ch := getc(KEYBOARD);
    if upper(ch) = 'Y' then
      displayHourAngleInDegrees := TRUE;
    elsif upper(ch) = 'N' then
      displayHourAngleInDegrees := FALSE;
    end if;
    writeln(scr, "Display RA in Degrees (Y/N)? ");
    ch := getc(KEYBOARD);
    if upper(ch) = 'Y' then
      displayRAinDegrees := TRUE;
    elsif upper(ch) = 'N' then
      displayRAinDegrees := FALSE;
    end if;
  end func; (* setup *)


const proc: plotInnerPlanets is func
  local
    var text: win is STD_NULL;
    var integer: julianDayNumber is 0;
    var integer: planetNumber is 0;
  begin
    clear(curr_win);
    win := openWindow(scr, 2, 135, 7, 32);
    showLocalTime(win);
    julianDayNumber := julianDayNumber(aTime);
    writeln(win, "julianDayNum =" <& julianDayNumber lpad 10);
    writeln(win);
    writeln(win, "Horizontal cursor keys:");
    writeln(win, "  Select time step");
    writeln(win, "Vertical cursor keys:");
    write(win,   "  Add or subtract a time step");
    fcircle(centerX, windowHeight - centerY, 3, light_cyan);
    setPosXY(scr, centerX + 4, windowHeight - centerY);
    write(scr, "Sun");
    for planetNumber range 1 to 4 do
      plotPlanet(planetNumber, scaleForInnerPlanets);
    end for;
  end func; (* plotInnerPlanets *)


const proc: plotOuterPlanets is func
  local
    var text: win is STD_NULL;
    var integer: julianDayNumber is 0;
    var integer: planetNumber is 0;
  begin
    clear(curr_win);
    win := openWindow(scr, 2, 135, 7, 32);
    showLocalTime(win);
    julianDayNumber := julianDayNumber(aTime);
    writeln(win, "julianDayNum =" <& julianDayNumber lpad 10);
    writeln(win);
    writeln(win, "Horizontal cursor keys:");
    writeln(win, "  Select time step");
    writeln(win, "Vertical cursor keys:");
    write(win,   "  Add or subtract a time step");
    point(centerX, windowHeight - centerY, light_cyan);
    plotPlanet(3, scaleForOuterPlanets);
    for planetNumber range 5 to 9 do
      plotPlanet(planetNumber, scaleForOuterPlanets);
    end for;
  end func; (* plotOuterPlanets *)


const proc: changeDateTime is func
  local
    var text: win is STD_NULL;
  begin
    clear(curr_win);
    showTime;
    setPos(scr, 1, 1);
    case currTimeChange of
      when {1}:
        write(scr, "Enter year: ");
        read(aTime.year);
      when {2}:
        write(scr, "Enter month: ");
        read(aTime.month);
      when {3}:
        write(scr, "Enter day: ");
        read(aTime.day);
      when {4}:
        write(scr, "Enter hour: ");
        read(aTime.hour);
      when {5}:
        write(scr, "Enter minute: ");
        read(aTime.minute);
      when {6}:
        write(scr, "Enter second: ");
        read(aTime.second);
    end case;
    showMenu;
  end func;


const proc: changeLongLat is func
  begin
    clear(curr_win);
    showTime;
    setPos(scr, 1, 1);
    writeln("Enter Longitude");
    readln(posAtEarth.longitude);
    writeln("Enter Latitude");
    readln(posAtEarth.latitude);
    showMenu;
  end func;


const func float: computeProjectedHeliocentricLongitude (in heliocentricPlanetType: planet,
    in orbitType: orbitDescr) is func
  result
    var float: projectedHeliocentricLongitude is 0.0;
  local
    var float: y is 0.0;
    var float: x is 0.0;
    var float: z is 0.0;
  begin
    y := sinD(planet.heliocentricLongitude - orbitDescr.longitudeOfTheAscendingNode) *
         cosD(orbitDescr.inclination);
    x := cosD(planet.heliocentricLongitude - orbitDescr.longitudeOfTheAscendingNode);
    z := atan2D(y, x);
    projectedHeliocentricLongitude := z + orbitDescr.longitudeOfTheAscendingNode;
  end func;


const proc: ProjectPlanetOntoEclipticalPlane (inout heliocentricPlanetType: planet,
    in orbitType: orbitDescr) is func
  begin
    planet.projectedHeliocentricLongitude := computeProjectedHeliocentricLongitude(planet, orbitDescr);
    degreesToRange(planet.projectedHeliocentricLongitude);
    planet.projectedRadius := planet.radiusVector * cosD(planet.heliocentricEclipticLatitude);
  end func;


const func eclipticalCoordinateType: calculateEclipticalCoordinates (in heliocentricPlanetType: planet,
    in heliocentricPlanetType: earth) is func
  result
    var eclipticalCoordinateType: eclipticalCoordinate is eclipticalCoordinateType.value;
  begin
    if planet.radiusVector < earth.radiusVector then
      eclipticalCoordinate.geocentricEclipticLongitude := earth.heliocentricLongitude + 180.0 +
          atanD(planet.projectedRadius * sinD(earth.heliocentricLongitude -
          planet.projectedHeliocentricLongitude) / (earth.radiusVector -
          planet.projectedRadius * cosD(earth.heliocentricLongitude -
          planet.projectedHeliocentricLongitude)));
    else
      eclipticalCoordinate.geocentricEclipticLongitude := planet.projectedHeliocentricLongitude +
          atanD(earth.radiusVector * sinD(planet.projectedHeliocentricLongitude -
          earth.heliocentricLongitude) / (planet.projectedRadius -
          earth.radiusVector * cosD(planet.projectedHeliocentricLongitude -
          earth.heliocentricLongitude)));
    end if;
    degreesToRange(eclipticalCoordinate.geocentricEclipticLongitude);
    eclipticalCoordinate.geocentricLatitude := atanD(planet.projectedRadius *
        tanD(planet.heliocentricEclipticLatitude) *
        sinD(eclipticalCoordinate.geocentricEclipticLongitude -
        planet.projectedHeliocentricLongitude) / (earth.radiusVector *
        sinD(planet.projectedHeliocentricLongitude - earth.heliocentricLongitude)));
  end func;


const func equatorialCoordinateType: calculateEquatorialCoordinates (
    in eclipticalCoordinateType: eclipticalCoordinate) is func
  result
    var equatorialCoordinateType: equatorialCoordinate is equatorialCoordinateType.value;
  begin
    equatorialCoordinate.rightAscension := eclipticalCoordinatesToRightAscension(
        eclipticalCoordinate.geocentricEclipticLongitude,
        eclipticalCoordinate.geocentricLatitude);
    equatorialCoordinate.declination := eclipticalCoordinatesToDeclination(
        eclipticalCoordinate.geocentricEclipticLongitude,
        eclipticalCoordinate.geocentricLatitude);
  end func;


const func float: calculateHourAngle (in equatorialCoordinateType: equatorialCoordinate,
    in float: localSiderealTime) is func
  result
    var float: hourAngle is 0.0;
  begin
    hourAngle := localSiderealTime - equatorialCoordinate.rightAscension;
    if hourAngle < 0.0 then
      hourAngle +:= 24.0;
    end if;
  end func;


const func horizontalCoordinateType: calculateHorizontalCoordinates (
    in equatorialCoordinateType: equatorialCoordinate,
    in earthCoordinateType: posAtEarth, in float: hourAngle) is func
  result
    var horizontalCoordinateType: horizontalCoordinate is horizontalCoordinateType.value;
  begin
    horizontalCoordinate.altitude := asinD(sinD(equatorialCoordinate.declination) * sinD(posAtEarth.latitude) +
        cosD(equatorialCoordinate.declination) * cosD(posAtEarth.latitude) *
        cosD(hourAngle * 15.0));
    horizontalCoordinate.azimuth := acosD((sinD(equatorialCoordinate.declination) -
        sinD(posAtEarth.latitude) * sinD(horizontalCoordinate.altitude)) /
        (cosD(posAtEarth.latitude) * cosD(horizontalCoordinate.altitude)));
    if sinD(hourAngle * 15.0) > 0.0 then
      horizontalCoordinate.azimuth := 360.0 - horizontalCoordinate.azimuth;
    end if;
  end func;


const func geocentricPlanetType: genGeocentricPos (in heliocentricPlanetType: planet,
    in heliocentricPlanetType: earth, in float: localSiderealTime,
    in earthCoordinateType: posAtEarth) is func
  result
    var geocentricPlanetType: geocentricPos is geocentricPlanetType.value;
  begin
    geocentricPos.eclipticalCoordinate := calculateEclipticalCoordinates(planet, earth);
    geocentricPos.equatorialCoordinate := calculateEquatorialCoordinates(geocentricPos.eclipticalCoordinate);
    geocentricPos.hourAngle := calculateHourAngle(geocentricPos.equatorialCoordinate, localSiderealTime);
    geocentricPos.horizontalCoordinate := calculateHorizontalCoordinates(geocentricPos.equatorialCoordinate,
        posAtEarth, geocentricPos.hourAngle);
  end func;


const func integer: panoramaXPos (in float: azimuth) is
  return panoramaXMin + round(azimuth / 360.0 * sizePanorama);


const func integer: panoramaYPos (in float: altitude) is
  return panoramaHorizont - round(altitude  / 360.0 * sizePanorama);


const proc: drawAllPlanets (in heliocentricPlanetType: earth) is func
  local
    var integer: planetNumber is 0;
    var heliocentricPlanetType: planet is heliocentricPlanetType.value;
    var geocentricPlanetType: geocentricPos is geocentricPlanetType.value;
  begin
    for planetNumber range 1 to 9 do
      if planetNumber <> 3 then
        planet := locatePositionOfPlanetInItsOrbitalPlane(orbitDescr[planetNumber]);
        ProjectPlanetOntoEclipticalPlane(planet, orbitDescr[planetNumber]);
        geocentricPos := genGeocentricPos(planet, earth, localSiderealTime, posAtEarth);
        if geocentricPos.horizontalCoordinate.altitude >= 0.0 then
          fcircle(panoramaXPos(geocentricPos.horizontalCoordinate.azimuth),
              panoramaYPos(geocentricPos.horizontalCoordinate.altitude), 2, light_cyan);
        end if;
      end if;
    end for;
  end func;


const proc: drawSun (in heliocentricPlanetType: earth) is func
  local
    var heliocentricPlanetType: sun is heliocentricPlanetType.value;
    var geocentricPlanetType: geocentricPos is geocentricPlanetType.value;
  begin
    geocentricPos := genGeocentricPos(sun, earth, localSiderealTime, posAtEarth);
    if geocentricPos.horizontalCoordinate.altitude >= 0.0 then
      fcircle(panoramaXPos(geocentricPos.horizontalCoordinate.azimuth),
          panoramaYPos(geocentricPos.horizontalCoordinate.altitude), 4, yellow);
    end if;
  end func;


const proc: drawHorizont is func
  local
    var integer: azimuth is 0;
    var integer: length is 0;
    const array string: name is [] ("North", "East", "South", "West", "North");
    var integer: pos is 1;
  begin
    lineTo(0, panoramaHorizont, 1023, panoramaHorizont, white);
    for azimuth range 0 to 360 do
      if azimuth rem 90 = 0 then
        length := 12;
        setPosXY(scr, panoramaXPos(flt(azimuth)) - 12, panoramaHorizont + 25);
        writeln(scr, name[pos]);
        incr(pos);
      elsif azimuth rem 45 = 0 then
        length := 10;
      elsif azimuth rem 10 = 0 then
        length := 8;
      elsif azimuth rem 5 = 0 then
        length := 6;
      elsif azimuth rem 5 = 0 then
        length := 4;
      else
        length := 2;
      end if;
      line(panoramaXPos(flt(azimuth)), panoramaHorizont, 0, length, white);
    end for;
  end func;


const proc: writePoint (in integer: x, in integer: y, in float: vis) is func
  begin
    if vis >= 5.0 then
      point(x, y, dark_gray);
    elsif vis >= 4.5 then
      point(x, y, Gray);
    elsif vis >= 4.0 then
      point(x, y, light_gray);
    elsif vis >= 3.5 then
      point(x, y, white);
    elsif vis >= 3.0 then
      point(x, y, white);
      point(succ(x), y, dark_gray);
      point(pred(x), y, dark_gray);
      point(x, succ(y), dark_gray);
      point(x, pred(y), dark_gray);
    elsif vis >= 2.5 then
      point(x, y, white);
      point(succ(x), y, Gray);
      point(pred(x), y, Gray);
      point(x, succ(y), Gray);
      point(x, pred(y), Gray);
    elsif vis >= 2.0 then
      point(x, y, white);
      point(succ(x), y, light_gray);
      point(pred(x), y, light_gray);
      point(x, succ(y), light_gray);
      point(x, pred(y), light_gray);
    elsif vis >= 1.5 then
      line(pred(x), y, 2, 0, white);
      point(x, succ(y), white);
      point(x, pred(y), white);
    elsif vis >= 1.0 then
      line(pred(x), y, 2, 0, white);
      point(x, succ(y), white);
      point(x, pred(y), white);
      point(succ(x), succ(y), dark_gray);
      point(succ(x), pred(y), dark_gray);
      point(pred(x), succ(y), dark_gray);
      point(pred(x), pred(y), dark_gray);
    elsif vis >= 0.5 then
      line(pred(x), y, 2, 0, white);
      point(x, succ(y), white);
      point(x, pred(y), white);
      point(succ(x), succ(y), Gray);
      point(succ(x), pred(y), Gray);
      point(pred(x), succ(y), Gray);
      point(pred(x), pred(y), Gray);
    elsif vis >= 0.0 then
      line(pred(x), y, 2, 0, white);
      point(x, succ(y), white);
      point(x, pred(y), white);
      point(succ(x), succ(y), light_gray);
      point(succ(x), pred(y), light_gray);
      point(pred(x), succ(y), light_gray);
      point(pred(x), pred(y), light_gray);
    elsif vis >= -0.5 then
      rect (pred(x), pred(y), 3, 3, white);
    else
      fcircle(x, y, 2, white);
    end if;
  end func;


const proc: drawStars is func
  local
    var integer: index is 0;
    var equatorialCoordinateType: equatorialCoordinate is equatorialCoordinateType.value;
    var float: hourAngle is 0.0;
    var horizontalCoordinateType: horizontalCoordinate is horizontalCoordinateType.value;
    var integer: xPos is 0;
    var integer: yPos is 0;
  begin
    for index range length(stars) downto 1 do
      equatorialCoordinate.rightAscension := stars[index].rightAscension;
      equatorialCoordinate.declination := stars[index].declination;
      hourAngle := calculateHourAngle(equatorialCoordinate, localSiderealTime);
      horizontalCoordinate := calculateHorizontalCoordinates(equatorialCoordinate,
          posAtEarth, hourAngle);
      if horizontalCoordinate.altitude >= 0.0 then
        xPos := panoramaXPos(horizontalCoordinate.azimuth);
        yPos := panoramaYPos(horizontalCoordinate.altitude);
        writePoint(xPos, yPos, stars[index].magnitude);
        if xPos + round(sizePanorama) < windowWidth then
          writePoint(xPos + round(sizePanorama), yPos, stars[index].magnitude);
        end if;
        if xPos - round(sizePanorama) >= 0 then
          writePoint(xPos - round(sizePanorama), yPos, stars[index].magnitude);
        end if;
      end if;
    end for;
  end func;


const proc: locatePlanet (in integer: planetNumber) is func
  local
    var text: win is STD_NULL;
    var heliocentricPlanetType: planet is heliocentricPlanetType.value;
    var heliocentricPlanetType: earth is heliocentricPlanetType.value;
    var geocentricPlanetType: geocentricPos is geocentricPlanetType.value;
    var float: azimuthOfRise is 0.0;
    var float: azimuthOfSet is 0.0;
    var float: localSiderealTime_Rise is 0.0;
    var float: localSiderealTime_Set  is 0.0;
    var float: localCivilTime_Rise is 0.0;
    var float: localCivilTime_Set  is 0.0;
    var float: phase is 0.0;
    var float: distanceFromEarth is 0.0;
    var float: diameter is 0.0;
    var float: magnitude is 0.0;
  begin
    clear(curr_win, black);
    showTime;
    planet := locatePositionOfPlanetInItsOrbitalPlane(orbitDescr[planetNumber]);

    win := openWindow(scr, 2, 3, 7, 51);
    setPos(win, 1, 1);
    write  (win, "                          ");
    writeln(win, orbitDescr[planetNumber].name);
    writeln(win, "mean longitude       = " <& angleAsDegrees(planet.meanLongitude));
    writeln(win, "mean anomaly         = " <& angleAsDegrees(planet.meanAnomaly));
    writeln(win, "true anomaly         = " <& angleAsDegrees(planet.trueAnomaly));
    writeln(win, "heliocentric long    = " <& angleAsDegrees(planet.heliocentricLongitude));
    writeln(win, "radius vector        = " <& planet.radiusVector digits 3 lpad 8 <& " AU");
    write  (win, "helio ecliptic lat.  = " <& angleAsDegrees(planet.heliocentricEclipticLatitude));

    if planetNumber <> 3 then

      earth := locatePositionOfPlanetInItsOrbitalPlane(orbitDescr[3]);
      win := openWindow(scr, 2, 40, 7, 14);
      setPos(win, 1, 1);
      writeln(win, "     Earth");
      writeln(win, angleAsDegrees(earth.meanLongitude));
      writeln(win, angleAsDegrees(earth.meanAnomaly));
      writeln(win, angleAsDegrees(earth.trueAnomaly));
      writeln(win, angleAsDegrees(earth.heliocentricLongitude));
      writeln(win, earth.radiusVector digits 3 lpad 8 <& " AU");
      write  (win, angleAsDegrees(0.0));


      (* Project Position of Planet onto plane of ecliptic *)
      ProjectPlanetOntoEclipticalPlane(planet, orbitDescr[planetNumber]);

      win := openWindow(scr, 11, 12, 2, 68);
      setPos(win, 1, 1);
      write(win, "Projected Longitude  = ");
      write(win, angleAsDegrees(planet.projectedHeliocentricLongitude));
      writeln(win);
      write(win, "Projected Radius     = ");
      write(win, planet.projectedRadius digits 3 lpad 8);
      write(win, " AU");


      (* Calculate Ecliptical Coordinates *)
      geocentricPos.eclipticalCoordinate := calculateEclipticalCoordinates(planet, earth);

      win := openWindow(scr, 19, 3, 5, 22);
      setPos(win, 1, 1);
      writeln(win, "      ECLIPTIC");
      writeln(win);
      writeln(win, " Long = " <& angleAsDegrees(
          geocentricPos.eclipticalCoordinate.geocentricEclipticLongitude));
      write(win, " Lat  = " <& angleAsDegrees(
          geocentricPos.eclipticalCoordinate.geocentricLatitude));


      (* Calculate Equatorial Coordinates *)
      geocentricPos.equatorialCoordinate := calculateEquatorialCoordinates(
          geocentricPos.eclipticalCoordinate);

      win := openWindow(scr, 19, 29, 5, 23);
      setPos(win, 1, 1);
      writeln(win, "     EQUATORIAL");
      writeln(win);
      write(win, " RA  = ");
      if displayRAinDegrees then
        writeln(win, angleAsDegrees(geocentricPos.equatorialCoordinate.rightAscension * 15.0));
      else
        writeln(win, angleAsHours(geocentricPos.equatorialCoordinate.rightAscension * 15.0));
      end if;
      writeln(win, " DEC = " <& angleAsDegrees(geocentricPos.equatorialCoordinate.declination));


      geocentricPos.hourAngle := calculateHourAngle(geocentricPos.equatorialCoordinate, localSiderealTime);
      write(win, " HA  = ");
      if displayHourAngleInDegrees then
        write(win, angleAsDegrees(geocentricPos.hourAngle * 15.0));
      else
        write(win, angleAsHours(geocentricPos.hourAngle * 15.0));
      end if;


      (* Calculate Horizontal Coordinates *)
      geocentricPos.horizontalCoordinate := calculateHorizontalCoordinates(
          geocentricPos.equatorialCoordinate,
          posAtEarth, geocentricPos.hourAngle);

      win := openWindow(scr, 19, 56, 5, 23);
      setPos(win, 1, 1);
      writeln(win, "      HORIZONTAL");
      writeln(win);
      writeln(win, " Alt  = " <& angleAsDegrees(geocentricPos.horizontalCoordinate.altitude));
      write(win, " Azim = " <& angleAsDegrees(geocentricPos.horizontalCoordinate.azimuth));


      (* Calculate Time and Position of Rise and Set *)

      win := openWindow(scr, 15, 4, 2, 50);
      setPos(win, 1, 1);
      azimuthOfRise := sinD(geocentricPos.equatorialCoordinate.declination) / cosD(posAtEarth.latitude);
      if azimuthOfRise < -1.0 or azimuthOfRise > 1.0 then
        if geocentricPos.horizontalCoordinate.altitude < 0.0 then
          writeln(win, "never rises above horizon");
        else
          writeln(win, "never sets below horizon");
        end if;
      else
        azimuthOfRise := acosD(azimuthOfRise);
        azimuthOfSet := 360.0 - azimuthOfRise;
        localSiderealTime_Rise := 24.0 + geocentricPos.equatorialCoordinate.rightAscension -
            (acosD(-tanD(posAtEarth.latitude) *
            tanD(geocentricPos.equatorialCoordinate.declination))) / 15.0;
        if localSiderealTime_Rise > 24.0 then
          localSiderealTime_Rise := localSiderealTime_Rise - 24.0;
        end if;
        localSiderealTime_Set := geocentricPos.equatorialCoordinate.rightAscension +
            (acosD(-tanD(posAtEarth.latitude) *
            tanD(geocentricPos.equatorialCoordinate.declination))) / 15.0;
        if localSiderealTime_Set > 24.0 then
          localSiderealTime_Set := localSiderealTime_Set - 24.0;
        end if;

        localCivilTime_Set := localSiderealToCivilTime(localSiderealTime_Set,
            posAtEarth.longitude, aTime);
        if localCivilTime_Set < 0.0 then
          localCivilTime_Set := localCivilTime_Set + 24.0;
        end if;

        localCivilTime_Rise := localSiderealToCivilTime(localSiderealTime_Rise,
            posAtEarth.longitude, aTime);
        if localCivilTime_Rise < 0.0 then
          localCivilTime_Rise := localCivilTime_Rise + 24.0;
        end if;

        write(win, "Rises at ");
        write(win, trunc(localCivilTime_Rise) lpad 2 <& ":" <&
            minutes(localCivilTime_Rise) lpad 2 <& " local time");
        write(win, "  Azimuth ");
        write(win, angleAsDegrees(azimuthOfRise));
        writeln(win);
        write(win, "Sets at  ");
        write(win, trunc(localCivilTime_Set) lpad 2 <& ":" <&
            minutes(localCivilTime_Set) lpad 2 <& " local time");
        write(win, "  Azimuth ");
        write(win, angleAsDegrees(azimuthOfSet));

        line(panoramaXPos(azimuthOfRise), panoramaHorizont - 4, 0, 4, light_red);
        line(panoramaXPos(azimuthOfSet), panoramaHorizont - 4, 0, 4, light_red);
      end if;


      (* Calculate Phase, Distance, Diameter, Magnitude *)

      win := openWindow(scr, 12, 57, 5, 23);
      setPos(win, 1, 1);

      phase := (1.0 + cosD(geocentricPos.eclipticalCoordinate.geocentricEclipticLongitude -
          planet.heliocentricLongitude)) / 2.0;
      write(win, "Phase     = ");
      write(win, 100.0 * phase digits 2 lpad 6);
      writeln(win, "%");
      distanceFromEarth := sqrt(earth.radiusVector * earth.radiusVector
                             + planet.radiusVector * planet.radiusVector
                             - 2.0 * earth.radiusVector * planet.radiusVector *
                             cosD(planet.heliocentricLongitude -
                             earth.heliocentricLongitude));
      write(win, "Distance  = ");
      write(win, distanceFromEarth digits 2 lpad 6);
      writeln(win, " AU");

      diameter := orbitDescr[planetNumber].size / distanceFromEarth;
      write(win, "Diameter  = ");
      write(win, diameter digits 2 lpad 6);
      writeln(win, "\"");

      write(win, "magnitude = ");
      if distanceFromEarth * planet.radiusVector /
                   orbitDescr[planetNumber].brightness * sqrt(phase) <> 0.0 then
        magnitude := 2.17147 * log(distanceFromEarth * planet.radiusVector /
                     orbitDescr[planetNumber].brightness * sqrt(phase)) - 26.7;
        write(win, magnitude digits 2 lpad 6);
      else
        write(win, "*****" lpad 6);
      end if;

      drawHorizont;
      drawStars;
      if geocentricPos.horizontalCoordinate.altitude >= 0.0 then
        fcircle(panoramaXPos(geocentricPos.horizontalCoordinate.azimuth),
            panoramaYPos(geocentricPos.horizontalCoordinate.altitude), 2, light_cyan);
      end if;
      drawSun(earth);
    else
      drawHorizont;
      drawStars;
      drawAllPlanets(planet);
      drawSun(planet);
    end if;
  end func; (* locatePlanet *)


const proc: showName is func
  local
    var integer: mouseXPos is 0;
    var integer: mouseYPos is 0;
    var integer: index is 0;
    var equatorialCoordinateType: equatorialCoordinate is equatorialCoordinateType.value;
    var float: hourAngle is 0.0;
    var horizontalCoordinateType: horizontalCoordinate is horizontalCoordinateType.value;
    var heliocentricPlanetType: sun is heliocentricPlanetType.value;
    var heliocentricPlanetType: earth is heliocentricPlanetType.value;
    var heliocentricPlanetType: planet is heliocentricPlanetType.value;
    var geocentricPlanetType: geocentricPos is geocentricPlanetType.value;
    var integer: xPos is 0;
    var integer: yPos is 0;
    var integer: delta is 0;
    var integer: minDelta is 0;
    var integer: selectedStar is 0;
    var integer: selectedPlanet is 0;
    var string: name is "";
  begin
    mouseXPos := getxpos(KEYBOARD);
    mouseYPos := getypos(KEYBOARD);
    for index range 1 to length(stars) do
      equatorialCoordinate.rightAscension := stars[index].rightAscension;
      equatorialCoordinate.declination := stars[index].declination;
      hourAngle := calculateHourAngle(equatorialCoordinate, localSiderealTime);
      horizontalCoordinate := calculateHorizontalCoordinates(equatorialCoordinate,
          posAtEarth, hourAngle);
      if horizontalCoordinate.altitude >= 0.0 then
        xPos := panoramaXPos(horizontalCoordinate.azimuth);
        yPos := panoramaYPos(horizontalCoordinate.altitude);
        delta := (mouseXPos - xPos) ** 2 + (mouseYPos - yPos) ** 2;
        if selectedStar = 0 or delta < minDelta then
          selectedStar := index;
          minDelta := delta;
        end if;
        if xPos + round(sizePanorama) < windowWidth then
          delta := (mouseXPos - xPos - round(sizePanorama)) ** 2 + (mouseYPos - yPos) ** 2;
          if selectedStar = 0 or delta < minDelta then
            selectedStar := index;
            minDelta := delta;
          end if;
        end if;
        if xPos - round(sizePanorama) >= 0 then
          delta := (mouseXPos - xPos + round(sizePanorama)) ** 2 + (mouseYPos - yPos) ** 2;
          if selectedStar = 0 or delta < minDelta then
            selectedStar := index;
            minDelta := delta;
          end if;
        end if;
      end if;
    end for;
    earth := locatePositionOfPlanetInItsOrbitalPlane(orbitDescr[3]);
    for index range 1 to 9 do
      if index <> 3 then
        planet := locatePositionOfPlanetInItsOrbitalPlane(orbitDescr[index]);
        ProjectPlanetOntoEclipticalPlane(planet, orbitDescr[index]);
        geocentricPos := genGeocentricPos(planet, earth, localSiderealTime, posAtEarth);
        if geocentricPos.horizontalCoordinate.altitude >= 0.0 then
          xPos := panoramaXPos(geocentricPos.horizontalCoordinate.azimuth);
          yPos := panoramaYPos(geocentricPos.horizontalCoordinate.altitude);
          delta := (mouseXPos - xPos) ** 2 + (mouseYPos - yPos) ** 2;
          if delta <= minDelta then
            selectedStar := 0;
            selectedPlanet := index;
            minDelta := delta;
          end if;
        end if;
      end if;
    end for;
    geocentricPos := genGeocentricPos(sun, earth, localSiderealTime, posAtEarth);
    if geocentricPos.horizontalCoordinate.altitude >= 0.0 then
      xPos := panoramaXPos(geocentricPos.horizontalCoordinate.azimuth);
      yPos := panoramaYPos(geocentricPos.horizontalCoordinate.altitude);
      delta := (mouseXPos - xPos) ** 2 + (mouseYPos - yPos) ** 2;
      if delta <= minDelta then
        selectedStar := 0;
        selectedPlanet := 0;
        minDelta := delta;
      end if;
    end if;
    if selectedStar <> 0 then
      equatorialCoordinate.rightAscension := stars[selectedStar].rightAscension;
      equatorialCoordinate.declination := stars[selectedStar].declination;
      hourAngle := calculateHourAngle(equatorialCoordinate, localSiderealTime);
      horizontalCoordinate := calculateHorizontalCoordinates(equatorialCoordinate,
          posAtEarth, hourAngle);
      xPos := panoramaXPos(horizontalCoordinate.azimuth);
      yPos := panoramaYPos(horizontalCoordinate.altitude);
      name := stars[selectedStar].name;
    elsif selectedPlanet <> 0 then
      planet := locatePositionOfPlanetInItsOrbitalPlane(orbitDescr[selectedPlanet]);
      ProjectPlanetOntoEclipticalPlane(planet, orbitDescr[selectedPlanet]);
      geocentricPos := genGeocentricPos(planet, earth, localSiderealTime, posAtEarth);
      xPos := panoramaXPos(geocentricPos.horizontalCoordinate.azimuth);
      yPos := panoramaYPos(geocentricPos.horizontalCoordinate.altitude);
      name := orbitDescr[selectedPlanet].name;
    else
      geocentricPos := genGeocentricPos(sun, earth, localSiderealTime, posAtEarth);
      xPos := panoramaXPos(geocentricPos.horizontalCoordinate.azimuth);
      yPos := panoramaYPos(geocentricPos.horizontalCoordinate.altitude);
      name := "Sun";
    end if;
    if minDelta <= 900 then
      setPosXY(scr, xPos + 4, yPos);
      write(name);
      if xPos + round(sizePanorama) < windowWidth then
        setPosXY(scr, xPos + round(sizePanorama) + 4, yPos);
        write(name);
      end if;
      if xPos - round(sizePanorama) >= 0 then
        setPosXY(scr, xPos - round(sizePanorama) + 4, yPos);
        write(name);
      end if;
    end if;
  end func;


const proc: initOrbitDescr is func
  begin
    orbitDescr[1].name :="Mercury";
    orbitDescr[1].orbitalPeriod := 87.96934; # 0.2408469 years
    orbitDescr[1].meanLongitudeAtEpoch := 252.25032350;
    orbitDescr[1].longitudeOfPerihelion := 77.45779628;
    orbitDescr[1].longitudeOfPerihelionPerCentury := 0.16047689;
    orbitDescr[1].eccentricity := 0.20563593;
    orbitDescr[1].eccentricityPerCentury := 0.00001906;
    orbitDescr[1].axis := 0.38709927;
    orbitDescr[1].axisPerCentury := 0.00000037;
    orbitDescr[1].inclination := 7.00497902;
    orbitDescr[1].longitudeOfTheAscendingNode := 48.33167;
    orbitDescr[1].argumentOfThePerihelion := 29.12478;
    orbitDescr[1].size := 6.74;
    orbitDescr[1].brightness := 0.000001918;

    orbitDescr[2].name :="Venus";
    orbitDescr[2].orbitalPeriod := 224.70069; # 0.6151970 years
    orbitDescr[2].meanLongitudeAtEpoch := 181.97909950;
    orbitDescr[2].longitudeOfPerihelion := 131.60246718;
    orbitDescr[2].longitudeOfPerihelionPerCentury := 0.00268329;
    orbitDescr[2].eccentricity := 0.00677672;
    orbitDescr[2].eccentricityPerCentury := -0.00004107;
    orbitDescr[2].axis := 0.72333566;
    orbitDescr[2].axisPerCentury := 0.00000390;
    orbitDescr[2].inclination := 3.39467605;
    orbitDescr[2].longitudeOfTheAscendingNode := 76.68069;
    orbitDescr[2].argumentOfThePerihelion := 54.85229;
    orbitDescr[2].size := 16.92;
    orbitDescr[2].brightness := 0.00001721;

    orbitDescr[3].name :="Earth";
    orbitDescr[3].orbitalPeriod := 365.256366; # 1.0000175 years
    orbitDescr[3].meanLongitudeAtEpoch := 100.46457166;
    orbitDescr[3].longitudeOfPerihelion := 102.93768193;
    orbitDescr[3].longitudeOfPerihelionPerCentury := 0.32327364;
    orbitDescr[3].eccentricity := 0.01671123;
    orbitDescr[3].eccentricityPerCentury := -0.00004392;
    orbitDescr[3].axis := 1.00000261;
    orbitDescr[3].axisPerCentury := 0.00000562;
    orbitDescr[3].inclination := 0.00001531;
    orbitDescr[3].longitudeOfTheAscendingNode := 348.73936;
    orbitDescr[3].argumentOfThePerihelion := 114.20783;
    orbitDescr[3].size := 17.0;
    orbitDescr[3].brightness := 0.0;

    orbitDescr[4].name :="Mars";
    orbitDescr[4].orbitalPeriod := 686.9600; # 1.8808 years
    orbitDescr[4].meanLongitudeAtEpoch := 355.44656795;
    orbitDescr[4].longitudeOfPerihelion := 336.05637041;
    orbitDescr[4].longitudeOfPerihelionPerCentury := 0.44441088;
    orbitDescr[4].eccentricity := 0.09339410;
    orbitDescr[4].eccentricityPerCentury := 0.00007882;
    orbitDescr[4].axis := 1.52371034;
    orbitDescr[4].axisPerCentury := 0.00001847;
    orbitDescr[4].inclination := 1.84969142;
    orbitDescr[4].longitudeOfTheAscendingNode := 49.57854;
    orbitDescr[4].argumentOfThePerihelion := 286.46230;
    orbitDescr[4].size := 9.36;
    orbitDescr[4].brightness := 0.00000454;

    orbitDescr[5].name :="Jupiter";
    orbitDescr[5].orbitalPeriod := 4333.2867; # 11.86 years
    orbitDescr[5].meanLongitudeAtEpoch := 34.39644051;
    orbitDescr[5].longitudeOfPerihelion := 14.72847983;
    orbitDescr[5].longitudeOfPerihelionPerCentury := 0.21252668;
    orbitDescr[5].eccentricity := 0.04838624;
    orbitDescr[5].eccentricityPerCentury := -0.00013253;
    orbitDescr[5].axis := 5.20288700;
    orbitDescr[5].axisPerCentury := -0.00011607;
    orbitDescr[5].inclination := 1.30439695;
    orbitDescr[5].longitudeOfTheAscendingNode := 100.55615;
    orbitDescr[5].argumentOfThePerihelion := 274.19770;
    orbitDescr[5].size := 196.74;
    orbitDescr[5].brightness := 0.0001994;

    orbitDescr[6].name :="Saturn";
    orbitDescr[6].orbitalPeriod := 10756.1995; # 29.45 years
    orbitDescr[6].meanLongitudeAtEpoch := 49.95424423;
    orbitDescr[6].longitudeOfPerihelion := 92.59887831;
    orbitDescr[6].longitudeOfPerihelionPerCentury := -0.41897216;
    orbitDescr[6].eccentricity := 0.05386179;
    orbitDescr[6].eccentricityPerCentury := -0.00050991;
    orbitDescr[6].axis := 9.53667594;
    orbitDescr[6].axisPerCentury := -0.00125060;
    orbitDescr[6].inclination := 2.48599187;
    orbitDescr[6].longitudeOfTheAscendingNode := 113.71504;
    orbitDescr[6].argumentOfThePerihelion := 338.71690;
    orbitDescr[6].size := 165.60;
    orbitDescr[6].brightness := 0.0001740;

    orbitDescr[7].name :="Uranus";
    orbitDescr[7].orbitalPeriod := 30707.4896; # 84.07 years
    orbitDescr[7].meanLongitudeAtEpoch := 313.23810451;
    orbitDescr[7].longitudeOfPerihelion := 170.95427630;
    orbitDescr[7].longitudeOfPerihelionPerCentury := 0.40805281;
    orbitDescr[7].eccentricity := 0.04725744;
    orbitDescr[7].eccentricityPerCentury := -0.00004397;
    orbitDescr[7].axis := 19.18916464;
    orbitDescr[7].axisPerCentury := -0.00196176;
    orbitDescr[7].inclination := 0.77263783;
    orbitDescr[7].longitudeOfTheAscendingNode := 74.22988;
    orbitDescr[7].argumentOfThePerihelion := 96.73436;
    orbitDescr[7].size := 65.80;
    orbitDescr[7].brightness := 0.00007768;

    orbitDescr[8].name :="Neptune";
    orbitDescr[8].orbitalPeriod := 60223.3528; # 164.88 years
    orbitDescr[8].meanLongitudeAtEpoch := 304.87997031;
    orbitDescr[8].longitudeOfPerihelion := 44.96476227;
    orbitDescr[8].longitudeOfPerihelionPerCentury := -0.32241464;
    orbitDescr[8].eccentricity := 0.00859048;
    orbitDescr[8].eccentricityPerCentury := 0.00005105;
    orbitDescr[8].axis := 30.06992276;
    orbitDescr[8].axisPerCentury := 0.00026291;
    orbitDescr[8].inclination := 1.77004347;
    orbitDescr[8].longitudeOfTheAscendingNode := 131.72169;
    orbitDescr[8].argumentOfThePerihelion := 273.24966;
    orbitDescr[8].size := 62.20;
    orbitDescr[8].brightness := 0.00007597;

    orbitDescr[9].name :="Pluto";
    orbitDescr[9].orbitalPeriod := 90613.3055; # 248.09 years
    orbitDescr[9].meanLongitudeAtEpoch := 238.92903833;
    orbitDescr[9].longitudeOfPerihelion := 224.06891629;
    orbitDescr[9].longitudeOfPerihelionPerCentury := -0.04062942;
    orbitDescr[9].eccentricity := 0.24882730;
    orbitDescr[9].eccentricityPerCentury := 0.00005170;
    orbitDescr[9].axis := 39.48211675;
    orbitDescr[9].axisPerCentury := -0.00031596;
    orbitDescr[9].inclination := 17.14001206;
    orbitDescr[9].longitudeOfTheAscendingNode := 110.30347;
    orbitDescr[9].argumentOfThePerihelion := 113.76329;
    orbitDescr[9].size := 8.2;
    orbitDescr[9].brightness := 0.000004073;
  end func; (* initOrbitDescr *)


const proc: main is func
  local
    var char: ch is ' ';
    var char: command is 'M';
    var char: old_command is 'M';
  begin
    screen(windowWidth, windowHeight);
    selectInput(curr_win, KEY_CLOSE, TRUE);
    clear(curr_win, black);
    KEYBOARD := GRAPH_KEYBOARD;
    scr := open(curr_win);
    IN := openEcho(KEYBOARD, scr);
    IN := openLine(IN);
    OUT := scr;
    aTime := truncToSecond(time(NOW));
    initOrbitDescr;

    showMenu;
    repeat
      ch := getc(KEYBOARD);
      if ch = KEY_RIGHT then
        if currTimeChange < 7 then
          incr(currTimeChange);
        end if;
      elsif ch = KEY_LEFT then
        if currTimeChange > 1 then
          decr(currTimeChange);
        end if;
      elsif ch = KEY_UP then
        case currTimeChange of
          when {1}: aTime -:= 1 . YEARS;
          when {2}: aTime -:= 1 . MONTHS;
          when {3}: aTime -:= 1 . DAYS;
          when {4}: aTime -:= 1 . HOURS;
          when {5}: aTime -:= 1 . MINUTES;
          when {6}: aTime -:= 1 . SECONDS;
          when {7}: if aTime.timeZone > -720 then
                      aTime.timeZone -:= 60;
                    end if;
        end case;
      elsif ch = KEY_DOWN then
        case currTimeChange of
          when {1}: aTime +:= 1 . YEARS;
          when {2}: aTime +:= 1 . MONTHS;
          when {3}: aTime +:= 1 . DAYS;
          when {4}: aTime +:= 1 . HOURS;
          when {5}: aTime +:= 1 . MINUTES;
          when {6}: aTime +:= 1 . SECONDS;
          when {7}: if aTime.timeZone < 840 then
                      aTime.timeZone +:= 60;
                    end if;
        end case;
      else
        old_command := command;
        command := ch;
      end if;
      case command of
        when {'M', 'm'}: showMenu;
        when {'I', 'i'}: plotInnerPlanets;
        when {'O', 'o'}: plotOuterPlanets;
        when {'D', 'd'}: changeDateTime;
                         command := 'M';
        when {'L', 'l'}: changeLongLat;
                         command := 'M';
        when {'S', 's'}: setup;
        when {'1', '2', '3', '4', '5', '6', '7', '8', '9'}:
          locatePlanet(ord(command) - ord('0'));
        when {KEY_MOUSE1}:
          if old_command in {'1', '2', '3', '4', '5', '6', '7', '8', '9'} then
            showName;
          end if;
          command := old_command;
        otherwise:
          command := old_command;
      end case;
      setPos(scr, 58, 1);
      write(scr, "(M)enu (I)nner (O)uter (D)ate (L)ong (Q)uit  (1-9) Planet");
    until ch = 'Q' or ch = 'q' or ch = KEY_CLOSE;
  end func;