# N-body problem

N-body problem is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

The ${\displaystyle N}$-body problem is an inherent feature of a physical description of mechanical system of ${\displaystyle N}$ interacting objects. The system is described by Newton's Law

${\displaystyle {\vec {F}}=m_{i}{\frac {d^{2}{\vec {r}}_{i}}{dt^{2}}},i=1..N}$

with continious combined force from other bodies

${\displaystyle F_{i}=\sum _{j\neq i}^{}F_{ij},i=1..N}$

Exact formulation of first mechanical problem is, given initial coordinates and velocities ${\displaystyle t=0}$, to find coordinates and velocities at ${\displaystyle t=t_{0}}$ with accuracy more then given relative value ${\displaystyle \epsilon ={\frac {{\vec {r}}-{\vec {r}}'}{|{\vec {r}}|}}}$

As well known from physical background, only the choice of ${\displaystyle N=2}$ can be analytically solved.

Write a simulation of three masses interacting under gravitation and show their evolution over a number of time-steps (at least 20).

## C

### A little history

A long, long time ago ( well, more than 10 years, which is a very long time when it comes to computing ), I implemented a BGI program in C for the Mitchell-Green Gravity Set. A chaotic dynamical system, the Gravity set, MGGS in short, is a plot of the paths of massive bodies as they move under the influence of gravity. Those were the days when I was crazy about Graphics and the MGGS orbits brought to life all of the high school Physics I had learnt. Orbits, oscillations, parabolic and hyperbolic trajectories suddenly came to life on my computer screen. I didn't stop there and even implemented a 3D version in OpenGL. This was on a Pentium-4 system with 512 MB RAM and I remember the machine groaning when the OpenGL simulation ran for even 10 bodies.

And then I got a job, and in case you are wondering, yes, I can't find either program anywhere. Sad, because I am sure I implemented far more Physics in them than in these ones. The only bit extra I have built into the two implementations below is collision detection. To keep things simple, all bodies are assumed to be point particles and collisions to be perfectly elastic.

### And Now

One time step is assumed to be unit time, t = 1, hence there is no time variable in the implementations. The TCL example uses a time step of 10^-12, which frankly is too small for the C Standard Library and needs specialized techniques or libraries to be handled. Input data is read from a file and hence you can plot as many bodies as you want, keep in mind that this is computationally intensive, so intensive that on an i5, 8 GB RAM system, I didn't have to use any delays in the Graphics implementation for animation. Ideally such tasks are well suited for parallel processing. That's why Supercomputers are required to study weather, galaxies, black holes and molecules. All systems, big or small, can be reduced to N-body problems or definitely play a big role in their study.

### File Format

For both implementations, the file format is as follows :

<Gravitational Constant> <Number of bodies(N)> <Time step>
<Mass of M1>
<Position of M1 in x,y,z co-ordinates>
<Initial velocity of M1 in x,,y,z components>
...
<And so on for N bodies>


### The Implementations

#### Textual or Console

 #include<stdlib.h>#include<stdio.h>#include<math.h> typedef struct{	double x,y,z;}vector; int bodies,timeSteps;double *masses,GravConstant;vector *positions,*velocities,*accelerations; vector addVectors(vector a,vector b){	vector c = {a.x+b.x,a.y+b.y,a.z+b.z}; 	return c;} vector scaleVector(double b,vector a){	vector c = {b*a.x,b*a.y,b*a.z}; 	return c;} vector subtractVectors(vector a,vector b){	vector c = {a.x-b.x,a.y-b.y,a.z-b.z}; 	return c;} double mod(vector a){	return sqrt(a.x*a.x + a.y*a.y + a.z*a.z);} void initiateSystem(char* fileName){	int i;	FILE* fp = fopen(fileName,"r"); 	fscanf(fp,"%lf%d%d",&GravConstant,&bodies,&timeSteps); 	masses = (double*)malloc(bodies*sizeof(double));	positions = (vector*)malloc(bodies*sizeof(vector));	velocities = (vector*)malloc(bodies*sizeof(vector));	accelerations = (vector*)malloc(bodies*sizeof(vector)); 	for(i=0;i<bodies;i++){		fscanf(fp,"%lf",&masses[i]);		fscanf(fp,"%lf%lf%lf",&positions[i].x,&positions[i].y,&positions[i].z);		fscanf(fp,"%lf%lf%lf",&velocities[i].x,&velocities[i].y,&velocities[i].z);	} 	fclose(fp);} void resolveCollisions(){	int i,j; 	for(i=0;i<bodies-1;i++)		for(j=i+1;j<bodies;j++){			if(positions[i].x==positions[j].x && positions[i].y==positions[j].y && positions[i].z==positions[j].z){				vector temp = velocities[i];				velocities[i] = velocities[j];				velocities[j] = temp;			}		}} void computeAccelerations(){	int i,j; 	for(i=0;i<bodies;i++){		accelerations[i].x = 0;		accelerations[i].y = 0;		accelerations[i].z = 0;		for(j=0;j<bodies;j++){			if(i!=j){				accelerations[i] = addVectors(accelerations[i],scaleVector(GravConstant*masses[j]/pow(mod(subtractVectors(positions[i],positions[j])),3),subtractVectors(positions[j],positions[i])));			}		}	}} void computeVelocities(){	int i; 	for(i=0;i<bodies;i++)		velocities[i] = addVectors(velocities[i],accelerations[i]);} void computePositions(){	int i; 	for(i=0;i<bodies;i++)		positions[i] = addVectors(positions[i],addVectors(velocities[i],scaleVector(0.5,accelerations[i])));} void simulate(){	computeAccelerations();	computePositions();	computeVelocities();	resolveCollisions();} int main(int argC,char* argV[]){	int i,j; 	if(argC!=2)		printf("Usage : %s <file name containing system configuration data>",argV[0]);	else{		initiateSystem(argV[1]);		printf("Body   :     x              y               z           |           vx              vy              vz   ");		for(i=0;i<timeSteps;i++){			printf("\nCycle %d\n",i+1);			simulate();			for(j=0;j<bodies;j++)				printf("Body %d : %lf\t%f\t%lf\t|\t%lf\t%lf\t%lf\n",j+1,positions[j].x,positions[j].y,positions[j].z,velocities[j].x,velocities[j].y,velocities[j].z);		}	}	return 0;}

Input file: same data as the TCL example.

0.01 3 20
1
0 0 0
0.01 0 0
0.1
1 1 0
0 0 0.02
0.001
0 1 1
0.01 -0.01 -0.01


Invocation and output (x,y,z denote position components in 3-space, vx,vy,vz denote velocity components) :

C:\rosettaCode>mggs.exe mggsData.txt
Body   :     x              y               z           |           vx              vy              vz
Cycle 1
Body 1 : 0.010177	0.000179	0.000002	|	0.010354	0.000357	0.000004
Body 2 : 0.998230	0.998232	0.020002	|	-0.003539	-0.003536	0.020004
Body 3 : 0.010177	0.988232	0.988055	|	0.010354	-0.013536	-0.013889

Cycle 2
Body 1 : 0.020709	0.000718	0.000011	|	0.010710	0.000721	0.000014
Body 2 : 0.992907	0.992896	0.039971	|	-0.007109	-0.007138	0.019935
Body 3 : 0.020717	0.972888	0.972173	|	0.010727	-0.017153	-0.017876

Cycle 3
Body 1 : 0.031600	0.001625	0.000034	|	0.011072	0.001094	0.000033
Body 2 : 0.983985	0.983910	0.059834	|	-0.010735	-0.010835	0.019790
Body 3 : 0.031643	0.953868	0.952235	|	0.011125	-0.020886	-0.021999

Cycle 4
Body 1 : 0.042858	0.002913	0.000081	|	0.011443	0.001481	0.000060
Body 2 : 0.971393	0.971163	0.079509	|	-0.014448	-0.014659	0.019561
Body 3 : 0.042981	0.931039	0.928087	|	0.011552	-0.024772	-0.026299

Cycle 5
Body 1 : 0.054492	0.004595	0.000160	|	0.011826	0.001884	0.000097
Body 2 : 0.955030	0.954509	0.098909	|	-0.018278	-0.018649	0.019238
Body 3 : 0.054766	0.904225	0.899522	|	0.012018	-0.028857	-0.030829

Cycle 6
Body 1 : 0.066517	0.006691	0.000281	|	0.012224	0.002308	0.000145
Body 2 : 0.934759	0.933760	0.117931	|	-0.022265	-0.022849	0.018806
Body 3 : 0.067040	0.873197	0.866280	|	0.012530	-0.033199	-0.035655

Cycle 7
Body 1 : 0.078950	0.009225	0.000456	|	0.012642	0.002759	0.000206
Body 2 : 0.910400	0.908677	0.136456	|	-0.026454	-0.027316	0.018244
Body 3 : 0.079856	0.837662	0.828023	|	0.013101	-0.037871	-0.040861

Cycle 8
Body 1 : 0.091815	0.012227	0.000702	|	0.013086	0.003245	0.000284
Body 2 : 0.881722	0.878958	0.154340	|	-0.030902	-0.032122	0.017523
Body 3 : 0.093281	0.797239	0.784313	|	0.013749	-0.042975	-0.046559

Cycle 9
Body 1 : 0.105140	0.015737	0.001035	|	0.013564	0.003775	0.000383
Body 2 : 0.848429	0.844216	0.171401	|	-0.035684	-0.037362	0.016600
Body 3 : 0.107405	0.751427	0.734579	|	0.014498	-0.048649	-0.052908

Cycle 10
Body 1 : 0.118964	0.019805	0.001481	|	0.014085	0.004362	0.000509
Body 2 : 0.810137	0.803953	0.187408	|	-0.040900	-0.043166	0.015414
Body 3 : 0.122346	0.699554	0.678056	|	0.015384	-0.055097	-0.060138

Cycle 11
Body 1 : 0.133337	0.024498	0.002071	|	0.014662	0.005025	0.000672
Body 2 : 0.766343	0.757509	0.202050	|	-0.046687	-0.049720	0.013868
Body 3 : 0.138268	0.640690	0.613685	|	0.016460	-0.062633	-0.068603

Cycle 12
Body 1 : 0.148327	0.029907	0.002851	|	0.015317	0.005792	0.000888
Body 2 : 0.716377	0.703998	0.214889	|	-0.053246	-0.057302	0.011810
Body 3 : 0.155406	0.573482	0.539941	|	0.017816	-0.071782	-0.078886

Cycle 13
Body 1 : 0.164025	0.036157	0.003887	|	0.016079	0.006709	0.001184
Body 2 : 0.659310	0.642172	0.225282	|	-0.060887	-0.066351	0.008976
Body 3 : 0.174112	0.495836	0.454475	|	0.019596	-0.083511	-0.092045

Cycle 14
Body 1 : 0.180564	0.043437	0.005286	|	0.017000	0.007852	0.001613
Body 2 : 0.593807	0.570186	0.232208	|	-0.070119	-0.077621	0.004875
Body 3 : 0.194929	0.404136	0.353320	|	0.022038	-0.099890	-0.110265

Cycle 15
Body 1 : 0.198150	0.052049	0.007234	|	0.018171	0.009372	0.002283
Body 2 : 0.517817	0.485100	0.233878	|	-0.081861	-0.092550	-0.001535
Body 3 : 0.218605	0.290860	0.228583	|	0.025314	-0.126661	-0.139210

Cycle 16
Body 1 : 0.217126	0.062542	0.010117	|	0.019781	0.011614	0.003484
Body 2 : 0.427899	0.381659	0.226654	|	-0.097974	-0.114332	-0.012913
Body 3 : 0.244268	0.131956	0.057562	|	0.026013	-0.191148	-0.202831

Cycle 17
Body 1 : 0.238346	0.076539	0.015221	|	0.022658	0.016380	0.006723
Body 2 : 0.317489	0.248502	0.200967	|	-0.122846	-0.151982	-0.038461
Body 3 : 0.075592	-0.559591	-0.487315	|	-0.363366	-1.191945	-0.886924

Cycle 18
Body 1 : 0.263123	0.097523	0.026918	|	0.026898	0.025587	0.016672
Body 2 : 0.173428	0.050424	0.112716	|	-0.165275	-0.244174	-0.138041
Body 3 : -0.286241	-1.745597	-1.369528	|	-0.360299	-1.180066	-0.877501

Cycle 19
Body 1 : 0.270854	0.113045	0.061923	|	-0.011436	0.005457	0.053339
Body 2 : 0.199821	-0.093105	-0.208666	|	0.218061	-0.042882	-0.504723
Body 3 : -0.646318	-2.924909	-2.246453	|	-0.359856	-1.178559	-0.876350

Cycle 20
Body 1 : 0.258572	0.116046	0.112038	|	-0.013129	0.000544	0.046890
Body 2 : 0.426346	-0.111425	-0.681150	|	0.234987	0.006241	-0.440245
Body 3 : -1.006089	-4.103186	-3.122591	|	-0.359686	-1.177995	-0.875924



#### Graphics

One of the golden rules of Physics : Draw a diagram. It doesn't have to be artistic, but a properly drawn diagram whose proportions are realistic simplifies the solution of almost every Physics problem. That was the case with this task as well, when I first ran the console version, I noticed that Body 1's position and velocity is varying wildly. I put it down to the chaotic nature of the system and got busy writing the Graphics part and it's when I found that nothing was being drawn on the screen, I took a good look at the implementation and discovered many glaring mistakes. Gravitation is an attractive force, it never repels ( we will talk about Dark Energy later). When I saw the particles moving away from each other, I knew the acceleration vector had the wrong sign.

So yes, always draw a diagram, always do a Graphics implementation of any simulation. Input file for this differs from the console one since the BGI screen is the positive or 1st quadrant of the XY-plane. Hence z components are 0, you can use non-zero values but you won't see any change ( or really weird ones ) on the screen.

0.01 3 20
1.0
100.0 100.0 0.0
0.0 0.0 0.0
2.0
600.0 100.0 0.0
0.0 0.0 0.0
3
500.0 500.0 0.0
0.0 0.0 0.0


The time step although present in the file, and read, is not used here, the loop goes on indefinitely. Requires the WinBGIm library.

 #include<graphics.h>#include<stdlib.h>#include<stdio.h>#include<math.h> typedef struct{	double x,y,z;}vector; int bodies,timeSteps;double *masses,GravConstant;vector *positions,*velocities,*accelerations; vector addVectors(vector a,vector b){	vector c = {a.x+b.x,a.y+b.y,a.z+b.z}; 	return c;} vector scaleVector(double b,vector a){	vector c = {b*a.x,b*a.y,b*a.z}; 	return c;} vector subtractVectors(vector a,vector b){	vector c = {a.x-b.x,a.y-b.y,a.z-b.z}; 	return c;} double mod(vector a){	return sqrt(a.x*a.x + a.y*a.y + a.z*a.z);} void initiateSystem(char* fileName){	int i;	FILE* fp = fopen(fileName,"r"); 	fscanf(fp,"%lf%d%d",&GravConstant,&bodies,&timeSteps); 	masses = (double*)malloc(bodies*sizeof(double));	positions = (vector*)malloc(bodies*sizeof(vector));	velocities = (vector*)malloc(bodies*sizeof(vector));	accelerations = (vector*)malloc(bodies*sizeof(vector)); 	for(i=0;i<bodies;i++){		fscanf(fp,"%lf",&masses[i]);		fscanf(fp,"%lf%lf%lf",&positions[i].x,&positions[i].y,&positions[i].z);		fscanf(fp,"%lf%lf%lf",&velocities[i].x,&velocities[i].y,&velocities[i].z);	} 	fclose(fp);} void resolveCollisions(){	int i,j; 	for(i=0;i<bodies-1;i++)		for(j=i+1;j<bodies;j++){			if(positions[i].x==positions[j].x && positions[i].y==positions[j].y && positions[i].z==positions[j].z){				vector temp = velocities[i];				velocities[i] = velocities[j];				velocities[j] = temp;			}		}} void computeAccelerations(){	int i,j; 	for(i=0;i<bodies;i++){		accelerations[i].x = 0;		accelerations[i].y = 0;		accelerations[i].z = 0;		for(j=0;j<bodies;j++){			if(i!=j){				accelerations[i] = addVectors(accelerations[i],scaleVector(GravConstant*masses[j]/pow(mod(subtractVectors(positions[i],positions[j])),3),subtractVectors(positions[j],positions[i])));			}		}	}} void computeVelocities(){	int i; 	for(i=0;i<bodies;i++)		velocities[i] = addVectors(velocities[i],accelerations[i]);} void computePositions(){	int i; 	for(i=0;i<bodies;i++)		positions[i] = addVectors(positions[i],addVectors(velocities[i],scaleVector(0.5,accelerations[i])));} void simulate(){	computeAccelerations();	computePositions();	computeVelocities();	resolveCollisions();} void plotOrbits(){	int i; 	for(i=0;i<bodies;i++)		putpixel(positions[i].x,positions[i].y,i%15 + 1);} int main(int argC,char* argV[]){	int i,j; 	if(argC!=4)		printf("Usage : %s <file name containing system configuration data and width and height of window>",argV[0]);	else{		initiateSystem(argV[1]);		initwindow(atoi(argV[2]),atoi(argV[3]),"N Body Simulation");		while(1){			simulate();			plotOrbits();		}	}	return 0;}

### Further work

And as in every experiment report/journal paper/whatever, time for some philosophizing, a more rigorous implementation will require the bodies to be realistic, this means they will not be point particles and hence the dynamics will be much more complex. Tidal effects, collisions, orbits will require much more computation. When I say realistic, I don't mean rigid bodies because even rigid bodies are an idealization of reality. And all of this, even with classical mechanics, employing the Newtonian version of Gravitation will be a worthy task for any competent Physicist, let alone a Programmer. Introduce the General Theory of Relativity and everything changes. I can go on and on, but this is a coding, not a rambling site. I will end by saying that the Inverse Square Law, a name often used for Newtonian Gravitation, is just one class of forces studied in Dynamics. The force relation can take any form, the N-body problem can therefore be solved in (countably/uncountably ?) infinite ways.

## C#

Translation of: D
using System;using System.IO; namespace NBodyProblem {    class Vector3D {        public Vector3D(double x, double y, double z) {            X = x;            Y = y;            Z = z;        }         public double X { get; }        public double Y { get; }        public double Z { get; }         public double Mod() {            return Math.Sqrt(X * X + Y * Y + Z * Z);        }         public static Vector3D operator +(Vector3D lhs, Vector3D rhs) {            return new Vector3D(lhs.X + rhs.X, lhs.Y + rhs.Y, lhs.Z + rhs.Z);        }         public static Vector3D operator -(Vector3D lhs, Vector3D rhs) {            return new Vector3D(lhs.X - rhs.X, lhs.Y - rhs.Y, lhs.Z - rhs.Z);        }         public static Vector3D operator *(Vector3D lhs, double rhs) {            return new Vector3D(lhs.X * rhs, lhs.Y * rhs, lhs.Z * rhs);        }    }     class NBody {        private readonly double gc;        private readonly int bodies;        private readonly int timeSteps;        private readonly double[] masses;        private readonly Vector3D[] positions;        private readonly Vector3D[] velocities;        private readonly Vector3D[] accelerations;         public NBody(string fileName) {            string[] lines = File.ReadAllLines(fileName);             string[] gbt = lines[0].Split();            gc = double.Parse(gbt[0]);            bodies = int.Parse(gbt[1]);            timeSteps = int.Parse(gbt[2]);             masses = new double[bodies];            positions = new Vector3D[bodies];            velocities = new Vector3D[bodies];            accelerations = new Vector3D[bodies];            for (int i = 0; i < bodies; ++i) {                masses[i] = double.Parse(lines[i * 3 + 1]);                positions[i] = Decompose(lines[i * 3 + 2]);                velocities[i] = Decompose(lines[i * 3 + 3]);            }             Console.WriteLine("Contents of {0}", fileName);            foreach (string line in lines) {                Console.WriteLine(line);            }            Console.WriteLine();            Console.Write("Body   :      x          y          z    |");            Console.WriteLine("     vx         vy         vz");        }         public int GetTimeSteps() {            return timeSteps;        }         private Vector3D Decompose(string line) {            string[] xyz = line.Split();            double x = double.Parse(xyz[0]);            double y = double.Parse(xyz[1]);            double z = double.Parse(xyz[2]);            return new Vector3D(x, y, z);        }         private void ComputeAccelerations() {            for (int i = 0; i < bodies; ++i) {                accelerations[i] = new Vector3D(0, 0, 0);                for (int j = 0; j < bodies; ++j) {                    if (i != j) {                        double temp = gc * masses[j] / Math.Pow((positions[i] - positions[j]).Mod(), 3);                        accelerations[i] = accelerations[i] + (positions[j] - positions[i]) * temp;                    }                }            }        }         private void ComputeVelocities() {            for (int i = 0; i < bodies; ++i) {                velocities[i] = velocities[i] + accelerations[i];            }        }         private void ComputePositions() {            for (int i = 0; i < bodies; ++i) {                positions[i] = positions[i] + velocities[i] + accelerations[i] * 0.5;            }        }         private void ResolveCollisions() {            for (int i = 0; i < bodies; ++i) {                for (int j = i + 1; j < bodies; ++j) {                    if (positions[i].X == positions[j].X                     && positions[i].Y == positions[j].Y                     && positions[i].Z == positions[j].Z) {                        Vector3D temp = velocities[i];                        velocities[i] = velocities[j];                        velocities[j] = temp;                    }                }            }        }         public void Simulate() {            ComputeAccelerations();            ComputePositions();            ComputeVelocities();            ResolveCollisions();        }         public void PrintResults() {            for (int i = 0; i < bodies; ++i) {                Console.WriteLine(                    "Body {0} : {1,9:F6}  {2,9:F6}  {3,9:F6} | {4,9:F6}  {5,9:F6}  {6,9:F6}",                    i + 1,                    positions[i].X, positions[i].Y, positions[i].Z,                    velocities[i].X, velocities[i].Y, velocities[i].Z                );            }        }    }     class Program {        static void Main(string[] args) {            NBody nb = new NBody("nbody.txt");             for (int i = 0; i < nb.GetTimeSteps(); ++i) {                Console.WriteLine();                Console.WriteLine("Cycle {0}", i + 1);                nb.Simulate();                nb.PrintResults();            }        }    }}
Output:
Contents of nbody.txt
0.01 3 20
1
0 0 0
0.01 0 0
0.1
1 1 0
0 0 0.02
0.001
0 1 1
0.01 -0.01 -0.01

Body   :      x          y          z    |     vx         vy         vz

Cycle 1
Body 1 :  0.010177   0.000179   0.000002 |  0.010354   0.000357   0.000004
Body 2 :  0.998230   0.998232   0.020002 | -0.003539  -0.003536   0.020004
Body 3 :  0.010177   0.988232   0.988055 |  0.010354  -0.013536  -0.013889

Cycle 2
Body 1 :  0.020709   0.000718   0.000011 |  0.010710   0.000721   0.000014
Body 2 :  0.992907   0.992896   0.039971 | -0.007109  -0.007138   0.019935
Body 3 :  0.020717   0.972888   0.972173 |  0.010727  -0.017153  -0.017876

Cycle 3
Body 1 :  0.031600   0.001625   0.000034 |  0.011072   0.001094   0.000033
Body 2 :  0.983985   0.983910   0.059834 | -0.010735  -0.010835   0.019790
Body 3 :  0.031643   0.953868   0.952235 |  0.011125  -0.020886  -0.021999

Cycle 4
Body 1 :  0.042858   0.002913   0.000081 |  0.011443   0.001481   0.000060
Body 2 :  0.971393   0.971163   0.079509 | -0.014448  -0.014659   0.019561
Body 3 :  0.042981   0.931039   0.928087 |  0.011552  -0.024772  -0.026299

Cycle 5
Body 1 :  0.054492   0.004595   0.000160 |  0.011826   0.001884   0.000097
Body 2 :  0.955030   0.954509   0.098909 | -0.018278  -0.018649   0.019238
Body 3 :  0.054766   0.904225   0.899522 |  0.012018  -0.028857  -0.030829

Cycle 6
Body 1 :  0.066517   0.006691   0.000281 |  0.012224   0.002308   0.000145
Body 2 :  0.934759   0.933760   0.117931 | -0.022265  -0.022849   0.018806
Body 3 :  0.067040   0.873197   0.866280 |  0.012530  -0.033199  -0.035655

Cycle 7
Body 1 :  0.078950   0.009225   0.000456 |  0.012642   0.002759   0.000206
Body 2 :  0.910400   0.908677   0.136456 | -0.026454  -0.027316   0.018244
Body 3 :  0.079856   0.837662   0.828023 |  0.013101  -0.037871  -0.040861

Cycle 8
Body 1 :  0.091815   0.012227   0.000702 |  0.013086   0.003245   0.000284
Body 2 :  0.881722   0.878958   0.154340 | -0.030902  -0.032122   0.017523
Body 3 :  0.093281   0.797239   0.784313 |  0.013749  -0.042975  -0.046559

Cycle 9
Body 1 :  0.105140   0.015737   0.001035 |  0.013564   0.003775   0.000383
Body 2 :  0.848429   0.844216   0.171401 | -0.035684  -0.037362   0.016600
Body 3 :  0.107405   0.751427   0.734579 |  0.014498  -0.048649  -0.052908

Cycle 10
Body 1 :  0.118964   0.019805   0.001481 |  0.014085   0.004362   0.000509
Body 2 :  0.810137   0.803953   0.187408 | -0.040900  -0.043166   0.015414
Body 3 :  0.122346   0.699554   0.678056 |  0.015384  -0.055097  -0.060138

Cycle 11
Body 1 :  0.133337   0.024498   0.002071 |  0.014662   0.005025   0.000672
Body 2 :  0.766343   0.757509   0.202050 | -0.046687  -0.049720   0.013868
Body 3 :  0.138268   0.640690   0.613685 |  0.016460  -0.062633  -0.068603

Cycle 12
Body 1 :  0.148327   0.029907   0.002851 |  0.015317   0.005792   0.000888
Body 2 :  0.716377   0.703998   0.214889 | -0.053246  -0.057302   0.011810
Body 3 :  0.155406   0.573482   0.539941 |  0.017816  -0.071782  -0.078886

Cycle 13
Body 1 :  0.164025   0.036157   0.003887 |  0.016079   0.006709   0.001184
Body 2 :  0.659310   0.642172   0.225282 | -0.060887  -0.066351   0.008976
Body 3 :  0.174112   0.495836   0.454475 |  0.019596  -0.083511  -0.092045

Cycle 14
Body 1 :  0.180564   0.043437   0.005286 |  0.017000   0.007852   0.001613
Body 2 :  0.593807   0.570186   0.232208 | -0.070119  -0.077621   0.004875
Body 3 :  0.194929   0.404136   0.353320 |  0.022038  -0.099890  -0.110265

Cycle 15
Body 1 :  0.198150   0.052049   0.007234 |  0.018171   0.009372   0.002283
Body 2 :  0.517817   0.485100   0.233878 | -0.081861  -0.092550  -0.001535
Body 3 :  0.218605   0.290860   0.228583 |  0.025314  -0.126661  -0.139210

Cycle 16
Body 1 :  0.217126   0.062542   0.010117 |  0.019781   0.011614   0.003484
Body 2 :  0.427899   0.381659   0.226654 | -0.097974  -0.114332  -0.012913
Body 3 :  0.244268   0.131956   0.057562 |  0.026013  -0.191148  -0.202831

Cycle 17
Body 1 :  0.238346   0.076539   0.015221 |  0.022658   0.016380   0.006723
Body 2 :  0.317489   0.248502   0.200967 | -0.122846  -0.151982  -0.038461
Body 3 :  0.075592  -0.559591  -0.487315 | -0.363366  -1.191945  -0.886924

Cycle 18
Body 1 :  0.263123   0.097523   0.026918 |  0.026898   0.025587   0.016672
Body 2 :  0.173428   0.050424   0.112716 | -0.165275  -0.244174  -0.138041
Body 3 : -0.286241  -1.745597  -1.369528 | -0.360299  -1.180066  -0.877501

Cycle 19
Body 1 :  0.270854   0.113045   0.061923 | -0.011436   0.005457   0.053339
Body 2 :  0.199821  -0.093105  -0.208666 |  0.218061  -0.042882  -0.504723
Body 3 : -0.646318  -2.924909  -2.246453 | -0.359856  -1.178559  -0.876350

Cycle 20
Body 1 :  0.258572   0.116046   0.112038 | -0.013129   0.000544   0.046890
Body 2 :  0.426346  -0.111425  -0.681150 |  0.234987   0.006241  -0.440245
Body 3 : -1.006089  -4.103186  -3.122591 | -0.359686  -1.177995  -0.875924

## D

Translation of: Kotlin

#### Configurations

##### Configuration generator

Generate a radially symmetrical configuration for a given number of bodies and an initial rotational velocity

require'trig'GEN =: 3 : 0   1 GEN y:   m =. 1   r =. 1   p =. r*(|:(2,y)$(cos,sin)(i.y)*2*pi%y),.0 v =. x*(|:(2,y)$(cos,sin)(pi%2)+(i.y)*2*pi%y),.0   (i.y),.(y#m),.p,.v)
##### Generate a few cases.

The 3-body static equilibrium case is proposed as the basis for validation, and for comparison of alternative methods.

#### 3D

Although the animation is limited to two dimensions, we can verify that the implementation really does support three, by changing the axis of rotation in the initialization function:

cs:{(0;cos 2*pi*x%y;sin 2*pi*x%y)}     / x-axissc:{(0;sin 2*pi*x%y;-cos 2*pi*x%y)}    / objects on screen move up and down cs:{(cos 2*pi*x%y;0;sin 2*pi*x%y)}     / y-axissc:{(sin 2*pi*x%y;0;-cos 2*pi*x%y)}    / objects on screen move left and right cs:{(cos 2*pi*x%y;sin 2*pi*x%y;0)}     / z-axissc:{(sin 2*pi*x%y;-cos 2*pi*x%y;0)}    / objects on screen go round and round

## Kotlin

Translation of: C
// version 1.2.0 import java.io.Fileimport kotlin.math.sqrtimport kotlin.math.pow class Vector3D(val x: Double, val y: Double, val z: Double) {     operator fun plus(v: Vector3D) = Vector3D(x + v.x, y + v.y, z + v.z)    operator fun minus(v: Vector3D) = Vector3D(x - v.x, y - v.y, z - v.z)    operator fun times(s: Double) = Vector3D(s * x, s * y, s * z)     val mod = sqrt(x * x + y * y + z * z)} val origin = Vector3D(0.0, 0.0, 0.0) class NBody(fileName: String) {     val gc: Double    val bodies: Int    val timeSteps: Int    val masses: DoubleArray    val positions: Array<Vector3D>    val velocities: Array<Vector3D>    val accelerations: Array<Vector3D>     init {        val f = File(fileName)        val lines = f.readLines()        val (g, b, t) = lines[0].split(' ')        gc = g.toDouble()        bodies = b.toInt()        timeSteps = t.toInt()        masses = DoubleArray(bodies)        positions = Array<Vector3D>(bodies) { origin }        velocities = Array<Vector3D>(bodies) { origin }        accelerations = Array<Vector3D>(bodies) { origin }        for (i in 0 until bodies) {            masses[i] = lines[i * 3 + 1].toDouble()            positions[i] = decompose(lines[i * 3 + 2])            velocities[i] = decompose(lines[i * 3 + 3])        }        println("Contents of $fileName") println(f.readText()) print("Body : x y z |") println(" vx vy vz") } private fun decompose(line: String): Vector3D { val (x, y, z) = line.split(' ').map { it.toDouble() } return Vector3D(x, y, z) } private fun resolveCollisions() { for (i in 0 until bodies) { for (j in i + 1 until bodies) { if (positions[i].x == positions[j].x && positions[i].y == positions[j].y && positions[i].z == positions[j].z) { val temp = velocities[i] velocities[i] = velocities[j] velocities[j] = temp } } } } private fun computeAccelerations() { for (i in 0 until bodies) { accelerations[i] = origin for (j in 0 until bodies) { if (i != j) { val temp = gc * masses[j] / (positions[i] - positions[j]).mod.pow(3) accelerations[i] += (positions[j] - positions[i]) * temp } } } } private fun computeVelocities() { for (i in 0 until bodies) velocities[i] += accelerations[i] } private fun computePositions() { for (i in 0 until bodies) { positions[i] += velocities[i] + accelerations[i] * 0.5 } } fun simulate() { computeAccelerations() computePositions() computeVelocities() resolveCollisions() } fun printResults() { val fmt = "Body %d : % 8.6f % 8.6f % 8.6f | % 8.6f % 8.6f % 8.6f" for (i in 0 until bodies) { println(fmt.format( i + 1, positions[i].x, positions[i].y, positions[i].z, velocities[i].x, velocities[i].y, velocities[i].z )) } }} fun main(args: Array<String>) { val fileName = "nbody.txt" val nb = NBody(fileName) for (i in 1..nb.timeSteps) { println("\nCycle$i")        nb.simulate()        nb.printResults()    }}
Output:
Contents of nbody.txt
0.01 3 20
1
0 0 0
0.01 0 0
0.1
1 1 0
0 0 0.02
0.001
0 1 1
0.01 -0.01 -0.01

Body   :      x          y          z    |     vx         vy         vz

Cycle 1
Body 1 :  0.010177   0.000179   0.000002 |  0.010354   0.000357   0.000004
Body 2 :  0.998230   0.998232   0.020002 | -0.003539  -0.003536   0.020004
Body 3 :  0.010177   0.988232   0.988055 |  0.010354  -0.013536  -0.013889

Cycle 2
Body 1 :  0.020709   0.000718   0.000011 |  0.010710   0.000721   0.000014
Body 2 :  0.992907   0.992896   0.039971 | -0.007109  -0.007138   0.019935
Body 3 :  0.020717   0.972888   0.972173 |  0.010727  -0.017153  -0.017876

Cycle 3
Body 1 :  0.031600   0.001625   0.000034 |  0.011072   0.001094   0.000033
Body 2 :  0.983985   0.983910   0.059834 | -0.010735  -0.010835   0.019790
Body 3 :  0.031643   0.953868   0.952235 |  0.011125  -0.020886  -0.021999

Cycle 4
Body 1 :  0.042858   0.002913   0.000081 |  0.011443   0.001481   0.000060
Body 2 :  0.971393   0.971163   0.079509 | -0.014448  -0.014659   0.019561
Body 3 :  0.042981   0.931039   0.928087 |  0.011552  -0.024772  -0.026299

Cycle 5
Body 1 :  0.054492   0.004595   0.000160 |  0.011826   0.001884   0.000097
Body 2 :  0.955030   0.954509   0.098909 | -0.018278  -0.018649   0.019238
Body 3 :  0.054766   0.904225   0.899522 |  0.012018  -0.028857  -0.030829

Cycle 6
Body 1 :  0.066517   0.006691   0.000281 |  0.012224   0.002308   0.000145
Body 2 :  0.934759   0.933760   0.117931 | -0.022265  -0.022849   0.018806
Body 3 :  0.067040   0.873197   0.866280 |  0.012530  -0.033199  -0.035655

Cycle 7
Body 1 :  0.078950   0.009225   0.000456 |  0.012642   0.002759   0.000206
Body 2 :  0.910400   0.908677   0.136456 | -0.026454  -0.027316   0.018244
Body 3 :  0.079856   0.837662   0.828023 |  0.013101  -0.037871  -0.040861

Cycle 8
Body 1 :  0.091815   0.012227   0.000702 |  0.013086   0.003245   0.000284
Body 2 :  0.881722   0.878958   0.154340 | -0.030902  -0.032122   0.017523
Body 3 :  0.093281   0.797239   0.784313 |  0.013749  -0.042975  -0.046559

Cycle 9
Body 1 :  0.105140   0.015737   0.001035 |  0.013564   0.003775   0.000383
Body 2 :  0.848429   0.844216   0.171401 | -0.035684  -0.037362   0.016600
Body 3 :  0.107405   0.751427   0.734579 |  0.014498  -0.048649  -0.052908

Cycle 10
Body 1 :  0.118964   0.019805   0.001481 |  0.014085   0.004362   0.000509
Body 2 :  0.810137   0.803953   0.187408 | -0.040900  -0.043166   0.015414
Body 3 :  0.122346   0.699554   0.678056 |  0.015384  -0.055097  -0.060138

Cycle 11
Body 1 :  0.133337   0.024498   0.002071 |  0.014662   0.005025   0.000672
Body 2 :  0.766343   0.757509   0.202050 | -0.046687  -0.049720   0.013868
Body 3 :  0.138268   0.640690   0.613685 |  0.016460  -0.062633  -0.068603

Cycle 12
Body 1 :  0.148327   0.029907   0.002851 |  0.015317   0.005792   0.000888
Body 2 :  0.716377   0.703998   0.214889 | -0.053246  -0.057302   0.011810
Body 3 :  0.155406   0.573482   0.539941 |  0.017816  -0.071782  -0.078886

Cycle 13
Body 1 :  0.164025   0.036157   0.003887 |  0.016079   0.006709   0.001184
Body 2 :  0.659310   0.642172   0.225282 | -0.060887  -0.066351   0.008976
Body 3 :  0.174112   0.495836   0.454475 |  0.019596  -0.083511  -0.092045

Cycle 14
Body 1 :  0.180564   0.043437   0.005286 |  0.017000   0.007852   0.001613
Body 2 :  0.593807   0.570186   0.232208 | -0.070119  -0.077621   0.004875
Body 3 :  0.194929   0.404136   0.353320 |  0.022038  -0.099890  -0.110265

Cycle 15
Body 1 :  0.198150   0.052049   0.007234 |  0.018171   0.009372   0.002283
Body 2 :  0.517817   0.485100   0.233878 | -0.081861  -0.092550  -0.001535
Body 3 :  0.218605   0.290860   0.228583 |  0.025314  -0.126661  -0.139210

Cycle 16
Body 1 :  0.217126   0.062542   0.010117 |  0.019781   0.011614   0.003484
Body 2 :  0.427899   0.381659   0.226654 | -0.097974  -0.114332  -0.012913
Body 3 :  0.244268   0.131956   0.057562 |  0.026013  -0.191148  -0.202831

Cycle 17
Body 1 :  0.238346   0.076539   0.015221 |  0.022658   0.016380   0.006723
Body 2 :  0.317489   0.248502   0.200967 | -0.122846  -0.151982  -0.038461
Body 3 :  0.075592  -0.559591  -0.487315 | -0.363366  -1.191945  -0.886924

Cycle 18
Body 1 :  0.263123   0.097523   0.026918 |  0.026898   0.025587   0.016672
Body 2 :  0.173428   0.050424   0.112716 | -0.165275  -0.244174  -0.138041
Body 3 : -0.286241  -1.745597  -1.369528 | -0.360299  -1.180066  -0.877501

Cycle 19
Body 1 :  0.270854   0.113045   0.061923 | -0.011436   0.005457   0.053339
Body 2 :  0.199821  -0.093105  -0.208666 |  0.218061  -0.042882  -0.504723
Body 3 : -0.646318  -2.924909  -2.246453 | -0.359856  -1.178559  -0.876350

Cycle 20
Body 1 :  0.258572   0.116046   0.112038 | -0.013129   0.000544   0.046890
Body 2 :  0.426346  -0.111425  -0.681150 |  0.234987   0.006241  -0.440245
Body 3 : -1.006089  -4.103186  -3.122591 | -0.359686  -1.177995  -0.875924


## Octave

Translation of: Perl 6

We'll show only the positions in the output.

global G = 6.674e-11;   # gravitational constantglobal au = 150e9;      # astronomical unit global ma = 2e30;       # mass of the Sunglobal mb = 6e24;       # mass of Earthglobal mc = 7.34e22;    # mass of the Moon function ret = ABCdot(ABC, t)global G au ma mb mc ret = ABC;ret(1:9) = ABC(10:18); a = ABC(1:3);b = ABC(4:6);c = ABC(7:9); ab = norm(a - b);bc = norm(b - c);ca = norm(c - a); ret(10:12) = G*(mb/ab^3*(b - a) + mc/ca^3*(c - a));ret(13:15) = G*(ma/ab^3*(a - b) + mc/bc^3*(c - b));ret(16:18) = G*(ma/ca^3*(a - c) + mb/bc^3*(b - c)); endfunction; seconds_in_a_month = 60*60*24*27;seconds_in_a_year  = 60*60*24*365.25;t = (0:0.05:1);ABC = vec([0 0 0                                                           # Sun positionau 0 0                                                          # Earth position0.998*au 0 0                                                    # Moon position0 0 0                                                           # Sun speed0 2*pi*au/seconds_in_a_year 0                                   # Earth speed0 2*pi*(au/seconds_in_a_year + 0.002*au/seconds_in_a_month) 0    # Moon speed]');  disp(lsode('ABCdot',ABC, t)(1:20,1:9));
Output:
   0.0000e+00   0.0000e+00   0.0000e+00   1.5000e+11   0.0000e+00   0.0000e+00   1.4970e+11   0.0000e+00   0.0000e+00
2.2520e-11   1.8822e-19   0.0000e+00   1.5000e+11   1.4933e+03   0.0000e+00   1.4970e+11   1.5337e+03   0.0000e+00
9.0080e-11   1.3338e-18   0.0000e+00   1.5000e+11   2.9865e+03   0.0000e+00   1.4970e+11   3.0673e+03   0.0000e+00
2.0268e-10   3.1665e-18   0.0000e+00   1.5000e+11   4.4798e+03   0.0000e+00   1.4970e+11   4.6010e+03   0.0000e+00
3.6032e-10   6.7638e-18   0.0000e+00   1.5000e+11   5.9731e+03   0.0000e+00   1.4970e+11   6.1347e+03   0.0000e+00
5.6300e-10   1.2128e-17   0.0000e+00   1.5000e+11   7.4663e+03   0.0000e+00   1.4970e+11   7.6683e+03   0.0000e+00
8.1072e-10   1.9549e-17   0.0000e+00   1.5000e+11   8.9596e+03   0.0000e+00   1.4970e+11   9.2020e+03   0.0000e+00
1.1035e-09   3.0116e-17   0.0000e+00   1.5000e+11   1.0453e+04   0.0000e+00   1.4970e+11   1.0736e+04   0.0000e+00
1.4413e-09   4.3249e-17   0.0000e+00   1.5000e+11   1.1946e+04   0.0000e+00   1.4970e+11   1.2269e+04   0.0000e+00
1.8241e-09   1.0985e-16   0.0000e+00   1.5000e+11   1.3439e+04   0.0000e+00   1.4970e+11   1.3803e+04   0.0000e+00
2.2520e-09   1.8655e-16   0.0000e+00   1.5000e+11   1.4933e+04   0.0000e+00   1.4970e+11   1.5337e+04   0.0000e+00
2.7249e-09   2.7013e-16   0.0000e+00   1.5000e+11   1.6426e+04   0.0000e+00   1.4970e+11   1.6870e+04   0.0000e+00
3.2429e-09   3.6059e-16   0.0000e+00   1.5000e+11   1.7919e+04   0.0000e+00   1.4970e+11   1.8404e+04   0.0000e+00
3.8059e-09   4.5794e-16   0.0000e+00   1.5000e+11   1.9412e+04   0.0000e+00   1.4970e+11   1.9938e+04   0.0000e+00
4.4139e-09   5.6217e-16   0.0000e+00   1.5000e+11   2.0906e+04   0.0000e+00   1.4970e+11   2.1471e+04   0.0000e+00
5.0670e-09   6.7328e-16   0.0000e+00   1.5000e+11   2.2399e+04   0.0000e+00   1.4970e+11   2.3005e+04   0.0000e+00
5.7651e-09   7.9128e-16   0.0000e+00   1.5000e+11   2.3892e+04   0.0000e+00   1.4970e+11   2.4539e+04   0.0000e+00
6.5083e-09   9.1616e-16   0.0000e+00   1.5000e+11   2.5386e+04   0.0000e+00   1.4970e+11   2.6072e+04   0.0000e+00
7.2965e-09   1.0479e-15   0.0000e+00   1.5000e+11   2.6879e+04   0.0000e+00   1.4970e+11   2.7606e+04   0.0000e+00
8.1297e-09   1.1866e-15   0.0000e+00   1.5000e+11   2.8372e+04   0.0000e+00   1.4970e+11   2.9140e+04   0.0000e+00

## Pascal

This was written in Turbo Pascal, prompted by an article in an astronomy magazine about the orbit of an asteroid about a planet about a sun. It expects to receive command-line parameters for the run, and wants to use the screen graphics to animate the result. Alas, the turbo .bgi protocol no longer works on modern computers that have screens with many more dots. It contains options to select different methods for computation. An interesting point is that if the central mass is fixed in position (rather than wobbling about the centre of mass, as it does), then the Trojan orbit positions are not stable.

 {$N+ Crunch only the better sort of numbers, thanks.}{$M 50000,24000,24000] {Stack,minheap,maxheap.}Program Swirl; Uses Graph, Crt;{   Calculates and displays the orbit of an asteroid around a sun and a planet, as described in the article by C.C. Foster, in Sky & Telescope, September 1994, which included the source code of a programme written in basic, so some changes here.     The situation: the asteroid is at X = (AX,AY) with velocity V = (VAX,VAY) and suffers acceleration from the central sun and the steadily orbiting planet. If the acceleration at the start of a time step is A = (acX,acY)  then X = (A.dt/2 + V)dt + X   at the end of a time step, A presumed constant.   and V = A.dt + V             (so update X before altering V) But that is a first-order method. The whole point is that A is not constant. Accordingly, use this X as an estimate of the position at the end of a time step, and calculate a fresh A at the new position. Use the average of this and the first A to perform the actual step.    The next step in the escalation is the classic fourth-order Runge-Kutta method that alas involves a lot of repetitious code, or else trickery with loops. Whichever, there are four evaluations per step: at the start, then using that to probe ahead a half step, then again a half step, then a full step probe, then a weighted average to perform the step. While this gives fifth-order accuracy for each step, and thus fourth order accuracy at the end of a sequence of steps, this is rather painful. This is to say that halving the step size would reduce the upper bound of the error in a single step by a factor of 2**5, and that for the result after a sequence of steps by a factor of 2**4. However, this requires an upper bound on the value of the higher order derivatives of the function over the entirety of the area of interest, and for this problem, it does not exist, for the acceleration is unbounded above as a particle approaches a mass. It varies by 1/r**2 and as r approaches zero, there can be surprises.    Unfortunately, the next stages, say Milne's Predictor-Corrector (which has the advantage of maintaining  the recent history, just what procedure Apastron needs), followed by the extrapolation to zero stepsize methods of Bulirsch and Stoer require heavy administration.    This programme uses a few tricks to enable a step size of 1 so that no actual multiplication need be made, and likewise takes G as being one. It is the work of moments to change the step size with these methods, but the difficulty lies in when to change the step size and by how much. Remember that with a large step, singularities can be stepped over without noticing them, whereas a small step involves not just more work but also accumulating round-off. The predictor- corrector methods offer a good check on the step size (if you bother to compare the difference between the predicted and corrected values) but require a lot of administrative effort when changing it. However, this is not all to the bad as you might as well add some checks that the step size is not being changed too enthusiastically while you're at it. The problem is that the step size is being changed according to the behaviour of the results, but the behaviour of the results depends on the step size as well as the current region.    Rather than pursue this matter (and make the effort) I have retained the fixed size step (but perhaps an option sometime later?) to allow comparison with the original article in Sky & Telescope (and also to make things easy for procedure Apastron). For close comparison you will need to fix the central sun (i.e. Wobble is false) and use Euler's method. The article however uses units convenient to its programme: the planet is at radius 200, the asteroid at 100, both starting on the x-axis. The planet has a mass of 5 to the sun's mass of 70 (which is more like that of a companion sun than a planet), with a year of 920 steps and the asteroid has an initial velocity in the y-direction of 0.7. (See why you shouldn't use full stops as decimal points!)    These conditions result in the asteroid having successive periods of 466, 452, 477, 437, etc steps in the original article, and the invocation     WANDER 920 0.0714285  0.5 0  0.0  1.1  r f     meaning      920 steps to the planet's orbit,      a planetary mass of 5/70 of the sun      asteroid's radial distance being half that of the planet      at zero degrees (i.e. on the x-axis)      zero radial velocity      circular velocity of 1.1 (an estimate from trials) times that of a circular         orbit at the asteroid's initial radius      Runge-Kutta method      Sun fixed at the centre of mass     Is about as close as I can get, giving an asteroid period of 390 or so.    But the orbit is very unstable. With Euler's method, the picture is even worse with the asteroid spiralling outwards. I think that this is due to too large a step size and a low-order method, as test runs with the planetary mass set zero behave better, provided that there are many steps to an asteroid orbit and that it doesn't pass close to a mass. StepSize control would help, but it would mean more programming effort. Alternatively, I could have blundered, but I've stared at this so long that if there is a blunder, I can't see it.    I prefer to believe that the step size is the problem. In a different programme on a different machine, a circular orbit of 90 mins about the earth (just above the atmosphere) with full units for G, etc, and a step size of one second (thus no multiply) so that there were 5,400 steps to an orbit, Euler's method after one orbit had the radial distance 2% high and increasing steadily, whereas the 2'nd order Euler's method with a two second step size (so the same number of evaluations, as two per step) and also no multiplication (i.e. (a + a')/2 *dt) gave an error of one part on 6,000,000 (the precision of the 32-bit floating point numbers on an IBM1130) and oscillating.    On the other hand, a simple arrangement such as the asteroid in a Trojan orbit with respect to the planet works well enough. Here the asteroid is in the same orbit as the planet, but leading (or lagging) by sixty degrees, thus:     WANDER 920 0.001 1 +60  0 1     Test runs give an unstable orbit unless the planet is much smaller than 5/70 of the sun's mass, and also require that the sun not be fixed to the centre of mass. Having lots more steps per orbit doesn't change this. Thus, even though the sun's wobble is small, it is necessary and anyway, physically correct. In the case of a more massive planet, the obvious extension is to three equal masses in mutual orbit and invoking symmetry, a nice circular orbit would have them spaced at a hundred and twenty degrees.    Incidentally, no attempt is made to spot occasions when the asteroid's position is exactly that of some attracting mass, so zero divisions can cause dismay. Except that they won't happen as a result of a computed step unless you are very unlucky, especially if the floating-point hardware is available, with eighteen digits that must all coincide. On the other hand, you can easily specify an initial position that will provoke a division by zero so if that is what you have asked for, then that is what you will receive.}  {Perpetrated by R.N.McLean (whom God preserve), Victoria University, Feb. VMMI.}  Type Grist = {$IFOPT N+} extended {$ELSE} real {$ENDIF}; Type Chaff = {$IFOPT N+} single   {$ELSE} real {$ENDIF}; const esc = #27; var colour,lastcolour,xmax,ymax,txtheight: integer; var stretch: grist; Var Asitwas: Word; Type AnyString = string[80];  Function Min(i,j: integer): integer; BEGIN if i <= j then min:=i else min:=j; END; Function Max(i,j: integer): integer; BEGIN if i >= j then max:=i else max:=j; END; Procedure Beep(f,t: integer); BEGIN Sound(f); delay(t); NoSound; END; Procedure Croak(Gasp: string);        {A lethal word.}  BEGIN   TextMode(Asitwas);   WriteLn;   Writeln(Gasp);   HALT;  END;  Function MemGrab(Var p: pointer;s: word): boolean;  BEGIN   p:=nil;   if s > 0 then GetMem(p,s);   MemGrab:=p<>nil;  END; Procedure MemDrop(Var p: pointer;s: word);  BEGIN   if (s > 0) and (p <> nil) then FreeMem(p,s);   p:=nil;                              {Oh that all usage would first check!}  END;  FUNCTION Trim(S : anystring) : anystring; var L1,L2 : integer; BEGIN   L1 := 1;   WHILE (L1 <= LENGTH(S)) AND (S[L1] = ' ') DO INC(L1);   L2 := LENGTH(S);   WHILE (S[L2] = ' ') AND (L2 > L1) DO DEC(L2);   IF L2 >= L1 THEN Trim := COPY(S,L1,L2 - L1 + 1) ELSE Trim := ''; END; {Of Trim.}  FUNCTION Ifmt(Digits : longint) : anystring;  var  S : anystring;  BEGIN   STR(Digits,S);   Ifmt := Trim(S);  END; { Ifmt }  FUNCTION Ffmt(Digits: grist; Width,Decimals: integer): anystring;  var   S : anystring;   L : integer; { a finger }  BEGIN   if digits = 0 then begin Ffmt:='0'; exit; end; {Mumble.}   IF Abs(Digits) < 0.0000001 THEN STR(Digits,S)    ELSE STR(Digits:Width:Decimals,S);   s:=trim(s);   if copy(s,1,2) = '0.' then s:=copy(s,2,Length(s) - 1);   if copy(s,1,3) = '-0.' then s:=concat('-',copy(s,3,Length(s) - 2));   l:=Length(s);   if pos('.',s) > 0 then while (l > 0) and (s[l] = '0') do l:=l - 1;   if s[l] = '.' then l:=l - 1;         {No non-zero fractional part!}   S := COPY(S,1,L);                    {Drop trailing boredom.}   l:=Pos('.',S);                       {Avoid sequences such as 3.1415.}   if l > 0 then s[l]:=chr(249);        {By using decimal points, please.}   Ffmt := S;  END; { Ffmt }  Function Atan2(x,y: grist): grist;  var a: grist;  BEGIN   if x = 0 then    if y > 0 then atan2:=pi/2     else if y < 0 then atan2:=-pi/2      else atan2:=0    else if x > 0 then if y >= 0 then atan2:=arctan(y/x) else atan2:=2*pi + arctan(y/x)     else atan2:=pi - arctan(-y/x);  END;  Procedure PrepareTheCanvas;  Var mode: integer;  Var SavedResultCode: integer;  Var ax,ay: word;  Function BGIFound(Here: Anystring): boolean;   Var driver: integer;   BEGIN    Driver:=Detect;    InitGraph(driver,mode,Here);    SavedResultCode:=GraphResult;      {Mumble!}    BGIFound:=SavedResultCode = 0;     {Missing file gives -3.}   END; {Of BGIFound.}  Var aplace: anystring;  Var Palette: paletteType;  BEGIN   Mode:=0;   if not BGIFound('') then    if not BGIFound('C:\TP7\BGI') then     if not BGIFound('D:\TP4') then      begin       WriteLn;       WriteLn('Trouble with Borland''s Graphic Interface!');       WriteLn('Message: ',GraphErrorMsg(SavedResultCode));       WriteLn('I''ve tried "", "C:\TP4", and "D:\TP4"...');       Repeat        Write('Another place to look? (e.g. a:\bgiset)');        ReadLn(Aplace);         {As opposed to Read. weird.}        If Aplace = '' then         begin          writeln;          write('No .bgi, no piccy.');          halt;         end;       until BGIFound(Aplace);      end;   SetViewPort(0,0,GetMaxX,GetMaxY,True); {Clipping on, thanks. No apparent slowdown.}   SetGraphMode(Mode);   GetPalette(Palette);   LastColour := Palette.Size - 1; {Colours are 0:15, 16 thereof.}   if LastColour <= 0 then LastColour:=1;   xmax:=GetMaxX; ymax:=GetMaxY;   txtheight:=TextHeight('I');   GetAspectRatio(ax,ay); stretch:=ay/ax; {I hope this handles the screen's physical dimensions}  END;                                    {as well as the number of dots in the x and y directions.}  Procedure Splot(x,y:integer; bumf: anystring);  var l: integer;  BEGIN   l:=TextWidth(Bumf);   if x + l > xmax then x:=xmax - l;   SetFillStyle(EmptyFill,Black); Bar(x,y,x + l,y + txtheight - 1);   OutTextXY(x,y,bumf);  END;  Procedure FlashHint(msg: string);  var i,h,ph,pl,y0,rt: integer;  var l1,l2,nl: integer;  var board: pointer; BoardSize: word;  BEGIN   nl:=0; pl:=0; l2:=1;   repeat    l1:=l2;    while (l2 <= Length(msg)) and (msg[l2] <> '%') do l2:=l2 + 1;    nl:=nl + 1; pl:=max(pl,TextWidth(Copy(Msg,l1,l2 - l1)));    l2:=l2 + 1;   until l2 > Length(msg);   h:=TextHeight('Idiotic') + 1;   ph:=1 + nl*h;   y0:=GetMaxY div 6 + h;   BoardSize:=ImageSize(0,0,pl,ph);   If BoardSize <= 0 then begin Beep(6666,200); exit; end;   if KeyPressed then exit;        {Last chance to cut and run.}   if not memgrab(board,BoardSize) then    begin     OutTextXY(0,y0,'Memory shortage...');     Beep(8888,200);     exit; {Thus save on too many "if"s below.}    end;   GetImage(0,y0,pl,y0 + ph,board^);    {Doesn't check that Board is not nil!&^$#@!} SetFillStyle(EmptyFill,Black); Bar(0,y0,pl,y0 + ph); SetColor(White); i:=0; l2:=1; repeat l1:=l2; while (l2 <= Length(msg)) and (msg[l2] <> '%') do l2:=l2 + 1; OutTextXY(0,y0 + 1 + i*h,Copy(Msg,l1,l2 - l1)); i:=i + 1; l2:=l2 + 1; until l2 > Length(msg); rt:=28 + Length(msg); {Read time, more for longer messages.} i:=0; while not KeyPressed and (i < rt) do begin delay(100); i:=i + 1; end; PutImage(0,y0,board^,NormalPut); MemDrop(board,BoardSize); END; PROCEDURE LurkWith(Msg: string); Var t,w: integer; BEGIN w:=111; repeat for t:=0 to w do begin if KeyPressed then exit; delay(100); end; FlashHint(msg); w:=222; until KeyPressed; END; Function KeyFondle: char; {Equivalent to ReadKey, except...} Var ticker: integer; { after a delay, it gives a hint.} BEGIN {Screen and keyboard are connected by a computer...} Ticker:=666; {A delay counter.} While not KeyPressed do {If nothing has happened, } begin { twiddle my thumbs.} ticker:=ticker - 1; {My patience is being exhausted.} if ticker > 0 then Delay(60) {Another irretrievable loss.} else FlashHint('Press a key!'); {Hullo sailor!} end; {Eventually, a key is pressed.} KeyFondle:=ReadKey {Yum.} END; Var Wobble,Dotty: boolean; Var Style: integer; Var Year,PM,SM,AR,AA,ARV,ACV: grist; Var n: longint; Procedure WanderAlong; const suncolour = Yellow; const planetcolour = LightGreen; const asteroidcolour = White; const anglecolour = cyan; const periodcolour = LightRed; Const radiuscolour = LightBlue; Const KeplerColour = Magenta; Const EnergyColour = White; Const dt = 1; var W,P,T, SX,SY,RS, PX,PY,RP, AX,AY,VAX,VAY,acX,acY, DSX,DSY,DS2,DS3, DPX,DPY,DP2,DP3: grist; var cxp,cyp,sxp,syp,pxp,pyp,axp,ayp: integer; var bxp,txp,byp,typ: integer; var nturn: longint; var T1,T2: grist; var pAX,pAY,ppAX,ppAY,D2,pD2,ppD2: grist; Var nstep: longint; var stepwise,bored: boolean; Function AsteroidEnergy: grist; {The asteroid mass is considered divided out.} BEGIN {Thus, 1/2 mv**2 becomes 1/2 v**2.} AsteroidEnergy:=(VAX*VAX + VAY*VAY)/2 - SM/sqrt((AX - SX)*(AX - SX) + (AY - SY)*(AY - SY)) - PM/sqrt((AX - PX)*(AX - PX) + (AY - PY)*(AY - PY)); END; {Likewise, potential Energy = -G*M1*M2/r12 } Procedure SplotStep; BEGIN SetColor(LightBlue); Splot(0,0,'Steps: ' + ifmt(nstep)); Splot(4*cxp div 3,0,'R='+ffmt(sqrt(AX*AX + AY*AY),8,2)); SetColor(EnergyColour); Splot(4*cxp div 3,Txtheight,'E/m: '+ ffmt(AsteroidEnergy,8,2) + ' '); END; const types = 5; Const tummy = 666; var v: array[1..types,1..tummy] of Chaff; var vmax,vmin,vs: array[1..types] of Chaff; var show,varies: array[1..types] of boolean; const vcolour: array[1..types] of byte = (anglecolour,periodcolour,radiuscolour,KeplerColour,EnergyColour); var nv: integer; Procedure Splash; var i,j: integer; BEGIN SetFillStyle(EmptyFill,Black); Bar(bxp,byp,txp,typ); SetColor(Green); Line(bxp,byp,txp,byp); line(txp,byp,txp,typ); Line(txp,typ,bxp,typ); Line(bxp,typ,bxp,byp); for i:=1 to types do if show[i] and varies[i] then if dotty then for j:=1 to nv-1 do putpixel(bxp + round((v[i,j] - vmin[i])*vs[i]),byp+j,vcolour[i]) else begin setcolor(vcolour[i]); MoveTo(bxp + round((v[i,1] - vmin[i])*vs[i]),byp+1); for j:=2 to nv-1 do LineTo(bxp + round((v[i,j] - vmin[i])*vs[i]),byp+j); end; END; Procedure DealWith(ch: char); var redraw: boolean; BEGIN Case upcase(ch) of ESC: bored:=true; 'J': Dotty:=not Dotty; 'A': show[1]:=not show[1]; 'P': show[2]:=not show[2]; 'R': show[3]:=not show[3]; 'K': show[4]:=not show[4]; 'E': show[5]:=not show[5]; 'S':stepwise:=not stepwise; '?':FlashHint('Pressing a key may help.%APRK or E selects dots to plot%S for Stepwise%ESC to quit.'); end; if pos(upcase(ch),'JAPRKE') > 0 then splash; END; Procedure Swallow(v1,v2,v3,v4,v5: grist); var i,j,l: integer; var hic,flux: boolean; BEGIN flux:=nv = 3; if (nv >= typ - byp) or (nv >= tummy) then begin flux:=true; l:=nv div 3; for i:=1 to types do for j:=l+1 to nv do v[i,j - l]:=v[i,j]; nv:=nv - l; for i:=1 to types do begin vmin[i]:=v[i,1]; vmax[i]:=v[i,1]; end; for i:=1 to types do for j:=2 to nv do begin if v[i,j] > vmax[i] then vmax[i]:=v[i,j]; if v[i,j] < vmin[i] then vmin[i]:=v[i,j]; end; for i:=1 to types do begin varies[i]:=vmax[i] <> vmin[i]; if varies[i] then vs[i]:=(txp - bxp)/(vmax[i] - vmin[i]); end; end; nv:=nv + 1; v[1,nv]:=v1; v[2,nv]:=v2; v[3,nv]:=v3; v[4,nv]:=v4; v[5,nv]:=v5; if nv = 1 then for i:=1 to types do begin varies[i]:=false; vmax[i]:=v[i,1]; vmin[i]:=v[i,1]; end else for i:=1 to types do begin hic:=false; if v[i,nv] > vmax[i] then begin vmax[i]:=v[i,nv]; hic:=true; end; if v[i,nv] < vmin[i] then begin vmin[i]:=v[i,nv]; hic:=true; end; if hic then vs[i]:=(txp - bxp)/(vmax[i] - vmin[i]); if not varies[i] then varies[i]:=vmax[i] <> vmin[i]; flux:=flux or hic; end; if nv > 3 then begin if flux then splash; for i:=1 to types do if show[i] and varies[i] then if dotty then putpixel(bxp + round((v[i,nv] - vmin[i])*vs[i]),byp+nv,vcolour[i]) else begin setcolor(vcolour[i]); Line(bxp + round((v[i,nv-1] - vmin[i])*vs[i]),byp+nv-1,bxp + round((v[i,nv] - vmin[i])*vs[i]),byp+nv); end; end; END; {of Swallow.} Procedure Apastron; {Given: at the last three time steps, DS2, pD2, and ppD2, being (squares of) the distances from the centre of mass, and further, DS2 < pD2, pD2 > ppD2. Wanted: the time of greatest distance from the sun. So, fit a parabola to the three equally spaced values ppD2, pD2, DS2, differentiate and set to zero for the extremum, solving for x. y = ax**2 + bx + c, and why not have x = -1, 0, +1 Thus c = y0, y-1 is ppD2 so c = pD2 b = (y1 - y-1)/2 y0 is pD2 b = (DS2 - ppD2)/2 a = (y1 + y-1 - 2y0)/2 y1 is DS2. a = (DS2 + ppD2 - 2pDS2)/2 The extremum is when 2ax + b = 0 so x = -b/2a when y is an extremum, call this h. and x = 0 corresponds to the middle time, which was one back, so T-dt Having found the time of greatest distance, a small adjustment from T gives the location of the asteroid then (only a small change from time T) and thus the direction of the apastron. Since the acceleration is at right angles to the velocity it won't have much effect on the angle, and being minimal, not much effect on the distance. But, I'm not doing the arithmetic.} Var X,Y,a,r,h,k,e: grist; BEGIN {splot(0,20,'ppD2:'+ffmt(ppD2,15,2)+' ');} {splot(0,30,' pD2:'+ffmt(pD2,15,3)+'> ');} {splot(0,40,' D2:'+ffmt(D2,15,3)+' ');} h:=(ppD2 - DS2)/(DS2 + ppD2 - 2*pD2) /2; {splot(0,50,' h:'+ffmt(h,15,3+' '));} t2:=T + (h - 1){*dt}; x:=((AX + ppAX - 2*pAX)/2*h + (AX - ppAX)/2)*h + pAX; y:=((AY + ppAY - 2*pAY)/2*h + (AY - ppAY)/2)*h + pAY; r:=sqrt(x*x + y*y); SetColor(AngleColour); Line(cxp,cyp, cxp + round(x),cyp - round(y/stretch)); a:=atan2(x,y)*180/pi; nturn:=nturn + 1; h:=(t2 - t1)/P*Year; {Use the planet's 'year'.} k:=8*r*r*r/h/h; {Kepler's law: G(m + m')/4pi = d**3/t**2} e:=(VAX*VAX + VAY*VAY)/2 - SM/sqrt((x - SX)*(x - SX) + (y - SY)*(y - SY)) - PM/sqrt((x - PX)*(x - PX) + (y - PY)*(y - PY)); Splot(0,0 + 1*TxtHeight,'Orbits: ' + ifmt(nturn)); SetColor(AngleColour); Splot(0,ymax - 1*txtheight,'Apastron: ' + ffmt(a,8,3) + 'ø '); SetColor(PeriodColour); Splot(0,ymax - 2*txtheight,'Period: ' + ffmt(h,12,3)+' '); SetColor(RadiusColour); Splot(cxp*4 div 3,ymax - 1*txtheight,'Rmax: ' + ffmt(r,12,4)+' '); SetColor(KeplerColour); Splot(cxp*4 div 3,ymax - 2*txtheight,'Kepler: ' + ffmt(k,12,3)+' '); Swallow(a,h,r,k,e); t1:=t2; END; {of Apastron.} Const Methodname: array[1..3] of string = ('Euler','Euler2','RK4'); var nSX,nSY,nPX,nPY,nAX,nAY: grist; Procedure EulerStep; {First-order.} BEGIN DSX:=SX - AX; DSY:=SY - AY; DS2:=DSX*DSX + DSY*DSY; DS3:=sqrt(DS2)*DS2; DPX:=PX - AX; DPY:=PY - AY; DP2:=DPX*DPX + DPY*DPY; DP3:=sqrt(DP2)*DP2; acX:=SM*DSX/DS3 + PM*DPX/DP3; acY:=SM*DSY/DS3 + PM*DPY/DP3; END; Procedure SecondOrder; {Modified Euler, 2nd order predictor-corrector, etc..} BEGIN EulerStep; nAX:=AX + (0.5*acX{*dt} + VAX){*dt}; nAY:=AY + (0.5*acY{*dt} + VAY){*dt}; DSX:=nSX - nAX; DSY:=nSY - nAY; DS2:=DSX*DSX + DSY*DSY; DS3:=sqrt(DS2)*DS2; DPX:=nPX - nAX; DPY:=nPY - nAY; DP2:=DPX*DPX + DPY*DPY; DP3:=sqrt(DP2)*DP2; acX:=(acX + SM*DSX/DS3 + PM*DPX/DP3)/2; acY:=(acY + SM*DSY/DS3 + PM*DPY/DP3)/2; END; Procedure RungeKutta4; {The classic.} var h2,hsx,hsy,hpx,hpy,tax,tay,k1x,k1y,k2x,k2y,k3x,k3y: grist; BEGIN EulerStep; {The story at 0.} h2:=dt/2; {Half a time step.} hSX:=(SX + nSX)/2; hSY:=(SY + nSY)/2; {I'm not calling for a fresh sin and cos.} hPX:=(PX + nPX)/2; hPY:=(PY + nPY)/2; {Halfway will do.} tAX:=AX + (0.5*acX*h2 + VAX)*h2; tAY:=AY + (0.5*acY*h2 + VAY)*h2; DSX:=hSX - tAX; DSY:=hSY - tAY; DS2:=DSX*DSX + DSY*DSY; DS3:=sqrt(DS2)*DS2; DPX:=hPX - tAX; DPY:=hPY - tAY; DP2:=DPX*DPX + DPY*DPY; DP3:=sqrt(DP2)*DP2; k1x:=SM*DSX/DS3 + PM*DPX/DP3; k1y:=SM*DSY/DS3 + PM*DPY/DP3; tAX:=AX + (0.5*k1X*h2 + VAX)*h2; tAY:=AY + (0.5*k1Y*h2 + VAY)*h2; DSX:=hSX - tAX; DSY:=hSY - tAY; DS2:=DSX*DSX + DSY*DSY; DS3:=sqrt(DS2)*DS2; DPX:=hPX - tAX; DPY:=hPY - tAY; DP2:=DPX*DPX + DPY*DPY; DP3:=sqrt(DP2)*DP2; k2x:=SM*DSX/DS3 + PM*DPX/DP3; k2y:=SM*DSY/DS3 + PM*DPY/DP3; tAX:=AX + (0.5*k2X{*dt} + VAX){*dt}; tAY:=AY + (0.5*k2Y{*dt} + VAY){*dt}; DSX:=nSX - tAX; DSY:=nSY - tAY; DS2:=DSX*DSX + DSY*DSY; DS3:=sqrt(DS2)*DS2; DPX:=nPX - tAX; DPY:=nPY - tAY; DP2:=DPX*DPX + DPY*DPY; DP3:=sqrt(DP2)*DP2; k3x:=SM*DSX/DS3 + PM*DPX/DP3; k3y:=SM*DSY/DS3 + PM*DPY/DP3; acX:=(acX + 2*(k1x + k2x) + k3x)/6; acY:=(acY + 2*(k1y + k2y) + k3y)/6; END; Const BlobSize = 8; var r,vcirc: grist; var d,m,distant2: grist; var coswt,sinwt: grist; var rsp,rpp: integer; {Radius of Sun in Pixels, and of Planet.} var pspix,pppix,pswas,ppwas: pointer;{Points to Sun pix, Planet pix.} var szspix,szppix: word; {Sizes of the pix map.} var i1,i2: integer; Var ch: char; Var text: string; BEGIN bored:=false; Dotty:=true; Stepwise:=false; for i1:=1 to types do show[i1]:=true; bxp:=xmax; pyp:=round(ymax*stretch); if pyp < bxp then bxp:=pyp; {Oh for square dots!##$%$#!} cxp:=bxp div 2; cyp:=round(bxp/stretch/2); {The centre of mass on the screen.} bxp:=bxp + 1; txp:=xmax; {Base and top of a side box.} byp:=0; typ:=ymax; {Wherein will appear a plot.} nv:=0; {No values saved for the box.} R:=cxp - blobsize; {Available space on the screen, with safety.} RP:=R; {All for the planet's radius.} RS:=0; {If the sun is motionless.} M:=SM + PM; {Total mass of the system.} if pm = 0 then wobble:=false; {The asteroid is as if massless.} if wobble then {But if the sun wobbles, } begin { then a re-adjustment.} RS:=R*PM/M; {Both orbit about the centre of mass.} RP:=R*SM/M; {RS + RP = R.} t:=RS; if RP > t then t:=RP; {Re-scale to fill the screen.} RS:=RS*R/t; RP:=RP*R/t; end; D:=RS + RP; {Distance between planet and sun.} Distant2:=6*D*D; {A long way out.} M:=4*pi*pi*D*D*D/({G}year*year) /M; {Use the proper orbital period, with G = 1.} SM:=SM*M; PM:=PM*M; M:=SM + PM; {Rescale so as to allow dt = 1.} P:=year; W:=2*Pi/P; {2Pi radians in one orbit, needing time P} {dt:=P/year; {'Year' Timesteps complete one P.} rsp:=Round(BlobSize*Sqrt(SM/M)); rpp:=Round(BlobSize*Sqrt(PM/M)); if rsp < 1 then rsp:=1; if rpp < 1 then rpp:=1; {Prepare a solar disc.} sxp:=cxp - round(RS); syp:=cyp - 0; szspix:=ImageSize(sxp-rsp,syp-rsp,sxp+rsp,syp+rsp); if not(MemGrab(pspix,szspix) and MemGrab(pswas,szspix)) then croak('Can''t grab memory for the solar disc.'); GetImage(sxp-rsp,syp-rsp,sxp+rsp,syp+rsp,pswas^); SetColor(SunColour); Circle(sxp,syp,rsp); SetFillStyle(SolidFill,SunColour); FloodFill(sxp,syp,SunColour); GetImage(sxp-rsp,syp-rsp,sxp+rsp,syp+rsp,pspix^); PutImage(sxp-rsp,syp-rsp,pswas^,0); {Remove the sun during preparation.} {Prepare a planetary disc.} pxp:=cxp + round(RP); pyp:=cyp + 0; szppix:=ImageSize(pxp-rpp,pyp-rpp,pxp+rpp,pyp+rpp); If not(MemGrab(pppix,szppix) and MemGrab(ppwas,szppix)) then Croak('Can''t grab memory for the planetary disc.'); GetImage(pxp-rpp,pyp-rpp,pxp+rpp,pyp+rpp,ppwas^); SetColor(PlanetColour); Circle(pxp,pyp,rpp); SetFillStyle(SolidFill,PlanetColour); {FloodFill(pxp,pyp,PlanetColour);} GetImage(pxp-rpp,pyp-rpp,pxp+rpp,pyp+rpp,pppix^); PutImage(pxp-rpp,pyp-rpp,ppwas^,0); {Place some geometry.} SetColor(1); Line(cxp - round(rp),cyp,cxp + round(rp),cyp); Line(cxp,cyp + round(rp/stretch),cxp,cyp - round(rp/stretch)); Circle(cxp,cyp,Round(RP)); SetColor(5); Line(cxp-1,cyp,cxp+1,cyp); Line(cxp,cyp-1,cxp,cyp+1); GetImage(sxp-rsp,syp-rsp,sxp+rsp,syp+rsp,pswas^); {+Geometry, no sun.} GetImage(pxp-rpp,pyp-rpp,pxp+rpp,pyp+rpp,ppwas^); {+Geometry, no planet.} PutImage(sxp-rsp,syp-rsp,pspix^,0); SetColor(LightBlue); Splot(0,0+2*txtheight,'Method: ' + methodname[style]); if wobble then text:='Wobble' else text:='Fixed'; setcolor(SunColour); Splot(0,0+3*txtheight,text); SX:=0; SY:=0; {Place the sun at the centre of mass.} PX:=RP; PY:=0; {Place the planet on the x-axis.} if wobble then SX:=-RS; {Opposite sides.} nSX:=SX; nSY:=SY; {In case it doesn't wobble.} axp:=0; ayp:=0; {Previous asteroid place.} D:=RP*AR; {Distance out from the C of M.} AX:=D*Cos(AA); AY:=D*Sin(AA); VCirc:=D*2*Pi/Sqrt(4*pi*pi*D*D*D/M); {Velocity of circular orbit.} VAX:=VCirc*(-ACV*Sin(AA) + ARV*Cos(AA)); VAY:=VCirc*(+ACV*Cos(AA) + ARV*Sin(AA)); pAX:=AX; ppAX:=AX; pAY:=AY;ppAY:=AY;{Shouldn't be needed, but fear uninitialisation.} D2:=AX*AX + AY*AY; {Distance from the centre of mass.} pD2:=D2;ppD2:=D2; {We start at the Apastron.} T1:=0; nturn:=0; nstep:=0; {The story at time T, the start of a time step dt: The sun's mass SM at (SX,SY) and (nSX,nSY) at the end of the time step. The planet PM (PX,PY) (nPX,nPY) The asteroid (AX,AY) velocity (VAX,VAY) To compute its acceleration (acX,acY) and therefrom a new value for (AX,AY) and (VAX,VAY).} repeat nstep:=nstep + 1; {Well, it is about to be made.} t:=nstep{*dt}; {This time is after the step has been made.} CosWT:=Cos(W*T); SinWT:=Sin(W*T); {The compiler probably does a poor job.} nPX:=RP*Coswt;nPY:=RP*sinwt; {And at the end of the step.} if wobble then {In proper centre-of-mass arrangements?} begin {Yes. The sun is wobbling.} nSX:=-RS*CosWT; nSY:=-RS*SinWT; {And at the end of the time step.} end; {Enough fooling around. Now for the real task.} Case style of 1:EulerStep; 2:SecondOrder; 3:RungeKutta4; end; AX:=AX + (0.5*acX{*dt} + VAX){*dt}; AY:=AY + (0.5*acY{*dt} + VAY){*dt}; VAX:=VAX + acX{*dt}; VAY:=VAY + acY{*dt}; PX:=nPX; PY:=nPY; if wobble then begin SX:=nSX; SY:=nSY; end; {Thus attain the end of the time step.} PutPixel(axp,ayp,0); {Scrub the old asteroidal spot.} axp:=cxp + round(AX); ayp:=cyp - Round(AY/stretch); PutPixel(axp,ayp,AsteroidColour); {Place the new asteroidal spot.} PutImage(pxp-rpp,pyp-rpp,ppwas^,0); {Scrub the old planetary spot.} pxp:=cxp + round(PX); pyp:=cyp - round(PY/stretch); GetImage(pxp-rpp,pyp-rpp,pxp+rpp,pyp+rpp,ppwas^); PutImage(pxp-rpp,pyp-rpp,pppix^,0); {Place the new planetary spot.} if wobble then {The sun might be mobile, too.} begin {So see if it has moved noticeably.} i1:=cxp + round(SX); i2:=cyp - round(SY/stretch); if (i1 <> sxp) or (i2 <> syp) then {At least a whole dot in some direction?} begin {Yep.} PutImage(sxp-rsp,syp-rsp,pswas^,0); {Restore what was there.} sxp:=i1; syp:=i2; {And jump to the new position.} GetImage(sxp-rsp,syp-rsp,sxp+rsp,syp+rsp,pswas^); PutImage(sxp-rsp,syp-rsp,pspix^,0); {Splot.} end; {So much for a noticeable move.} end; {So much for sun movement.} D2:=AX*AX + AY*AY; {The new position.} if (D2 < pD2) and (pD2 > ppD2) then Apastron; ppAX:=pAX; ppAY:=pAY; pAX:=AX; pAY:=AY; ppD2:=pD2; pD2:=D2; if stepwise or (nstep mod 100 = 0) then splotstep; if D2 > Distant2 then {Too distant for close interactions?} if AsteroidEnergy > 0 then {Yeah, perhaps a getaway is in progress.} begin SplotStep; Splot(cxp,cyp div 3,'Escaping!'); Bored:=true; end; if stepwise or keypressed then dealwith(KeyFondle); until bored; repeat Text:=''; if nv > 3 then text:='Adjust dots, or%'; text:=text + 'ESC to quit.'; LurkWith(Text); ch:=KeyFondle; if ch <> esc then dealwith(ch); until ch = esc; MemDrop(pspix,szspix); MemDrop(pswas,szspix); MemDrop(pppix,szppix); MemDrop(ppwas,szppix); END; {Of WanderAlong.} Procedure Grunt(n: integer); Var i: integer; Var Unflushed: boolean; Procedure Z(Text: string); {Roll some text.} BEGIN {With screen pauses.} if unflushed then ClrEol; {Perhaps bumf lurks on this line.} WriteLn(Text); i:=i + 1; {Writes only to the end of text, not eol.} if i >= Hi(WindMax) then {Have we reached the bottom?} begin {Yes, the display would soon scroll up.} if unflushed then clreol; {A last remnant.} unflushed:=false; {Once scrolling starts, new lines are blank.} Write('(Press a key)'); {A hint, offering out-by-one possibilities.} if ReadKey = #0 then if ReadKey = esc then; {Ignore a key.} GoToXY(1,wherey); ClrEol; {Scrub the hint.} i:=0; {Restart the count.} end; {So much for a screen full.} END; {So much for that line.} BEGIN i:=n; Unflushed:=true; Z(' A massive planet is in a fixed circular orbit about its sun.'); Z('Between it and the sun is an asteroid that wanders as it can.'); Z('Although this problem is (just) within the scope of analytic methods,'); Z('the asteroid''s movement is computed by numerical integration, with one'); Z('step per ''day'' of the planet''s orbital period. As the simulation proceeds,'); Z('various attributes of its orbit are plotted in a box to the side of the screen'); Z('for each completed orbit. This is taken to be when it reaches a maximum'); Z('distance from the centre of mass, or in other words, starts to fall back.'); Z('No attempt is made to identify other orbital centres, such as the planet'); Z('or the Trojan points of its orbit or whatever other stabilities exist, nor'); Z('more baroque orbits such as about both the sun and the planet.'); Z(' The attributes plotted (radius, angle of apastron, etc) can be selected'); Z('or unselected by pressing the key for the first letter of the name, thus the'); Z('plot of the successive values for the Apastron may be suppressed by pressing'); Z('the A key, and so on. The plots appear as separate dots for each value, but if'); Z('you prefer Joined-up plots, pressing the J key will switch states.'); Z(' The S key switches between Stepwise and continual execution; pressing some'); Z('other key (e.g. Return) then advances the calculation one step each time.'); Z(' Everything depends on the initial state, and a distressingly long list of'); Z('parameters only begins to explore the opportunities, as follows:'); Z(''); Z(' YearLength: days in the planet''s year, e.g. 1000'); Z(' PlanetMass: as a fraction of the sun''s mass, e.g. 0.0125'); Z(' AsteroidR: as a fraction of the planet''s distance from the C of M'); Z(' AsteroidA: angle with the x-axis (in degrees)'); Z(' AsteroidRV: radial velocity as a factor of the circular orbit velocity'); Z(' AsteroidCV: circular orbit velocity factor'); Z(' Style: E = Euler, M = Euler 2''nd, R = Runge-Kutta.'); Z(' SunMotion: W = wobble about C of M, F = (improperly) fixed at C of M.'); Z(''); Z(' You need not supply all of these, but as they are not named, there can be'); Z('no gap in the list. Thus, if you supply the third parameter then you must'); Z('also supply the first and second, but need not supply those past the third.'); Z('If no parameters are offered, the default is equivalent to'); Z(''); Z(' WANDER 1000 0.15 0.5 0 0.0 -0.75 R W'); Z(''); Z(' This will use the classic fourth-order Runge-Kutta scheme for numerical'); Z('integration of a differential equation for the case where a massive planet'); Z('(the Earth''s mass is 1/300000 of the mass of the Sun) orbits in a fixed'); Z('circle having a thousand one-day steps, and a (massless) asteroid starts off'); Z('on the x-axis halfway between the planet and the centre of mass with an'); Z('initial velocity straight down of three quarters of the speed a circular orbit'); Z('of that radius would require in the absence of planetary perturbations.'); Z(''); Z(' Don''t forget the leading 0 in 0.15 that turbo pascal demands...'); Z(' ESC to quit.'); END; Procedure Reject(Gripe: anystring); BEGIN Writeln; Writeln('Unsavoury: ',Gripe); Writeln; Grunt(3); Halt; END; Const Usage: array[1..8] of string = ( 'Steps in a planetary year', 'Planetary mass ratio', 'Asteroid''s initial radial distance', 'Asteroid''s initial angular position', 'Asteroid''s initial outwards velocity', 'Asteroid''s initial orbital velocity', 'Numerical Integration Method', 'Fixed or Wobbling Sun position'); Procedure Scrog(i: integer; var v: grist); var hic: integer; BEGIN hic:=-666; if paramstr(i) <> '' then Val(Paramstr(i),v,hic); if hic <> 0 then Reject(paramstr(i) + ' for ' + usage[i]); END; Function ScrogT(i: integer; var n: integer; List: string): boolean; var c: char; var stupid: string[1]; BEGIN if paramstr(i) = '' then scrogt:=false else begin stupid:=copy(Paramstr(i),1,1); c:=char(hi(integer(stupid))); n:=Pos(upcase(c),List); if n <= 0 then Reject(paramstr(i) + ' not one of ' + list + ' for ' + usage[i]); scrogt:=true; end; END; {$F+} Function HeapFull(Size: word): integer; {$F-} Begin HeapFull:=1; End; {Sez "If full, return a null pointer" to GetMem.}{Damnit, Turbo pascal's pointer using procedures don't check for null pointers!} var i: Integer; BEGIN {The main programme at last.} HeapError:=@HeapFull; {Memory shortages?} AsItWas:=LastMode; Writeln(' A Wanderer in an unstable orbit.'); Style:=3; Wobble:=true; Year:=1000; SM:=1; PM:=0.15; AR:=0.5; AA:=0; ARV:=0; ACV:=-0.75; if paramcount > 0 then begin if paramstr(1) = '?' then begin Grunt(1); exit; end; Scrog(1,year); year:=round(year); {Whole steps only!} Scrog(2,PM); Scrog(3,AR); Scrog(4,AA); Scrog(5,ARV); Scrog(6,ACV); if ScrogT(7,Style,'EMR') then style:=max(1,min(style,3)); if ScrogT(8,i,'FW') then Wobble:=i = 2; end; AA:=Pi*AA/180; PrepareTheCanvas; WanderAlong; TextMode(AsItWas); END.  ### Relativistic orbits Another magazine article described a calculation for orbits twisted by the gravity around black holes. Only a two-body problem, and only just, because the orbiting body's mass is ignored.  {$N+ Crunch only the better sort of numbers, thanks.}Program Swirl; Uses Graph, Crt;{   Calculates and displays orbits strongly twisted by gravity, as described in the article by J. Gallmeier et al, in Sky & Telescope, October 1995, which included the source code of a programme written in basic, so some changes here.    An actual pulsar in orbit about another compact body has parameters EC = .61713 and SG = .000025, resulting in a precession of .0037 degrees/orbit.    More fierce relativistic effects involve close passages by a black hole. A dark disc appears corresponding to the event horizon diameter for a non-rotating black hole (the screen is normally 'black' so the black hole is not depicted as black) and a circle at twice its diameter is drawn; particles venturing within that bound will not repeat their orbit as they are either on the way down the drain, or else on a hyperbolic flypast.    This calculation is for orbits that continue indefinitely, and the authors have used the classic Runge-Kutta fourth-order method to numerically integrate the relativistic equations twisting the axes of the orbit.} {Perpetrated by R.N.McLean (whom God preserve), Victoria University, Nov. VMM.}  Type Grist = {$IFOPT N+} extended {$ELSE} real {$ENDIF}; const esc = #27; var colour,lastcolour,xmax,ymax,txtheight: integer; var stretch: grist; Type AnyString = string[80]; FUNCTION Trim(S : anystring) : anystring; var L1,L2 : integer; BEGIN L1 := 1; WHILE (L1 <= LENGTH(S)) AND (S[L1] = ' ') DO INC(L1); L2 := LENGTH(S); WHILE (S[L2] = ' ') AND (L2 > L1) DO DEC(L2); IF L2 >= L1 THEN Trim := COPY(S,L1,L2 - L1 + 1) ELSE Trim := ''; END; {Of Trim.} FUNCTION Ifmt(Digits : longint) : anystring; var S : anystring; BEGIN STR(Digits,S); Ifmt := Trim(S); END; { Ifmt } FUNCTION Ffmt(Digits: grist; Width,Decimals: integer): anystring; var S : anystring; L : integer; { a finger } BEGIN if digits = 0 then begin Ffmt:='0'; exit; end; {Mumble.} IF Abs(Digits) < 0.0000001 THEN STR(Digits,S) ELSE STR(Digits:Width:Decimals,S); s:=trim(s); if copy(s,1,2) = '0.' then s:=copy(s,2,Length(s) - 1); if copy(s,1,3) = '-0.' then s:=concat('-',copy(s,3,Length(s) - 2)); l:=Length(s); if pos('.',s) > 0 then while (l > 0) and (s[l] = '0') do l:=l - 1; if s[l] = '.' then l:=l - 1; {No non-zero fractional part!} S := COPY(S,1,L); {Drop trailing boredom.} l:=Pos('.',S); {Avoid sequences such as 3.1415.} if l > 0 then s[l]:=chr(249); {By using decimal points, please.} Ffmt := S; END; { Ffmt } Procedure Rogbiv; BEGIN; SetPalette( 0,Black); SetPalette( 1,DarkGray); {Looks navy blue to me.} SetPalette( 2,Blue); SetPalette( 3,LightBlue); SetPalette( 4,LightCyan); SetPalette( 5,Cyan); SetPalette( 6,LightGreen); SetPalette( 7,Green); SetPalette( 8,LightGray); {Looks light brown to me.} SetPalette( 9,Yellow); SetPalette(10,Brown); {Looks yellow to me.} SetPalette(11,Red); SetPalette(12,LightRed); SetPalette(13,Magenta); SetPalette(14,LightMagenta); SetPalette(15,White); END; {Of Rogbiv.} Procedure PrepareTheCanvas; Var mode: integer; Var SavedResultCode: integer; Var ax,ay: word; Function BGIFound(Here: Anystring): boolean; Var driver: integer; BEGIN Driver:=Detect; InitGraph(driver,mode,Here); SavedResultCode:=GraphResult; {Mumble!} BGIFound:=SavedResultCode = 0; {Missing file gives -3.} END; {Of BGIFound.} Var aplace: anystring; Var Palette: paletteType; BEGIN Mode:=0; if not BGIFound('') then if not BGIFound('C:\TP7\BGI') then if not BGIFound('D:\TP4') then begin WriteLn; WriteLn('Trouble with Borland''s Graphic Interface!'); WriteLn('Message: ',GraphErrorMsg(SavedResultCode)); WriteLn('I''ve tried "", "C:\TP4", and "D:\TP4"...'); Repeat Write('Another place to look? (e.g. a:\bgiset)'); ReadLn(Aplace); {As opposed to Read. weird.} If Aplace = '' then begin writeln; write('No .bgi, no piccy.'); halt; end; until BGIFound(Aplace); end; SetViewPort(0,0,GetMaxX,GetMaxY,True); {Clipping on, thanks. No apparent slowdown.} SetGraphMode(Mode); GetPalette(Palette); LastColour := Palette.Size - 1; {Colours are 0:15, 16 thereof.} if LastColour <= 0 then LastColour:=1; if LastColour = 15 then Rogbiv; xmax:=GetMaxX; ymax:=GetMaxY; txtheight:=TextHeight('I'); GetAspectRatio(ax,ay); stretch:=ay/ax; {I hope this handles the screen's physical dimensions} END; {as well as the number of dots in the x and y directions.} Procedure Splot(x,y:integer; bumf: anystring); var l: integer; BEGIN l:=TextWidth(Bumf); if x + l > xmax then x:=xmax - l; SetFillStyle(EmptyFill,Black); Bar(x,y,x + l,y + txtheight - 1); OutTextXY(x,y,bumf); END; Var nstep: longint; Procedure SplotStep; BEGIN Splot(xmax,0 + TxtHeight,'Stepcount: ' + ifmt(nstep)); END; Var EC,SG,SS,SA: grist; Var n: longint; Procedure PrecessRelativistically; var SX,SY, C0,C2,RH, Q,SN,QN,DN,HN,DP,SP, K1,K2,K3,K4,L1,L2,L3,L4: grist; var xc,yc,px,py: integer; Procedure Fullturn; BEGIN Line(xc,yc, round(PX),round(PY)); n:=n + 1; Splot(xmax,0,'Completed orbits ' + ifmt(n)); SplotStep; if n = 1 then {First time round?} begin {Yep. DP and DN straddle zero. Find S where D = 0.} SA:=(SP - DP*(SN - SP)/(DN - DP))*180/Pi - 360; Splot(xmax,ymax - txtheight,'Precession per orbit ' + ffmt(SA,15,5) + 'ø'); end; {y:x=0 given (x1,y1) and (x2,y2): y = y1 - x1*(y2 - y1)/(x2 - x1)} colour:=colour + 1; if colour > Lastcolour then colour:=3; setcolor(colour); END; Var ch: char; BEGIN Splot(0,0,'Eccentricity ' + ffmt(EC,9,6)); Splot(0,ymax - txtheight,'Relativistic Factor ' + ffmt(SG,9,6)); Splot(xmax,0 + 2*txtheight,'StepSize ' + ffmt(SS,8,3)); xc:=xmax div 2; yc:=ymax div 2; n:=0; nstep:=0; SS:=0.0009*(SS - 0.9)*(1 - 0.9*exp(10*EC - 9)); {Step size stuff...} SY:=yc/(1 + EC); SX:=SY*stretch; C0:=(1 - SG*(3 + EC*EC)/(6 + 2*EC))/(1 - EC*EC); RH:=SG*(1 - EC*EC)/(3 + EC); C2:=1.5*RH; SA:=0; Colour:=1; SetColor(Colour); Circle(xc,yc,round(SX*RH)); SetFillStyle(SolidFill,Colour); FloodFill(xc,yc,Colour); Circle(xc,yc,round(SX*RH*2)); Colour:=2; SetColor(Colour); Line(0,yc, xmax,yc); Line(xc,0, xc,ymax); Colour:=3; SetColor(Colour); SN:=0; QN:=1/(1 + EC); DN:=0; repeat DP:=DN; SP:=SN; HN:=SS*QN*QN; Q:=QN; K1:=(C0 - Q + C2*Q*Q)*HN; L1:=HN*DN; Q:=QN + L1/2; K2:=(C0 - Q + C2*Q*Q)*HN; L2:=HN*(DN + K1/2); Q:=QN + L2/2; K3:=(C0 - Q + C2*Q*Q)*HN; L3:=HN*(DN + K2/2); Q:=QN + L3; K4:=(C0 - Q + C2*Q*Q)*HN; L4:=HN*(DN + K3); QN:=QN + (L1 + (L2 + L3)*2 + L4)/6; DN:=DN + (K1 + (K2 + K3)*2 + K4)/6; SN:=SN + HN; px:=xc + Round(SX*cos(SN)/QN); py:=yc - Round(SY*sin(SN)/QN); PutPixel(px,py,Colour); nstep:=nstep + 1; if (DN > 0) and (DP < 0) then fullturn else if nstep mod 250 = 0 then splotstep; until keypressed; if readkey <> esc then repeat until keypressed; END; {Of PrecessRelativistically.} Procedure Grunt; BEGIN WriteLn('Draws orbits under relativistic conditions.'); WriteLn('Activate with two parameters:'); WriteLn(' Orbital Eccentricity (0 to 0ù9),'); WriteLn(' Relativistic factor (0 to 0ù999)'); WriteLn('An optional third parameter is the step size (1 to 10, e.g. 5).'); WriteLn('(If no parameters are supplied, an example run results)'); WriteLn('Eg. Swirl 0.5 0.5 5'); WriteLn('Don''t forget the leading zeroes that turbo pascal demands...'); WriteLn('ESC to quit.'); END; Procedure Reject(Gripe: anystring); BEGIN Writeln; Writeln('Unsavoury: ',Gripe); Writeln; Grunt; Halt; END; Var AsItWas: word; Var i: integer; Var Hic: array[1..3] of integer; BEGIN AsItWas:=LastMode; Writeln(' Precession of Elliptical Orbits.'); EC:=0.8; SG:=0.3826057; SS:=5; {These values result in 72 degrees/orbit...} if paramcount > 0 then begin if paramstr(1) = '?' then begin Grunt; exit; end; Val(ParamStr(1),EC,hic[1]); Val(ParamStr(2),SG,hic[2]); hic[3]:=0; if paramcount >= 3 then Val(paramstr(3),SS,Hic[3]); for i:=1 to 3 do if hic[i] > 0 then Reject('Parameter ' + ifmt(i) + ': ' + paramstr(i)); end; if (ec < 0) or (ec >= 1) then Reject('an eccentricity of ' + ffmt(ec,15,5)); if (sg < 0) or (sg >= 1) then Reject('a relativistic factor of ' + ffmt(sg,15,5)); if ss <= 0 then Reject('a step size of ' + ffmt(ss,15,5)); PrepareTheCanvas; PrecessRelativistically; TextMode(AsItWas); WriteLn('Orbital eccentricity:',EC:15:5); WriteLn('Relativistic factor :',SG:15:5); WriteLn('Precession per orbit:',SA:15:5); WriteLn('Orbits completed: ',n); WriteLn('Calculation steps:',nstep); END.  ## Perl 6 We'll try to simulate the Sun+Earth+Moon system, with plausible astronomical values. We use a 18-dimension vector ${\displaystyle ABC}$. The first nine dimensions are the positions of the three bodies. The other nine are the velocities. This allows us to write the dynamics as a first-temporal derivative equation, since ${\displaystyle {\frac {d\mathrm {position} }{dt}}=\mathrm {speed} }$ and thus ${\displaystyle {\frac {dABC_{1..9}}{dt}}=ABC_{10..18}}$ To keep things compact, we'll only display the first 20 lines of output. # Simple Vector implementationmulti infix:<+>(@a, @b) { @a Z+ @b }multi infix:<->(@a, @b) { @a Z- @b }multi infix:<*>($r, @a) { $r X* @a }multi infix:</>(@a,$r) { @a X/ $r }sub norm { sqrt [+] @_ X** 2 } # Runge-Kutta stuffsub runge-kutta(&yp) { return -> \t, \y, \δt { my$a = δt * yp( t, y );        my $b = δt * yp( t + δt/2, y +$a/2 );        my $c = δt * yp( t + δt/2, y +$b/2 );        my $d = δt * yp( t + δt, y +$c );        ($a + 2*($b + $c) +$d) / 6;    }} # gravitational constantconstant G = 6.674e-11;# astronomical unitconstant au = 150e9; # time constants in secondsconstant year = 365.25*24*60*60;constant month = 21*24*60*60; # masses in kgconstant $ma = 2e30; # Sunconstant$mb = 6e24;     # Earthconstant $mc = 7.34e22; # Moon my &dABC = runge-kutta my &f = sub ($t, @ABC ) {    my @a = @ABC[0..2];    my @b = @ABC[3..5];    my @c = @ABC[6..8];     my $ab = norm(@a - @b); my$ac = norm(@a - @c);    my $bc = norm(@b - @c); return [ flat @ABC[@(9..17)], map G * *,$mb/$ab**3 * (@b - @a) +$mc/$ac**3 * (@c - @a),$ma/$ab**3 * (@a - @b) +$mc/$bc**3 * (@c - @b),$ma/$ac**3 * (@a - @c) +$mb/$bc**3 * (@b - @c); ];} loop ( my ($t, @ABC) = 0,        0, 0, 0,                                 # Sun position        au, 0, 0,                                # Earth position        0.998*au, 0, 0,                          # Moon position        0, 0, 0,                                 # Sun speed        0, 2*pi*au/year, 0,                      # Earth speed        0, 2*pi*(au/year + 0.002*au/month), 0    # Moon speed    ;    $t < .2; ($t, @ABC) »+=« (.01, dABC($t, @ABC, .01))) { printf "t = %.02f : %s\n",$t, @ABC.fmt("%+.3e");}
Output:
t = 0.00 : +0.000e+00 +0.000e+00 +0.000e+00 +1.500e+11 +0.000e+00 +0.000e+00 +1.497e+11 +0.000e+00 +0.000e+00 +0.000e+00 +0.000e+00 +0.000e+00 +0.000e+00 +2.987e+04 +0.000e+00 +0.000e+00 +3.090e+04 +0.000e+00
t = 0.01 : +9.008e-13 +5.981e-22 +0.000e+00 +1.500e+11 +2.987e+02 +0.000e+00 +1.497e+11 +3.090e+02 +0.000e+00 +1.802e-10 +1.794e-19 +0.000e+00 -5.987e-05 +2.987e+04 +0.000e+00 -1.507e-05 +3.090e+04 +0.000e+00
t = 0.02 : +3.603e-12 +4.785e-21 +0.000e+00 +1.500e+11 +5.973e+02 +0.000e+00 +1.497e+11 +6.181e+02 +0.000e+00 +3.603e-10 +7.177e-19 +0.000e+00 -1.197e-04 +2.987e+04 +0.000e+00 -3.014e-05 +3.090e+04 +0.000e+00
t = 0.03 : +8.107e-12 +1.615e-20 +0.000e+00 +1.500e+11 +8.960e+02 +0.000e+00 +1.497e+11 +9.271e+02 +0.000e+00 +5.405e-10 +1.615e-18 +0.000e+00 -1.796e-04 +2.987e+04 +0.000e+00 -4.521e-05 +3.090e+04 +0.000e+00
t = 0.04 : +1.441e-11 +3.828e-20 +0.000e+00 +1.500e+11 +1.195e+03 +0.000e+00 +1.497e+11 +1.236e+03 +0.000e+00 +7.206e-10 +2.871e-18 +0.000e+00 -2.395e-04 +2.987e+04 +0.000e+00 -6.028e-05 +3.090e+04 +0.000e+00
t = 0.05 : +2.252e-11 +7.476e-20 +0.000e+00 +1.500e+11 +1.493e+03 +0.000e+00 +1.497e+11 +1.545e+03 +0.000e+00 +9.008e-10 +4.486e-18 +0.000e+00 -2.993e-04 +2.987e+04 +0.000e+00 -7.535e-05 +3.090e+04 +0.000e+00
t = 0.06 : +3.243e-11 +1.292e-19 +0.000e+00 +1.500e+11 +1.792e+03 +0.000e+00 +1.497e+11 +1.854e+03 +0.000e+00 +1.081e-09 +6.460e-18 +0.000e+00 -3.592e-04 +2.987e+04 +0.000e+00 -9.041e-05 +3.090e+04 +0.000e+00
t = 0.07 : +4.414e-11 +2.051e-19 +0.000e+00 +1.500e+11 +2.091e+03 +0.000e+00 +1.497e+11 +2.163e+03 +0.000e+00 +1.261e-09 +8.792e-18 +0.000e+00 -4.191e-04 +2.987e+04 +0.000e+00 -1.055e-04 +3.090e+04 +0.000e+00
t = 0.08 : +5.765e-11 +3.062e-19 +0.000e+00 +1.500e+11 +2.389e+03 +0.000e+00 +1.497e+11 +2.472e+03 +0.000e+00 +1.441e-09 +1.148e-17 +0.000e+00 -4.789e-04 +2.987e+04 +0.000e+00 -1.206e-04 +3.090e+04 +0.000e+00
t = 0.09 : +7.296e-11 +4.360e-19 +0.000e+00 +1.500e+11 +2.688e+03 +0.000e+00 +1.497e+11 +2.781e+03 +0.000e+00 +1.621e-09 +1.453e-17 +0.000e+00 -5.388e-04 +2.987e+04 +0.000e+00 -1.356e-04 +3.090e+04 +0.000e+00
t = 0.10 : +9.008e-11 +5.981e-19 +0.000e+00 +1.500e+11 +2.987e+03 +0.000e+00 +1.497e+11 +3.090e+03 +0.000e+00 +1.802e-09 +1.794e-17 +0.000e+00 -5.987e-04 +2.987e+04 +0.000e+00 -1.507e-04 +3.090e+04 +0.000e+00
t = 0.11 : +1.090e-10 +7.961e-19 +0.000e+00 +1.500e+11 +3.285e+03 +0.000e+00 +1.497e+11 +3.399e+03 +0.000e+00 +1.982e-09 +2.171e-17 +0.000e+00 -6.586e-04 +2.987e+04 +0.000e+00 -1.658e-04 +3.090e+04 +0.000e+00
t = 0.12 : +1.297e-10 +1.034e-18 +0.000e+00 +1.500e+11 +3.584e+03 +0.000e+00 +1.497e+11 +3.709e+03 +0.000e+00 +2.162e-09 +2.584e-17 +0.000e+00 -7.184e-04 +2.987e+04 +0.000e+00 -1.808e-04 +3.090e+04 +0.000e+00
t = 0.13 : +1.522e-10 +1.314e-18 +0.000e+00 +1.500e+11 +3.882e+03 +0.000e+00 +1.497e+11 +4.018e+03 +0.000e+00 +2.342e-09 +3.032e-17 +0.000e+00 -7.783e-04 +2.987e+04 +0.000e+00 -1.959e-04 +3.090e+04 +0.000e+00
t = 0.14 : +1.766e-10 +1.641e-18 +0.000e+00 +1.500e+11 +4.181e+03 +0.000e+00 +1.497e+11 +4.327e+03 +0.000e+00 +2.522e-09 +3.517e-17 +0.000e+00 -8.382e-04 +2.987e+04 +0.000e+00 -2.110e-04 +3.090e+04 +0.000e+00
t = 0.15 : +2.027e-10 +2.019e-18 +0.000e+00 +1.500e+11 +4.480e+03 +0.000e+00 +1.497e+11 +4.636e+03 +0.000e+00 +2.702e-09 +4.037e-17 +0.000e+00 -8.980e-04 +2.987e+04 +0.000e+00 -2.260e-04 +3.090e+04 +0.000e+00
t = 0.16 : +2.306e-10 +2.450e-18 +0.000e+00 +1.500e+11 +4.778e+03 +0.000e+00 +1.497e+11 +4.945e+03 +0.000e+00 +2.883e-09 +4.593e-17 +0.000e+00 -9.579e-04 +2.987e+04 +0.000e+00 -2.411e-04 +3.090e+04 +0.000e+00
t = 0.17 : +2.603e-10 +2.938e-18 +0.000e+00 +1.500e+11 +5.077e+03 +0.000e+00 +1.497e+11 +5.254e+03 +0.000e+00 +3.063e-09 +5.186e-17 +0.000e+00 -1.018e-03 +2.987e+04 +0.000e+00 -2.562e-04 +3.090e+04 +0.000e+00
t = 0.18 : +2.919e-10 +3.488e-18 +0.000e+00 +1.500e+11 +5.376e+03 +0.000e+00 +1.497e+11 +5.563e+03 +0.000e+00 +3.243e-09 +5.814e-17 +0.000e+00 -1.078e-03 +2.987e+04 +0.000e+00 -2.712e-04 +3.090e+04 +0.000e+00
t = 0.19 : +3.252e-10 +4.102e-18 +0.000e+00 +1.500e+11 +5.674e+03 +0.000e+00 +1.497e+11 +5.872e+03 +0.000e+00 +3.423e-09 +6.477e-17 +0.000e+00 -1.138e-03 +2.987e+04 +0.000e+00 -2.863e-04 +3.090e+04 +0.000e+00

## Python

import math class Vector:    def __init__(self, x, y, z):        self.x = x        self.y = y        self.z = z     def __add__(self, other):        return Vector(self.x + other.x, self.y + other.y, self.z + other.z)     def __sub__(self, other):        return Vector(self.x - other.x, self.y - other.y, self.z - other.z)     def __mul__(self, other):        return Vector(self.x * other, self.y * other, self.z * other)     def __div__(self, other):        return Vector(self.x / other, self.y / other, self.z / other)     def __eq__(self, other):        if isinstance(other, Vector):            return self.x == other.x and self.y == other.y and self.z == other.z        return False     def __ne__(self, other):        return not self.__eq__(other)     def __str__(self):        return '({x}, {y}, {z})'.format(x=self.x, y=self.y, z=self.z)     def abs(self):        return math.sqrt(self.x*self.x + self.y*self.y + self.z*self.z) origin = Vector(0, 0, 0) class NBody:    def __init__(self, fileName):        with open(fileName, "r") as fh:            lines = fh.readlines()            gbt = lines[0].split()            self.gc = float(gbt[0])            self.bodies = int(gbt[1])            self.timeSteps = int(gbt[2])            self.masses = [0.0 for i in range(self.bodies)]            self.positions = [origin for i in range(self.bodies)]            self.velocities = [origin for i in range(self.bodies)]            self.accelerations = [origin for i in range(self.bodies)]            for i in range(self.bodies):                self.masses[i] = float(lines[i*3 + 1])                self.positions[i] = self.__decompose(lines[i*3 + 2])                self.velocities[i] = self.__decompose(lines[i*3 + 3])             print "Contents of", fileName            for line in lines:                print line.rstrip()            print            print "Body   :      x          y          z    |",            print "     vx         vy         vz"     def __decompose(self, line):        xyz = line.split()        x = float(xyz[0])        y = float(xyz[1])        z = float(xyz[2])        return Vector(x, y, z)     def __computeAccelerations(self):        for i in xrange(self.bodies):            self.accelerations[i] = origin            for j in xrange(self.bodies):                if i != j:                    temp = self.gc * self.masses[j] / math.pow((self.positions[i] - self.positions[j]).abs(), 3)                    self.accelerations[i] += (self.positions[j] - self.positions[i]) * temp        return None     def __computePositions(self):        for i in xrange(self.bodies):            self.positions[i] += self.velocities[i] + self.accelerations[i] * 0.5        return None     def __computeVelocities(self):        for i in xrange(self.bodies):            self.velocities[i] += self.accelerations[i]        return None     def __resolveCollisions(self):        for i in xrange(self.bodies):            for j in xrange(self.bodies):                if self.positions[i] == self.positions[j]:                    (self.velocities[i], self.velocities[j]) = (self.velocities[j], self.velocities[i])        return None     def simulate(self):        self.__computeAccelerations()        self.__computePositions()        self.__computeVelocities()        self.__resolveCollisions()        return None     def printResults(self):        fmt = "Body %d : % 8.6f  % 8.6f  % 8.6f | % 8.6f  % 8.6f  % 8.6f"        for i in xrange(self.bodies):            print fmt % (i+1, self.positions[i].x, self.positions[i].y, self.positions[i].z, self.velocities[i].x, self.velocities[i].y, self.velocities[i].z)        return None nb = NBody("nbody.txt")for i in xrange(nb.timeSteps):    print "\nCycle %d" % (i + 1)    nb.simulate()    nb.printResults()
Output:
Contents of nbody.txt
0.01 3 20
1
0 0 0
0.01 0 0
0.1
1 1 0
0 0 0.02
0.001
0 1 1
0.01 -0.01 -0.01

Body   :      x          y          z    |      vx         vy         vz

Cycle 1
Body 1 :  0.010177   0.000179   0.000002 |  0.010354   0.000357   0.000004
Body 2 :  0.998230   0.998232   0.020002 | -0.003539  -0.003536   0.020004
Body 3 :  0.010177   0.988232   0.988055 |  0.010354  -0.013536  -0.013889

Cycle 2
Body 1 :  0.020709   0.000718   0.000011 |  0.010710   0.000721   0.000014
Body 2 :  0.992907   0.992896   0.039971 | -0.007109  -0.007138   0.019935
Body 3 :  0.020717   0.972888   0.972173 |  0.010727  -0.017153  -0.017876

Cycle 3
Body 1 :  0.031600   0.001625   0.000034 |  0.011072   0.001094   0.000033
Body 2 :  0.983985   0.983910   0.059834 | -0.010735  -0.010835   0.019790
Body 3 :  0.031643   0.953868   0.952235 |  0.011125  -0.020886  -0.021999

Cycle 4
Body 1 :  0.042858   0.002913   0.000081 |  0.011443   0.001481   0.000060
Body 2 :  0.971393   0.971163   0.079509 | -0.014448  -0.014659   0.019561
Body 3 :  0.042981   0.931039   0.928087 |  0.011552  -0.024772  -0.026299

Cycle 5
Body 1 :  0.054492   0.004595   0.000160 |  0.011826   0.001884   0.000097
Body 2 :  0.955030   0.954509   0.098909 | -0.018278  -0.018649   0.019238
Body 3 :  0.054766   0.904225   0.899522 |  0.012018  -0.028857  -0.030829

Cycle 6
Body 1 :  0.066517   0.006691   0.000281 |  0.012224   0.002308   0.000145
Body 2 :  0.934759   0.933760   0.117931 | -0.022265  -0.022849   0.018806
Body 3 :  0.067040   0.873197   0.866280 |  0.012530  -0.033199  -0.035655

Cycle 7
Body 1 :  0.078950   0.009225   0.000456 |  0.012642   0.002759   0.000206
Body 2 :  0.910400   0.908677   0.136456 | -0.026454  -0.027316   0.018244
Body 3 :  0.079856   0.837662   0.828023 |  0.013101  -0.037871  -0.040861

Cycle 8
Body 1 :  0.091815   0.012227   0.000702 |  0.013086   0.003245   0.000284
Body 2 :  0.881722   0.878958   0.154340 | -0.030902  -0.032122   0.017523
Body 3 :  0.093281   0.797239   0.784313 |  0.013749  -0.042975  -0.046559

Cycle 9
Body 1 :  0.105140   0.015737   0.001035 |  0.013564   0.003775   0.000383
Body 2 :  0.848429   0.844216   0.171401 | -0.035684  -0.037362   0.016600
Body 3 :  0.107405   0.751427   0.734579 |  0.014498  -0.048649  -0.052908

Cycle 10
Body 1 :  0.118964   0.019805   0.001481 |  0.014085   0.004362   0.000509
Body 2 :  0.810137   0.803953   0.187408 | -0.040900  -0.043166   0.015414
Body 3 :  0.122346   0.699554   0.678056 |  0.015384  -0.055097  -0.060138

Cycle 11
Body 1 :  0.133337   0.024498   0.002071 |  0.014662   0.005025   0.000672
Body 2 :  0.766343   0.757509   0.202050 | -0.046687  -0.049720   0.013868
Body 3 :  0.138268   0.640690   0.613685 |  0.016460  -0.062633  -0.068603

Cycle 12
Body 1 :  0.148327   0.029907   0.002851 |  0.015317   0.005792   0.000888
Body 2 :  0.716377   0.703998   0.214889 | -0.053246  -0.057302   0.011810
Body 3 :  0.155406   0.573482   0.539941 |  0.017816  -0.071782  -0.078886

Cycle 13
Body 1 :  0.164025   0.036157   0.003887 |  0.016079   0.006709   0.001184
Body 2 :  0.659310   0.642172   0.225282 | -0.060887  -0.066351   0.008976
Body 3 :  0.174112   0.495836   0.454475 |  0.019596  -0.083511  -0.092045

Cycle 14
Body 1 :  0.180564   0.043437   0.005286 |  0.017000   0.007852   0.001613
Body 2 :  0.593807   0.570186   0.232208 | -0.070119  -0.077621   0.004875
Body 3 :  0.194929   0.404136   0.353320 |  0.022038  -0.099890  -0.110265

Cycle 15
Body 1 :  0.198150   0.052049   0.007234 |  0.018171   0.009372   0.002283
Body 2 :  0.517817   0.485100   0.233878 | -0.081861  -0.092550  -0.001535
Body 3 :  0.218605   0.290860   0.228583 |  0.025314  -0.126661  -0.139210

Cycle 16
Body 1 :  0.217126   0.062542   0.010117 |  0.019781   0.011614   0.003484
Body 2 :  0.427899   0.381659   0.226654 | -0.097974  -0.114332  -0.012913
Body 3 :  0.244268   0.131956   0.057562 |  0.026013  -0.191148  -0.202831

Cycle 17
Body 1 :  0.238346   0.076539   0.015221 |  0.022658   0.016380   0.006723
Body 2 :  0.317489   0.248502   0.200967 | -0.122846  -0.151982  -0.038461
Body 3 :  0.075592  -0.559591  -0.487315 | -0.363366  -1.191945  -0.886924

Cycle 18
Body 1 :  0.263123   0.097523   0.026918 |  0.026898   0.025587   0.016672
Body 2 :  0.173428   0.050424   0.112716 | -0.165275  -0.244174  -0.138041
Body 3 : -0.286241  -1.745597  -1.369528 | -0.360299  -1.180066  -0.877501

Cycle 19
Body 1 :  0.270854   0.113045   0.061923 | -0.011436   0.005457   0.053339
Body 2 :  0.199821  -0.093105  -0.208666 |  0.218061  -0.042882  -0.504723
Body 3 : -0.646318  -2.924909  -2.246453 | -0.359856  -1.178559  -0.876350

Cycle 20
Body 1 :  0.258572   0.116046   0.112038 | -0.013129   0.000544   0.046890
Body 2 :  0.426346  -0.111425  -0.681150 |  0.234987   0.006241  -0.440245
Body 3 : -1.006089  -4.103186  -3.122591 | -0.359686  -1.177995  -0.875924

## Tcl

Works with: Tcl version 8.6
package require Tcl 8.6 set G 0.01set epsilon 1e-12 proc acceleration.gravity {positions masses} {    global G epsilon    set i -1    lmap position $positions mass$masses {	incr i	set dp2 [lmap p $position {expr 0.0}] set j -1 foreach pj$positions mj $masses { if {[incr j] ==$i} continue	    set dp [lmap p1 $position p2$pj {expr {$p1-$p2}}]	    set d3 [expr {		sqrt(		    [tcl::mathop::+ {*}[lmap p $dp {expr {$p ** 2}}]]		    + $epsilon ) ** 3 }] # Add epsilon here? set dp2 [lmap a$dp2 b $dp {expr {$a - $G*$mj*$b/$d3}}]	}	set dp2    }} # The rest of the system; simple numeric solution of differential equationsproc velocity {velocities accelerations} {    lmap v $velocities a$accelerations {	lmap vi $v ai$a {expr {$vi +$ai}}    }}proc position {positions velocities} {    lmap p $positions v$velocities {	lmap pi $p vi$v {expr {$pi +$vi}}    }}proc timestep {masses positions velocities} {    set accelerations [acceleration.gravity $positions$masses]    set velocities [velocity $velocities$accelerations]    set positions [position $positions$velocities]    list $positions$velocities} # Combine to make a simulation engineproc simulate {masses positions velocities {steps 10}} {    set p $positions set v$velocities    for {set i 0} {$i <$steps} {incr i} {	lassign [timestep $masses$p $v] p v puts [lmap pos$p {format (%.5f,%.5f,%.5f) {*}$pos}] }} Demonstrating for 20 steps: set m {1 .1 .001}set p {{0 0 0} {1 1 0} {0 1 1}}set v {{0.01 0 0} {0 0 0.02} {0.01 -0.01 -0.01}}simulate$m $p$v 20
Output:
(0.01035,0.00036,0.00000) (0.99646,0.99646,0.02000) (0.01035,0.98646,0.98611)
(0.02107,0.00108,0.00002) (0.98934,0.98931,0.03994) (0.02108,0.96930,0.96822)
(0.03214,0.00218,0.00005) (0.97856,0.97843,0.05973) (0.03221,0.94837,0.94617)
(0.04359,0.00367,0.00011) (0.96402,0.96368,0.07928) (0.04377,0.92350,0.91977)
(0.05544,0.00557,0.00021) (0.94559,0.94487,0.09851) (0.05581,0.89447,0.88875)
(0.06768,0.00790,0.00036) (0.92308,0.92176,0.11729) (0.06837,0.86100,0.85279)
(0.08036,0.01070,0.00057) (0.89626,0.89406,0.13548) (0.08152,0.82271,0.81146)
(0.09350,0.01400,0.00086) (0.86482,0.86136,0.15292) (0.09535,0.77910,0.76420)
(0.10714,0.01786,0.00126) (0.82839,0.82318,0.16939) (0.10996,0.72951,0.71025)
(0.12133,0.02234,0.00179) (0.78643,0.77885,0.18457) (0.12552,0.67304,0.64858)
(0.13614,0.02754,0.00251) (0.73827,0.72746,0.19806) (0.14224,0.60833,0.57768)
(0.15167,0.03357,0.00346) (0.68293,0.66775,0.20924) (0.16048,0.53332,0.49521)
(0.16805,0.04065,0.00476) (0.61899,0.59779,0.21710) (0.18075,0.44439,0.39722)
(0.18551,0.04908,0.00659) (0.54422,0.51449,0.21990) (0.20386,0.33411,0.27581)
(0.20445,0.05946,0.00934) (0.45470,0.41210,0.21398) (0.22971,0.17975,0.10892)
(0.22573,0.07337,0.01421) (0.34228,0.27744,0.18934) (0.19746,-0.27259,-0.30685)
(0.25165,0.09535,0.02601) (0.18355,0.06171,0.09512) (0.16823,-0.69092,-0.69109)
(0.21467,0.08626,0.10162) (0.65371,0.15666,-0.63735) (0.13969,-1.10220,-1.06883)
(0.17838,0.07727,0.17608) (1.11702,0.25051,-1.35831) (0.11149,-1.51049,-1.44390)
(0.14224,0.06831,0.25028) (1.57874,0.34406,-2.07666) (0.08347,-1.91722,-1.81757)