NnBody

N-Body

See http://openeuphoria.org/forum/119406.wc

N-Body Euphoria

-- 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())  
 

N-Body Euphoria (Alternate)

-- The Computer Language Shootout Benchmarks 
-- http://benchmarksgame.alioth.debian.org 
--  
--  Converted to Euphoria by Jason Gade 
--  Version 6: vectorized position and velocity operations 
--  run: exu nbody.ex N 
 
-- with profile 
without warning 
without type_check 
 
include get.e 
 
constant    PI = 3.141592653589793, 
            SOLAR_MASS = 4 * PI * PI, 
            DAYS_PER_YEAR = 365.24 
 
-- point fields 
constant    X = 1, 
            Y = 2, 
            Z = 3 
 
-- body fields 
constant    POS = 1, 
            VELOCITY = 2, 
            MASS = 3 
 
sequence sun = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, SOLAR_MASS} 
sequence 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} 
 
sequence bodies = {sun, jupiter, saturn, uranus, neptune} 
 
procedure offset_momentum() 
    for i = 1 to length(bodies) do 
        bodies[1][VELOCITY] -= bodies[i][VELOCITY] * bodies[i][MASS] / SOLAR_MASS 
    end for 
end procedure 
 
function energy() 
    atom e = 0 
    sequence dx, vx 
 
    for i = 1 to length(bodies) do 
        vx = bodies[i][VELOCITY] * bodies[i][VELOCITY] 
        e += 0.5 * bodies[i][MASS] * (vx[X] + vx[Y] + vx[Z]) 
        for j = i+1 to length(bodies) do 
            dx = bodies[i][POS] - bodies[j][POS] 
            dx *= dx 
            e -= bodies[i][MASS] * bodies[j][MASS] / sqrt(dx[X] + dx[Y] + dx[Z]) 
        end for 
    end for 
 
    return e 
 
end function 
 
procedure advance(atom dt) 
    sequence dx, dx2 
    atom dx2sum, distance, mag 
 
    for i = 1 to length(bodies) do 
        for j = i+1 to length(bodies) do 
            dx = bodies[i][POS] - bodies[j][POS] 
            dx2 = dx * dx 
            dx2sum = dx2[X] + dx2[Y] + dx2[Z] 
            distance = sqrt(dx2sum) 
            mag = dt / (distance * dx2sum) 
 
            bodies[i][VELOCITY] -= dx * bodies[j][MASS] * mag 
            bodies[j][VELOCITY] += dx * bodies[i][MASS] * mag             
        end for 
        bodies[i][POS] += dt * bodies[i][VELOCITY] 
    end for 
             
end procedure 
 
procedure main(sequence argv) 
    object n  
     
    if length(argv) > 2 then 
        n = value(argv[3]) 
    else 
        n = {GET_FAIL, 0} 
    end if 
     
    if n[1] = GET_SUCCESS then 
            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(command_line()) 
 
 

N-Body Python

# The Computer Language Benchmarks Game 
# http://shootout.alioth.debian.org/ 
# 
# originally by Kevin Carson 
# modified by Tupteq, Fredrik Johansson, and Daniel Nanz 
# modified by Maciej Fijalkowski 
# 2to3 
 
import sys 
 
def combinations(l): 
    result = [] 
    for x in range(len(l) - 1): 
        ls = l[x+1:] 
        for y in ls: 
            result.append((l[x],y)) 
    return result 
 
PI = 3.14159265358979323 
SOLAR_MASS = 4 * PI * PI 
DAYS_PER_YEAR = 365.24 
 
