# How to plot multiple orbits with more than two functions

Using a program given to me (which was designed to have a particle orbiting another of infinite mass), I tried to extend it to the motion of two particles orbiting a common centre of mass. I have a few functions of multiple variables and I am unsure as to how to make sure the functions are all calling at the same rate. Here is my attempt at an extension:

Clear["Global"];
dt = 0.04;
Gm1 = 1;
Gm2 = 10*Gm1;
v1i = {0, 1.2};
v2i = {0, 1.0};
r1i = {1, 0};
r2i = {-1, 0};

f1[r1_, r2_] := -Gm2*(r1 - r2)/((r1 - r2).(r1 - r2))^(3/2);
f2[r1_, r2_] := -Gm1*(r2 - r1)/((r1 - r2).(r1 - r2))^(3/2);
v1[{v1pre_, r1pre_, r2pre_}] := v1pre + f1[r1pre, r2pre]*dt
v2[{v2pre_, r2pre_, r1pre_}] := v2pre + f2[r1pre, r2pre]*dt
r1[{v1pre_, r1pre_}] := r1pre + v1[r1pre]*dt
r2[{v2pre_, r2pre_}] := r2pre + v2[r2pre]*dt


Originally, there was only 1 velocity, force and position function. Now there is two of each. Essentially all I want is a function that will plot this. If there is any.

• @RunnyKine - how do you make it so that the actual program shows up like that? – Karl Sep 25 '14 at 2:27
• At the top of the input box, there are tools for editing your post. You can highlight your code and click on the curly braces button. – RunnyKine Sep 25 '14 at 2:43
• en.wikipedia.org/wiki/Two-body_problem – Dr. belisarius Sep 25 '14 at 3:07
• @belisarius - I sent you an email. I hope you received it. – Karl Sep 26 '14 at 3:50

dt=0.04 is too large to plot accurate orbits.

Arguments for v1,v2,r1,r2 are adjusted.

Clear["Global"];
dt = 0.0001;
Gm1 = 1;
Gm2 = 10*Gm1;
v1i = {0, 1.2};
v2i = {0, 1.0};
r1i = {1, 0};
r2i = {-1, 0};

f1[r1_, r2_] := -Gm2*(r1 - r2)/((r1 - r2).(r1 - r2))^(3/2);
f2[r1_, r2_] := -Gm1*(r2 - r1)/((r1 - r2).(r1 - r2))^(3/2);
v1[{v1pre_, r1pre_, r2pre_}] := v1pre + f1[r1pre, r2pre]*dt
v2[{v2pre_, r1pre_, r2pre_}] := v2pre + f2[r1pre, r2pre]*dt
r1[{v1pre_, r1pre_, r2pre_}] := r1pre + v1[{v1pre, r1pre, r2pre}]*dt
r2[{v2pre_, r1pre_, r2pre_}] := r2pre + v2[{v2pre, r1pre, r2pre}]*dt

ListLinePlot@
Transpose[
Take[#, -2] & /@
NestList[{v1[Drop[#, {2}]], v2[Drop[#, {1}]], r1[Drop[#, {2}]],
r2[Drop[#, {1}]]} &, {v1i, v2i, r1i, r2i}, 50000], {2, 1}] 