Turning a 2D 3 body problem into 3D not so easy?

Question

I have programmed a 3 body problem, that can in theory work for more than 3 bodies. I apologize in advance for the length of the code that I posted, I don't really know how much one would need to help.

My problem in a nutshell is that when I switched the program from 2D to 3D by adding a component to my initial vectors and making my ListPlot 3D, it didn't work any more and I couldn't see what was even happening. (I just want a way for me to see the path each of my particles takes in 3D)

On a side note... how do I create my Orbit1, Orbit2 and Orbit3 in a loop? so that I don't have to just keep adding terms when I add a particle.

Clear["Global`"];  
Off[General::spell, General::spell1, General::infy];
dt = 0.01;
bodies = 3;
rnum := RandomReal[{-1.00, 1.00}];
rnum2 := RandomReal[{-0.07, 0.07}];

PosVel = Table[0, {0}, {j, 2}];
Do[{
   Pos = {rnum, rnum, rnum},
   Vel = {rnum2, rnum2, rnum2},
   AppendTo[PosVel, {Pos, Vel}]
   }, {j, 1, bodies}];

r = Table[0, {i, bodies}, {j, bodies}];
rmag = Table[0, {i, bodies}, {j, bodies}];
rsq = Table[0, {i, bodies}, {j, bodies}];

Do[{
   r[[i, j]] = (PosVel[[i, 1]] - PosVel[[j, 1]]),
   rmag[[i, j]] = Norm[(PosVel[[i, 1]] - PosVel[[j, 1]])],
   rsq[[i, j]] = (Norm[(PosVel[[i, 1]] - PosVel[[j, 1]])])^2
   }, {j, 1, bodies}, {i, 1, bodies}]

f = -1/rmag^3;
Do[{
  f[[i, i]] = 0
  }, {i, 1, bodies}]

Force = Table[0, {i, bodies}];
Do[{
  Force[[i]] = Sum[f[[i, j]]*r[[i, j]], {j, 1, bodies}]
  }, {i, 1, bodies}]

n = 3000;
Clear[Orbit];
Orbit = Table[0, {i, n}, {j, bodies}];

Do[{
  PosVel[[i, 2]] = PosVel[[i, 2]] + Force[[i]]*dt,
  PosVel[[i, 1]] = PosVel[[i, 1]] + PosVel[[i, 2]]*dt,

  Do[{
    r[[i, j]] = (PosVel[[i, 1]] - PosVel[[j, 1]]),
    rmag[[i, j]] = Norm[(PosVel[[i, 1]] - PosVel[[j, 1]])],
    rsq[[i, j]] = (Norm[(PosVel[[i, 1]] - PosVel[[j, 1]])])^2
    }, {j, 1, bodies}, {i, 1, bodies}],
  f = -1/rmag^3,
  Do[{
    f[[i, i]] = 0
    }, {i, 1, bodies}],
  Do[{
    Force[[i]] = Sum[f[[i, j]]*r[[i, j]], {j, 1, bodies}]
    }, {i, 1, bodies}],
  Orbit[[j, i]] = PosVel[[i, 1]]
  }, {i, 1, bodies}, {j, 1, n}]

Orbit1 = Table[Orbit[[i, 1]], {i, 1, n}];
Orbit2 = Table[Orbit[[i, 2]], {i, 1, n}];
Orbit3 = Table[Orbit[[i, 3]], {i, 1, n}];

Animate[
 ListPlot3D[{Orbit1[[1 ;; j]], Orbit2[[1 ;; j]], Orbit3[[1 ;; j]]},
  Mesh -> None,
  Epilog -> {PointSize -> 0.05, Blue, 
    Point[{Orbit1[[j]], Orbit2[[j]], Orbit3[[j]]}]}], {j, 1, n, 1}]

Answer 1

in response to your comment, I wrote a small example using map:

First create some data:

bodies = 4;
createBody[size_] := Module[
  {rp := RandomReal[{-1.00 size, 1.00 size}],
   rv := RandomReal[{-0.07 size, 0.07 size}]},
  {{rp, rp, rp}, {rv, rv, rv}, 1./bodies}]

data = Map[createBody, Range[bodies]]

Output:

{{{0.353151, -0.362113, -0.177183}, {-0.033645, -0.0480111, -0.0448151}, 0.25},
 {{0.377195, 0.53147, 1.16556}, {0.00865239, -0.0489291, -0.0692621}, 0.25}, 
 {{0.123244, 2.74798, 1.83201}, {0.142131, -0.0563642, -0.185028}, 0.25}, 
 {{-1.05087, 3.08883, -2.36998}, {-0.0721776, 0.187565, -0.0197556}, 0.25}}

Now we tell mathematica where it can find positions, velocities etc. I'm not sure how much this slows the program down, but it makes the code easier to read.

positions := Table[data[[n, 1]], {n, 1, Length[data]}]
velocities := Table[data[[n, 2]], {n, 1, Length[data]}]
masses := Table[data[[n, 3]], {n, 1, Length[data]}]

Now we define the force on one body:

force[body_] := G *
  Sum[If[body != others, 
    masses[[others]]/
     Norm[positions[[others]] - 
       positions[[body]]]^3*(positions[[others]] - positions[[body]]),
     0], {others, Length[data]}]

And map this computation over all the bodies:

Map[force, Range[bodies]]

Output:

{{-0.00403823 G, 0.0791318 G, 0.0830302 G}, {-0.0101503 G, -0.00272245 G, -0.0759115 G}, 
 {0.00265276 G, -0.058108 G, -0.0355709 G}, {0.0115358 G, -0.0183013 G, 0.0284522 G}}

I'll be gone for a few hours now, but I'll gladly expand this some more later on if you like.