BODIES = { 
    'sun': ([0.0, 0.0, 0.0], [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) } 
 
 
SYSTEM = list(BODIES.values()) 
PAIRS = combinations(SYSTEM) 
 
 
def advance(dt, n, bodies=SYSTEM, pairs=PAIRS): 
 
    for i in range(n): 
        for (([x1, y1, z1], v1, m1), 
             ([x2, y2, z2], v2, m2)) in pairs: 
            dx = x1 - x2 
            dy = y1 - y2 
            dz = z1 - z2 
            mag = dt * ((dx * dx + dy * dy + dz * dz) ** (-1.5)) 
            b1m = m1 * mag 
            b2m = m2 * mag 
            v1[0] -= dx * b2m 
            v1[1] -= dy * b2m 
            v1[2] -= dz * b2m 
            v2[0] += dx * b1m 
            v2[1] += dy * b1m 
            v2[2] += dz * b1m 
        for (r, [vx, vy, vz], m) in bodies: 
            r[0] += dt * vx 
            r[1] += dt * vy 
            r[2] += dt * vz 
 
 
def report_energy(bodies=SYSTEM, pairs=PAIRS, e=0.0): 
 
    for (((x1, y1, z1), v1, m1), 
         ((x2, y2, z2), v2, m2)) in pairs: 
        dx = x1 - x2 
        dy = y1 - y2 
        dz = z1 - z2 
        e -= (m1 * m2) / ((dx * dx + dy * dy + dz * dz) ** 0.5) 
    for (r, [vx, vy, vz], m) in bodies: 
        e += m * (vx * vx + vy * vy + vz * vz) / 2. 
    print("%.9f" % e) 
 
def offset_momentum(ref, bodies=SYSTEM, px=0.0, py=0.0, pz=0.0): 
 
    for (r, [vx, vy, vz], m) in bodies: 
        px -= vx * m 
        py -= vy * m 
        pz -= vz * m 
    (r, v, m) = ref 
    v[0] = px / m 
    v[1] = py / m 
    v[2] = pz / m 
 
def main(n, ref='sun'): 
    offset_momentum(BODIES[ref]) 
    report_energy() 
    advance(0.01, n) 
    report_energy() 
 
if __name__ == '__main__': 
    main(int(sys.argv[1])) 
 
 
	 

N-Body Perl

# The Computer Language Shootout 
# http://shootout.alioth.debian.org/ 
# 
# contributed by Christoph Bauer 
# converted into Perl by M?rton Papp 
# fixed and cleaned up by Danny Sauer 
# optimized by Jesse Millikan 
 
use constant PI            => 3.141592653589793; 
use constant SOLAR_MASS    => (4 * PI * PI); 
use constant DAYS_PER_YEAR => 365.24; 
 
#  Globals for arrays... Oh well. 
#  Almost every iteration is a range, so I keep the last index rather than a count. 
my (@xs, @ys, @zs, @vxs, @vys, @vzs, @mass, $last); 
 
sub advance($) 
{ 
  my ($dt) = @_; 
  my ($mm, $mm2, $j, $dx, $dy, $dz, $distance, $mag); 
 
#  This is faster in the outer loop... 
  for (0..$last) { 
#  But not in the inner loop. Strange. 
    for ($j = $_ + 1; $j < $last + 1; $j++) { 
      $dx = $xs[$_] - $xs[$j]; 
      $dy = $ys[$_] - $ys[$j]; 
      $dz = $zs[$_] - $zs[$j]; 
      $distance = sqrt($dx * $dx + $dy * $dy + $dz * $dz); 
      $mag = $dt / ($distance * $distance * $distance); 
      $mm = $mass[$_] * $mag; 
      $mm2 = $mass[$j] * $mag; 
      $vxs[$_] -= $dx * $mm2; 
      $vxs[$j] += $dx * $mm; 
      $vys[$_] -= $dy * $mm2; 
      $vys[$j] += $dy * $mm; 
      $vzs[$_] -= $dz * $mm2; 
      $vzs[$j] += $dz * $mm; 
    } 
 
# We're done with planet $_ at this point 
# This could be done in a seperate loop, but it's slower 
    $xs[$_] += $dt * $vxs[$_]; 
    $ys[$_] += $dt * $vys[$_]; 
    $zs[$_] += $dt * $vzs[$_]; 
  } 
} 
 
sub energy 
{ 
  my ($e, $i, $dx, $dy, $dz, $distance); 
 
  $e = 0.0; 
  for $i (0..$last) { 
    $e += 0.5 * $mass[$i] * 
          ($vxs[$i] * $vxs[$i] + $vys[$i] * $vys[$i] + $vzs[$i] * $vzs[$i]); 
    for ($i + 1..$last) { 
      $dx = $xs[$i] - $xs[$_]; 
      $dy = $ys[$i] - $ys[$_]; 
      $dz = $zs[$i] - $zs[$_]; 
      $distance = sqrt($dx * $dx + $dy * $dy + $dz * $dz); 
      $e -= ($mass[$i] * $mass[$_]) / $distance; 
    } 
  } 
  return $e; 
} 
 
sub offset_momentum 
{ 
  my ($px, $py, $pz) = (0.0, 0.0, 0.0); 
 
  for (0..$last) { 
    $px += $vxs[$_] * $mass[$_]; 
    $py += $vys[$_] * $mass[$_]; 
    $pz += $vzs[$_] * $mass[$_]; 
  } 
  $vxs[0] = - $px / SOLAR_MASS; 
  $vys[0] = - $py / SOLAR_MASS; 
  $vzs[0] = - $pz / SOLAR_MASS; 
} 
 
# @ns = ( sun, jupiter, saturn, uranus, neptune ) 
@xs = (0, 4.84143144246472090e+00, 8.34336671824457987e+00, 1.28943695621391310e+01, 1.53796971148509165e+01); 
@ys = (0, -1.16032004402742839e+00, 4.12479856412430479e+00, -1.51111514016986312e+01, -2.59193146099879641e+01); 
@zs = (0, -1.03622044471123109e-01, -4.03523417114321381e-01, -2.23307578892655734e-01, 1.79258772950371181e-01); 
@vxs = map {$_ * DAYS_PER_YEAR} 
  (0, 1.66007664274403694e-03, -2.76742510726862411e-03, 2.96460137564761618e-03, 2.68067772490389322e-03); 
@vys = map {$_ * DAYS_PER_YEAR} 
  (0, 7.69901118419740425e-03, 4.99852801234917238e-03, 2.37847173959480950e-03, 1.62824170038242295e-03); 
@vzs = map {$_ * DAYS_PER_YEAR} 
  (0, -6.90460016972063023e-05, 2.30417297573763929e-05, -2.96589568540237556e-05, -9.51592254519715870e-05); 
@mass = map {$_ * SOLAR_MASS} 
  (1, 9.54791938424326609e-04, 2.85885980666130812e-04, 4.36624404335156298e-05, 5.15138902046611451e-05); 
 
$last = @xs - 1; 
 
offset_momentum(); 
printf ("%.9f\n", energy()); 
 
my $n = $ARGV[0]; 
 
# This does not, in fact, consume N*4 bytes of memory 
for (1..$n){ 
  advance(0.01); 
} 
 
printf ("%.9f\n", energy()); 
 
 
	 

Search



Quick Links

User menu

Not signed in.

Misc Menu