File:Galaxy collision toy model.gif

mcm1 = 1;
mcm2 = 1;
soft = 0.01;
poscm = {{0, 0}, {-7, 3}};
vcm = {{0, 0}, {0.5, 0}};
acccm = {-mcm2 (poscm[[1]] - poscm[[2]])/(Norm[poscm[[1]] - poscm[[2]]]^3 + soft), -mcm1 (poscm[[2]] - poscm[[1]])/(Norm[poscm[[1]] - poscm[[2]]]^3 + soft)};
dt = 0.01;
tmax = 2 10^3;
poscmold = poscm;
poscm = poscmold + vcm dt + acccm/2 dt^2;
evocm = Reap[Do[
     acccm = {-mcm2 (poscm[[1]] - poscm[[2]])/(Norm[poscm[[1]] - poscm[[2]]]^3 + soft), -mcm1 (poscm[[2]] - poscm[[1]])/(Norm[poscm[[1]] - poscm[[2]]]^3 + soft)};
     poscmoldold = poscmold;
     poscmold = poscm;
     poscm = 2 poscmold - poscmoldold + acccm dt^2;
     Sow[poscm];
     , tmax]][[2, 1]];

posstmp = RandomPoint[Disk[], 500];
poss = Join[posstmp, Map[(evocm[[1, 2]] + #) &, posstmp]  ];
ns = Dimensions[poss][[1]];
accs = Table[-mcm1 (poss[[j]] - evocm[[1, 1]])/(Norm[evocm[[1, 1]] - poss[[j]]]^3 + soft) - mcm2 (poss[[j]] - evocm[[1, 2]])/(Norm[evocm[[1, 2]] - poss[[j]]]^3 + soft), {j, 1, ns}];
vstmp = Table[
   Sqrt[Norm[accs[[j]] ] *Norm[evocm[[1, 1]] - poss[[j]]] ]/Norm[evocm[[1, 1]] - poss[[j]]] {poss[[j, 2]], -poss[[j, 1]]}, {j,1, ns/2}];
vs = Join[vstmp, Map[(vcm[[2]] + #) &, vstmp] ];
vs = vs + RandomReal[{-0.1, 0.1}, {ns, 2}];

possold = poss;
poss = possold + vs dt + accs/2 dt^2;
evos = Reap[For[k = 1, k <= tmax, k++,
     accs = 
      Table[-mcm1 (poss[[j]] - evocm[[k, 1]])/(Norm[evocm[[k, 1]] - poss[[j]]]^3 + soft) - mcm2 (poss[[j]] - evocm[[k, 2]])/(Norm[evocm[[k, 2]] - poss[[j]]]^3 + soft), {j, 1, ns}];
     possoldold = possold;
     possold = poss;
     poss = 2 possold - possoldold + accs dt^2;
     Sow[poss];
     ]][[2, 1]];

p1 = Table[Graphics[{Black, Disk[evocm[[j, 1]], 0.2], Disk[evocm[[j, 2]], 0.2], Purple, Table[Disk[evos[[j, k]], 0.05 ], {k, 1, ns/2}], Orange, Table[Disk[evos[[j, k]], 0.05 ], {k, ns/2 + 1, ns}]
     }, PlotRange -> {{-10, 6}, {-4, 6}}], {j, 1, tmax, 10}];
ListAnimate[p1]