N-Body

#!/usr/bin/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, bodies=SYSTEM, pairs=PAIRS):

    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():
    offset_momentum(BODIES['sun'])
    for t in range(10):
	report_energy()
	for n in range(100000):
	    advance(0.01)

    report_energy()

if __name__ == '__main__':
    main()

#!/usr/bin/perl -w

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

for (1..10)
{
  printf ("%.9f\n", energy());

  # This does not, in fact, consume N*4 bytes of memory
  for (1..100000){
    advance(0.01);
  }

}

printf ("%.9f\n", energy());

#!/usr/bin/env gbs3

Class Body

  Static Private SOLAR_MASS As Float = 4 * Pi * Pi
  Private Const DAYS_PER_YEAR As Float = 365.24

  Public X As Float
  Public Y As Float
  Public Z As Float
  Public VX As Float
  Public VY As Float
  Public VZ As Float
  Public Mass As Float

  Static Public Sub Jupiter() As Body

    Dim P As New Body
    p.x = 4.84143144246472090e+00
    p.y = -1.16032004402742839e+00
    p.z = -1.03622044471123109e-01
    p.vx = 1.66007664274403694e-03 * DAYS_PER_YEAR
    p.vy = 7.69901118419740425e-03 * DAYS_PER_YEAR
    p.vz = -6.90460016972063023e-05 * DAYS_PER_YEAR
    p.mass = 9.54791938424326609e-04 * SOLAR_MASS
    return p

  End

  Static Public Sub Saturn() As Body

    Dim P As New Body
    p.x = 8.34336671824457987e+00
    p.y = 4.12479856412430479e+00
    p.z = -4.03523417114321381e-01
    p.vx = -2.76742510726862411e-03 * DAYS_PER_YEAR
    p.vy = 4.99852801234917238e-03 * DAYS_PER_YEAR
    p.vz = 2.30417297573763929e-05 * DAYS_PER_YEAR
    p.mass = 2.85885980666130812e-04 * SOLAR_MASS
    return p

  End

  Static Public Sub Uranus() As Body

    Dim P As New Body

    p.x = 1.28943695621391310e+01
    p.y = -1.51111514016986312e+01
    p.z = -2.23307578892655734e-01
    p.vx = 2.96460137564761618e-03 * DAYS_PER_YEAR
    p.vy = 2.37847173959480950e-03 * DAYS_PER_YEAR
    p.vz = -2.96589568540237556e-05 * DAYS_PER_YEAR
    p.mass = 4.36624404335156298e-05 * SOLAR_MASS
    return p

  End

  Static Public Sub Neptune() As Body

    Dim P As New Body

    p.x = 1.53796971148509165e+01
    p.y = -2.59193146099879641e+01
    p.z = 1.79258772950371181e-01
    p.vx = 2.68067772490389322e-03 * DAYS_PER_YEAR
    p.vy = 1.62824170038242295e-03 * DAYS_PER_YEAR
    p.vz = -9.51592254519715870e-05 * DAYS_PER_YEAR
    p.mass = 5.15138902046611451e-05 * SOLAR_MASS
    return p

  End

  Static Public Sub Sun() As Body

    Dim P As New Body

    p.mass = SOLAR_MASS
    return p

  End

  Public Sub OffsetMomentum(px As Float, py As Float, pz As Float) As Body

    vx = -px / SOLAR_MASS
    vy = -py / SOLAR_MASS
    vz = -pz / SOLAR_MASS

    Return Me

  End

End Class


Class NBodySystem

  Private Bodies As Body[]

  Public Sub _new()

    Dim PX, PY, PZ As Float
    Dim I As Integer

    Bodies = [ Body.Sun(), Body.Jupiter(), Body.Saturn(), Body.Uranus(), Body.Neptune() ]

    For I = 0 To Bodies.Max
      PX += Bodies[I].vx * Bodies[I].mass
      PY += Bodies[I].vy * Bodies[I].mass
      PZ += Bodies[I].vz * Bodies[I].mass
    Next

    Bodies[0].offsetMomentum(PX, PY, PZ)

  End

  Public Sub Advance(dt As Float)

    Dim I, J As Integer
    Dim iBody, jBody As Body
    Dim dx, dy, dz As Float
    Dim dSquared, fMag As Float
    Dim iMass, jMass, iMag, jMag As Float

    For I = 0 To Bodies.Max
      iBody = Bodies[I]
      iMass = iBody.mass

      For J = I + 1 To Bodies.Max
        jBody = Bodies[J]
	jMass = jBody.mass

        dx = iBody.x - jBody.x
        dy = iBody.y - jBody.y
        dz = iBody.z - jBody.z

        dSquared = dx * dx + dy * dy + dz * dz
        fMag = dt / (dSquared * Sqr(dSquared))
	iMag = iMass * fMag
	jMag = jMass * fMag

        iBody.vx -= dx * jMag
	iBody.vy -= dy * jMag
	iBody.vz -= dz * jMag

	jBody.vx += dx * iMag
	jBody.vy += dy * iMag
	jBody.vz += dz * iMag
      Next
    Next

    For Each iBody in Bodies
      iBody.x += dt * iBody.vx
      iBody.y += dt * iBody.vy
      iBody.z += dt * iBody.vz
    Next

  End

  Public Sub Energy() As Float

    Dim dx, dy, dz, distance, E As Float
    Dim iBody, jBody As Body
    Dim I, J As Integer

    For I = 0 To Bodies.Max

      iBody = bodies[i]
      E += 0.5 * iBody.mass * (iBody.vx * iBody.vx + iBody.vy * iBody.vy + iBody.vz * iBody.vz)

      For J = I + 1 To Bodies.Max

	jBody = Bodies[J]
        dx = iBody.x - jBody.x
        dy = iBody.y - jBody.y
        dz = iBody.z - jBody.z

        distance = Sqr(dx*dx + dy*dy + dz*dz)
        E -= (iBody.mass * jBody.mass) / distance

      Next
    Next

    return E

  End

End Class


Dim S As New NBodySystem
Dim I, N As Integer

For N = 1 To 10

  Print S.Energy()

  For I = 1 To 100000
    S.Advance(0.01)
  Next

Next

Print S.Energy()