nbody
1 Using Class Type | Top |
# The Great Computer Language Shootout
# http://shootout.alioth.debian.org/
# contributed by Isaac Gouy, tuned by Mike Pall
# Ported to Dao by Limin Fu
global const PI = 3.1415926535897932300;
global const SOLAR_MASS = 4 * PI * PI;
global const DAYS_PER_YEAR = 365.2400;
final class OneBody
{
my x, y, z = 0.00;
my vx, vy, vz = 0.00;
my mass = 0.00;
routine Print(){ stdio.println( x, y, z, vx, vy, vz ); }
}
Jupiter = OneBody
{
x => 4.84143144246472090E+00 ;
y => -1.16032004402742839E+00 ;
z => -1.03622044471123109E-01 ;
vx => 1.66007664274403694E-03 * DAYS_PER_YEAR ;
vy => 7.69901118419740425E-03 * DAYS_PER_YEAR ;
vz => -6.90460016972063023E-05 * DAYS_PER_YEAR ;
mass => 9.54791938424326609E-04 * SOLAR_MASS
}
Saturn = OneBody
{
x => 8.34336671824457987E+00 ;
y => 4.12479856412430479E+00 ;
z => -4.03523417114321381E-01 ;
vx => -2.76742510726862411E-03 * DAYS_PER_YEAR ;
vy => 4.99852801234917238E-03 * DAYS_PER_YEAR ;
vz => 2.30417297573763929E-05 * DAYS_PER_YEAR ;
mass => 2.85885980666130812E-04 * SOLAR_MASS
}
Uranus = OneBody{
x => 1.28943695621391310E+01 ;
y => -1.51111514016986312E+01 ;
z => -2.23307578892655734E-01 ;
vx => 2.96460137564761618E-03 * DAYS_PER_YEAR ;
vy => 2.37847173959480950E-03 * DAYS_PER_YEAR ;
vz => -2.96589568540237556E-05 * DAYS_PER_YEAR ;
mass => 4.36624404335156298E-05 * SOLAR_MASS
}
Neptune = OneBody
{
x => 1.53796971148509165E+01 ;
y => -2.59193146099879641E+01 ;
z => 1.79258772950371181E-01 ;
vx => 2.68067772490389322E-03 * DAYS_PER_YEAR ;
vy => 1.62824170038242295E-03 * DAYS_PER_YEAR ;
vz => -9.51592254519715870E-05 * DAYS_PER_YEAR ;
mass => 5.15138902046611451E-05 * SOLAR_MASS
}
Sun = OneBody{ x => 0.00, y => 0.00, z => 0.00,
vx => 0.00, vy => 0.00, vz => 0.00, mass => SOLAR_MASS }
routine advance( bodies : list<OneBody>, nbody : int, dt : double )
{
for( i = 0 : nbody-1 ){
bi = bodies[i];
bix = bi.x; biy = bi.y; biz = bi.z; bimass = bi.mass;
bivx = bi.vx; bivy = bi.vy; bivz = bi.vz;
for( j = i+1 : nbody-1 ){
bj = bodies[j];
dx = bix-bj.x;
dy = biy-bj.y;
dz = biz-bj.z;
distance = (dx*dx + dy*dy + dz*dz)**0.500;
mag = dt / (distance * distance * distance);
bim = bimass * mag;
bjm = bj.mass * mag;
bivx = bivx - (dx * bjm);
bivy = bivy - (dy * bjm);
bivz = bivz - (dz * bjm);
bj.vx = bj.vx + (dx * bim);
bj.vy = bj.vy + (dy * bim);
bj.vz = bj.vz + (dz * bim);
#{
#}
}
bi.vx = bivx;
bi.vy = bivy;
bi.vz = bivz;
}
for( i = 0 : nbody-1 ){
bi = bodies[i];
bi.x = bi.x + (dt * bi.vx);
bi.y = bi.y + (dt * bi.vy);
bi.z = bi.z + (dt * bi.vz);
}
}
routine energy( bodies : list<OneBody> , nbody : int )
{
e = 0.00;
for( i = 0 : nbody-1 ){
bi = bodies[i];
vx = bi.vx; vy = bi.vy; vz = bi.vz; bim = bi.mass;
e = e + (0.500 * bim * (vx*vx + vy*vy + vz*vz));
for( j = i+1 : nbody-1 ){
bj = bodies[j];
dx = bi.x-bj.x; dy = bi.y-bj.y; dz = bi.z-bj.z;
distance = (dx*dx + dy*dy + dz*dz)**0.500;
e = e - ((bim * bj.mass) / distance);
}
}
return e
}
routine offsetMomentum( b : list<OneBody>, nbody : int )
{
px = 0.00; py = 0.00; pz = 0.00
for( i = 0 : nbody-1 ){
bi = b[i];
bim = bi.mass;
px = px + (bi.vx * bim);
py = py + (bi.vy * bim);
pz = pz + (bi.vz * bim);
}
b[0].vx = -px / SOLAR_MASS
b[0].vy = -py / SOLAR_MASS
b[0].vz = -pz / SOLAR_MASS
}
global bodies = { Sun, Jupiter, Saturn, Uranus, Neptune };
routine main( N=1000 )
{
#{
N = 1;
#}
nbody = bodies.size();
offsetMomentum( bodies, nbody );
stdio.printf( "%0.9f\n", energy(bodies, nbody) );
for( i = 1 : N ) advance(bodies, nbody, 0.0100);
stdio.printf( "%0.9f\n", energy(bodies, nbody) );
}
# http://shootout.alioth.debian.org/
# contributed by Isaac Gouy, tuned by Mike Pall
# Ported to Dao by Limin Fu
global const PI = 3.1415926535897932300;
global const SOLAR_MASS = 4 * PI * PI;
global const DAYS_PER_YEAR = 365.2400;
final class OneBody
{
my x, y, z = 0.00;
my vx, vy, vz = 0.00;
my mass = 0.00;
routine Print(){ stdio.println( x, y, z, vx, vy, vz ); }
}
Jupiter = OneBody
{
x => 4.84143144246472090E+00 ;
y => -1.16032004402742839E+00 ;
z => -1.03622044471123109E-01 ;
vx => 1.66007664274403694E-03 * DAYS_PER_YEAR ;
vy => 7.69901118419740425E-03 * DAYS_PER_YEAR ;
vz => -6.90460016972063023E-05 * DAYS_PER_YEAR ;
mass => 9.54791938424326609E-04 * SOLAR_MASS
}
Saturn = OneBody
{
x => 8.34336671824457987E+00 ;
y => 4.12479856412430479E+00 ;
z => -4.03523417114321381E-01 ;
vx => -2.76742510726862411E-03 * DAYS_PER_YEAR ;
vy => 4.99852801234917238E-03 * DAYS_PER_YEAR ;
vz => 2.30417297573763929E-05 * DAYS_PER_YEAR ;
mass => 2.85885980666130812E-04 * SOLAR_MASS
}
Uranus = OneBody{
x => 1.28943695621391310E+01 ;
y => -1.51111514016986312E+01 ;
z => -2.23307578892655734E-01 ;
vx => 2.96460137564761618E-03 * DAYS_PER_YEAR ;
vy => 2.37847173959480950E-03 * DAYS_PER_YEAR ;
vz => -2.96589568540237556E-05 * DAYS_PER_YEAR ;
mass => 4.36624404335156298E-05 * SOLAR_MASS
}
Neptune = OneBody
{
x => 1.53796971148509165E+01 ;
y => -2.59193146099879641E+01 ;
z => 1.79258772950371181E-01 ;
vx => 2.68067772490389322E-03 * DAYS_PER_YEAR ;
vy => 1.62824170038242295E-03 * DAYS_PER_YEAR ;
vz => -9.51592254519715870E-05 * DAYS_PER_YEAR ;
mass => 5.15138902046611451E-05 * SOLAR_MASS
}
Sun = OneBody{ x => 0.00, y => 0.00, z => 0.00,
vx => 0.00, vy => 0.00, vz => 0.00, mass => SOLAR_MASS }
routine advance( bodies : list<OneBody>, nbody : int, dt : double )
{
for( i = 0 : nbody-1 ){
bi = bodies[i];
bix = bi.x; biy = bi.y; biz = bi.z; bimass = bi.mass;
bivx = bi.vx; bivy = bi.vy; bivz = bi.vz;
for( j = i+1 : nbody-1 ){
bj = bodies[j];
dx = bix-bj.x;
dy = biy-bj.y;
dz = biz-bj.z;
distance = (dx*dx + dy*dy + dz*dz)**0.500;
mag = dt / (distance * distance * distance);
bim = bimass * mag;
bjm = bj.mass * mag;
bivx = bivx - (dx * bjm);
bivy = bivy - (dy * bjm);
bivz = bivz - (dz * bjm);
bj.vx = bj.vx + (dx * bim);
bj.vy = bj.vy + (dy * bim);
bj.vz = bj.vz + (dz * bim);
#{
#}
}
bi.vx = bivx;
bi.vy = bivy;
bi.vz = bivz;
}
for( i = 0 : nbody-1 ){
bi = bodies[i];
bi.x = bi.x + (dt * bi.vx);
bi.y = bi.y + (dt * bi.vy);
bi.z = bi.z + (dt * bi.vz);
}
}
routine energy( bodies : list<OneBody> , nbody : int )
{
e = 0.00;
for( i = 0 : nbody-1 ){
bi = bodies[i];
vx = bi.vx; vy = bi.vy; vz = bi.vz; bim = bi.mass;
e = e + (0.500 * bim * (vx*vx + vy*vy + vz*vz));
for( j = i+1 : nbody-1 ){
bj = bodies[j];
dx = bi.x-bj.x; dy = bi.y-bj.y; dz = bi.z-bj.z;
distance = (dx*dx + dy*dy + dz*dz)**0.500;
e = e - ((bim * bj.mass) / distance);
}
}
return e
}
routine offsetMomentum( b : list<OneBody>, nbody : int )
{
px = 0.00; py = 0.00; pz = 0.00
for( i = 0 : nbody-1 ){
bi = b[i];
bim = bi.mass;
px = px + (bi.vx * bim);
py = py + (bi.vy * bim);
pz = pz + (bi.vz * bim);
}
b[0].vx = -px / SOLAR_MASS
b[0].vy = -py / SOLAR_MASS
b[0].vz = -pz / SOLAR_MASS
}
global bodies = { Sun, Jupiter, Saturn, Uranus, Neptune };
routine main( N=1000 )
{
#{
N = 1;
#}
nbody = bodies.size();
offsetMomentum( bodies, nbody );
stdio.printf( "%0.9f\n", energy(bodies, nbody) );
for( i = 1 : N ) advance(bodies, nbody, 0.0100);
stdio.printf( "%0.9f\n", energy(bodies, nbody) );
}
2 Using Tuple Type | Top |
# The Great Computer Language Shootout
# http://shootout.alioth.debian.org/
# contributed by Isaac Gouy, tuned by Mike Pall
# Ported to Dao by Limin Fu
use math
global const PI = 3.1415926535897932300;
global const SOLAR_MASS = 4 * PI * PI;
global const DAYS_PER_YEAR = 365.2400;
typedef tuple<x:double,y:double,z:double,vx:double,vy:double,vz:double,mass:double> OneBody;
# define by field enumeration:
Jupiter = OneBody
{
x => 4.84143144246472090E+00 ;
y => -1.16032004402742839E+00 ;
z => -1.03622044471123109E-01 ;
vx => 1.66007664274403694E-03 * DAYS_PER_YEAR ;
vy => 7.69901118419740425E-03 * DAYS_PER_YEAR ;
vz => -6.90460016972063023E-05 * DAYS_PER_YEAR ;
mass => 9.54791938424326609E-04 * SOLAR_MASS
}
#{
It can also be defined as a tuple of the same type as OneBody,
namely, a tuple with the same names and item types as OneBody.
Jupiter =
(
x => 4.84143144246472090E+00 ,
y => -1.16032004402742839E+00 ,
z => -1.03622044471123109E-01 ,
vx => 1.66007664274403694E-03 * DAYS_PER_YEAR ,
vy => 7.69901118419740425E-03 * DAYS_PER_YEAR ,
vz => -6.90460016972063023E-05 * DAYS_PER_YEAR ,
mass => 9.54791938424326609E-04 * SOLAR_MASS
)
#}
#{
It can also be created from a tuple without item names:
Jupiter = (
4.84143144246472090E+00 ,
-1.16032004402742839E+00 ,
-1.03622044471123109E-01 ,
1.66007664274403694E-03 * DAYS_PER_YEAR ,
7.69901118419740425E-03 * DAYS_PER_YEAR ,
-6.90460016972063023E-05 * DAYS_PER_YEAR ,
mass => 9.54791938424326609E-04 * SOLAR_MASS )
#}
Saturn = OneBody
{
x => 8.34336671824457987E+00 ;
y => 4.12479856412430479E+00 ;
z => -4.03523417114321381E-01 ;
vx => -2.76742510726862411E-03 * DAYS_PER_YEAR ;
vy => 4.99852801234917238E-03 * DAYS_PER_YEAR ;
vz => 2.30417297573763929E-05 * DAYS_PER_YEAR ;
mass => 2.85885980666130812E-04 * SOLAR_MASS
}
Uranus = OneBody{
x => 1.28943695621391310E+01 ;
y => -1.51111514016986312E+01 ;
z => -2.23307578892655734E-01 ;
vx => 2.96460137564761618E-03 * DAYS_PER_YEAR ;
vy => 2.37847173959480950E-03 * DAYS_PER_YEAR ;
vz => -2.96589568540237556E-05 * DAYS_PER_YEAR ;
mass => 4.36624404335156298E-05 * SOLAR_MASS
}
Neptune = OneBody
{
x => 1.53796971148509165E+01 ;
y => -2.59193146099879641E+01 ;
z => 1.79258772950371181E-01 ;
vx => 2.68067772490389322E-03 * DAYS_PER_YEAR ;
vy => 1.62824170038242295E-03 * DAYS_PER_YEAR ;
vz => -9.51592254519715870E-05 * DAYS_PER_YEAR ;
mass => 5.15138902046611451E-05 * SOLAR_MASS
}
Sun = OneBody{ x => 0.00, y => 0.00, z => 0.00,
vx => 0.00, vy => 0.00, vz => 0.00, mass => SOLAR_MASS }
routine advance( bodies : list<OneBody>, nbody : int, dt : double )
{
for( i = 0 : nbody-1 ){
bi = bodies[i];
bix = bi.x; biy = bi.y; biz = bi.z; bimass = bi.mass;
bivx = bi.vx; bivy = bi.vy; bivz = bi.vz;
for( j = i+1 : nbody-1 ){
bj = bodies[j];
dx = bix-bj.x;
dy = biy-bj.y;
dz = biz-bj.z;
distance = sqrt(dx*dx + dy*dy + dz*dz);
mag = dt / (distance * distance * distance);
bim = bimass * mag;
bjm = bj.mass * mag;
bivx = bivx - (dx * bjm);
bivy = bivy - (dy * bjm);
bivz = bivz - (dz * bjm);
bj.vx = bj.vx + (dx * bim);
bj.vy = bj.vy + (dy * bim);
bj.vz = bj.vz + (dz * bim);
}
bi.vx = bivx;
bi.vy = bivy;
bi.vz = bivz;
}
for( i = 0 : nbody-1 ){
bi = bodies[i];
bi.x = bi.x + (dt * bi.vx);
bi.y = bi.y + (dt * bi.vy);
bi.z = bi.z + (dt * bi.vz);
}
}
routine energy( bodies : list<OneBody> , nbody : int )
{
e = 0.00;
for( i = 0 : nbody-1 ){
bi = bodies[i];
vx = bi.vx; vy = bi.vy; vz = bi.vz; bim = bi.mass;
e = e + (0.500 * bim * (vx*vx + vy*vy + vz*vz));
for( j = i+1 : nbody-1 ){
bj = bodies[j];
dx = bi.x-bj.x; dy = bi.y-bj.y; dz = bi.z-bj.z;
distance = sqrt(dx*dx + dy*dy + dz*dz);
e = e - ((bim * bj.mass) / distance);
}
}
return e
}
routine offsetMomentum( b : list<OneBody>, nbody : int )
{
px = 0.00; py = 0.00; pz = 0.00
for( i = 0 : nbody-1 ){
bi = b[i];
bim = bi.mass;
px = px + (bi.vx * bim);
py = py + (bi.vy * bim);
pz = pz + (bi.vz * bim);
}
b[0].vx = -px / SOLAR_MASS
b[0].vy = -py / SOLAR_MASS
b[0].vz = -pz / SOLAR_MASS
}
global bodies = { Sun, Jupiter, Saturn, Uranus, Neptune };
routine main( N=1000 )
{
nbody = bodies.size();
offsetMomentum( bodies, nbody );
stdio.printf( "%0.9f\n", energy(bodies, nbody) );
for( i = 1 : N ) advance(bodies, nbody, 0.0100);
stdio.printf( "%0.9f\n", energy(bodies, nbody) );
}
# http://shootout.alioth.debian.org/
# contributed by Isaac Gouy, tuned by Mike Pall
# Ported to Dao by Limin Fu
use math
global const PI = 3.1415926535897932300;
global const SOLAR_MASS = 4 * PI * PI;
global const DAYS_PER_YEAR = 365.2400;
typedef tuple<x:double,y:double,z:double,vx:double,vy:double,vz:double,mass:double> OneBody;
# define by field enumeration:
Jupiter = OneBody
{
x => 4.84143144246472090E+00 ;
y => -1.16032004402742839E+00 ;
z => -1.03622044471123109E-01 ;
vx => 1.66007664274403694E-03 * DAYS_PER_YEAR ;
vy => 7.69901118419740425E-03 * DAYS_PER_YEAR ;
vz => -6.90460016972063023E-05 * DAYS_PER_YEAR ;
mass => 9.54791938424326609E-04 * SOLAR_MASS
}
#{
It can also be defined as a tuple of the same type as OneBody,
namely, a tuple with the same names and item types as OneBody.
Jupiter =
(
x => 4.84143144246472090E+00 ,
y => -1.16032004402742839E+00 ,
z => -1.03622044471123109E-01 ,
vx => 1.66007664274403694E-03 * DAYS_PER_YEAR ,
vy => 7.69901118419740425E-03 * DAYS_PER_YEAR ,
vz => -6.90460016972063023E-05 * DAYS_PER_YEAR ,
mass => 9.54791938424326609E-04 * SOLAR_MASS
)
#}
#{
It can also be created from a tuple without item names:
Jupiter = (
4.84143144246472090E+00 ,
-1.16032004402742839E+00 ,
-1.03622044471123109E-01 ,
1.66007664274403694E-03 * DAYS_PER_YEAR ,
7.69901118419740425E-03 * DAYS_PER_YEAR ,
-6.90460016972063023E-05 * DAYS_PER_YEAR ,
mass => 9.54791938424326609E-04 * SOLAR_MASS )
#}
Saturn = OneBody
{
x => 8.34336671824457987E+00 ;
y => 4.12479856412430479E+00 ;
z => -4.03523417114321381E-01 ;
vx => -2.76742510726862411E-03 * DAYS_PER_YEAR ;
vy => 4.99852801234917238E-03 * DAYS_PER_YEAR ;
vz => 2.30417297573763929E-05 * DAYS_PER_YEAR ;
mass => 2.85885980666130812E-04 * SOLAR_MASS
}
Uranus = OneBody{
x => 1.28943695621391310E+01 ;
y => -1.51111514016986312E+01 ;
z => -2.23307578892655734E-01 ;
vx => 2.96460137564761618E-03 * DAYS_PER_YEAR ;
vy => 2.37847173959480950E-03 * DAYS_PER_YEAR ;
vz => -2.96589568540237556E-05 * DAYS_PER_YEAR ;
mass => 4.36624404335156298E-05 * SOLAR_MASS
}
Neptune = OneBody
{
x => 1.53796971148509165E+01 ;
y => -2.59193146099879641E+01 ;
z => 1.79258772950371181E-01 ;
vx => 2.68067772490389322E-03 * DAYS_PER_YEAR ;
vy => 1.62824170038242295E-03 * DAYS_PER_YEAR ;
vz => -9.51592254519715870E-05 * DAYS_PER_YEAR ;
mass => 5.15138902046611451E-05 * SOLAR_MASS
}
Sun = OneBody{ x => 0.00, y => 0.00, z => 0.00,
vx => 0.00, vy => 0.00, vz => 0.00, mass => SOLAR_MASS }
routine advance( bodies : list<OneBody>, nbody : int, dt : double )
{
for( i = 0 : nbody-1 ){
bi = bodies[i];
bix = bi.x; biy = bi.y; biz = bi.z; bimass = bi.mass;
bivx = bi.vx; bivy = bi.vy; bivz = bi.vz;
for( j = i+1 : nbody-1 ){
bj = bodies[j];
dx = bix-bj.x;
dy = biy-bj.y;
dz = biz-bj.z;
distance = sqrt(dx*dx + dy*dy + dz*dz);
mag = dt / (distance * distance * distance);
bim = bimass * mag;
bjm = bj.mass * mag;
bivx = bivx - (dx * bjm);
bivy = bivy - (dy * bjm);
bivz = bivz - (dz * bjm);
bj.vx = bj.vx + (dx * bim);
bj.vy = bj.vy + (dy * bim);
bj.vz = bj.vz + (dz * bim);
}
bi.vx = bivx;
bi.vy = bivy;
bi.vz = bivz;
}
for( i = 0 : nbody-1 ){
bi = bodies[i];
bi.x = bi.x + (dt * bi.vx);
bi.y = bi.y + (dt * bi.vy);
bi.z = bi.z + (dt * bi.vz);
}
}
routine energy( bodies : list<OneBody> , nbody : int )
{
e = 0.00;
for( i = 0 : nbody-1 ){
bi = bodies[i];
vx = bi.vx; vy = bi.vy; vz = bi.vz; bim = bi.mass;
e = e + (0.500 * bim * (vx*vx + vy*vy + vz*vz));
for( j = i+1 : nbody-1 ){
bj = bodies[j];
dx = bi.x-bj.x; dy = bi.y-bj.y; dz = bi.z-bj.z;
distance = sqrt(dx*dx + dy*dy + dz*dz);
e = e - ((bim * bj.mass) / distance);
}
}
return e
}
routine offsetMomentum( b : list<OneBody>, nbody : int )
{
px = 0.00; py = 0.00; pz = 0.00
for( i = 0 : nbody-1 ){
bi = b[i];
bim = bi.mass;
px = px + (bi.vx * bim);
py = py + (bi.vy * bim);
pz = pz + (bi.vz * bim);
}
b[0].vx = -px / SOLAR_MASS
b[0].vy = -py / SOLAR_MASS
b[0].vz = -pz / SOLAR_MASS
}
global bodies = { Sun, Jupiter, Saturn, Uranus, Neptune };
routine main( N=1000 )
{
nbody = bodies.size();
offsetMomentum( bodies, nbody );
stdio.printf( "%0.9f\n", energy(bodies, nbody) );
for( i = 1 : N ) advance(bodies, nbody, 0.0100);
stdio.printf( "%0.9f\n", energy(bodies, nbody) );
}
view count 494 times
created at 2009-02-20, 16:09 GMT
modified at 2009-02-20, 16:20 GMT
created at 2009-02-20, 16:09 GMT
modified at 2009-02-20, 16:20 GMT
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | 29 | 30 | 31 |
fu: Many thanks (Jul.04,04:29)
klabim: fixed Hi, great, now my test works now :- ). (Jun.30,17:51)
Nightwalker: Few suggestions (Jul.03,14:37)