Hello David,
This is my code for the trajectory of the Pioneer 10 spacecraft using a small solar system, (Sun, Earth, and Jupiter) and only accounting for gravitational forces. I'm working in AU units and trying to use the correct initial conditions for position and velocity at launch to the location of the craft in 1987. That's why the time of integration stops at 15 years after launch. So for the project I want to make sure I get accurate coordinates in the (x,y,z) plane from the Sun and compare my calculation with measured data.
In[1]:=
Clear[r1,r2,r3,r4,m1,m2,m3,m4,G,M,x,y,z,t]
r1={x1[t],y1[t],z1[t]};
r2={x2[t],y2[t],z2[t]};
r3={x3[t],y3[t],z3[t]};
r4={x4[t],y4[t],z4[t]};
r12=r2-r1;
r13=r3-r1;
r14=r4-r1;
r23=r3-r2;
r21=r1-r2;
r24=r4-r2;
r31=r1-r3;
r32=r2-r3;
r34=r4-r3;
r41=r1-r4;
r42=r2-r4;
r43=r3-r4;
equation1= D[r1,{t,2}]- G m2 r12/(r12.r12)^(3/2+\[Delta])-
m3 G r13/(r13.r13)^(3/2+\[Delta])- G m4 r14/(r14.r14)^(3/2+\[Delta]);
equation2= D[r2,{t,2}]-m1 G r21/(r21.r21)^(3/2+\[Delta])-
m3 G r23/(r23.r23)^(3/2+\[Delta])- G m4 r24/(r24.r24)^(3/2+\[Delta]);
equation3= D[r3,{t,2}]- m2 G r32/(r32.r32)^(3/2+\[Delta])-
m1 G r31/(r31.r31)^(3/2+\[Delta])- G m4 r34/(r34.r34)^(3/2+\[Delta]);
equation4= D[r4,{t,2}]- G m2 r42/(r42.r42)^(3/2+\[Delta])-
m3 G r43/(r43.r43)^(3/2+\[Delta])- G m1 r41/(r41.r41)^(3/2+\[Delta]);
In[26]:=
Clear[M,m1,m2,m3,m4,G]
list1=Thread[(r1/.t->0)-{0,0,0}==0];
list2=Thread[(D[r1,t]/.t->0)-{0,0,0}==0];
list3=Thread[(r2/.t->0)-{1,0,0}==0];
list4=Thread[(D[r2,t]/.t->0)-{0,Sqrt[G M/x2[0]],0}==0];
list5=Thread[(r3/.t->0)-{4.95,0,0}==0];
list6=Thread[(D[r3,t]/.t->0)-{0,-Sqrt[G M/x3[0]],0}==0];
list7=Thread[(r4/.t->0)-{0.999940891,1.000041014,1.000000001}==0];
list8=Thread[(D[r4,t]/.t->0)-{0,(2.4512*10^-7),0}==0];
EquationList=Join[Thread[equation1==0],
Thread[equation2==0],
Thread[equation3==0],Thread[equation4==0],list1,list2,list3,list4,list5,list6,list7,list8];
In[36]:=
G=(4 \[Pi]^2)/M;
M=m1+m2+m3+m4;
m1=1;
m2=0.000003003;
m3=0.001;
m4=0;
tmax=15;
\[Delta]=0;
In[44]:=
Clear[x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4]
Planets=NDSolve[Flatten[{EquationList,WhenEvent[Sqrt[{x4[t],y4[t],z4[t]}.{x4[t],y4[t],z4[t]}]==40.01,{"StopIntegration",tstop=t}]}],Join[r1,r2,r3,r4],{t,0,tmax},
MaxSteps->100000,AccuracyGoal->10,WorkingPrecision->20]//Quiet;
Planet1[t_]=r1/.Planets;
Planet2[t_]=r2/.Planets;
Planet3[t_]=r3/.Planets;
Planet4[t_]=r4/.Planets
In[50]:= Planet4[0.000340921065169852074016373895628471032722346106772892051`20]