Pastey nbody update
- Posted by mattlewis
(admin)
Oct 05, 2012
-- The Computer Language Shootout Benchmarks
-- http://shootout.alioth.debian.org/
--
-- Converted to Euphoria by Jason Gade
-- Optimized by Matt Lewis
-- run: eui nbody.ex N
without warning
without type_check
include get.e
constant PI = 3.141592653589793,
SOLAR_MASS = 4 * PI * PI,
DAYS_PER_YEAR = 365.24
-- struct planet
constant name = 1,
x = 2,
y = 3,
z = 4,
vx = 5,
vy = 6,
vz = 7,
mass = 8
-- end struct
sequence bodies
bodies = {
{ "Sun", 0, 0, 0, 0, 0, 0, SOLAR_MASS},
{ "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,
9.54791938424326609e-04 * SOLAR_MASS },
{ "Saturn",
8.34336671824457987e+00,
4.12479856412430479e+00,
-4.03523417114321381e-01,
-2.76742510726862411e-03 * DAYS_PER_YEAR,
4.99852801234917238e-03 * DAYS_PER_YEAR,
2.30417297573763929e-05 * DAYS_PER_YEAR,
2.85885980666130812e-04 * SOLAR_MASS },
{ "Uranus",
1.28943695621391310e+01,
-1.51111514016986312e+01,
-2.23307578892655734e-01,
2.96460137564761618e-03 * DAYS_PER_YEAR,
2.37847173959480950e-03 * DAYS_PER_YEAR,
-2.96589568540237556e-05 * DAYS_PER_YEAR,
4.36624404335156298e-05 * SOLAR_MASS },
{ "Neptune",
1.53796971148509165e+01,
-2.59193146099879641e+01,
1.79258772950371181e-01,
2.68067772490389322e-03 * DAYS_PER_YEAR,
1.62824170038242295e-03 * DAYS_PER_YEAR,
-9.51592254519715870e-05 * DAYS_PER_YEAR,
5.15138902046611451e-05 * SOLAR_MASS }
}
constant NBODIES = length(bodies)
sequence NAME = {}, X = {}, Y = {}, Z = {}, VX = {}, VY = {}, VZ = {}, MASS = {}
for i = 1 to NBODIES do
for f = name to mass do
switch f do
case name then NAME = append( NAME, bodies[i][f] )
case x then X &= bodies[i][f]
case y then Y &= bodies[i][f]
case z then Z &= bodies[i][f]
case vx then VX &= bodies[i][f]
case vy then VY &= bodies[i][f]
case vz then VZ &= bodies[i][f]
case mass then MASS &= bodies[i][f]
end switch
end for
end for
procedure advance(atom dt)
atom dx, dy, dz, distance, mag, mass, distance_2, mass_mag
for i = 1 to NBODIES do
for j = i + 1 to NBODIES do
dx = X[i] - X[j]
dy = Y[i] - Y[j]
dz = Z[i] - Z[j]
distance_2 = dx*dx
distance_2 += dy*dy
distance_2 += dz*dz
distance = sqrt(distance_2)
distance *= distance_2
mag = dt / distance
mass = MASS[j]
mass_mag = mass * mag
VX[i] -= dx * mass_mag
VY[i] -= dy * mass_mag
VZ[i] -= dz * mass_mag
mass = MASS[i]
mass_mag = mass * mag
VX[j] += dx * mass_mag
VY[j] += dy * mass_mag
VZ[j] += dz * mass_mag
end for
X[i] += dt * VX[i]
Y[i] += dt * VY[i]
Z[i] += dt * VZ[i]
end for
end procedure -- advance
function energy()
atom e, dx, dy, dz, distance, mass_i
e = 0.0
for i = 1 to NBODIES do
mass_i = MASS[i]
e += 0.5 * mass_i * (VX[i]*VX[i] +
VY[i]*VY[i] +
VZ[i]*VZ[i])
for j = i + 1 to NBODIES do
dx = X[i] - X[j]
dy = Y[i] - Y[j]
dz = Z[i] - Z[j]
distance = sqrt(dx*dx + dy*dy + dz*dz)
e -= (mass_i*MASS[j])/distance
end for
end for
return e
end function -- energy
procedure offset_momentum()
atom px, py, pz
px = 0.0
py = 0.0
pz = 0.0
for i = 1 to NBODIES do
px += VX[i] * MASS[i]
py += VY[i] * MASS[i]
pz += VZ[i] * MASS[i]
end for
VX[1] = - px / SOLAR_MASS
VY[1] = - py / SOLAR_MASS
VZ[1] = - pz / SOLAR_MASS
end procedure -- offset_momentum
procedure main(sequence argv)
object n
if length(argv) > 2 then
n = value(argv[3])
n = n[2]
else
n = 1000
end if
offset_momentum()
printf(1, "%.9f\n", energy())
for i = 1 to n do
advance(0.01)
end for
printf(1, "%.9f\n", energy())
end procedure -- main
main(command_line())