nbody2.jl

const solMass = 4 * pi^2
const daysPYear = 365.24

immutable ThreeD
   x::Float64
   y::Float64
   z::Float64
end

# all these functions modify a global array, a better idea is to pass them the array
function offsetMomentum!(velocities::Array{ThreeD,1},masses::Array{Float64,1})
   px = 0.0
   py = 0.0
   pz = 0.0

   for i = 1:length(velocities)
      px += velocities[i].x * masses[i]
      py += velocities[i].y * masses[i]
      pz += velocities[i].z * masses[i]
   end

   velocities[1] = ThreeD(-px/solMass,-py/solMass,-pz/solMass);
   return;
end

function energy(positions::Array{ThreeD,1},velocities::Array{ThreeD,1},masses::Array{Float64,1})
   egy = 0.0
   for i = 1:length(masses)
      egy += (0.5 * masses[i] * (velocities[i].x^2 + velocities[i].y^2 + velocities[i].z^2))
      for j = (i + 1):length(masses)
         dist = sqrt((positions[i].x - positions[j].x)^2 + (positions[i].y - positions[j].y)^2 +(positions[i].z - positions[j].z)^2)
         egy -= ((masses[i] * masses[j]) / dist)
      end
   end

   return egy
end

function advance!(dt::Float64,positions::Array{ThreeD,1},velocities::Array{ThreeD,1},masses::Array{Float64,1})
   for i = 1:(length(positions) - 1)
      for j = (i+1):length(positions)
         @inbounds dx = positions[i].x - positions[j].x
         @inbounds dy = positions[i].y - positions[j].y
         @inbounds dz = positions[i].z - positions[j].z

         dist = sqrt(dx^2 + dy^2 + dz^2)
         mag = dt / dist^3
      
         @inbounds velocities[i] = ThreeD(velocities[i].x - dx*masses[j]*mag,velocities[i].y - dy*masses[j]*mag,velocities[i].z - dz*masses[j]*mag);
         @inbounds velocities[j] = ThreeD(velocities[j].x + dx*masses[i]*mag,velocities[j].y + dy*masses[i]*mag,velocities[j].z + dz*masses[i]*mag);
      end
   end

   for i = 1:length(positions)
      @inbounds positions[i] = ThreeD(positions[i].x + dt*velocities[i].x,positions[i].y + dt*velocities[i].y,positions[i].z + dt*velocities[i].z)
   end
end

function main()
   masses = Float64[
      solMass,
      9.54791938424326609e-04 * solMass,
      2.85885980666130812e-04 * solMass,
      4.36624404335156298e-05 * solMass,
      5.15138902046611451e-05 * solMass,
   ]

   positions = ThreeD[
      ThreeD(0,0,0),
      ThreeD(4.84143144246472090e+00, -1.16032004402742839e+00, -1.03622044471123109e-01),
      ThreeD(8.34336671824457987e+00,  4.12479856412430479e+00, -4.03523417114321381e-01),
      ThreeD(1.28943695621391310e+01, -1.51111514016986312e+01, -2.23307578892655734e-01),
      ThreeD(1.53796971148509165e+01, -2.59193146099879641e+01,  1.79258772950371181e-01)
   ]

   velocities = ThreeD[
      ThreeD(0,0,0),
      ThreeD( 1.66007664274403694e-03 * daysPYear, 7.69901118419740425e-03 * daysPYear, -6.90460016972063023e-05 * daysPYear),
      ThreeD(-2.76742510726862411e-03 * daysPYear, 4.99852801234917238e-03 * daysPYear,  2.30417297573763929e-05 * daysPYear),
      ThreeD( 2.96460137564761618e-03 * daysPYear, 2.37847173959480950e-03 * daysPYear, -2.96589568540237556e-05 * daysPYear),
      ThreeD( 2.68067772490389322e-03 * daysPYear, 1.62824170038242295e-03 * daysPYear, -9.51592254519715870e-05 * daysPYear)
   ]

   @time offsetMomentum!(velocities,masses)
   @time en = energy(positions,velocities,masses)
   println(en)
   @time for i = 1:50_000_000
      advance!(0.01,positions,velocities,masses)
   end
   @time en = energy(positions,velocities,masses)
   println(en)
end