N-body problem: Difference between revisions
Added Algol 68
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
(Added Algol 68) |
||
| (5 intermediate revisions by 4 users not shown) | |||
Line 27:
F (fileName)
V lines = File(fileName
V gbt = lines[0].split(‘ ’)
.gc = Float(gbt[0])
Line 108:
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
Line 3,164 ⟶ 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}}==
Line 5,316 ⟶ 5,769:
(formerly Perl 6)
We'll try to simulate the Sun+Earth+Moon system, with plausible
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 5,807 ⟶ 6,260:
{{libheader|Wren-ioutil}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="
import "./fmt" for Fmt
class Vector3D {
| |||