nbody.p6

#
# The Great Computer Language Shootout
# http://shootout.alioth.debian.org/
#
# contributed by Christoph Bauer
# converted into Perl by Márton Papp
#

my num $pi = 3.141592653589793e0;
my num $solar_mass =(4 * $pi * $pi);
my num $days_per_year =365.24e0;

class Body is rw {
    has num $.mass;
    has num ($.x, $.y, $.z);
    has num ($.vx, $.vy, $.vz);
}

sub advance (@bodies, num $dt)
{
    my int $i;
    my int $j;

    my int $nbodies = @bodies.elems;

    loop ($i = 0; $i < $nbodies; $i++) {
        my $b := @bodies[$i];
        my $b2;
        my num $dx;
        my num $dy;
        my num $dz;
        my num $distance;
        my num $mag;
        loop ($j = $i + 1; $j < $nbodies; $j++) {
            $b2 := @bodies[$j];
            $dx = $b.x - $b2.x;
            $dy = $b.y - $b2.y;
            $dz = $b.z - $b2.z;

            $distance = sqrt($dx * $dx + $dy * $dy + $dz * $dz);
            $mag = $dt / ($distance * $distance * $distance);

            my num $bmm = $b.mass * $mag;
            my num $b2mm = $b2.mass * $mag;

            $b.vx = $b.vx - $dx * $b2mm;
            $b.vy = $b.vy - $dy * $b2mm;
            $b.vz = $b.vz - $dz * $b2mm;
            $b2.vx = $b2.vx + $dx * $bmm;
            $b2.vy = $b2.vy + $dy * $bmm;
            $b2.vz = $b2.vz + $dz * $bmm;
        }
    }

    loop ($i = 0; $i < $nbodies; $i++) {
        my $b := @bodies[$i];
        $b.x = $b.x + $dt * $b.vx;
        $b.y = $b.y + $dt * $b.vy;
        $b.z = $b.z + $dt * $b.vz;
    }
}

sub energy(@bodies)
{
    my num $e;
    my int $i;
    my num $j;
    my int $nbodies = @bodies.elems;
    $e = 0e0;
    loop ($i = 0; $i < $nbodies; $i++) {
        my $b = @bodies[$i];
        $e += 0.5 * $b.mass * ($b.vx * $b.vx + $b.vy * $b.vy + $b.vz * $b.vz);
        my ($dx,$dy,$dz,$distance);
        loop ($j = $i + 1; $j < $nbodies; $j++) {
            my $b2 = @bodies[$j];
            $dx = $b.x - $b2.x;
            $dy = $b.y - $b2.y;
            $dz = $b.z - $b2.z;
            $distance = sqrt($dx * $dx + $dy * $dy + $dz * $dz);
            $e -= ($b.mass * $b2.mass) / $distance;
        }
    }
    return $e;
}

sub offset_momentum(@bodies)
{
    my num $px = 0e0;
    my num $py = 0e0;
    my num $pz = 0e0;
    my int $i;
    my int $nbodies = @bodies.elems;
    loop ($i = 0; $i < $nbodies; $i++) {
        my num $bim = @bodies[$i].mass;
        $px += @bodies[$i].vx * $bim;
        $py += @bodies[$i].vy * $bim;
        $pz += @bodies[$i].vz * $bim;
    }
    @bodies[0].vx = - $px / $solar_mass;
    @bodies[0].vy = - $py / $solar_mass;
    @bodies[0].vz = - $pz / $solar_mass;
}


my @bodies = (
    # sun
    Body.new( x => 0e0, y => 0e0,z => 0e0,vx => 0e0, vy => 0e0,vz => 0e0, mass=> $solar_mass ),
    # jupiter
    Body.new(
        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
    Body.new(
        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
    Body.new(
        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
    Body.new(
        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
    )
);


my int $n = (@*ARGS[0] // 10_000).Int;

offset_momentum(@bodies);

printf("%.9f\n", energy(@bodies));
for ^$n
{
    advance(@bodies, 0.01e0);
}
#printf("%.9f\n", energy(@bodies));
say energy(@bodies);