Astronomy App for motion of planets around sun.

[Alchemisst]

Below is a astronomy program which depicts motion of planets around the sun.It also includes Halley comet.Also,there is also included a simulation of a rocket which can travel from one planet to another.

EXPORT SolarSystem()
BEGIN
// gravitational constant
LOCAL G := 6.6743E-11;

// masses of the sun and the planets (in kg)
LOCAL M := {1.989E30, 3.301E23, 4.867E24, 5.972E24, 6.39E23, 1.898E27, 5.683E26, 8.681E25};

// radii of the sun and the planets (in meters)
LOCAL R := {696340000, 2439700, 6051800, 6371000, 3396200, 69911000, 58232000, 25362000};

// initial positions and velocities of the sun and the planets (in meters and meters per second)
LOCAL X0 := {0, 57910000000, 108200000000, 149600000000, 227900000000, 778500000000, 1427000000000, 2871000000000};
LOCAL Y0 := {0, 0, 0, 0, 0, 0, 0, 0};
LOCAL VX0 := {0, 0, 0, 0, 0, 0, 0, 0};
LOCAL VY0 := {0, 47890, 35020, 29780, 24130, 13060, 9690, 6830};

// initial position and velocity of Halley's Comet (in meters and meters per second)
LOCAL XH := 8.783E10;
LOCAL YH := 5.774E10;
LOCAL VXH := -54120;
LOCAL VYH := 34140;

// initialize plot
PICTSTART();
RECT_P(0, 0, 320, 240);
PICTEND();
PICTSTART();
AXIS(-3E12, -3E12, 3E12, 3E12);
PICTEND();

// loop for simulating motion
LOCAL X := X0; // current x positions of the sun and the planets
LOCAL Y := Y0; // current y positions of the sun and the planets
LOCAL VX := VX0; // current x velocities of the sun and the planets
LOCAL VY := VY0; // current y velocities of the sun and the planets
LOCAL XHC := XH; // current x position of Halley's Comet
LOCAL YHC := YH; // current y position of Halley's Comet
LOCAL VXHC := VXH; // current x velocity of Halley's Comet
LOCAL VYHC := VYH; // current y velocity of Halley's Comet
LOCAL DT := 86400; // time step (1 day)
LOCAL T := 0; // current time
LOCAL N := SIZE(M); // number of celestial bodies
LOCAL D := {0, 0, 0, 0, 0, 0, 0, 0}; // distances between the sun and the planets
LOCAL Fx := {0, 0, 0, 0, 0, 0, 0, 0}; // x components of gravitational forces on the sun and the planets
LOCAL Fy := {0, 0, 0, 0, 0, 0, 0, 0}; // y components of gravitational forces on the sun and the planets

// loop for simulating motion
REPEAT
// calculate distances and gravitational forces
FOR i FROM 2 TO N DO
D(i) := SQRT((X(i)-X(1))^2 + (Y(i)-Y(1))^2);
Fx(i) := -G*M(1)*M(i)*(X(i)-X(1))/D(i)^3;
Fy(i) := -G*M(1)*M(i)*(Y(i)-Y(1))/D(i)^3;
Fx(1) := Fx(1) - Fx(i);
Fy(1) := Fy(1) - Fy(i);
END;

// calculate distances and gravitational forces for Halley's Comet
LOCAL DH := SQRT((XHC-X(1))^2 + (YHC-Y(1))^2);
LOCAL FXH := -G*M(1)*M(N+1)*(XHC-X(1))/DH^3;
LOCAL FYH := -G*M(1)*M(N+1)*(YHC-Y(1))/DH^3;
Fx(1) := Fx(1) + FXH;
Fy(1) := Fy(1) + FYH;

// calculate accelerations, velocities, and positions
FOR i FROM 1 TO N+1 DO
VX(i) := VX(i) + Fx(i)/M(i)*DT;
VY(i) := VY(i) + Fy(i)/M(i)*DT;
IF i == N+1 THEN
VXHC := VXHC + FXH/M(N+1)*DT;
VYHC := VYHC + FYH/M(N+1)*DT;
XHC := XHC + VXHC*DT;
YHC := YHC + VYHC*DT;
ELSE
X(i) := X(i) + VX(i)*DT;
Y(i) := Y(i) + VY(i)*DT;
END;
END;

// plot positions and velocities
PICTSTART();
FOR i FROM 1 TO N+1 DO
POINT_P(X(i), Y(i));
LINE_P(X(i), Y(i), X(i)+VX(i)/1E5, Y(i)+VY(i)/1E5);
END;
POINT_P(XHC, YHC);
LINE_P(XHC, YHC, XHC+VXHC/1E5, YHC+VYHC/1E5);
PICTEND();

// update time
T := T + DT;
UNTIL false;
END;

EXPORT LaunchRocket(PlanetIndex)
BEGIN
// gravitational constant
LOCAL G := 6.6743E-11;

// masses of the sun and the planets (in kg)
LOCAL M := {1.989E30, 3.301E23, 4.867E24, 5.972E24, 6.39E23, 1.898E27, 5.683E26, 8.681E25};

// radii of the sun and the planets (in meters)
LOCAL R := {696340000, 2439700, 6051800, 6371000, 3396200, 69911000, 58232000, 25362000};

// initial positions and velocities of the sun and the planets (in meters and meters per second)
LOCAL X0 := {0, 57910000000, 108200000000, 149600000000, 227900000000, 778500000000, 1427000000000, 2871000000000};
LOCAL Y0 := {0, 0, 0, 0, 0, 0, 0, 0};
LOCAL VX0 := {0, 0, 0, 0, 0, 0, 0, 0};
LOCAL VY0 := {0, 47890, 35020, 29780, 24130, 13060, 9690, 6830};

// gravitational force on the rocket
LOCAL FR := G*M(PlanetIndex)*1E3/R(PlanetIndex)^2;

// initial position and velocity of the rocket (on Earth)
LOCAL XR := X0(4) + R(4);
LOCAL YR := Y0(4);
LOCAL VXR := 0;
LOCAL VYR := SQRT(2*FR);

// initialize plot
PICTSTART();
RECT_P(0, 0, 320, 240);
PICTEND();
PICTSTART();
AXIS(-3E12, -3E12, 3E12, 3E12);
PICTEND();

// loop for simulating motion
LOCAL X := X0; // current x positions of the sun and the planets
LOCAL Y := Y0; // current y positions of the sun and the planets
LOCAL VX := VX0; // current x velocities of the sun and the planets
LOCAL VY := VY0; // current y velocities of the sun and the planets
LOCAL DT := 86400; // time step (1 day)
LOCAL T := 0; // current time
LOCAL N := SIZE(M); // number of celestial bodies
LOCAL D := {0, 0, 0, 0, 0, 0, 0, 0}; // distances between the sun and the planets
LOCAL Fx := {0, 0, 0, 0, 0, 0, 0

[Alchemisst]

Advanced App includes rocket simulation travelling from one planet to another and I will soon release an app showing live time map of rockets and asteroids around it.