N-body problem: Difference between revisions
Added Algol 68
(Added Go) |
(Added Algol 68) |
||
| (26 intermediate revisions by 14 users not shown) | |||
Line 4:
The system is described by Newton's Law
: <math>\vec F = m_i \frac{d^2 \vec r_{i}}{dt^2}, i=1..N</math>
with
: <math>F_i = \sum^{}_{j\neq i} F_{ij}, i=1..N</math>
Line 13:
;Task
Write a simulation of three masses interacting under gravitation and show their evolution over a number of time-steps (at least 20).
=={{header|11l}}==
{{trans|Python}}
<syntaxhighlight lang="11l">-V origin = (0.0, 0.0, 0.0)
T NBody
Float gc
Int bodies
Int timeSteps
[Float] masses
[(Float, Float, Float)] positions, velocities, accelerations
F (fileName)
V lines = File(fileName).read().split("\n")
V gbt = lines[0].split(‘ ’)
.gc = Float(gbt[0])
.bodies = Int(gbt[1])
.timeSteps = Int(gbt[2])
.masses = [0.0] * .bodies
.positions = [:origin] * .bodies
.velocities = [:origin] * .bodies
.accelerations = [:origin] * .bodies
L(i) 0 .< .bodies
.masses[i] = Float(lines[i * 3 + 1])
.positions[i] = .__decompose(lines[i * 3 + 2])
.velocities[i] = .__decompose(lines[i * 3 + 3])
print(‘Contents of ’fileName)
L(line) lines
print(line)
print()
print(‘Body : x y z |’, end' ‘ ’)
print(‘ vx vy vz’)
F __decompose(line)
V xyz = line.split(‘ ’)
V x = Float(xyz[0])
V y = Float(xyz[1])
V z = Float(xyz[2])
R (x, y, z)
F __computeAccelerations()
L(i) 0 .< .bodies
.accelerations[i] = :origin
L(j) 0 .< .bodies
I i != j
V temp = .gc * .masses[j] / (length(.positions[i] - .positions[j]) ^ 3)
.accelerations[i] += (.positions[j] - .positions[i]) * temp
F __computePositions()
L(i) 0 .< .bodies
.positions[i] += .velocities[i] + .accelerations[i] * 0.5
F __computeVelocities()
L(i) 0 .< .bodies
.velocities[i] += .accelerations[i]
F __resolveCollisions()
L(i) 0 .< .bodies
L(j) 0 .< .bodies
I .positions[i] == .positions[j]
swap(&.velocities[i], &.velocities[j])
F simulate()
.__computeAccelerations()
.__computePositions()
.__computeVelocities()
.__resolveCollisions()
F printResults()
L(i) 0 .< .bodies
print(‘Body #. : #2.6 #2.6 #2.6 | #2.6 #2.6 #2.6’.format(i + 1, .positions[i].x, .positions[i].y, .positions[i].z, .velocities[i].x, .velocities[i].y, .velocities[i].z))
V nb = NBody(‘nbody.txt’)
L(i) 0 .< nb.timeSteps
print("\nCycle #.".format(i + 1))
nb.simulate()
nb.printResults()</syntaxhighlight>
{{out}}
<pre>
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 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
</pre>
=={{header|ALGOL 68}}==
{{Trans|Python}}
<syntaxhighlight lang="algol68">
BEGIN # N Body Problem - translated from the Python sample #
MODE VECTOR = STRUCT( REAL x, y, z );
OP + = ( VECTOR self, other )VECTOR:
( x OF self + x OF other, y OF self + y OF other, z OF self + z OF other);
OP - = ( VECTOR self, other )VECTOR:
( x OF self - x OF other, y OF self - y OF other, z OF self - z OF other );
OP * = ( VECTOR self, REAL other )VECTOR:
( x OF self * other, y OF self * other, z OF self * other );
OP / = ( VECTOR self, REAL other )VECTOR:
( x OF self / other, y OF self / other, z OF self / other );
OP = = ( VECTOR self, other )BOOL:
x OF self = x OF other AND y OF self = y OF other AND z OF self = z OF other;
OP /= = ( VECTOR self, other )BOOL: NOT ( self = other );
OP TOSTRING = ( INT v )STRING: whole( v, 0 );
OP TOSTRING = ( REAL v )STRING: fixed( v, -9, 6 );
OP TOSTRING = ( VECTOR self )STRING:
TOSTRING x OF self + " " + TOSTRING y OF self + " " + TOSTRING z OF self;
OP ABS = ( VECTOR self )REAL:
sqrt( ( x OF self * x OF self ) + ( y OF self * y OF self ) + ( z OF self * z OF self ) );
OP +:= = ( REF VECTOR self, VECTOR other )REF VECTOR: self := self + other;
VECTOR origin = ( 0, 0, 0 );
MODE NBODY = STRUCT( REAL gc
, INT bodies
, INT time steps
, REF[]REAL masses
, REF[]VECTOR positions, velocities, accelerations
);
PRIO INIT = 1;
OP INIT = ( REF NBODY self, STRING filename )VOID:
IF FILE data;
open( data, filename, stand back channel ) /= 0
THEN
print( ( "Unable to open ", filename, newline ) );
stop
ELSE
get( data, ( gc OF self, bodies OF self, time steps OF self, newline ) );
masses OF self := HEAP[ 1 : bodies OF self ]REAL;
positions OF self := HEAP[ 1 : bodies OF self ]VECTOR;
velocities OF self := HEAP[ 1 : bodies OF self ]VECTOR;
accelerations OF self := HEAP[ 1 : bodies OF self ]VECTOR;
FOR i TO bodies OF self DO
get( data, ( ( masses OF self )[ i ], newline ) );
get( data, ( x OF ( positions OF self )[ i ]
, y OF ( positions OF self )[ i ]
, z OF ( positions OF self )[ i ]
, newline
)
);
get( data, ( x OF ( velocities OF self )[ i ]
, y OF ( velocities OF self )[ i ]
, z OF ( velocities OF self )[ i ]
, newline
)
);
( accelerations OF self )[ i ] := origin
OD;
close( data );
print( ( "Contents of ", filename, newline ) );
print( ( TOSTRING gc OF self, " ", TOSTRING bodies OF self, " " ) );
print( ( TOSTRING time steps OF self, newline ) );
FOR i TO bodies OF self DO
print( ( TOSTRING ( masses OF self )[ i ], newline ) );
print( ( TOSTRING ( positions OF self )[ i ], newline ) );
print( ( TOSTRING ( velocities OF self )[ i ], newline ) )
OD;
print( ( newline ) );
print( ( "Body : x y z |" ) );
print( ( " vx vy vz", newline ) )
FI # INIT # ;
OP COMPUTEACCELERATIONS = ( REF NBODY self )VOID:
FOR i TO bodies OF self DO
( accelerations OF self )[ i ] := origin;
FOR j TO bodies OF self DO
IF i /= j THEN
REAL temp = ( gc OF self * ( masses OF self )[ j ] )
/ ( ABS ( ( positions OF self )[ i ] - ( positions OF self )[ j ] ) ^ 3 );
( accelerations OF self )[ i ]
+:= ( ( positions OF self )[ j ] - ( positions OF self )[ i ] ) * temp
FI
OD
OD # COMPUTEACCELERATIONS # ;
OP COMPUTEPOSITIONS = ( REF NBODY self )VOID:
FOR i TO bodies OF self DO
( positions OF self )[ i ] +:= ( velocities OF self )[ i ] + ( accelerations OF self )[ i ] * 0.5
OD # COMPUTEPOSITIONS # ;
OP COMPUTEVELOCITIES = ( REF NBODY self )VOID:
FOR i TO bodies OF self DO
( velocities OF self )[ i ] +:= ( accelerations OF self )[ i ]
OD # COMPUTEVELOCITES # ;
OP RESOLVECOLLISIONS = ( REF NBODY self )VOID:
FOR i TO bodies OF self DO
FOR j TO bodies OF self DO
IF ( positions OF self )[ i ] = ( positions OF self )[ j ] THEN
VECTOR vj = ( velocities OF self )[ j ];
( velocities OF self )[ j ] := ( velocities OF self )[ i ];
( velocities OF self )[ i ] := vj
FI
OD
OD # RESOLVECOLLISIONS # ;
OP SIMULATE = ( REF NBODY self )VOID:
BEGIN
COMPUTEACCELERATIONS self;
COMPUTEPOSITIONS self;
COMPUTEVELOCITIES self;
RESOLVECOLLISIONS self
END # SIMULATE # ;
NBODY nb;
nb INIT "nbody.txt";
FOR i TO time steps OF nb DO
print( ( newline, "Cycle ", TOSTRING i, newline ) );
SIMULATE nb;
FOR i TO bodies OF nb DO
print( ( "Body ", TOSTRING i
, " : ", TOSTRING ( positions OF nb )[ i ]
, " | ", TOSTRING ( velocities OF nb )[ i ]
, newline
)
)
OD
OD
END
</syntaxhighlight>
{{out}}
<pre>
Contents of nbody.txt
0.010000 3 20
1.000000
0.000000 0.000000 0.000000
0.010000 0.000000 0.000000
0.100000
1.000000 1.000000 0.000000
0.000000 0.000000 0.020000
0.001000
0.000000 1.000000 1.000000
0.010000 -0.010000 -0.010000
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 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
</pre>
=={{header|C}}==
Line 36 ⟶ 338:
===The Implementations===
====Textual or Console====
<syntaxhighlight lang="c">
#include<stdlib.h>
#include<stdio.h>
Line 160 ⟶ 460:
return 0;
}
</syntaxhighlight>
Input file: same data as the TCL example.
<pre>
Line 281 ⟶ 581:
====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
----
Line 298 ⟶ 598:
</pre>
The time step although present in the file, and read, is not used here, the loop goes on indefinitely. Requires the [http://www.cs.colorado.edu/~main/bgi/cs1300/ WinBGIm] library.
<syntaxhighlight lang="c">
#include<graphics.h>
#include<stdlib.h>
Line 428 ⟶ 726:
return 0;
}
</syntaxhighlight>
===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.
=={{header|C sharp|C#}}==
{{trans|D}}
<syntaxhighlight lang="csharp">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();
}
}
}
}</syntaxhighlight>
{{out}}
<pre>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</pre>
=={{header|C++}}==
{{trans|Java}}
<syntaxhighlight lang="cpp">#include <algorithm>
#include <iomanip>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
class Vector {
private:
double px, py, pz;
public:
Vector() : px(0.0), py(0.0), pz(0.0) {
// empty
}
Vector(double x, double y, double z) : px(x), py(y), pz(z) {
// empty
}
double mod() const {
return sqrt(px*px + py * py + pz * pz);
}
Vector operator+(const Vector& rhs) const {
return Vector(px + rhs.px, py + rhs.py, pz + rhs.pz);
}
Vector operator-(const Vector& rhs) const {
return Vector(px - rhs.px, py - rhs.py, pz - rhs.pz);
}
Vector operator*(double s) const {
return Vector(px*s, py*s, pz*s);
}
bool operator==(const Vector& rhs) const {
return px == rhs.px
&& py == rhs.py
&& pz == rhs.pz;
}
friend std::istream& operator>>(std::istream&, Vector&);
friend std::ostream& operator<<(std::ostream&, Vector&);
};
std::istream& operator>>(std::istream& in, Vector& v) {
return in >> v.px >> v.py >> v.pz;
}
std::ostream& operator<<(std::ostream& out, Vector& v) {
auto precision = out.precision();
auto width = out.width();
out << std::fixed << std::setw(width) << std::setprecision(precision) << v.px << " ";
out << std::fixed << std::setw(width) << std::setprecision(precision) << v.py << " ";
out << std::fixed << std::setw(width) << std::setprecision(precision) << v.pz;
return out;
}
const Vector ORIGIN{ 0.0, 0.0, 0.0 };
class NBody {
private:
double gc;
int bodies;
int timeSteps;
std::vector<double> masses;
std::vector<Vector> positions;
std::vector<Vector> velocities;
std::vector<Vector> accelerations;
void resolveCollisions() {
for (int i = 0; i < bodies; ++i) {
for (int j = i + 1; j < bodies; ++j) {
if (positions[i] == positions[j]) {
std::swap(velocities[i], velocities[j]);
}
}
}
}
void computeAccelerations() {
for (int i = 0; i < bodies; ++i) {
accelerations[i] = ORIGIN;
for (int j = 0; j < bodies; ++j) {
if (i != j) {
double temp = gc * masses[j] / pow((positions[i] - positions[j]).mod(), 3);
accelerations[i] = accelerations[i] + (positions[j] - positions[i]) * temp;
}
}
}
}
void computeVelocities() {
for (int i = 0; i < bodies; ++i) {
velocities[i] = velocities[i] + accelerations[i];
}
}
void computePositions() {
for (int i = 0; i < bodies; ++i) {
positions[i] = positions[i] + velocities[i] + accelerations[i] * 0.5;
}
}
public:
NBody(std::string& fileName) {
using namespace std;
ifstream ifs(fileName);
if (!ifs.is_open()) {
throw runtime_error("Could not open file.");
}
ifs >> gc >> bodies >> timeSteps;
masses.resize(bodies);
positions.resize(bodies);
fill(positions.begin(), positions.end(), ORIGIN);
velocities.resize(bodies);
fill(velocities.begin(), velocities.end(), ORIGIN);
accelerations.resize(bodies);
fill(accelerations.begin(), accelerations.end(), ORIGIN);
for (int i = 0; i < bodies; ++i) {
ifs >> masses[i] >> positions[i] >> velocities[i];
}
cout << "Contents of " << fileName << '\n';
cout << gc << ' ' << bodies << ' ' << timeSteps << '\n';
for (int i = 0; i < bodies; ++i) {
cout << masses[i] << '\n';
cout << positions[i] << '\n';
cout << velocities[i] << '\n';
}
cout << "\nBody : x y z | vx vy vz\n";
}
int getTimeSteps() {
return timeSteps;
}
void simulate() {
computeAccelerations();
computePositions();
computeVelocities();
resolveCollisions();
}
friend std::ostream& operator<<(std::ostream&, NBody&);
};
std::ostream& operator<<(std::ostream& out, NBody& nb) {
for (int i = 0; i < nb.bodies; ++i) {
out << "Body " << i + 1 << " : ";
out << std::setprecision(6) << std::setw(9) << nb.positions[i];
out << " | ";
out << std::setprecision(6) << std::setw(9) << nb.velocities[i];
out << '\n';
}
return out;
}
int main() {
std::string fileName = "nbody.txt";
NBody nb(fileName);
for (int i = 0; i < nb.getTimeSteps(); ++i) {
std::cout << "\nCycle " << i + 1 << '\n';
nb.simulate();
std::cout << nb;
}
return 0;
}</syntaxhighlight>
{{out}}
<pre>Contents of nbody.txt
0.01 3 20
1
0.000000 0.000000 0.000000
0.010000 0.000000 0.000000
0.100000
1.000000 1.000000 0.000000
0.000000 0.000000 0.020000
0.001000
0.000000 1.000000 1.000000
0.010000 -0.010000 -0.010000
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</pre>
=={{header|D}}==
{{trans|Kotlin}}
<
import std.array;
import std.conv;
Line 611 ⟶ 1,468:
nb.printResults();
}
}</
{{out}}
<pre>Contents of nbody.txt
Line 779 ⟶ 1,636:
More seriously, later Fortran has additional routines, and confusion can arise: thus subroutine FREE was renamed to SFREE, and array INDEX renamed to INDEXS, though calling it SINDEX was the first idea - but that would require additional declarations so as to keep it integer. The habit of not bothering with minor declarations lingers on with exploratory programmes...
<syntaxhighlight lang="fortran">
SUBROUTINE PLOTS(DX,DY,DEVICE)
INTEGER DEVICE(10)
Line 1,398 ⟶ 2,255:
9001 FORMAT (' Closest approach:',E15.6)
END
</syntaxhighlight>
===Results===
Line 1,517 ⟶ 2,374:
The idea was to start with an initial state having the Sun, Earth and Moon all in a straight line along the x-axis, but after some ad-hoc messing about, confusions escalated to the degree that an ad-hoc calculation prog. was in order. This turns out to need a cube root function, and a proper solution for this raises again the utility of palindromic functions. Later Fortran supplies intrinsic functions such as EXPONENT(x) which returns the exponent part of the floating-point number, ''x'' and extracting this would be useful in devising an initial value for an iterative refinement calculation. Clearly, one wants power/3 for this, and so being able to write something like <code>EXPONENT(t) = EXPONENT(x)/3</code> would help, along with similar usage of the FRACTION(x) function - omitting details such as the remainder when the power is divided by three. The same ideas arise with the square root function, though here, SQRT is already supplied. Alas, only the SUBSTR intrinsic function of pl/i has palindromic usage.
<syntaxhighlight lang="fortran">
Calculate some parameters for the solar system of Sun, Earth and Moon.
IMPLICIT REAL*8 (A-Z) !No integers need apply.
Line 1,587 ⟶ 2,444:
WRITE (6,*) 2*PI,"2Pi: distance per year for R = 1."
END
</syntaxhighlight>
Notably, wanting a circular orbit to ease inspection of the results meant that the actual semi-major axis of the earth's elliptical orbit could not be used as in the first trials. Instead, the orbital period was taken as a year (wrongly so, as this includes the precession of the earth's axis of revolution, with a period of about 26,000 years) and the radius of a circular orbit having that period determined. This however is affected by whether the calculation is for the sun alone (so, a massless Earth), or, the sun and the earth together, or the sun, earth and moon together. As well, there is no z-action: the moon's orbital plane is deemed the same as that of the earth and it isn't. Results:
<pre>
Line 1,762 ⟶ 2,619:
=={{header|Go}}==
{{trans|C}}
<
import (
Line 1,883 ⟶ 2,740:
}
}
}</
Contents of nbody.txt:
Line 2,024 ⟶ 2,881:
<li>Velocity is { vx, vy, vz }.</li>
</ul>
<
I =: 0&{"1
M =: 1&{"1
Line 2,057 ⟶ 2,914:
out =: out,D z
z
)</
<h4>Integrator</h4>
<p>Iterate over a time interval. Plot the results (x and y coordinates over time). </p>
<
maxn =: 20000
require'plot'
Line 2,069 ⟶ 2,926:
dt NEXT^:x y
plot 0 1 { |: out
)</
<h4>Configurations</h4>
<h5>Configuration generator</h5>
<p>Generate a radially symmetrical configuration for a given number of bodies and an initial rotational velocity</p>
<
GEN =: 3 : 0
1 GEN y
Line 2,083 ⟶ 2,940:
(i.y),.(y#m),.p,.v
)
</syntaxhighlight>
<h5>Generate a few cases.</h5>
<p>The 3-body static equilibrium case is proposed as the basis for validation, and for comparison of alternative methods.</p>
<
case =. 3 : 'y;".y'
r =. 0 2 $ a:
Line 2,108 ⟶ 2,965:
¦ ¦2 1 _0.5 _0.866025 0 0.658037 _0.379918 0 ¦
+------------------------------------------------------------------------------+
</syntaxhighlight>
<h5>Static equilibrium case. </h5>
<p> Bodies orbit with constant kinetic energy.</p>
<
[<em>sinusoidal curves of constant amplitude</em>[https://commons.wikimedia.org/wiki/File:Nbody_j_2.jpg]]
<p>Determine orbital period, and compare configuration after two rotations vs. initial position. </p>
<syntaxhighlight lang="j">
v =: (1%3)^(1%4)
dt =: 0.001
Line 2,128 ⟶ 2,985:
¦ 0.080507 _0.0400152 0¦
+------------------------+
</syntaxhighlight>
<h5>Dynamic equilibrium case. </h5>
<p> Within a narrow range of initial velocities, the bodies orbit, maintaining symetry and swapping kinetic and potential energy back and forth. The ideal system is stable, but an integrating solver will not be, due to numerical precision. </p>
<
[<em>sinusoidal curves of varying amplitudes</em> [https://commons.wikimedia.org/wiki/File:Nbody_j_1.jpg]]
<h5>Perturbed case. </h5>
<p>Any perturbation from symmetry, however small, will eventually cause the bodies to escape. Note that the effect demonstrated here is due to the physics of the problem, independently of the numerical precision and increment size of the solver. </p>
<
[<em>curves start out nice then get jumbled and diverge</em> [https://commons.wikimedia.org/wiki/File:Nbody_j_3.jpg]]
<h5>n>3</h5>
<p>We can generate and execute tests for radially symmetrical configurations of any number of bodies, run the simulations, and compare the results with an explicitly determined result:</p>
<
f0 =. 4 : '1%((x r y)^2)' NB. force due to nth body
f =. 4 : '(x f0 y)*(1-cos 2*pi*y%x)%x r y'"0 NB. centripetal component of force
Line 2,191 ⟶ 3,048:
+--+--------+-------+---------¦
¦19¦3.04673 ¦2.06227¦0.0916438¦
+-----------------------------+</
<h5>Sensitivity to discretization</h5>
<p>Although the static-equilibrium configuration is stable, its
Line 2,197 ⟶ 3,054:
can examine the sensitivity to time increment, taking position error after
one rotation as a stability metric:</p>
<
dt =: y
GENDT0 3
Line 2,223 ⟶ 3,080:
+------+----------¦
¦1 ¦6.0354 ¦
+-----------------+</
=={{header|Java}}==
{{trans|Kotlin}}
<
import java.nio.file.Files;
import java.nio.file.Path;
Line 2,376 ⟶ 3,233:
}
}
}</
{{out}}
<pre>Contents of nbody.txt
Line 2,491 ⟶ 3,348:
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</pre>
=={{header|jq}}==
'''Works with jq and gojq, that is, the C and Go implementations of jq.'''
'''Adapted from [[#Java|Java]]'''
Warning: this adaptation from the Java entry does not address the inadequacy of the physics behind `resolveCollisions`; fortunately, in the simulation,
there are no collisions anyway.
'''Preliminaries'''
<syntaxhighlight lang=jq>
def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;
def rpad($len): tostring | ($len - length) as $l | .+ ("0" * $l)[:$l];
# Input: a non-negative number or a string representation of a non-negative number.
# Output: if tostring has no "e" or "E", adjust the number of
# decimal digits to be $n, otherwise adjust $n to allow room for the exponent.
def align_decimal($n):
tostring
| if test("e|E")
then capture( "^(?<x>[^eE]*)[eE](?<e>.*)$" )
| (.e | length + 1) as $e
| (.x | align_decimal($n - $e)) + "e" + .e
elif index(".")
then capture("(?<i>[0-9]*[.])(?<j>[0-9]{0," + ($n|tostring) + "})")
| .i + (.j|rpad($n))
else . + "." + ("0" * $n)
end ;
def Vector3D($x; $y; $z):
{$x, $y, $z};
def Vector3D_add($v):
.x += $v.x | .y += $v.y | .z += $v.z;
def Vector3D_add:
. as $in
| reduce range(1; length) as $i (.[0]; Vector3D_add($in[$i])) ;
# $a - $v
def Vector3D_minus($a; $v):
$a
| .x -= $v.x | .y -= $v.y | .z -= $v.z ;
# $v should be a number
def Vector3D_mult($v):
.x *= $v | .y *= $v | .z *= $v ;
# input: Vector3D
def norm: # i.e. mod
(.x * .x + .y * .y + .z * .z) | sqrt;
def Origin: Vector3D(0; 0; 0);
# input: a string representing three numbers
# aka decompose
def toVector3D:
[splits(" *") | tonumber]
| Vector3D(.[0]; .[1]; .[2]);
</syntaxhighlight>
'''The Task'''
<syntaxhighlight lang=jq>
# Read the parameters assuming invocation of jq includes the -n option
# Emit a JSON object {masses, positions, velocities, accelerations, lines, ...}
def NBodyProblem:
[inputs] as $lines
| [$lines[0] | splits(" *") | tonumber] as $gbt
| $gbt[1] as $bodies
| {gc: $gbt[0],
$bodies,
timeSteps: $gbt[2],
masses: null,
positions: null,
velocities: null,
accelerations: null,
$lines}
| reduce range(0;$bodies) as $i (.;
.masses[$i] = ($lines[$i * 3 + 1] | tonumber)
| .positions[$i] = ($lines[$i * 3 + 2] | toVector3D)
| .velocities[$i] = ($lines[$i * 3 + 3] | toVector3D) ) ;
# {bodies, positions, velocities}
def resolveCollisions:
reduce range(0; .bodies) as $i (.;
reduce range($i + 1; .bodies) as $j (.;
if .positions[$i] == .positions[$j]
then .velocities[$i] as $temp
| .velocities[$i] = .velocities[$j]
| .velocities[$j] = $temp
else .
end) );
# input {bodies, positions, accelarations}
def computeAccelerations:
.bodies as $bodies
| reduce range(0; $bodies) as $i (.;
.accelerations[$i] = Origin
| reduce range(0; $bodies) as $j (.;
if $i != $j
then (.gc * .masses[$j] /
pow( Vector3D_minus(.positions[$i]; .positions[$j]) | norm; 3) ) as $temp
| .accelerations[$i] = (
[ .accelerations[$i],
(Vector3D_minus(.positions[$j]; .positions[$i]) | Vector3D_mult($temp)) ]
| Vector3D_add )
else .
end ));
def computeVelocities:
reduce range(0; .bodies) as $i (.;
.velocities[$i] = ([.velocities[$i], .accelerations[$i] ] | Vector3D_add) );
def computePositions:
reduce range(0; .bodies) as $i (.;
.positions[$i] = ([.positions[$i], .velocities[$i], (.accelerations[$i] | Vector3D_mult(0.5) )]
| Vector3D_add) ) ;
def simulate:
computeAccelerations
| computePositions
| computeVelocities
| resolveCollisions;
def printResults:
def p:
tostring
| if startswith("-")
then "-" + (.[1:] | align_decimal(6))
else " " + align_decimal(6)
end
| lpad(8);
range(0; .bodies) as $i
| "Body \($i + 1) : \(.positions[$i].x |p) \(.positions[$i].y |p) \(.positions[$i].z |p) | "
+ "\(.velocities[$i].x |p) \(.velocities[$i].y |p) \(.velocities[$i].z|p)" ;
def prelude:
"Contents of input file",
.lines[],
"",
"Body : x y z | vx vy vz";
def task:
NBodyProblem
| (prelude,
foreach range (1; 1 + .timeSteps) as $i (.;
simulate;
"\nCycle \($i)", printResults) ) ;
task
</syntaxhighlight>
'''Invocation: gojq -nrRf n-body-problem.jq nbody.txt'''
{{output}}
<pre>
Contents of input file
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.010176 0.000178 0.000001 | 0.010353 0.000357 0.000003
Body 2 : 0.998230 0.998232 0.020001 | -0.003539 -0.003535 0.020003
Body 3 : 0.010176 0.988232 0.988055 | 0.010353 -0.013535 -0.013889
Cycle 2
Body 1 : 0.020708 0.000717 0.000010 | 0.010710 0.000720 0.000014
Body 2 : 0.992906 0.992895 0.039971 | -0.007108 -0.007137 0.019935
Body 3 : 0.020716 0.972887 0.972172 | 0.010726 -0.017153 -0.017876
Cycle 3
Body 1 : 0.031599 0.001625 0.000034 | 0.011072 0.001094 0.000033
Body 2 : 0.983984 0.983909 0.059833 | -0.010735 -0.010834 0.019789
Body 3 : 0.031642 0.953868 0.952235 | 0.011124 -0.020886 -0.021998
Cycle 4
Body 1 : 0.042857 0.002912 0.000081 | 0.011443 0.001480 0.000060
Body 2 : 0.971393 0.971162 0.079509 | -0.014447 -0.014659 0.019561
Body 3 : 0.042981 0.931039 0.928086 | 0.011552 -0.024771 -0.026298
Cycle 5
Body 1 : 0.054492 0.004594 0.000159 | 0.011825 0.001883 0.000097
Body 2 : 0.955030 0.954508 0.098908 | -0.018278 -0.018648 0.019238
Body 3 : 0.054766 0.904224 0.899522 | 0.012017 -0.028856 -0.030829
Cycle 6
Body 1 : 0.066517 0.006690 0.000280 | 0.012223 0.002308 0.000145
Body 2 : 0.934759 0.933759 0.117930 | -0.022264 -0.022849 0.018806
Body 3 : 0.067040 0.873197 0.866280 | 0.012529 -0.033198 -0.035654
Cycle 7
Body 1 : 0.078950 0.009224 0.000456 | 0.012642 0.002759 0.000206
Body 2 : 0.910399 0.908677 0.136455 | -0.026453 -0.027315 0.018244
Body 3 : 0.079855 0.837662 0.828022 | 0.013101 -0.037871 -0.040860
Cycle 8
Body 1 : 0.091814 0.012226 0.000701 | 0.013086 0.003245 0.000284
Body 2 : 0.881722 0.878958 0.154339 | -0.030902 -0.032121 0.017523
Body 3 : 0.093281 0.797239 0.784312 | 0.013749 -0.042975 -0.046559
Cycle 9
Body 1 : 0.105139 0.015736 0.001035 | 0.013563 0.003774 0.000382
Body 2 : 0.848428 0.844216 0.171401 | -0.035684 -0.037361 0.016600
Body 3 : 0.107404 0.751426 0.734579 | 0.014498 -0.048649 -0.052908
Cycle 10
Body 1 : 0.118963 0.019805 0.001481 | 0.014084 0.004361 0.000508
Body 2 : 0.810136 0.803952 0.187408 | -0.040899 -0.043166 0.015413
Body 3 : 0.122345 0.699554 0.678055 | 0.015383 -0.055096 -0.060138
Cycle 11
Body 1 : 0.133337 0.024498 0.002071 | 0.014662 0.005024 0.000671
Body 2 : 0.766343 0.757509 0.202049 | -0.046686 -0.049720 0.013868
Body 3 : 0.138267 0.640689 0.613685 | 0.016460 -0.062632 -0.068602
Cycle 12
Body 1 : 0.148326 0.029906 0.002851 | 0.015316 0.005791 0.000887
Body 2 : 0.716376 0.703998 0.214888 | -0.053245 -0.057301 0.011809
Body 3 : 0.155405 0.573482 0.539940 | 0.017815 -0.071781 -0.078886
Cycle 13
Body 1 : 0.164024 0.036156 0.003887 | 0.016079 0.006708 0.001184
Body 2 : 0.659310 0.642172 0.225281 | -0.060886 -0.066350 0.008976
Body 3 : 0.174111 0.495835 0.454475 | 0.019596 -0.083510 -0.092044
Cycle 14
Body 1 : 0.180564 0.043437 0.005285 | 0.016999 0.007852 0.001612
Body 2 : 0.593807 0.570186 0.232207 | -0.070119 -0.077621 0.004875
Body 3 : 0.194928 0.404135 0.353320 | 0.022038 -0.099889 -0.110265
Cycle 15
Body 1 : 0.198149 0.052049 0.007233 | 0.018170 0.009371 0.002282
Body 2 : 0.517817 0.485100 0.233877 | -0.081861 -0.092550 -0.001534
Body 3 : 0.218604 0.290860 0.228582 | 0.025313 -0.126660 -0.139210
Cycle 16
Body 1 : 0.217125 0.062542 0.010117 | 0.019781 0.011614 0.003484
Body 2 : 0.427899 0.381659 0.226653 | -0.097974 -0.114332 -0.012913
Body 3 : 0.244268 0.131955 0.057562 | 0.026012 -0.191148 -0.202830
Cycle 17
Body 1 : 0.238345 0.076539 0.015220 | 0.022657 0.016380 0.006723
Body 2 : 0.317489 0.248501 0.200966 | -0.122846 -0.151982 -0.038460
Body 3 : 0.075591 -0.559590 -0.487315 | -0.363365 -1.191945 -0.886924
Cycle 18
Body 1 : 0.263123 0.097523 0.026917 | 0.026897 0.025587 0.016671
Body 2 : 0.173428 0.050423 0.112715 | -0.165275 -0.244174 -0.138041
Body 3 : -0.286240 -1.745596 -1.369527 | -0.360299 -1.180066 -0.877501
Cycle 19
Body 1 : 0.270854 0.113045 0.061923 | -0.011436 0.005456 0.053338
Body 2 : 0.199821 -0.093104 -0.208666 | 0.218061 -0.042881 -0.504723
Body 3 : -0.646318 -2.924909 -2.246453 | -0.359855 -1.178559 -0.876350
Cycle 20
Body 1 : 0.258571 0.116045 0.112037 | -0.013129 0.000543 0.046890
Body 2 : 0.426345 -0.111425 -0.681150 | 0.234987 0.006240 -0.440245
Body 3 : -1.006088 -4.103186 -3.122590 | -0.359685 -1.177995 -0.875924
</pre>
=={{header|Julia}}==
Uses the NBodySimulator module.
<syntaxhighlight lang="julia">using StaticArrays, Plots, NBodySimulator
const stablebodies = [MassBody(SVector(0.0, 1.0, 0.0), SVector( 5.775e-6, 0.0, 0.0), 2.0),
MassBody(SVector(0.0,-1.0, 0.0), SVector(-5.775e-6, 0.0, 0.0), 2.0)]
const bodies = [
MassBody(SVector(0.0, 1.0, 0.0), SVector( 5.775e-6, 0.0, 0.0), 2.0),
MassBody(SVector(0.0,-1.0, 0.0), SVector(-5.775e-6, 0.0, 0.0), 2.0),
MassBody(SVector(0.0, 4.5, 0.0), SVector(-2.5e-6, 0.0, 0.0), 1.0)
]
const G = 6.673e-11 # m^3/kg/s^2
const timespan = (0.0, 1111150.0)
function nbodysim(nbodies, tspan)
system = GravitationalSystem(nbodies, G)
simulation = NBodySimulation(system, tspan)
run_simulation(simulation)
end
simresult = nbodysim(bodies, timespan)
plot(simresult)
</syntaxhighlight>
=={{header|K}}==
Line 2,496 ⟶ 3,647:
<h4>configuration generator</h4>
<p>The configurations shown here are based on those in the J task.</p>
<
cs:{(cos 2*pi*x%y;sin 2*pi*x%y)} / positions
sc:{(sin 2*pi*x%y;-cos 2*pi*x%y)} / velocities</
<h4>static equilibrium configuration</h4>
<p>The system is stable for an initial velocity for which centrifugal and gravitational forces are in balance</p>
<
g::1;vs::pow[1%3;1%4];v*::vs / gravitational constant, initial velocity
t::0;dt::0.01;rot:2 / time, time increment, rotational period</
<h4>dynamic equilibrium</h4>
<p>Within a narrow range of initial velocities the system oscillates between more and less kinetic and potential energy.</p>
<
/ v*::0.7;rot::2</
<h4>unstable</h4>
<p>The slightest perturbation from symmetry causes the system to become unstable.</p>
<
/ v+::(0 0;0 0;0 0.0001);rot:4 / wait for it ... wait for it ...</
<h4>physics</h4>
<
A1:{$[x=y;0;g*m[y]*DV[x;y]%pow[DS[x;y];3]]} / nth body centripetal acceleration component
A:{+/{A1/:[x;!n]}'!n} / total centripetal acceleration
Line 2,520 ⟶ 3,671:
v0:v;v+::dt*A[];p+::dt*0.5*v+v0 / update velocities and positions
t+::dt;$[t>tmax;do[msg::,"[";tick::{}];0] / update time, detect end of time
}</
<h4>events and graphics</h4>
<
sc:1.25;P::{w*0.5*1+(1%sc)*x} / display units
wh::(w;h);C::wh%2 / viewbox
Line 2,532 ⟶ 3,683:
r:(,:'P[p]-s%2),'`RGB,',:'(2#'s)#'!n / plot the bodies
r,:,(P[1 0];BW;~,/'+text@`i$msg) / plot the initial position marker
}</
<h4>3D</h4>
<p>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:</p>
<
sc:{(0;sin 2*pi*x%y;-cos 2*pi*x%y)} / objects on screen move up and down
Line 2,544 ⟶ 3,695:
cs:{(cos 2*pi*x%y;sin 2*pi*x%y;0)} / z-axis
sc:{(sin 2*pi*x%y;-cos 2*pi*x%y;0)} / objects on screen go round and round</
=={{header|Kotlin}}==
{{trans|C}}
<
import java.io.File
Line 2,669 ⟶ 3,820:
nb.printResults()
}
}</
{{out}}
Line 2,787 ⟶ 3,938:
Body 3 : -1.006089 -4.103186 -3.122591 | -0.359686 -1.177995 -0.875924
</pre>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">data = NBodySimulation[
"InverseSquare", {<|"Mass" -> 1, "Position" -> {0, 0},
"Velocity" -> {0, .5}|>,
<|"Mass" -> 1, "Position" -> {1, 1}, "Velocity" -> {0, -.5}|>,
<|"Mass" -> 1, "Position" -> {0, 1}, "Velocity" -> {0, 0}|>}, 5];
{"Time:", N[#], "Positions: ", data[All, "Position", #], "Velocities: ", data[All, "Velocity", #]} & /@ Subdivide[0, 5, 20] // Grid
ParametricPlot[Evaluate[data[All, "Position", t]], {t, 0, 4}]</syntaxhighlight>
{{out}}
<pre>Time: 0. Positions: {{0.,0.},{1.,1.},{0.,1.}} Velocities: {{0.,0.5},{0.,-0.5},{0.,0.}}
Time: 0.25 Positions: {{0.0130044,0.171909},{0.955177,0.864576},{0.0318183,0.963515}} Velocities: {{0.115009,0.902202},{-0.374278,-0.581613},{0.259269,-0.32059}}
Time: 0.5 Positions: {{0.0717548,0.488983},{0.79434,0.709843},{0.133905,0.801174}} Velocities: {{0.446796,1.88926},{-0.995361,-0.652841},{0.548564,-1.23642}}
Time: 0.75 Positions: {{0.349847,0.612073},{0.459794,0.409715},{0.190359,0.978212}} Velocities: {{-2.08701,0.957204},{1.09683,-1.90053},{0.99018,0.943324}}
Time: 1. Positions: {{-0.0171517,0.909988},{0.611433,0.132047},{0.405719,0.957965}} Velocities: {{-0.702007,1.50372},{0.386666,-0.728162},{0.315341,-0.775562}}
Time: 1.25 Positions: {{-0.0541207,1.23874},{0.68204,0.011408},{0.372081,0.749849}} Velocities: {{0.222905,1.0525},{0.18142,-0.246531},{-0.404325,-0.805968}}
Time: 1.5 Positions: {{0.039104,1.43964},{0.696163,0.00833041},{0.264733,0.552029}} Velocities: {{0.468249,0.584491},{-0.0929974,0.227301},{-0.375251,-0.811791}}
Time: 1.75 Positions: {{0.166265,1.54363},{0.608954,0.133385},{0.224782,0.322983}} Velocities: {{0.537312,0.264809},{-0.734718,0.808666},{0.197406,-1.07347}}
Time: 2. Positions: {{0.304798,1.57616},{0.297257,0.0619309},{0.397945,0.361908}} Velocities: {{0.565787,-0.00463666},{0.181379,-1.31263},{-0.747166,1.31727}}
Time: 2.25 Positions: {{0.445248,1.53798},{0.360772,-0.0914083},{0.19398,0.553427}} Velocities: {{0.544779,-0.309614},{0.220404,-0.167661},{-0.765182,0.477274}}
Time: 2.5 Positions: {{0.568905,1.41861},{0.397179,-0.0644824},{0.0339158,0.645877}} Velocities: {{0.426233,-0.647816},{0.0663481,0.354581},{-0.492581,0.293235}}
Time: 2.75 Positions: {{0.648579,1.21279},{0.391808,0.0824137},{-0.0403866,0.704793}} Velocities: {{0.188043,-1.00485},{-0.115952,0.83295},{-0.0720915,0.171896}}
Time: 3. Positions: {{0.64316,0.905955},{0.331764,0.376278},{0.0250761,0.717767}} Velocities: {{-0.318522,-1.50449},{-0.400444,1.65218},{0.718967,-0.147697}}
Time: 3.25 Positions: {{0.409627,0.623846},{0.255008,0.415972},{0.335365,0.960183}} Velocities: {{1.00799,0.824653},{-1.95929,-2.06579},{0.951299,1.24114}}
Time: 3.5 Positions: {{0.606343,0.868804},{-0.122419,0.069204},{0.516076,1.06199}} Velocities: {{0.865322,-0.797808},{-1.29972,-1.0888},{0.4344,1.88661}}
Time: 3.75 Positions: {{0.706052,0.996203},{-0.423763,-0.174935},{0.717711,1.17873}} Velocities: {{0.19751,2.02643},{-1.1312,-0.89257},{0.933693,-1.13386}}
Time: 4. Positions: {{0.869387,1.03249},{-0.694252,-0.384878},{0.824865,1.35239}} Velocities: {{0.193378,0.789463},{-1.03984,-0.796168},{0.846463,0.00670532}}
Time: 4.25 Positions: {{1.01635,1.12428},{-0.946242,-0.575969},{0.929894,1.45169}} Velocities: {{0.295969,0.0882097},{-0.979766,-0.736675},{0.683797,0.648465}}
Time: 4.5 Positions: {{1.14204,1.26545},{-1.18549,-0.754722},{1.04345,1.48927}} Velocities: {{0.565039,-0.766966},{-0.936353,-0.695536},{0.371314,1.4625}}
Time: 4.75 Positions: {{1.19611,1.39678},{-1.41526,-0.924622},{1.21915,1.52785}} Velocities: {{0.152075,2.47452},{-0.903116,-0.665026},{0.75104,-1.8095}}
Time: 5. Positions: {{1.33622,1.3904},{-1.63761,-1.08779},{1.3014,1.69739}} Velocities: {{0.0891981,0.848018},{-0.876639,-0.641268},{0.787441,-0.20675}}
[outputs a graphics object showing the trajectories of the masses]</pre>
=={{header|Nim}}==
{{trans|Go}}
<syntaxhighlight lang="nim">import math, os, strformat, strscans, strutils
type Vector = tuple[x, y, z: float]
func`+`(v1, v2: Vector): Vector = (v1.x + v2.x, v1.y + v2.y, v1.z + v2.z)
func `+=`(v1: var Vector; v2: Vector) =
v1.x += v2.x
v1.y += v2.y
v1.z += v2.z
func`-`(v1, v2: Vector): Vector = (v1.x - v2.x, v1.y - v2.y, v1.z - v2.z)
func `*`(v: Vector; m: float): Vector = (v.x * m, v.y * m, v.z * m)
func abs(v: Vector): float = sqrt(v.x * v.x + v.y * v.y + v.z * v.z)
type Simulation = object
bodies: int
timeSteps: int
masses: seq[float]
gc: float
positions: seq[Vector]
velocities: seq[Vector]
accelerations: seq[Vector]
proc emitError(linenum: Positive) =
raise newException(ValueError, "wrong data at line $#.".format(linenum))
proc initSimulation(fileName: string): Simulation =
let infile = fileName.open()
var line = infile.readLine()
if not line.scanf("$f $i $i", result.gc, result.bodies, result.timeSteps):
emitError(1)
result.masses.setlen(result.bodies)
result.positions.setLen(result.bodies)
result.velocities.setLen(result.bodies)
result.accelerations.setLen(result.bodies)
var linenum = 1
for i in 0..<result.bodies:
inc linenum
line = infile.readLine()
if not line.scanf("$f", result.masses[i]):
emitError(linenum)
inc linenum
line = infile.readLine()
if not line.scanf("$f $f $f", result.positions[i].x, result.positions[i].y, result.positions[i].z):
emitError(linenum)
inc linenum
line = infile.readLine()
if not line.scanf("$f $f $f", result.velocities[i].x, result.velocities[i].y, result.velocities[i].z):
emitError(linenum)
func resolveCollisions(sim: var Simulation) =
for i in 0..sim.bodies-2:
for j in i+1..sim.bodies-1:
if sim.positions[i] == sim.positions[j]:
swap sim.velocities[i], sim.velocities[j]
func computeAccelerations(sim: var Simulation) =
for i in 0..<sim.bodies:
sim.accelerations[i] = (0.0, 0.0, 0.0)
for j in 0..<sim.bodies:
if i != j:
let m = sim.gc * sim.masses[j] / abs(sim.positions[i] - sim.positions[j])^3
sim.accelerations[i] += (sim.positions[j] - sim.positions[i]) * m
func computeVelocities(sim: var Simulation) =
for i in 0..<sim.bodies:
sim.velocities[i] += sim.accelerations[i]
func computePositions(sim: var Simulation) =
for i in 0..<sim.bodies:
sim.positions[i] += sim.velocities[i] + sim.accelerations[i] * 0.5
func step(sim: var Simulation) =
sim.computeAccelerations()
sim.computePositions()
sim.computeVelocities()
sim.resolveCollisions()
proc printResults(sim: Simulation) =
for i in 0..<sim.bodies:
echo &"Body {i + 1}: ",
&"{sim.positions[i].x: 8.6f} {sim.positions[i].y: 8.6f} {sim.positions[i].z: 8.6f} | ",
&"{sim.velocities[i].x: 8.6f} {sim.velocities[i].y: 8.6f} {sim.velocities[i].z: 8.6f}"
if paramCount() != 1:
echo "Usage: $# <file name containing system configuration data>" % getAppFilename().lastPathPart
else:
var sim = initSimulation(paramStr(1))
echo "Body x y z | vx vy vz"
for step in 1..sim.timeSteps:
echo "\nCycle ", step
sim.step()
sim.printResults()</syntaxhighlight>
{{out}}
Using same data as Go.
<pre>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</pre>
=={{header|Octave}}==
{{trans|
We'll show only the positions in the output.
<
global au = 150e9; # astronomical unit
Line 2,834 ⟶ 4,228:
disp(lsode('ABCdot',ABC, t)(1:20,1:9));
</syntaxhighlight>
{{out}}
<pre> 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
Line 2,859 ⟶ 4,253:
=={{header|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.
<syntaxhighlight lang="pascal">
{$N+ Crunch only the better sort of numbers, thanks.}
{$M 50000,24000,24000] {Stack,minheap,maxheap.}
Line 3,707 ⟶ 5,101:
TextMode(AsItWas);
END.
</syntaxhighlight>
===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.
<syntaxhighlight lang="pascal">
{$N+ Crunch only the better sort of numbers, thanks.}
Program Swirl; Uses Graph, Crt;
Line 3,967 ⟶ 5,361:
WriteLn('Calculation steps:',nstep);
END.
</syntaxhighlight>
=={{header|Perl
Cycle through one Neptunian year.
<syntaxhighlight lang="perl">use strict;
use warnings;
use constant PI => 3.141592653589793;
use constant SOLAR_MASS => (4 * PI * PI);
use constant YEAR => 365.24;
sub solar_offset {
my($vxs, $vys, $vzs, $mass) = @_;
my($px, $py, $pz);
for (0..@$mass-1) {
$px += @$vxs[$_] * @$mass[$_];
$py += @$vys[$_] * @$mass[$_];
$pz += @$vzs[$_] * @$mass[$_];
}
(@$vxs[0], @$vys[0], @$vzs[0]) = ((-$px/SOLAR_MASS), (-$py/SOLAR_MASS), (-$pz/SOLAR_MASS));
}
sub orbit {
my ($dt, $xs, $ys, $zs, $vxs, $vys, $vzs, $mass) = @_;
for (0..@$mass-1) {
for (my $j = $_ + 1; $j < @$mass; $j++) {
my($dx, $dy, $dz) = ((@$xs[$_] - @$xs[$j]), (@$ys[$_] - @$ys[$j]), (@$zs[$_] - @$zs[$j]));
my $mag = $dt / ($dx**2 + $dy**2 + $dz**2)**(3/2);
my($mm, $mm2) = ((@$mass[$_] * $mag), (@$mass[$j] * $mag));
@$vxs[$_] -= $dx * $mm2;
@$vxs[$j] += $dx * $mm;
@$vys[$_] -= $dy * $mm2;
@$vys[$j] += $dy * $mm;
@$vzs[$_] -= $dz * $mm2;
@$vzs[$j] += $dz * $mm;
}
@$xs[$_] += $dt * @$vxs[$_];
@$ys[$_] += $dt * @$vys[$_];
@$zs[$_] += $dt * @$vzs[$_];
}
}
sub display {
my($t,$xs,$ys) = @_;
printf '%6d', $t;
printf '%7.2f' x 2, @$xs[$_],@$ys[$_] for 1..4;
print "\n";
}
my @ns = <Sun Jupiter Saturn Uranus Neptune >;
my @xs = <0 4.84143144e+00 8.34336671e+00 1.28943695e+01 1.53796971e+01>;
my @ys = <0 -1.16032004e+00 4.12479856e+00 -1.51111514e+01 -2.59193146e+01>;
my @zs = <0 -1.03622044e-01 -4.03523417e-01 -2.23307578e-01 1.79258772e-01>;
my @vxs = map {$_ * YEAR} <0 1.66007664e-03 -2.76742510e-03 2.96460137e-03 2.68067772e-03>;
my @vys = map {$_ * YEAR} <0 7.69901118e-03 4.99852801e-03 2.37847173e-03 1.62824170e-03>;
my @vzs = map {$_ * YEAR} <0 -6.90460016e-05 2.30417297e-05 -2.96589568e-05 -9.51592254e-05>;
my @mass = map {$_ * SOLAR_MASS} <1 9.54791938e-04 2.85885980e-04 4.36624404e-05 5.15138902e-05>;
solar_offset(\@vxs, \@vys, \@vzs, \@mass);
printf ' t' . '%14s'x4 . "\n", @ns[1..4];
display(0, \@xs, \@ys);
my $steps = 16567;
for my $t (1 .. $steps) {
orbit(0.01, \@xs, \@ys, \@zs, \@vxs, \@vys, \@vzs, \@mass);
display($t,\@xs, \@ys) unless $t % 1000;
}
display($steps, \@xs, \@ys);</syntaxhighlight>
{{out}}
<pre> Jupiter Saturn Uranus Neptune
t x y x y x y x y
0 4.84 -1.16 8.34 4.12 12.89 -15.11 15.38 -25.92
1000 1.42 -5.00 -8.75 3.29 19.78 -3.52 23.83 -18.26
2000 -3.44 -4.20 0.84 -10.06 17.37 9.73 28.84 -7.96
3000 -5.46 0.16 7.88 4.86 6.60 18.23 29.65 3.50
4000 -2.98 4.35 -9.06 2.44 -7.48 17.22 26.11 14.44
5000 2.11 4.56 1.68 -9.94 -17.23 6.52 18.74 23.26
6000 4.95 0.15 7.30 5.63 -16.69 -8.04 8.63 28.68
7000 2.57 -4.46 -9.31 1.56 -6.34 -17.96 -2.73 29.94
8000 -2.40 -4.82 2.49 -9.78 7.59 -18.13 -13.71 26.90
9000 -5.37 -1.03 6.71 6.26 17.81 -9.20 -22.74 20.04
10000 -3.94 3.59 -9.45 0.70 19.61 4.12 -28.58 10.36
11000 0.89 5.00 3.29 -9.51 12.12 15.46 -30.46 -0.75
12000 4.73 1.45 6.00 6.91 -1.36 18.98 -28.13 -11.76
13000 3.55 -3.63 -9.53 -0.20 -13.97 12.16 -21.93 -21.16
14000 -1.24 -5.19 4.06 -9.21 -18.34 -1.74 -12.70 -27.63
15000 -5.02 -2.17 5.29 7.42 -11.77 -14.64 -1.71 -30.27
16000 -4.70 2.63 -9.52 -1.06 1.54 -19.34 9.51 -28.67
16567 3.45 -3.74 -3.24 -9.52 9.24 -17.42 15.38 -25.92</pre>
=={{header|Phix}}==
{{trans|Kotlin}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">origin</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">},</span>
<span style="color: #000000;">nbody</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">"""
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
"""</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">lines</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">split</span><span style="color: #0000FF;">(</span><span style="color: #000000;">nbody</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">),</span>
<span style="color: #0000FF;">{{</span><span style="color: #000000;">grav_constant</span><span style="color: #0000FF;">,</span><span style="color: #000000;">bodies</span><span style="color: #0000FF;">,</span><span style="color: #000000;">timeSteps</span><span style="color: #0000FF;">}}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">scanf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lines</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span><span style="color: #008000;">"%f %f %f"</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">masses</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">bodies</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">positions</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">bodies</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">velocities</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">bodies</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">accelerations</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">bodies</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">bodies</span> <span style="color: #008080;">do</span>
<span style="color: #0000FF;">{{</span><span style="color: #000000;">masses</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]}}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">scanf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lines</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">*</span><span style="color: #000000;">3</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span><span style="color: #008000;">"%f"</span><span style="color: #0000FF;">)</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">positions</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">scanf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lines</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">*</span><span style="color: #000000;">3</span><span style="color: #0000FF;">],</span><span style="color: #008000;">"%f %f %f"</span><span style="color: #0000FF;">)</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">velocities</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">scanf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lines</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">*</span><span style="color: #000000;">3</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span><span style="color: #008000;">"%f %f %f"</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Body : x y z |"</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" vx vy vz\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">vmod</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">computeAccelerations</span><span style="color: #0000FF;">()</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">bodies</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">ai</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">origin</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">bodies</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">j</span> <span style="color: #008080;">then</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">positions</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">],</span><span style="color: #000000;">positions</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">temp</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">grav_constant</span><span style="color: #0000FF;">*</span><span style="color: #000000;">masses</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]/</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">vmod</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">),</span><span style="color: #000000;">3</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">ai</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ai</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">temp</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">accelerations</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">ai</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">computePositions</span><span style="color: #0000FF;">()</span>
<span style="color: #000000;">positions</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">positions</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">velocities</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">accelerations</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">)))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">computeVelocities</span><span style="color: #0000FF;">()</span>
<span style="color: #000000;">velocities</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">velocities</span><span style="color: #0000FF;">,</span><span style="color: #000000;">accelerations</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">resolveCollisions</span><span style="color: #0000FF;">()</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">bodies</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">i</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">bodies</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">positions</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]=</span><span style="color: #000000;">positions</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">then</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">velocities</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span><span style="color: #000000;">velocities</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">velocities</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">],</span><span style="color: #000000;">velocities</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">simulate</span><span style="color: #0000FF;">()</span>
<span style="color: #000000;">computeAccelerations</span><span style="color: #0000FF;">()</span>
<span style="color: #000000;">computePositions</span><span style="color: #0000FF;">()</span>
<span style="color: #000000;">computeVelocities</span><span style="color: #0000FF;">()</span>
<span style="color: #000000;">resolveCollisions</span><span style="color: #0000FF;">()</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">printResults</span><span style="color: #0000FF;">()</span>
<span style="color: #004080;">string</span> <span style="color: #000000;">fmt</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">"Body %d : %9.6f %9.6f %9.6f | %9.6f %9.6f %9.6f\n"</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">bodies</span> <span style="color: #008080;">do</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">i</span><span style="color: #0000FF;">}&</span><span style="color: #000000;">positions</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]&</span><span style="color: #000000;">velocities</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">timeSteps</span> <span style="color: #008080;">do</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\nCycle %d\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">i</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">simulate</span><span style="color: #0000FF;">()</span>
<span style="color: #000000;">printResults</span><span style="color: #0000FF;">()</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<!--</syntaxhighlight>-->
Output same as Python (and C and C# and D and Kotlin and Go)
=={{header|Python}}==
<syntaxhighlight lang="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()</syntaxhighlight>
{{out}}
<pre>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</pre>
=={{header|Raku}}==
(formerly Perl 6)
We'll try to simulate the Sun+Earth+Moon system, with plausible astrometric values.
We use a 18-dimension vector <math>ABC</math>. 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
Line 3,983 ⟶ 5,781:
To keep things compact, we'll only display the first 20 lines of output.
<syntaxhighlight lang="raku"
multi infix:<+>(@a, @b) { @a Z+ @b }
multi infix:<->(@a, @b) { @a Z- @b }
Line 4,047 ⟶ 5,845:
) {
printf "t = %.02f : %s\n", $t, @ABC.fmt("%+.3e");
}</
{{out}}
<pre>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
Line 4,069 ⟶ 5,867:
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</pre>
=={{header|Swift}}==
<syntaxhighlight lang="swift">import Foundation
public struct Vector {
public var px = 0.0
public var py = 0.0
public var pz = 0.0
public init(px: Double, py: Double, pz: Double) {
(self.px, self.py, self.pz) = (px, py, pz)
}
public init?(array: [Double]) {
guard array.count == 3 else {
return nil
}
(self.px, self.py, self.pz) = (array[0], array[1], array[2])
}
public func mod() -> Double {
(px * px + py * py + pz * pz).squareRoot()
}
static func + (lhs: Vector, rhs: Vector) -> Vector {
return Vector(
px: lhs.px + rhs.px,
py: lhs.py + rhs.py,
pz: lhs.pz + rhs.pz
)
}
static func - (lhs: Vector, rhs: Vector) -> Vector {
return Vector(
px: lhs.px - rhs.px,
py: lhs.py - rhs.py,
pz: lhs.pz - rhs.pz
)
}
static func * (lhs: Vector, rhs: Double) -> Vector {
return Vector(
px: lhs.px * rhs,
py: lhs.py * rhs,
pz: lhs.pz * rhs
)
}
}
extension Vector {
public static let origin = Vector(px: 0, py: 0, pz: 0)
}
extension Vector: Equatable {
public static func == (lhs: Vector, rhs: Vector) -> Bool {
return lhs.px == rhs.px && lhs.py == rhs.py && lhs.pz == rhs.pz
}
}
extension Vector: CustomStringConvertible {
public var description: String {
return String(format: "%.6f\t%.6f\t%.6f", px, py, pz)
}
}
public class NBody {
public let gravitationalConstant: Double
public let numBodies: Int
public let timeSteps: Int
public private(set) var masses: [Double]
public private(set) var positions: [Vector]
public private(set) var velocities: [Vector]
public private(set) var accelerations: [Vector]
public init?(file: String) {
guard let data = try? String(contentsOfFile: file) else {
return nil
}
print("Input file:\n\(data)")
let lines = data.components(separatedBy: "\n").map({ $0.components(separatedBy: " ") })
let worldData = lines.first!
guard worldData.count == 3,
let gc = Double(worldData[0]),
let bodies = Int(worldData[1]),
let timeSteps = Int(worldData[2]) else {
return nil
}
let defaultState = Array(repeating: Vector.origin, count: bodies)
self.gravitationalConstant = gc
self.numBodies = bodies
self.timeSteps = timeSteps
self.masses = Array(repeating: 0, count: bodies)
self.positions = defaultState
self.accelerations = defaultState
self.velocities = defaultState
let bodyData = lines.dropFirst().map({ $0.compactMap(Double.init) })
guard bodyData.count == bodies * 3 else {
return nil
}
for n in 0..<bodies {
masses[n] = bodyData[0 + n * 3][0]
guard let position = Vector(array: bodyData[1 + n * 3]),
let velocity = Vector(array: bodyData[2 + n * 3]) else {
return nil
}
positions[n] = position
velocities[n] = velocity
}
}
private func computeAccelerations() {
for i in 0..<numBodies {
accelerations[i] = .origin
for j in 0..<numBodies where i != j {
let t = gravitationalConstant * masses[j] / pow((positions[i] - positions[j]).mod(), 3)
accelerations[i] = accelerations[i] + (positions[j] - positions[i]) * t
}
}
}
private func resolveCollisions() {
for i in 0..<numBodies {
for j in 0..<numBodies where positions[i] == positions[j] {
velocities.swapAt(i, j)
}
}
}
private func computeVelocities() {
for i in 0..<numBodies {
velocities[i] = velocities[i] + accelerations[i]
}
}
private func computePositions() {
for i in 0..<numBodies {
positions[i] = positions[i] + velocities[i] + accelerations[i] * 0.5
}
}
public func printState() {
for i in 0..<numBodies {
print("Body \(i + 1): \(positions[i]) | \(velocities[i])")
}
}
public func simulate() {
computeAccelerations()
computePositions()
computeVelocities()
resolveCollisions()
}
}
guard let sim = NBody(file: "input.txt") else {
fatalError()
}
print()
print("Body : x y z | vx vy vz")
for i in 0..<sim.timeSteps {
print("Step \(i + 1)")
sim.simulate()
sim.printState()
print()
}
</syntaxhighlight>
{{out}}
<pre style="height: 20em; overflow: auto;">Input file:
0.01 3 20
1
0.000000 0.000000 0.000000
0.010000 0.000000 0.000000
0.100000
1.000000 1.000000 0.000000
0.000000 0.000000 0.020000
0.001000
0.000000 1.000000 1.000000
0.010000 -0.010000 -0.010000
Body : x y z | vx vy vz
Step 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
Step 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
Step 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
Step 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
Step 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
Step 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
Step 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
Step 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
Step 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
Step 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
Step 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
Step 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
Step 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
Step 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
Step 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
Step 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
Step 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
Step 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
Step 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
Step 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
Process finished with exit code 0
</pre>
=={{header|Tcl}}==
{{works with|Tcl|8.6}}
<
set G 0.01
Line 4,126 ⟶ 6,226:
puts [lmap pos $p {format (%.5f,%.5f,%.5f) {*}$pos}]
}
}</
Demonstrating for 20 steps:
<
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</
{{out}}
<pre>
Line 4,154 ⟶ 6,254:
(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)
</pre>
=={{header|Wren}}==
{{trans|Kotlin}}
{{libheader|Wren-ioutil}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="wren">import "./ioutil" for FileUtil
import "./fmt" for Fmt
class Vector3D {
construct new(x, y, z) {
_x = x
_y = y
_z = z
}
x { _x }
y { _y }
z { _z }
+(v) { Vector3D.new(_x + v.x, _y + v.y, _z + v.z) }
-(v) { Vector3D.new(_x - v.x, _y - v.y, _z - v.z) }
*(s) { Vector3D.new(s * _x, s * _y, s * _z) }
mod { (_x * _x + _y * _y + _z * _z).sqrt }
}
var Origin = Vector3D.new(0, 0, 0)
class NBody {
construct new(fileName) {
var lines = FileUtil.readLines(fileName)
var gbt = lines[0].split(" ")
_gc = Num.fromString(gbt[0])
_bodies = Num.fromString(gbt[1])
_timeSteps = Num.fromString(gbt[2])
_masses = List.filled(_bodies, 0)
_positions = List.filled(_bodies, null)
_velocities = List.filled(_bodies, null)
_accelerations = List.filled(_bodies, null)
for (i in 0..._bodies) {
_masses[i] = Num.fromString(lines[i * 3 + 1])
_positions[i] = decompose_(lines[i * 3 + 2])
_velocities[i] = decompose_(lines[i * 3 + 3])
}
System.print("Contents of %(fileName)")
System.print(lines.join("\n"))
System.write("Body : x y z |")
System.print(" vx vy vz")
}
timeSteps { _timeSteps }
decompose_(line) {
var xyz = line.split(" ").map { |e| Num.fromString(e) }.toList
return Vector3D.new(xyz[0], xyz[1], xyz[2])
}
resolveCollisions_() {
for (i in 0..._bodies) {
for (j in i + 1..._bodies) {
if (_positions[i].x == _positions[j].x &&
_positions[i].y == _positions[j].y &&
_positions[i].z == _positions[j].z) {
var temp = _velocities[i]
_velocities[i] = _velocities[j]
_velocities[j] = temp
}
}
}
}
computeAccelerations_() {
for (i in 0..._bodies) {
_accelerations[i] = Origin
for (j in 0..._bodies) {
if (i != j) {
var temp = _gc * _masses[j] / (_positions[i] - _positions[j]).mod.pow(3)
_accelerations[i] = _accelerations[i] + (_positions[j] - _positions[i]) * temp
}
}
}
}
computeVelocities_() {
for (i in 0..._bodies) _velocities[i] = _velocities[i] + _accelerations[i]
}
computePositions_() {
for (i in 0..._bodies) {
_positions[i] = _positions[i] + _velocities[i] + _accelerations[i] * 0.5
}
}
simulate() {
computeAccelerations_()
computePositions_()
computeVelocities_()
resolveCollisions_()
}
printResults() {
var fmt = "Body $d : $9.6f $9.6f $9.6f | $9.6f $9.6f $9.6f"
for (i in 0..._bodies) {
Fmt.lprint(fmt,
[
i + 1,
_positions[i].x,
_positions[i].y,
_positions[i].z,
_velocities[i].x,
_velocities[i].y,
_velocities[i].z
]
)
}
}
}
var fileName = "nbody.txt"
var nb = NBody.new(fileName)
for (i in 1..nb.timeSteps) {
System.print("\nCycle %(i)")
nb.simulate()
nb.printResults()
}</syntaxhighlight>
{{out}}
<pre>
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
</pre>
| |||