(* Define some constants *)
Gg=6.674 (10^(-11));(*gravity constant*)
OneYear=31557600;(*one year: 365.25 days in seconds*)
(*Planetary data, mostly from https://solarsystem.nasa.gov/planets/index.cfm and some from Wikipedia*)
Msun=1.989 (10^30);
(*For each planet: mass, semimajor axis, perihelion, eccentricity, period in days, inclination, ascending node, argument of periapsis*)
Mmercury=3.3010 (10^23);RAmercury=5.7909227 (10^10);RPmercury=4.600 (10^10);Emercury=0.20563593;Pmercury=87.97;ιmercury=7.005;Ωmercury=48.331;ωmercury=29.124;
Mvenus=4.8673 (10^24);RAvenus=1.0820948 (10^11);RPvenus=1.07476 (10^11);Evenus=0.00677672;Pvenus=224.70;ιvenus=3.39458;Ωvenus=76.678;ωvenus=55.186;
Mearth=5.9722 (10^24)+7.3477 (10^22);RAearth=1.4959826 (10^11);RPearth=1.47098 (10^11);Eearth=0.01671123;Pearth=365.26;ιearth=0.00005;Ωearth=−11.26064;ωearth=102.94719;
Mmars=6.4169 (10^23);RAmars=2.2794382 (10^11);RPmars=2.06655 (10^11);Emars=0.0933941;Pmars=686.98;ιmars=1.850;Ωmars=49.562;ωmars=286.537;
Mjupiter=1.8981 (10^27);RAjupiter=7.7834082 (10^11);RPjupiter=7.40680 (10^11);Ejupiter=0.04838624;Pjupiter=4332.82;ιjupiter=1.305;Ωjupiter=100.492;ωjupiter=275.066;
Msaturn=5.6832 (10^26);RAsaturn=1.4266664 (10^12);RPsaturn=1.34982 (10^12);Esaturn=0.05386179;Psaturn=10755.70;ιsaturn=2.485240;Ωsaturn=113.642811;ωsaturn=336.013862;
Muranus=8.6810 (10^25);RAuranus=2.8706582 (10^12);RPuranus=2.73500 (10^12);Euranus=0.04725744;Puranus=30687.15;ιuranus=0.772556;Ωuranus=73.999342;ωuranus=96.998857;
Mneptune=1.0241 (10^26);RAneptune=4.4983964 (10^12);RPneptune=4.45975 (10^12);Eneptune=0.00859048;Pneptune=60190.03;ιneptune=1.767975;Ωneptune=131.782974;ωneptune=273.219414;(*put the above information into arrays*)
Pnames={"Mercury","Venus  ","Earth  ","Mars   ","Jupiter","Saturn ","Uranus ","Neptune"};Pmass={Mmercury,Mvenus,Mearth,Mmars,Mjupiter,Msaturn,Muranus,Mneptune};Psemimaj={RAmercury,RAvenus,RAearth,RAmars,RAjupiter,RAsaturn,RAuranus,RAneptune};Pperihel={RPmercury,RPvenus,RPearth,RPmars,RPjupiter,RPsaturn,RPuranus,RPneptune};Peccen={Emercury,Evenus,Eearth,Emars,Ejupiter,Esaturn,Euranus,Eneptune};Pperdays={Pmercury,Pvenus,Pearth,Pmars,Pjupiter,Psaturn,Puranus,Pneptune};Pinclin={ιmercury,ιvenus,ιearth,ιmars,ιjupiter,ιsaturn,ιuranus,ιneptune};Pascenode={Ωmercury,Ωvenus,Ωearth,Ωmars,Ωjupiter,Ωsaturn,Ωuranus,Ωneptune};
Pargperi={ωmercury,ωvenus,ωearth,ωmars,ωjupiter,ωsaturn,ωuranus,ωneptune};
(*when the Sun is floating, boost the planet's intitial velocity by this factor to get the correct orbital period*)
(*adjustment for perihelion starting point*)
PveladjP={1.00003630,1.00004,1.00005,1.00003,1.000524,1.000186,1,1};
(*adjustment for semimajor axis starting point*)PveladjS={1.000044,1.00004,1.000045,1.000036,1.000529,1.00019,0.999898,0.999867};
(*adjustment for aphelion starting point*)PveladjA={1.0000835,1.000041,1.000048,1.000044,1.000534,1.000197,0.999885,0.999865};

(*This script will determine the perihelion advance of Mercury using the simple "Ring Planet" model. The calculations follow equations 12 and 15 in this paper: http://www.math.toronto.edu/~colliand/426_03/Papers03/C_Pollock.pdf*)
G=Gg;sum=0;F0=-((G π Msun Mmercury)/(RAmercury^2));
For[i=2,i≤8,i++,
Fa=G π Mmercury((Pmass[[i]] RAmercury)/(2π Psemimaj[[i]](((Psemimaj[[i]])^2)-(RAmercury^2))));Fda=G π Mmercury RAmercury ((Pmass[[i]] (((Psemimaj[[i]])^2)+(RAmercury^2)))/(2π Psemimaj[[i]](((((Psemimaj[[i]])^2)-(RAmercury^2)))^2)));
psi=π/(Sqrt[3+((-2 F0+Fda)/(F0+Fa))])-π;
arcseccent=(2π(1+psi)-2π)(180/π 3600) (100 Pearth/Pmercury);
Print[Pnames[[i]],"   ",arcseccent];
sum+=arcseccent;
];Print["Total = ",sum];

(* Orbit simulation equations - 2D.  This will be used for the test simulations involving "sun", "mercury", "venus" etc. *)
rdist[ptx_,pty_,psx_,psy_]:=√((((ptx-psx))^2)+(((pty-psy))^2));
(*total force*)Ftot[ptm_,ptx_,pty_,psm_,psx_,psy_]:=-G ptm psm/((rdist[ptx,pty,psx,psy])^2) ;
(*x/y force components*)Fxx[ptm_,ptx_,pty_,psm_,psx_,psy_]:=Ftot[ptm,ptx,pty,psm,psx,psy](ptx-psx) /rdist[ptx,pty,psx,psy];
Fyy[ptm_,ptx_,pty_,psm_,psx_,psy_]:=Ftot[ptm,ptx,pty,psm,psx,psy](pty-psy) /rdist[ptx,pty,psx,psy];
ma=1;xa=2;ya=3;
Fx[pT_,pS_]:=
Fxx[pT[[ma]],pT[[xa]],pT[[ya]],pS[[ma]],pS[[xa]],pS[[ya]]];
Fy[pT_,pS_]:=
Fyy[pT[[ma]],pT[[xa]],pT[[ya]],pS[[ma]],pS[[xa]],pS[[ya]]];
(*x/y acceleration components, as vectors*)Ax[pT_,pS_]:=Fx[pT,pS]/pT[[ma]];
Ay[pT_,pS_]:=Fy[pT,pS]/pT[[ma]];

(* Orbit simulation equations - 3D.  This will be used for the test proper simulations involving actual Sun and planets. *)
rdist[ptx_,pty_,ptz_,psx_,psy_,psz_]:=√((((ptx-psx))^2)+(((pty-psy))^2)+(((ptz-psz))^2));
(*total force*)Ftot[ptm_,ptx_,pty_,ptz_,psm_,psx_,psy_,psz_]:=-G ptm psm/((rdist[ptx,pty,ptz,psx,psy,psz])^2) ;
(*x/y/z force components*)Fxx[ptm_,ptx_,pty_,ptz_,psm_,psx_,psy_,psz_]:=Ftot[ptm,ptx,pty,ptz,psm,psx,psy,psz](ptx-psx) /rdist[ptx,pty,ptz,psx,psy,psz];
Fyy[ptm_,ptx_,pty_,ptz_,psm_,psx_,psy_,psz_]:=Ftot[ptm,ptx,pty,ptz,psm,psx,psy,psz](pty-psy) /rdist[ptx,pty,ptz,psx,psy,psz];
Fzz[ptm_,ptx_,pty_,ptz_,psm_,psx_,psy_,psz_]:=Ftot[ptm,ptx,pty,ptz,psm,psx,psy,psz](ptz-psz) /rdist[ptx,pty,ptz,psx,psy,psz];
ma=1;xa=2;ya=3;za=4;
Fx[pT_,pS_]:=
Fxx[pT[[ma]],pT[[xa]],pT[[ya]],pT[[za]],pS[[ma]],pS[[xa]],pS[[ya]],pS[[za]]];
Fy[pT_,pS_]:=
Fyy[pT[[ma]],pT[[xa]],pT[[ya]],pT[[za]],pS[[ma]],pS[[xa]],pS[[ya]],pS[[za]]];
Fz[pT_,pS_]:=
Fzz[pT[[ma]],pT[[xa]],pT[[ya]],pT[[za]],pS[[ma]],pS[[xa]],pS[[ya]],pS[[za]]];
(*x/y/z acceleration components, as vectors*)Ax[pT_,pS_]:=Fx[pT,pS]/pT[[ma]];
Ay[pT_,pS_]:=Fy[pT,pS]/pT[[ma]];
Az[pT_,pS_]:=Fz[pT,pS]/pT[[ma]];

(* Dummy Orbit simulation equations - 3D.  This is used to force Mathematica to use higher accuracy *)
Qrdist[pdx_,pdy_,pdz_]:=√((pdx^2)+(pdy^2)+(pdz^2));
(*total force*)
QFtot[ptm_,psm_,pdx_,pdy_,pdz_]:=-G ptm psm/((Qrdist[pdx,pdy,pdz])^2) ;
(*x/y/z force components*)
QFxx[ptm_,psm_,pdx_,pdy_,pdz_]:=QFtot[ptm,psm,pdx,pdy,pdz]pdx/Qrdist[pdx,pdy,pdz];QFyy[ptm_,psm_,pdx_,pdy_,pdz_]:=QFtot[ptm,psm,pdx,pdy,pdz]pdy/Qrdist[pdx,pdy,pdz];QFzz[ptm_,psm_,pdx_,pdy_,pdz_]:=QFtot[ptm,psm,pdx,pdy,pdz]pdz/Qrdist[pdx,pdy,pdz];
(*x/y/z force components, as vectors*)
QFx[pT_,pS_]:=
QFxx[pT[[ma]],pS[[ma]],pT[[xa]]-pS[[xa]],pT[[ya]]-pS[[ya]],pT[[za]]-pS[[za]]];
QFy[pT_,pS_]:=
QFyy[pT[[ma]],pS[[ma]],pT[[xa]]-pS[[xa]],pT[[ya]]-pS[[ya]],pT[[za]]-pS[[za]]];
QFz[pT_,pS_]:=
QFzz[pT[[ma]],pS[[ma]],pT[[xa]]-pS[[xa]],pT[[ya]]-pS[[ya]],pT[[za]]-pS[[za]]];
(*x/y/z acceleration components, as vectors*)
QAx[pT_,pS_]:=QFx[pT,pS]/pT[[ma]];
QAy[pT_,pS_]:=QFy[pT,pS]/pT[[ma]];
QAz[pT_,pS_]:=QFz[pT,pS]/pT[[ma]];

(*Transformation function.  This takes an x/y/z coordinate and translates it using ι/Ω/ω (inclination/ascending node/arg of periapsis). It outputs a vector(array) containing new x/y/z coordinates*)OrbTransXYZ[x_,y_,z_,ι_,Ω_,ω_]:={x (Cos[ι] Cos[ω] Cos[Ω]-Sin[ω] Sin[Ω])+y(-Cos[ι] Cos[Ω] Sin[ω]-Cos[ω] Sin[Ω])+z(-Cos[Ω] Sin[ι]),x (Cos[Ω] Sin[ω]+Cos[ι] Cos[ω] Sin[Ω])+y (Cos[ω] Cos[Ω]-Cos[ι] Sin[ω] Sin[Ω])+z(-Sin[ι] Sin[Ω]),x(Cos[ω] Sin[ι])+y(-Sin[ι] Sin[ω])+z Cos[ι]};
(*This does the same but using ι/Ω/ω in degrees rather than radians.*)OrbTransXYZdeg[x_,y_,z_,ι_,Ω_,ω_]:=OrbTransXYZ[x,y,z ,ι π/180,Ω π/180,ω π/180];

(*Inverse transformation function.  This takes a translated x/y/z coordinate and undoes it using ι/Ω/ω (inclination/ascending node/arg of periapsis). It outputs a vector(array) containing the original x/y/z coordinates*)
OrbUntran[x_,y_,z_,ι_,Ω_,ω_]:={x(Cos[ι] Cos[ω] Cos[Ω]-Sin[ω] Sin[Ω])+y(Cos[Ω] Sin[ω]+Cos[ι] Cos[ω] Sin[Ω])+z(Cos[ω] Sin[ι]),x(-Cos[ι] Cos[Ω] Sin[ω]-Cos[ω] Sin[Ω])+y(Cos[ω] Cos[Ω]-Cos[ι] Sin[ω] Sin[Ω])+z(-Sin[ι] Sin[ω]),x(-Cos[Ω] Sin[ι])+y(-Sin[ι] Sin[Ω])+z(Cos[ι])};(*This does the same but using ι/Ω/ω in degrees rather than radians.*)OrbUntranDeg[x_,y_,z_,ι_,Ω_,ω_]:=OrbUntran[x,y,z,ι π/180,Ω π/180,ω π/180];

(*Functions to describe "away velocity" (velocity of Mercury away from Sun) and the angle of Mercury around Sun in arcseconds.  These will be used in the Perihelion Finder scripts below.  There are several versions of these functions.  The one to use depends on the number of dimensions (2 or 3), whether the Sun is fixed or floating, and whether starting point of Mercury has been rotated using "regular" or HORIZONS coordinates.  The functions will be labelled VALANG1 VELANG2 etc.*)

(*VELANG1: fixed sun, 2 dimensions*)VelAwy[t_]:=(x1[t]x1'[t]+y1[t]y1'[t])/(Sqrt[((x1[t])^2)+((y1[t])^2)])/.result[[1]];AngArcsec[t_]:=ArcTan[x1[t],y1[t]]/π 180 3600/.result[[1]];

(*VELANG2: floating sun, 2 dimensions*)VelAwy[t_]:=((x1[t]-x0[t])(x1'[t]-x0'[t])+(y1[t]-y0[t])(y1'[t]-y0'[t]))/(Sqrt[(((x1[t]-x0[t]))^2)+(((y1[t]-y0[t]))^2)])/.result[[1]];AngArcsec[t_]:=ArcTan[(x1[t]-x0[t]),(y1[t]-y0[t])]/π 180 3600/.result[[1]];

(*VELANG3: floating sun, 2 dimensions with orbits rotated but not inclined*)VelAwy[t_]:=((x1[t]-x0[t])(x1'[t]-x0'[t])+(y1[t]-y0[t])(y1'[t]-y0'[t]))/(Sqrt[(((x1[t]-x0[t]))^2)+(((y1[t]-y0[t]))^2)])/.result[[1]];AngArcsec[t_]:=ArcTan[OrbUntranDeg[(x1[t]-x0[t]),(y1[t]-y0[t]),0,0ιmercury,Ωmercury,ωmercury][[1]]
,OrbUntranDeg[(x1[t]-x0[t]),(y1[t]-y0[t]),0,0ιmercury,Ωmercury,ωmercury][[2]]]/π 180 3600/.result[[1]];

(*VELANG4: fixed sun, 3 dimensions*)VelAwy[t_]:=(x1[t]x1'[t]+y1[t]y1'[t]+z1[t]z1'[t])/(Sqrt[((x1[t])^2)+((y1[t])^2)+((z1[t])^2)])/.result[[1]];AngArcsec[t_]:=ArcTan[OrbUntranDeg[x1[t],y1[t],z1[t],ιmercury,Ωmercury,ωmercury][[1]]
,OrbUntranDeg[x1[t],y1[t],z1[t],ιmercury,Ωmercury,ωmercury][[2]]]/π 180 3600/.result[[1]];

(*VELANG5: floating sun, 3 dimensions*)VelAwy[t_]:=((x1[t]-x0[t])(x1'[t]-x0'[t])+(y1[t]-y0[t])(y1'[t]-y0'[t])+(z1[t]-z0[t])(z1'[t]-z0'[t]))/(Sqrt[(((x1[t]-x0[t]))^2)+(((y1[t]-y0[t]))^2)+(((z1[t]-z0[t]))^2)])/.result[[1]];AngArcsec[t_]:=ArcTan[OrbUntranDeg[x1[t]-x0[t],y1[t]-y0[t],z1[t]-z0[t],ιmercury,Ωmercury,ωmercury][[1]]
,OrbUntranDeg[x1[t]-x0[t],y1[t]-y0[t],z1[t]-z0[t],ιmercury,Ωmercury,ωmercury][[2]]]/π 180 3600/.result[[1]];

(*VELANG6: floating sun, 3 dimensions, HORIZONS coords*)VelAwy[t_]:=((x1[t]-x0[t])(x1'[t]-x0'[t])+(y1[t]-y0[t])(y1'[t]-y0'[t])+(z1[t]-z0[t])(z1'[t]-z0'[t]))/(Sqrt[(((x1[t]-x0[t]))^2)+(((y1[t]-y0[t]))^2)+(((z1[t]-z0[t]))^2)])/.result[[1]];AngArcsec[t_]:=ArcTan[OrbUntranDeg[y1[t]-y0[t],-(x1[t]-x0[t]),z1[t]-z0[t],ιmercury,Ωmercury,ωmercury-90][[1]]
,OrbUntranDeg[y1[t]-y0[t],-(x1[t]-x0[t]),z1[t]-z0[t],ιmercury,Ωmercury,ωmercury-90][[2]]]/π 180 3600/.result[[1]];

(*Perihelion Finder Script 1: Find perihelions and create chart in the test simulations (containing "sun" and "mercury")*)
PeriAng={};StdPeriod=2π(Sqrt[(((R1/(1-e1)))^3)/(G M0)]);
prevMidT=0;firstAng=True;prevAng=0;startAng=1/3600AngArcsec[0];tt=0;divn=0.001;
While[tt<Per-divn ,
While[tt>divn && VelAwy[tt-divn]>0,tt-=divn];
While[tt<Per-divn && VelAwy[tt+divn]<0,tt+=divn];
If[firstAng==False,lowerT=tt-divn;upperT=tt+divn,lowerT=0;upperT=2 divn];
middleT=lowerT+(upperT-lowerT)(VelAwy[lowerT]/(VelAwy[lowerT]-VelAwy[upperT]));
If[middleT≥0 && middleT<Per,
thisAng=AngArcsec[middleT];  ,
thisAng=AngArcsec[lowerT]+(AngArcsec[upperT]-AngArcsec[lowerT])(middleT-lowerT)/(upperT-lowerT);
];
thisAng/=3600;
AppendTo[PeriAng,{middleT/60,thisAng-startAng}];
tt=middleT;
If[firstAng==False,StdPeriod=middleT-prevMidT];
prevMidT=middleT;
firstAng=False;
tt+=StdPeriod;
]

(*Perihelion Finder Script 2: Find perihelions and create chart in the regular simulations (containing Sun and Mercury)*)
PeriAng={};StdPeriod=87.969 86400;prevMidT=0;firstAng=True;prevAng=0;numdiffAng=0;totdiffAng=0;avgdiffAng=0;startAng=AngArcsec[0];tt=0;
(*get average angle move over first 10 years*)
While[tt<10 OneYear && tt<Per-100 ,
If[firstAng==False,lowerT=tt-100;upperT=tt+100,lowerT=0;upperT=200];
middleT=lowerT+(upperT-lowerT)(VelAwy[lowerT]/(VelAwy[lowerT]-VelAwy[upperT]));
If[middleT≥0 && middleT<Per,
thisAng=AngArcsec[middleT];  ,
thisAng=AngArcsec[lowerT]+(AngArcsec[upperT]-AngArcsec[lowerT])(middleT-lowerT)/(upperT-lowerT);
];
If[firstAng==True,startAng=thisAng];
diffAng=Abs[thisAng-prevAng];
If[firstAng==False && (numdiffAng==0 || diffAng<5totdiffAng/numdiffAng ),
totdiffAng+=diffAng;
numdiffAng++
];
tt=middleT;
If[firstAng==False,StdPeriod=middleT-prevMidT];
prevMidT=middleT;
prevAng=thisAng;
firstAng=False;
tt+=StdPeriod;
];
If[numdiffAng>0,avgdiffAng=totdiffAng/numdiffAng];
(*now accumulate the full period*)
StdPeriod=87.969 86400;prevMidT=0;firstAng=True;prevAng=0;missAng=1;tt=0;
While[tt<Per-100 ,
If[firstAng==False,lowerT=tt-100;upperT=tt+100,lowerT=0;upperT=200];
middleT=lowerT+(upperT-lowerT)(VelAwy[lowerT]/(VelAwy[lowerT]-VelAwy[upperT]));
If[middleT≥0 && middleT<Per,
thisAng=AngArcsec[middleT];  ,
thisAng=AngArcsec[lowerT]+(AngArcsec[upperT]-AngArcsec[lowerT])(middleT-lowerT)/(upperT-lowerT);
];
diffAng=Abs[thisAng-prevAng];
If[diffAng<2 missAng avgdiffAng || avgdiffAng==0,
AppendTo[PeriAng,{middleT/OneYear,thisAng-startAng}];
prevAng=thisAng;
missAng=1,
missAng++;
];
tt=middleT;
If[firstAng==False,StdPeriod=middleT-prevMidT];
prevMidT=middleT;
prevAng=thisAng;
firstAng=False;
tt+=StdPeriod;
]

(*Test perihelion advance with a fixed "sun" orbited by "mercury" and a fixed outer planet. This requires the 2D orbit simulation equations. Do this with AA=1,2,1/2,1/10,-1. For AA=1, this takes 5050sec to fully rotate the perihelion*)
AA=1;M0=50;M1=1/100;M2=2;
Pc0={M0,0,0};
Pc1={M1,x1[t],y1[t]};
G=1;R1=1;R2=7;e1=0.4;V1=(Sqrt[(1+e1)((G M0)/R1)]);V2=AA(Sqrt[((G M0)/R2)]);Pc2={M2,R2 Cos[(V2/R2)t],R2 Sin[(V2/R2)t]};Per=600;result=NDSolve[{ x1''[t]==Ax[Pc1,Pc0]+Ax[Pc1,Pc2], y1''[t]==Ay[Pc1,Pc0]+Ay[Pc1,Pc2],x1[0]==R1,y1[0]==0,x1'[0]==0,y1'[0]==V1},{x1,y1},{t,0,Per},MaxSteps->400000];ParametricPlot[{Evaluate[{x1[t],y1[t]}/.result],Evaluate[{Pc2[[2]],Pc2[[3]]}/.result]},{t,(Per-17),Per},AspectRatio->Automatic,PlotPoints->500,PlotRange->Automatic]

(*After running the above simulation, do the following: 
1. Run the function VELANG1
2. Run the Perihelion Finder Script 1
3. Plot the result using ListPlot *)

ListPlot[PeriAng,Joined->True]

(*Get line of best fit - in degrees per minute*)Normal[LinearModelFit[PeriAng,x,x]][[2]]/x

(*Simulate Mercury orbiting a fixed Sun and an outer planet in a fixed circular orbit.  The outer planet can be selected by changing the variable pno, e.g. pno=2 is Venus and pno=8 is Neptune.  
The outer planet has a circular orbit with radius equal to its semimajor axis.  To boost precision this will be done in 3D - even though only 2 dimensions are used.  Also a dummy inner planet is added inside Mercury's orbit.  This dummy planet forces Mathematica to use a smaller timestep but does not effect the other bodies in any way.
This simulation will be repeated further down with super-high accuracy using a different method*)
pno=2;M0=Msun;M1=Mmercury;M2=Pmass[[pno]];
(*Pc0=Sun, Pc1=Mercury, Pc2=outer planet, PcQ=dummy*)
Pc0={M0,0,0,0};
Pc1={M1,x1[t],y1[t],z1[t]};
G=Gg;R1=RPmercury;R2=Psemimaj[[pno]];e1= Emercury;V1=(Sqrt[(1+e1)((G M0)/R1)]);V2=(Sqrt[((G M0)/R2)]);Pc2={M2,R2 Cos[(V2/R2)t],R2 Sin[(V2/R2)t],0};PcQ={M1,xQ[t],yQ[t],zQ[t]};RQ=RPmercury*0.45;VQ=(Sqrt[((G M0)/RQ)]);Per=100 OneYear;result=NDSolve[{ x1''[t]==Ax[Pc1,Pc0]+Ax[Pc1,Pc2], y1''[t]==Ay[Pc1,Pc0]+Ay[Pc1,Pc2], z1''[t]==Az[Pc1,Pc0]+Az[Pc1,Pc2],
xQ''[t]==QAx[PcQ,Pc0],yQ''[t]==QAy[PcQ,Pc0],zQ''[t]==QAz[PcQ,Pc0],x1[0]==R1,y1[0]==0,z1[0]==0,x1'[0]==0,y1'[0]==V1,z1'[0]==0,xQ[0]==RQ,yQ[0]==0,zQ[0]==0,xQ'[0]==0,yQ'[0]==VQ,zQ'[0]==0},{x1,y1,z1,xQ,yQ,zQ},{t,0,Per},MaxSteps->2000000,AccuracyGoal->10];

(*After running the above simulation, do the following: 
1. Run the function VELANG1
2. Run the Perihelion Finder Script 2
3. Plot the result using ListPlot
4. Rerun the simulation and steps 2 & 3 for pno=3 to 8*)

ListPlot[PeriAng,Joined->True]

(*Get line of best fit - in arcseconds per 100 years*)Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Simulate test fixed "sun", floating "mercury" and fixed outer "jupiter" orbit.*)
M0=50;M1=1/100;M2=4;
Pc0={M0,0,0};
Pc1={M1,x1[t],y1[t]};
G=1;R1=1;R2=14;e1=0.4;V1=(Sqrt[(1+e1)((G M0)/R1)]);V2=(Sqrt[((G M0)/R2)]);Pc2={M2,R2 Cos[(V2/R2)t],R2 Sin[(V2/R2)t]};Per=1188.5;result=NDSolve[{ x1''[t]==Ax[Pc1,Pc0]+Ax[Pc1,Pc2], y1''[t]==Ay[Pc1,Pc0]+Ay[Pc1,Pc2],x1[0]==R1,y1[0]==0,x1'[0]==0,y1'[0]==V1},{x1,y1},{t,0,Per},MaxSteps->400000];

(*After running the above simulation, do the following: 
1. Run the function VELANG1 
2. Run the Perihelion Finder Script 1
3. Plot the result using ListPlot *)

ListPlot[PeriAng,Joined->True]

(*Get line of best fit - in degrees per minute*)Normal[LinearModelFit[PeriAng,x,x]][[2]]/x

(*Simulate test fixed "sun", floating "mercury" and fixed outer "venus"+"jupiter" orbit.*)
M0=50;M1=1/100;M2=2;M3=4;
Pc0={M0,0,0};
Pc1={M1,x1[t],y1[t]};
G=1;R1=1;R2=7;R3=14;e1=0.4;V1=(Sqrt[(1+e1)((G M0)/R1)]);V2=(Sqrt[((G M0)/R2)]);V3=(Sqrt[((G M0)/R3)]);Pc2={M2,R2 Cos[(V2/R2)t],R2 Sin[(V2/R2)t]};Pc3={M3,R3 Cos[(V3/R3)t],R3 Sin[(V3/R3)t]};Per=1200;result=NDSolve[{ x1''[t]==Ax[Pc1,Pc0]+Ax[Pc1,Pc2]+Ax[Pc1,Pc3], y1''[t]==Ay[Pc1,Pc0]+Ay[Pc1,Pc2]+Ay[Pc1,Pc3],x1[0]==R1,y1[0]==0,x1'[0]==0,y1'[0]==V1},{x1,y1},{t,0,Per},MaxSteps->400000];

(*After running the above simulation, do the following: 
1. Run the function VELANG1
2. Run the Perihelion Finder Script 1
3. Plot the result using ListPlot *)

ListPlot[PeriAng,Joined->True]

(*Get line of best fit - in degrees per minute*)Normal[LinearModelFit[PeriAng,x,x]][[2]]/x

(*Simulate Mercury orbiting a fixed Sun and surrounded by Jupiter and Saturn.  The outer planets will have a circular orbit with a radius equal to their semimajor axis.*){M0,M1,M5,M6}={Msun,Mmercury,Mjupiter,Msaturn};
Pc0={M0,0,0,0};
Pc1={M1,x1[t],y1[t],z1[t]};
G=Gg;{R1,R5,R6}={RPmercury,RAjupiter,RAsaturn};e1= Emercury;{V1,V5,V6}={(Sqrt[(1+e1)G M0/R1]),(Sqrt[(G M0)/R5]),(Sqrt[(G M0)/R6])};Pc5={M5,R5 Cos[(V5/R5)t],R5 Sin[(V5/R5)t],0};Pc6={M6,R6 Cos[(V6/R6)t],R6 Sin[(V6/R6)t],0};PcQ={M1,xQ[t],yQ[t],zQ[t]};RQ=RPmercury*0.45;VQ=(Sqrt[((G M0)/RQ)]);Per=113 OneYear;result=NDSolve[{ x1''[t]==Ax[Pc1,Pc0]+Ax[Pc1,Pc5]+Ax[Pc1,Pc6], y1''[t]==Ay[Pc1,Pc0]+Ay[Pc1,Pc5]+Ay[Pc1,Pc6], z1''[t]==Az[Pc1,Pc0]+Az[Pc1,Pc5]+Az[Pc1,Pc6],
xQ''[t]==QAx[PcQ,Pc0],yQ''[t]==QAy[PcQ,Pc0],zQ''[t]==QAz[PcQ,Pc0],x1[0]==R1,y1[0]==0,z1[0]==0,x1'[0]==0,y1'[0]==V1,z1'[0]==0,xQ[0]==RQ,yQ[0]==0,zQ[0]==0,xQ'[0]==0,yQ'[0]==VQ,zQ'[0]==0},{x1,y1,z1,xQ,yQ,zQ},{t,0,Per},MaxSteps->2000000,AccuracyGoal->10];

(*After running the above simulation, do the following: 
1. Run the function VELANG1
2. Run the Perihelion Finder Script 2
3. Plot the result using ListPlot
4. Rerun the simulation and steps 2 & 3 for pno=3 to 8*)

ListPlot[PeriAng,Joined->True]

(*Get line of best fit - in arcseconds per 100 years*)Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Simulate Mercury orbiting a fixed Sun and surrounded by all the outer planets in fixed orbits.  The outer planets will have a circular orbit with a radius equal to their semimajor axis.*){M0,M1,M2,M3,M4,M5,M6,M7,M8}={Msun,Mmercury,Mvenus,Mearth,Mmars,Mjupiter,Msaturn,Muranus,Mneptune};
Pc0={M0,0,0,0};
Pc1={M1,x1[t],y1[t],z1[t]};
G=Gg;{R1,R2,R3,R4,R5,R6,R7,R8}={RPmercury,RAvenus,RAearth,RAmars,RAjupiter,RAsaturn,RAuranus,RAneptune};e1= Emercury;{V1,V2,V3,V4,V5,V6,V7,V8}={(Sqrt[(1+e1)G M0/R1]),(Sqrt[(G M0)/R2]),(Sqrt[(G M0)/R3]),(Sqrt[(G M0)/R4]),(Sqrt[(G M0)/R5]),(Sqrt[(G M0)/R6]),(Sqrt[(G M0)/R7]),(Sqrt[(G M0)/R8])};Pc2={M2,R2 Cos[(V2/R2)t],R2 Sin[(V2/R2)t],0};Pc3={M3,R3 Cos[(V3/R3)t],R3 Sin[(V3/R3)t],0};Pc4={M4,R4 Cos[(V4/R4)t],R4 Sin[(V4/R4)t],0};Pc5={M5,R5 Cos[(V5/R5)t],R5 Sin[(V5/R5)t],0};Pc6={M6,R6 Cos[(V6/R6)t],R6 Sin[(V6/R6)t],0};Pc7={M7,R7 Cos[(V7/R7)t],R7 Sin[(V7/R7)t],0};Pc8={M8,R8 Cos[(V8/R8)t],R8 Sin[(V8/R8)t],0};PcQ={M1,xQ[t],yQ[t],zQ[t]};RQ=RPmercury*0.45;VQ=(Sqrt[((G M0)/RQ)]);Per=149 OneYear;result=NDSolve[{ x1''[t]==Ax[Pc1,Pc0]+Ax[Pc1,Pc2]+Ax[Pc1,Pc3]+Ax[Pc1,Pc4]+Ax[Pc1,Pc5]+Ax[Pc1,Pc6]+Ax[Pc1,Pc7]+Ax[Pc1,Pc8], y1''[t]==Ay[Pc1,Pc0]+Ay[Pc1,Pc2]+Ay[Pc1,Pc3]+Ay[Pc1,Pc4]+Ay[Pc1,Pc5]+Ay[Pc1,Pc6]+Ay[Pc1,Pc7]+Ay[Pc1,Pc8], z1''[t]==Az[Pc1,Pc0]+Az[Pc1,Pc2]+Az[Pc1,Pc3]+Az[Pc1,Pc4]+Az[Pc1,Pc5]+Az[Pc1,Pc6]+Az[Pc1,Pc7]+Az[Pc1,Pc8],
xQ''[t]==QAx[PcQ,Pc0],yQ''[t]==QAy[PcQ,Pc0],zQ''[t]==QAz[PcQ,Pc0],x1[0]==R1,y1[0]==0,z1[0]==0,x1'[0]==0,y1'[0]==V1,z1'[0]==0,xQ[0]==RQ,yQ[0]==0,zQ[0]==0,xQ'[0]==0,yQ'[0]==VQ,zQ'[0]==0},{x1,y1,z1,xQ,yQ,zQ},{t,0,Per},MaxSteps->2000000,AccuracyGoal->10];

(*After running the above simulation, do the following: 
1. Run the function VELANG1
2. Run the Perihelion Finder Script 2
3. Plot the result using ListPlot *)

ListPlot[PeriAng,Joined->True]

(*Get line of best fit - in arcseconds per 100 years*)Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Simulate test floating "sun", floating "mercury" and floating outer "venus" orbit.*)
M0=50;M1=1/100;M2=2;
Pc0={M0,x0[t],y0[t]};
Pc1={M1,x1[t],y1[t]};
Pc2={M2,x2[t],y2[t]};
G=1;R1=1;R2=7;e1= 0.4;V1=(Sqrt[(1+e1)((G M0)/R1)]);V2=2.7432;(*V2=(Sqrt[((G M0)/R2)]);*)Per=600;result=NDSolve[{ x0''[t]==Ax[Pc0,Pc2], y0''[t]==Ay[Pc0,Pc2], x1''[t]==Ax[Pc1,Pc0]+Ax[Pc1,Pc2], y1''[t]==Ay[Pc1,Pc0]+Ay[Pc1,Pc2], x2''[t]==Ax[Pc2,Pc0], y2''[t]==Ay[Pc2,Pc0],x0[0]==0,y0[0]==0,x0'[0]==0,y0'[0]==0,x1[0]==R1,y1[0]==0,x1'[0]==0,y1'[0]==V1,x2[0]==R2,y2[0]==0,x2'[0]==0,y2'[0]==V2},{x0,y0,x1,y1,x2,y2},{t,0,Per},MaxSteps->400000];(* For the above;
identical period (circle)  slope=4.78 deg/min;identical velocity (ellipse)  slope=5.684;
*)

(*After running the above simulation, do the following: 
1. Run the function VELANG2
2. Run the Perihelion Finder Script 1
3. Plot the result using ListPlot *)

ListPlot[PeriAng,Joined->True]

(*Get line of best fit - in degrees per minute*)Normal[LinearModelFit[PeriAng,x,x]][[2]]/x

(*Simulate Mercury orbiting a floating Sun and an outer planet in a circular orbit.  The outer planet can be selected by changing the variable pno, e.g. pno=2 is Venus and pno=8 is Neptune.  
The outer planet has a circular orbit with radius equal to its semimajor axis.  It's initial velocity will be boosted to ensure it is circular.*)
pno=2;M0=Msun;M1=Mmercury;M2=Pmass[[pno]];
(*Pc0=Sun, Pc1=Mercury, Pc2=outer planet, PcQ=dummy*)
Pc0={M0,x0[t],y0[t],z0[t]};
Pc1={M1,x1[t],y1[t],z1[t]};
Pc2={M2,x2[t],y2[t],z2[t]};
G=Gg;R1=RPmercury;R2=Psemimaj[[pno]];e1= Emercury;V1=(Sqrt[(1+e1)((G M0)/R1)]);V2=PveladjS[[pno]](Sqrt[((G M0)/R2)]);PcQ={M1,xQ[t],yQ[t],zQ[t]};RQ=RPmercury*0.45;VQ=(Sqrt[((G M0)/RQ)]);Per=100 OneYear;result=NDSolve[{x0''[t]==Ax[Pc0,Pc1]+Ax[Pc0,Pc2],y0''[t]==Ay[Pc0,Pc1]+Ay[Pc0,Pc2],z0''[t]==Az[Pc0,Pc1]+Az[Pc0,Pc2], x1''[t]==Ax[Pc1,Pc0]+Ax[Pc1,Pc2], y1''[t]==Ay[Pc1,Pc0]+Ay[Pc1,Pc2], z1''[t]==Az[Pc1,Pc0]+Az[Pc1,Pc2],x2''[t]==Ax[Pc2,Pc0]+Ax[Pc2,Pc1],y2''[t]==Ay[Pc2,Pc0]+Ay[Pc2,Pc1],z2''[t]==Az[Pc2,Pc0]+Az[Pc2,Pc1],
xQ''[t]==QAx[PcQ,Pc0],yQ''[t]==QAy[PcQ,Pc0],zQ''[t]==QAz[PcQ,Pc0],x0[0]==0,y0[0]==0,z0[0]==0,x0'[0]==0,y0'[0]==0,z0'[0]==0,x1[0]==R1,y1[0]==0,z1[0]==0,x1'[0]==0,y1'[0]==V1,z1'[0]==0,x2[0]==R2,y2[0]==0,z2[0]==0,x2'[0]==0,y2'[0]==V2,z2'[0]==0,xQ[0]==RQ,yQ[0]==0,zQ[0]==0,xQ'[0]==0,yQ'[0]==VQ,zQ'[0]==0},{x0,y0,z0,x1,y1,z1,x2,y2,z2,xQ,yQ,zQ},{t,0,Per},MaxSteps->2000000,AccuracyGoal->10];

(*After running the above simulation, do the following: 
1. Run the function VELANG2
2. Run the Perihelion Finder Script 2
3. Plot the result using ListPlot
4. Rerun the simulation and steps 2 & 3 for pno=3 to 8*)

ListPlot[PeriAng,Joined->True]

(*Get line of best fit - in arcseconds per 100 years*)Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Simulate Mercury orbiting a floating Sun and an outer planet in an eccentric orbit.  The outer planet can be selected by changing the variable pno, e.g. pno=2 is Venus and pno=8 is Neptune.  
The outer planet starts at its perihelion point.  It's initial velocity will be boosted to ensure the period is correct.*)
pno=2;M0=Msun;M1=Mmercury;M2=Pmass[[pno]];
(*Pc0=Sun, Pc1=Mercury, Pc2=outer planet, PcQ=dummy*)
Pc0={M0,x0[t],y0[t],z0[t]};
Pc1={M1,x1[t],y1[t],z1[t]};
Pc2={M2,x2[t],y2[t],z2[t]};
G=Gg;R1=RPmercury;R2=Pperihel[[pno]];e1= Emercury;e2=Peccen[[pno]];V1=(Sqrt[(1+e1)((G M0)/R1)]);V2=PveladjP[[pno]](Sqrt[(1+e2)(G M0)/R2]);PcQ={M1,xQ[t],yQ[t],zQ[t]};RQ=RPmercury*0.45;VQ=(Sqrt[((G M0)/RQ)]);Per=100 OneYear;result=NDSolve[{x0''[t]==Ax[Pc0,Pc1]+Ax[Pc0,Pc2],y0''[t]==Ay[Pc0,Pc1]+Ay[Pc0,Pc2],z0''[t]==Az[Pc0,Pc1]+Az[Pc0,Pc2], x1''[t]==Ax[Pc1,Pc0]+Ax[Pc1,Pc2], y1''[t]==Ay[Pc1,Pc0]+Ay[Pc1,Pc2], z1''[t]==Az[Pc1,Pc0]+Az[Pc1,Pc2],x2''[t]==Ax[Pc2,Pc0]+Ax[Pc2,Pc1],y2''[t]==Ay[Pc2,Pc0]+Ay[Pc2,Pc1],z2''[t]==Az[Pc2,Pc0]+Az[Pc2,Pc1],
xQ''[t]==QAx[PcQ,Pc0],yQ''[t]==QAy[PcQ,Pc0],zQ''[t]==QAz[PcQ,Pc0],x0[0]==0,y0[0]==0,z0[0]==0,x0'[0]==0,y0'[0]==0,z0'[0]==0,x1[0]==R1,y1[0]==0,z1[0]==0,x1'[0]==0,y1'[0]==V1,z1'[0]==0,x2[0]==R2,y2[0]==0,z2[0]==0,x2'[0]==0,y2'[0]==V2,z2'[0]==0,xQ[0]==RQ,yQ[0]==0,zQ[0]==0,xQ'[0]==0,yQ'[0]==VQ,zQ'[0]==0},{x0,y0,z0,x1,y1,z1,x2,y2,z2,xQ,yQ,zQ},{t,0,Per},MaxSteps->2000000,AccuracyGoal->10];

(*After running the above simulation, do the following: 
1. Run the function VELANG2
2. Run the Perihelion Finder Script 2
3. Plot the result using ListPlot
4. Rerun the simulation and steps 2 & 3 for pno=3 to 8*)

ListPlot[PeriAng,Joined->True]

(*Get line of best fit - in arcseconds per 100 years*)Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Simulate Mercury orbiting a floating Sun and all outer planets from Venus to Neptune with eccentric orbits.  The orbits will be rotated to match the actual planet's orbits but won't be inclined.*)
{M0,M1,M2,M3,M4,M5,M6,M7,M8}={Msun,Mmercury,Mvenus,Mearth,Mmars,Mjupiter,Msaturn,Muranus,Mneptune};
Pc0={M0,x0[t],y0[t],z0[t]};
Pc1={M1,x1[t],y1[t],z1[t]};
Pc2={M2,x2[t],y2[t],z2[t]};
Pc3={M3,x3[t],y3[t],z3[t]};
Pc4={M4,x4[t],y4[t],z4[t]};
Pc5={M5,x5[t],y5[t],z5[t]};
Pc6={M6,x6[t],y6[t],z6[t]};
Pc7={M7,x7[t],y7[t],z7[t]};
Pc8={M8,x8[t],y8[t],z8[t]};
PcQ={M1,xQ[t],yQ[t],zQ[t]};
G=Gg;{R1,R2,R3,R4,R5,R6,R7,R8}={RPmercury,RPvenus,RPearth,RPmars,RPjupiter,RPsaturn,RPuranus,RPneptune};velint=Array[0&,8];
For[i=1,i≤8,i++,velint[[i]]=PveladjP[[i]](Sqrt[(1+ Peccen[[i]])((G Msun)/(Pperihel[[i]]))])];{V1,V2,V3,V4,V5,V6,V7,V8}=velint;{Px1,Py1,Pz1}=OrbTransXYZdeg[R1,0,0,0ιmercury,Ωmercury,ωmercury];{Vx1,Vy1,Vz1}=OrbTransXYZdeg[0,V1,0,0ιmercury,Ωmercury,ωmercury];
{Px2,Py2,Pz2}=OrbTransXYZdeg[R2,0,0,0ιvenus,Ωvenus,ωvenus];{Vx2,Vy2,Vz2}=OrbTransXYZdeg[0,V2,0,0ιvenus,Ωvenus,ωvenus];
{Px3,Py3,Pz3}=OrbTransXYZdeg[R3,0,0,0ιearth,Ωearth,ωearth];{Vx3,Vy3,Vz3}=OrbTransXYZdeg[0,V3,0,0ιearth,Ωearth,ωearth];
{Px4,Py4,Pz4}=OrbTransXYZdeg[R4,0,0,0ιmars,Ωmars,ωmars];
{Vx4,Vy4,Vz4}=OrbTransXYZdeg[0,V4,0,0ιmars,Ωmars,ωmars];
{Px5,Py5,Pz5}=OrbTransXYZdeg[R5,0,0,0ιjupiter,Ωjupiter,ωjupiter];{Vx5,Vy5,Vz5}=OrbTransXYZdeg[0,V5,0,0ιjupiter,Ωjupiter,ωjupiter];
{Px6,Py6,Pz6}=OrbTransXYZdeg[R6,0,0,0ιsaturn,Ωsaturn,ωsaturn];{Vx6,Vy6,Vz6}=OrbTransXYZdeg[0,V6,0,0ιsaturn,Ωsaturn,ωsaturn];
{Px7,Py7,Pz7}=OrbTransXYZdeg[R7,0,0,0ιuranus,Ωuranus,ωuranus];{Vx7,Vy7,Vz7}=OrbTransXYZdeg[0,V7,0,0ιuranus,Ωuranus,ωuranus];
{Px8,Py8,Pz8}=OrbTransXYZdeg[R8,0,0,0ιneptune,Ωneptune,ωneptune];{Vx8,Vy8,Vz8}=OrbTransXYZdeg[0,V8,0,0ιneptune,Ωneptune,ωneptune];RQ=RPmercury*0.45;VQ=(Sqrt[((G M0)/RQ)]);Per=100 OneYear;result=NDSolve[{ x0''[t]==Ax[Pc0,Pc1]+Ax[Pc0,Pc2]+Ax[Pc0,Pc3]+Ax[Pc0,Pc4]+Ax[Pc0,Pc5]+Ax[Pc0,Pc6]+Ax[Pc0,Pc7]+Ax[Pc0,Pc8], y0''[t]==Ay[Pc0,Pc1]+Ay[Pc0,Pc2]+Ay[Pc0,Pc3]+Ay[Pc0,Pc4]+Ay[Pc0,Pc5]+Ay[Pc0,Pc6]+Ay[Pc0,Pc7]+Ay[Pc0,Pc8], z0''[t]==Az[Pc0,Pc1]+Az[Pc0,Pc2]+Az[Pc0,Pc3]+Az[Pc0,Pc4]+Az[Pc0,Pc5]+Az[Pc0,Pc6]+Az[Pc0,Pc7]+Az[Pc0,Pc8], x1''[t]==Ax[Pc1,Pc0]+Ax[Pc1,Pc2]+Ax[Pc1,Pc3]+Ax[Pc1,Pc4]+Ax[Pc1,Pc5]+Ax[Pc1,Pc6]+Ax[Pc1,Pc7]+Ax[Pc1,Pc8], y1''[t]==Ay[Pc1,Pc0]+Ay[Pc1,Pc2]+Ay[Pc1,Pc3]+Ay[Pc1,Pc4]+Ay[Pc1,Pc5]+Ay[Pc1,Pc6]+Ay[Pc1,Pc7]+Ay[Pc1,Pc8], z1''[t]==Az[Pc1,Pc0]+Az[Pc1,Pc2]+Az[Pc1,Pc3]+Az[Pc1,Pc4]+Az[Pc1,Pc5]+Az[Pc1,Pc6]+Az[Pc1,Pc7]+Az[Pc1,Pc8], x2''[t]==Ax[Pc2,Pc0]+Ax[Pc2,Pc1]+Ax[Pc2,Pc3]+Ax[Pc2,Pc4]+Ax[Pc2,Pc5]+Ax[Pc2,Pc6]+Ax[Pc2,Pc7]+Ax[Pc2,Pc8], y2''[t]==Ay[Pc2,Pc0]+Ay[Pc2,Pc1]+Ay[Pc2,Pc3]+Ay[Pc2,Pc4]+Ay[Pc2,Pc5]+Ay[Pc2,Pc6]+Ay[Pc2,Pc7]+Ay[Pc2,Pc8],
z2''[t]==Az[Pc2,Pc0]+Az[Pc2,Pc1]+Az[Pc2,Pc3]+Az[Pc2,Pc4]+Az[Pc2,Pc5]+Az[Pc2,Pc6]+Az[Pc2,Pc7]+Az[Pc2,Pc8],
x3''[t]==Ax[Pc3,Pc0]+Ax[Pc3,Pc1]+Ax[Pc3,Pc2]+Ax[Pc3,Pc4]+Ax[Pc3,Pc5]+Ax[Pc3,Pc6]+Ax[Pc3,Pc7]+Ax[Pc3,Pc8], y3''[t]==Ay[Pc3,Pc0]+Ay[Pc3,Pc1]+Ay[Pc3,Pc2]+Ay[Pc3,Pc4]+Ay[Pc3,Pc5]+Ay[Pc3,Pc6]+Ay[Pc3,Pc7]+Ay[Pc3,Pc8], z3''[t]==Az[Pc3,Pc0]+Az[Pc3,Pc1]+Az[Pc3,Pc2]+Az[Pc3,Pc4]+Az[Pc3,Pc5]+Az[Pc3,Pc6]+Az[Pc3,Pc7]+Az[Pc3,Pc8], x4''[t]==Ax[Pc4,Pc0]+Ax[Pc4,Pc1]+Ax[Pc4,Pc2]+Ax[Pc4,Pc3]+Ax[Pc4,Pc5]+Ax[Pc4,Pc6]+Ax[Pc4,Pc7]+Ax[Pc4,Pc8], y4''[t]==Ay[Pc4,Pc0]+Ay[Pc4,Pc1]+Ay[Pc4,Pc2]+Ay[Pc4,Pc3]+Ay[Pc4,Pc5]+Ay[Pc4,Pc6]+Ay[Pc4,Pc7]+Ay[Pc4,Pc8], z4''[t]==Az[Pc4,Pc0]+Az[Pc4,Pc1]+Az[Pc4,Pc2]+Az[Pc4,Pc3]+Az[Pc4,Pc5]+Az[Pc4,Pc6]+Az[Pc4,Pc7]+Az[Pc4,Pc8], x5''[t]==Ax[Pc5,Pc0]+Ax[Pc5,Pc1]+Ax[Pc5,Pc2]+Ax[Pc5,Pc3]+Ax[Pc5,Pc4]+Ax[Pc5,Pc6]+Ax[Pc5,Pc7]+Ax[Pc5,Pc8], y5''[t]==Ay[Pc5,Pc0]+Ay[Pc5,Pc1]+Ay[Pc5,Pc2]+Ay[Pc5,Pc3]+Ay[Pc5,Pc4]+Ay[Pc5,Pc6]+Ay[Pc5,Pc7]+Ay[Pc5,Pc8],
z5''[t]==Az[Pc5,Pc0]+Az[Pc5,Pc1]+Az[Pc5,Pc2]+Az[Pc5,Pc3]+Az[Pc5,Pc4]+Az[Pc5,Pc6]+Az[Pc5,Pc7]+Az[Pc5,Pc8],
x6''[t]==Ax[Pc6,Pc0]+Ax[Pc6,Pc1]+Ax[Pc6,Pc2]+Ax[Pc6,Pc3]+Ax[Pc6,Pc4]+Ax[Pc6,Pc5]+Ax[Pc6,Pc7]+Ax[Pc6,Pc8], 
y6''[t]==Ay[Pc6,Pc0]+Ay[Pc6,Pc1]+Ay[Pc6,Pc2]+Ay[Pc6,Pc3]+Ay[Pc6,Pc4]+Ay[Pc6,Pc5]+Ay[Pc6,Pc7]+Ay[Pc6,Pc8], 
z6''[t]==Az[Pc6,Pc0]+Az[Pc6,Pc1]+Az[Pc6,Pc2]+Az[Pc6,Pc3]+Az[Pc6,Pc4]+Az[Pc6,Pc5]+Az[Pc6,Pc7]+Az[Pc6,Pc8], 
x7''[t]==Ax[Pc7,Pc0]+Ax[Pc7,Pc1]+Ax[Pc7,Pc2]+Ax[Pc7,Pc3]+Ax[Pc7,Pc4]+Ax[Pc7,Pc5]+Ax[Pc7,Pc6]+Ax[Pc7,Pc8], 
y7''[t]==Ay[Pc7,Pc0]+Ay[Pc7,Pc1]+Ay[Pc7,Pc2]+Ay[Pc7,Pc3]+Ay[Pc7,Pc4]+Ay[Pc7,Pc5]+Ay[Pc7,Pc6]+Ay[Pc7,Pc8], 
z7''[t]==Az[Pc7,Pc0]+Az[Pc7,Pc1]+Az[Pc7,Pc2]+Az[Pc7,Pc3]+Az[Pc7,Pc4]+Az[Pc7,Pc5]+Az[Pc7,Pc6]+Az[Pc7,Pc8], 
x8''[t]==Ax[Pc8,Pc0]+Ax[Pc8,Pc1]+Ax[Pc8,Pc2]+Ax[Pc8,Pc3]+Ax[Pc8,Pc4]+Ax[Pc8,Pc5]+Ax[Pc8,Pc6]+Ax[Pc8,Pc7], 
y8''[t]==Ay[Pc8,Pc0]+Ay[Pc8,Pc1]+Ay[Pc8,Pc2]+Ay[Pc8,Pc3]+Ay[Pc8,Pc4]+Ay[Pc8,Pc5]+Ay[Pc8,Pc6]+Ay[Pc8,Pc7], 
z8''[t]==Az[Pc8,Pc0]+Az[Pc8,Pc1]+Az[Pc8,Pc2]+Az[Pc8,Pc3]+Az[Pc8,Pc4]+Az[Pc8,Pc5]+Az[Pc8,Pc6]+Az[Pc8,Pc7], 
xQ''[t]==QAx[PcQ,Pc0],yQ''[t]==QAy[PcQ,Pc0],zQ''[t]==QAz[PcQ,Pc0],
x0[0]==0,y0[0]==0,z0[0]==0,x0'[0]==0,y0'[0]==0,z0'[0]==0,
x1[0]==Px1,y1[0]==Py1,z1[0]==Pz1,x1'[0]==Vx1,y1'[0]==Vy1,z1'[0]==Vz1,
x2[0]==Px2,y2[0]==Py2,z2[0]==Pz2,x2'[0]==Vx2,y2'[0]==Vy2,z2'[0]==Vz2,
x3[0]==Px3,y3[0]==Py3,z3[0]==Pz3,x3'[0]==Vx3,y3'[0]==Vy3,z3'[0]==Vz3,
x4[0]==Px4,y4[0]==Py4,z4[0]==Pz4,x4'[0]==Vx4,y4'[0]==Vy4,z4'[0]==Vz4,
x5[0]==Px5,y5[0]==Py5,z5[0]==Pz5,x5'[0]==Vx5,y5'[0]==Vy5,z5'[0]==Vz5,
x6[0]==Px6,y6[0]==Py6,z6[0]==Pz6,x6'[0]==Vx6,y6'[0]==Vy6,z6'[0]==Vz6,
x7[0]==Px7,y7[0]==Py7,z7[0]==Pz7,x7'[0]==Vx7,y7'[0]==Vy7,z7'[0]==Vz7,
x8[0]==Px8,y8[0]==Py8,z8[0]==Pz8,x8'[0]==Vx8,y8'[0]==Vy8,z8'[0]==Vz8,xQ[0]==RQ,yQ[0]==0,zQ[0]==0,xQ'[0]==0,yQ'[0]==VQ,zQ'[0]==0},{x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8,xQ,yQ,zQ},{t,0,Per},MaxSteps->2000000,AccuracyGoal->10];

(*After running the above simulation, do the following: 
1. Run the function VELANG3
2. Run the Perihelion Finder Script 2
3. Plot the result using ListPlot *)

ListPlot[PeriAng,Joined->True]

(*Get line of best fit - in arcseconds per 100 years*)Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Simulate Mercury orbiting a floating Sun and all outer planets from Venus to Neptune with eccentric orbits.  The orbits will be rotated and inclined to match the actual planets' orbits.  This is a 3D situation.*){M0,M1,M2,M3,M4,M5,M6,M7,M8}={Msun,Mmercury,Mvenus,Mearth,Mmars,Mjupiter,Msaturn,Muranus,Mneptune};
Pc0={M0,x0[t],y0[t],z0[t]};
Pc1={M1,x1[t],y1[t],z1[t]};
Pc2={M2,x2[t],y2[t],z2[t]};
Pc3={M3,x3[t],y3[t],z3[t]};
Pc4={M4,x4[t],y4[t],z4[t]};
Pc5={M5,x5[t],y5[t],z5[t]};
Pc6={M6,x6[t],y6[t],z6[t]};
Pc7={M7,x7[t],y7[t],z7[t]};
Pc8={M8,x8[t],y8[t],z8[t]};
PcQ={M1,xQ[t],yQ[t],zQ[t]};
G=Gg;{R1,R2,R3,R4,R5,R6,R7,R8}={RPmercury,RPvenus,RPearth,RPmars,RPjupiter,RPsaturn,RPuranus,RPneptune};velint=Array[0&,8];
For[i=1,i≤8,i++,velint[[i]]=PveladjP[[i]](Sqrt[(1+ Peccen[[i]])((G Msun)/(Pperihel[[i]]))])];{V1,V2,V3,V4,V5,V6,V7,V8}=velint;
{Px1,Py1,Pz1}=OrbTransXYZdeg[R1,0,0,ιmercury,Ωmercury,ωmercury];{Vx1,Vy1,Vz1}=OrbTransXYZdeg[0,V1,0,ιmercury,Ωmercury,ωmercury];
{Px2,Py2,Pz2}=OrbTransXYZdeg[R2,0,0,ιvenus,Ωvenus,ωvenus];{Vx2,Vy2,Vz2}=OrbTransXYZdeg[0,V2,0,ιvenus,Ωvenus,ωvenus];
{Px3,Py3,Pz3}=OrbTransXYZdeg[R3,0,0,ιearth,Ωearth,ωearth];{Vx3,Vy3,Vz3}=OrbTransXYZdeg[0,V3,0,ιearth,Ωearth,ωearth];
{Px4,Py4,Pz4}=OrbTransXYZdeg[R4,0,0,ιmars,Ωmars,ωmars];
{Vx4,Vy4,Vz4}=OrbTransXYZdeg[0,V4,0,ιmars,Ωmars,ωmars];
{Px5,Py5,Pz5}=OrbTransXYZdeg[R5,0,0,ιjupiter,Ωjupiter,ωjupiter];{Vx5,Vy5,Vz5}=OrbTransXYZdeg[0,V5,0,ιjupiter,Ωjupiter,ωjupiter];
{Px6,Py6,Pz6}=OrbTransXYZdeg[R6,0,0,ιsaturn,Ωsaturn,ωsaturn];{Vx6,Vy6,Vz6}=OrbTransXYZdeg[0,V6,0,ιsaturn,Ωsaturn,ωsaturn];
{Px7,Py7,Pz7}=OrbTransXYZdeg[R7,0,0,ιuranus,Ωuranus,ωuranus];{Vx7,Vy7,Vz7}=OrbTransXYZdeg[0,V7,0,ιuranus,Ωuranus,ωuranus];
{Px8,Py8,Pz8}=OrbTransXYZdeg[R8,0,0,ιneptune,Ωneptune,ωneptune];{Vx8,Vy8,Vz8}=OrbTransXYZdeg[0,V8,0,ιneptune,Ωneptune,ωneptune];
RQ=RPmercury*0.45;VQ=(Sqrt[((G M0)/RQ)]);
{PxQ,PyQ,PzQ}={RQ,0,0};{VxQ,VyQ,VzQ}={0,VQ,0};
Per=200 OneYear;result=NDSolve[{
 x0''[t]==Ax[Pc0,Pc1]+Ax[Pc0,Pc2]+Ax[Pc0,Pc3]+Ax[Pc0,Pc4]+Ax[Pc0,Pc5]+Ax[Pc0,Pc6]+Ax[Pc0,Pc7]+Ax[Pc0,Pc8], y0''[t]==Ay[Pc0,Pc1]+Ay[Pc0,Pc2]+Ay[Pc0,Pc3]+Ay[Pc0,Pc4]+Ay[Pc0,Pc5]+Ay[Pc0,Pc6]+Ay[Pc0,Pc7]+Ay[Pc0,Pc8], z0''[t]==Az[Pc0,Pc1]+Az[Pc0,Pc2]+Az[Pc0,Pc3]+Az[Pc0,Pc4]+Az[Pc0,Pc5]+Az[Pc0,Pc6]+Az[Pc0,Pc7]+Az[Pc0,Pc8], x1''[t]==Ax[Pc1,Pc0]+Ax[Pc1,Pc2]+Ax[Pc1,Pc3]+Ax[Pc1,Pc4]+Ax[Pc1,Pc5]+Ax[Pc1,Pc6]+Ax[Pc1,Pc7]+Ax[Pc1,Pc8], y1''[t]==Ay[Pc1,Pc0]+Ay[Pc1,Pc2]+Ay[Pc1,Pc3]+Ay[Pc1,Pc4]+Ay[Pc1,Pc5]+Ay[Pc1,Pc6]+Ay[Pc1,Pc7]+Ay[Pc1,Pc8], z1''[t]==Az[Pc1,Pc0]+Az[Pc1,Pc2]+Az[Pc1,Pc3]+Az[Pc1,Pc4]+Az[Pc1,Pc5]+Az[Pc1,Pc6]+Az[Pc1,Pc7]+Az[Pc1,Pc8], x2''[t]==Ax[Pc2,Pc0]+Ax[Pc2,Pc1]+Ax[Pc2,Pc3]+Ax[Pc2,Pc4]+Ax[Pc2,Pc5]+Ax[Pc2,Pc6]+Ax[Pc2,Pc7]+Ax[Pc2,Pc8], y2''[t]==Ay[Pc2,Pc0]+Ay[Pc2,Pc1]+Ay[Pc2,Pc3]+Ay[Pc2,Pc4]+Ay[Pc2,Pc5]+Ay[Pc2,Pc6]+Ay[Pc2,Pc7]+Ay[Pc2,Pc8],
z2''[t]==Az[Pc2,Pc0]+Az[Pc2,Pc1]+Az[Pc2,Pc3]+Az[Pc2,Pc4]+Az[Pc2,Pc5]+Az[Pc2,Pc6]+Az[Pc2,Pc7]+Az[Pc2,Pc8],
x3''[t]==Ax[Pc3,Pc0]+Ax[Pc3,Pc1]+Ax[Pc3,Pc2]+Ax[Pc3,Pc4]+Ax[Pc3,Pc5]+Ax[Pc3,Pc6]+Ax[Pc3,Pc7]+Ax[Pc3,Pc8], y3''[t]==Ay[Pc3,Pc0]+Ay[Pc3,Pc1]+Ay[Pc3,Pc2]+Ay[Pc3,Pc4]+Ay[Pc3,Pc5]+Ay[Pc3,Pc6]+Ay[Pc3,Pc7]+Ay[Pc3,Pc8], z3''[t]==Az[Pc3,Pc0]+Az[Pc3,Pc1]+Az[Pc3,Pc2]+Az[Pc3,Pc4]+Az[Pc3,Pc5]+Az[Pc3,Pc6]+Az[Pc3,Pc7]+Az[Pc3,Pc8], x4''[t]==Ax[Pc4,Pc0]+Ax[Pc4,Pc1]+Ax[Pc4,Pc2]+Ax[Pc4,Pc3]+Ax[Pc4,Pc5]+Ax[Pc4,Pc6]+Ax[Pc4,Pc7]+Ax[Pc4,Pc8], y4''[t]==Ay[Pc4,Pc0]+Ay[Pc4,Pc1]+Ay[Pc4,Pc2]+Ay[Pc4,Pc3]+Ay[Pc4,Pc5]+Ay[Pc4,Pc6]+Ay[Pc4,Pc7]+Ay[Pc4,Pc8], z4''[t]==Az[Pc4,Pc0]+Az[Pc4,Pc1]+Az[Pc4,Pc2]+Az[Pc4,Pc3]+Az[Pc4,Pc5]+Az[Pc4,Pc6]+Az[Pc4,Pc7]+Az[Pc4,Pc8], x5''[t]==Ax[Pc5,Pc0]+Ax[Pc5,Pc1]+Ax[Pc5,Pc2]+Ax[Pc5,Pc3]+Ax[Pc5,Pc4]+Ax[Pc5,Pc6]+Ax[Pc5,Pc7]+Ax[Pc5,Pc8], y5''[t]==Ay[Pc5,Pc0]+Ay[Pc5,Pc1]+Ay[Pc5,Pc2]+Ay[Pc5,Pc3]+Ay[Pc5,Pc4]+Ay[Pc5,Pc6]+Ay[Pc5,Pc7]+Ay[Pc5,Pc8],
z5''[t]==Az[Pc5,Pc0]+Az[Pc5,Pc1]+Az[Pc5,Pc2]+Az[Pc5,Pc3]+Az[Pc5,Pc4]+Az[Pc5,Pc6]+Az[Pc5,Pc7]+Az[Pc5,Pc8],
x6''[t]==Ax[Pc6,Pc0]+Ax[Pc6,Pc1]+Ax[Pc6,Pc2]+Ax[Pc6,Pc3]+Ax[Pc6,Pc4]+Ax[Pc6,Pc5]+Ax[Pc6,Pc7]+Ax[Pc6,Pc8], 
y6''[t]==Ay[Pc6,Pc0]+Ay[Pc6,Pc1]+Ay[Pc6,Pc2]+Ay[Pc6,Pc3]+Ay[Pc6,Pc4]+Ay[Pc6,Pc5]+Ay[Pc6,Pc7]+Ay[Pc6,Pc8], 
z6''[t]==Az[Pc6,Pc0]+Az[Pc6,Pc1]+Az[Pc6,Pc2]+Az[Pc6,Pc3]+Az[Pc6,Pc4]+Az[Pc6,Pc5]+Az[Pc6,Pc7]+Az[Pc6,Pc8], 
x7''[t]==Ax[Pc7,Pc0]+Ax[Pc7,Pc1]+Ax[Pc7,Pc2]+Ax[Pc7,Pc3]+Ax[Pc7,Pc4]+Ax[Pc7,Pc5]+Ax[Pc7,Pc6]+Ax[Pc7,Pc8], 
y7''[t]==Ay[Pc7,Pc0]+Ay[Pc7,Pc1]+Ay[Pc7,Pc2]+Ay[Pc7,Pc3]+Ay[Pc7,Pc4]+Ay[Pc7,Pc5]+Ay[Pc7,Pc6]+Ay[Pc7,Pc8], 
z7''[t]==Az[Pc7,Pc0]+Az[Pc7,Pc1]+Az[Pc7,Pc2]+Az[Pc7,Pc3]+Az[Pc7,Pc4]+Az[Pc7,Pc5]+Az[Pc7,Pc6]+Az[Pc7,Pc8], 
x8''[t]==Ax[Pc8,Pc0]+Ax[Pc8,Pc1]+Ax[Pc8,Pc2]+Ax[Pc8,Pc3]+Ax[Pc8,Pc4]+Ax[Pc8,Pc5]+Ax[Pc8,Pc6]+Ax[Pc8,Pc7], 
y8''[t]==Ay[Pc8,Pc0]+Ay[Pc8,Pc1]+Ay[Pc8,Pc2]+Ay[Pc8,Pc3]+Ay[Pc8,Pc4]+Ay[Pc8,Pc5]+Ay[Pc8,Pc6]+Ay[Pc8,Pc7], 
z8''[t]==Az[Pc8,Pc0]+Az[Pc8,Pc1]+Az[Pc8,Pc2]+Az[Pc8,Pc3]+Az[Pc8,Pc4]+Az[Pc8,Pc5]+Az[Pc8,Pc6]+Az[Pc8,Pc7], 
xQ''[t]==QAx[PcQ,Pc0],yQ''[t]==QAy[PcQ,Pc0],zQ''[t]==QAz[PcQ,Pc0],
x0[0]==0,y0[0]==0,z0[0]==0,x0'[0]==0,y0'[0]==0,z0'[0]==0,
x1[0]==Px1,y1[0]==Py1,z1[0]==Pz1,x1'[0]==Vx1,y1'[0]==Vy1,z1'[0]==Vz1,
x2[0]==Px2,y2[0]==Py2,z2[0]==Pz2,x2'[0]==Vx2,y2'[0]==Vy2,z2'[0]==Vz2,
x3[0]==Px3,y3[0]==Py3,z3[0]==Pz3,x3'[0]==Vx3,y3'[0]==Vy3,z3'[0]==Vz3,
x4[0]==Px4,y4[0]==Py4,z4[0]==Pz4,x4'[0]==Vx4,y4'[0]==Vy4,z4'[0]==Vz4,
x5[0]==Px5,y5[0]==Py5,z5[0]==Pz5,x5'[0]==Vx5,y5'[0]==Vy5,z5'[0]==Vz5,
x6[0]==Px6,y6[0]==Py6,z6[0]==Pz6,x6'[0]==Vx6,y6'[0]==Vy6,z6'[0]==Vz6,
x7[0]==Px7,y7[0]==Py7,z7[0]==Pz7,x7'[0]==Vx7,y7'[0]==Vy7,z7'[0]==Vz7,
x8[0]==Px8,y8[0]==Py8,z8[0]==Pz8,x8'[0]==Vx8,y8'[0]==Vy8,z8'[0]==Vz8,xQ[0]==PxQ,yQ[0]==PyQ,zQ[0]==PzQ,xQ'[0]==VxQ,yQ'[0]==VyQ,zQ'[0]==VzQ},
{x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8,xQ,yQ,zQ},{t,0,Per},MaxSteps->2000000,AccuracyGoal->10];

(*After running the above simulation, do the following: 
1. Run the function VELANG5
2. Run the Perihelion Finder Script 2
3. Plot the result using ListPlot *)

ListPlot[PeriAng,Joined->True]

(*Get line of best fit - in arcseconds per 100 years*)Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Simulate Mercury orbiting a floating Sun and all outer planets from Venus to Neptune with eccentric orbits.  The orbits will be rotated and inclined to match the actual planets' orbits.  This is the same as above but this time the planets start at their aphelion position.*)
{M0,M1,M2,M3,M4,M5,M6,M7,M8}={Msun,Mmercury,Mvenus,Mearth,Mmars,Mjupiter,Msaturn,Muranus,Mneptune};Pc0={M0,x0[t],y0[t],z0[t]};
Pc1={M1,x1[t],y1[t],z1[t]};
Pc2={M2,x2[t],y2[t],z2[t]};
Pc3={M3,x3[t],y3[t],z3[t]};
Pc4={M4,x4[t],y4[t],z4[t]};
Pc5={M5,x5[t],y5[t],z5[t]};
Pc6={M6,x6[t],y6[t],z6[t]};
Pc7={M7,x7[t],y7[t],z7[t]};
Pc8={M8,x8[t],y8[t],z8[t]};
PcQ={M1,xQ[t],yQ[t],zQ[t]};
G=Gg;posint=Array[0&,8];
For[i=1,i≤8,i++,posint[[i]]=((1+Peccen[[i]])/(1-Peccen[[i]]))Pperihel[[i]]];{R1,R2,R3,R4,R5,R6,R7,R8}=posint;velint=Array[0&,8];
For[i=1,i≤8,i++,velint[[i]]=PveladjA[[i]](Sqrt[(1- Peccen[[i]])((G Msun)/(posint[[i]]))])];{V1,V2,V3,V4,V5,V6,V7,V8}=velint;R1=RPmercury;V1=PveladjP[[1]](Sqrt[(1+ Emercury)((G Msun)/R1)]);
{Px1,Py1,Pz1}=OrbTransXYZdeg[R1,0,0,ιmercury,Ωmercury,ωmercury];{Vx1,Vy1,Vz1}=OrbTransXYZdeg[0,V1,0,ιmercury,Ωmercury,ωmercury];
{Px2,Py2,Pz2}=OrbTransXYZdeg[R2,0,0,ιvenus,Ωvenus,ωvenus+180];{Vx2,Vy2,Vz2}=OrbTransXYZdeg[0,V2,0,ιvenus,Ωvenus,ωvenus+180];
{Px3,Py3,Pz3}=OrbTransXYZdeg[R3,0,0,ιearth,Ωearth,ωearth+180];{Vx3,Vy3,Vz3}=OrbTransXYZdeg[0,V3,0,ιearth,Ωearth,ωearth+180];
{Px4,Py4,Pz4}=OrbTransXYZdeg[R4,0,0,ιmars,Ωmars,ωmars+180];
{Vx4,Vy4,Vz4}=OrbTransXYZdeg[0,V4,0,ιmars,Ωmars,ωmars+180];
{Px5,Py5,Pz5}=OrbTransXYZdeg[R5,0,0,ιjupiter,Ωjupiter,ωjupiter+180];{Vx5,Vy5,Vz5}=OrbTransXYZdeg[0,V5,0,ιjupiter,Ωjupiter,ωjupiter+180];
{Px6,Py6,Pz6}=OrbTransXYZdeg[R6,0,0,ιsaturn,Ωsaturn,ωsaturn+180];{Vx6,Vy6,Vz6}=OrbTransXYZdeg[0,V6,0,ιsaturn,Ωsaturn,ωsaturn+180];
{Px7,Py7,Pz7}=OrbTransXYZdeg[R7,0,0,ιuranus,Ωuranus,ωuranus+180];{Vx7,Vy7,Vz7}=OrbTransXYZdeg[0,V7,0,ιuranus,Ωuranus,ωuranus+180];
{Px8,Py8,Pz8}=OrbTransXYZdeg[R8,0,0,ιneptune,Ωneptune,ωneptune+180];{Vx8,Vy8,Vz8}=OrbTransXYZdeg[0,V8,0,ιneptune,Ωneptune,ωneptune+180];
RQ=RPmercury*0.45;VQ=(Sqrt[((G M0)/RQ)]);
{PxQ,PyQ,PzQ}={RQ,0,0};{VxQ,VyQ,VzQ}={0,VQ,0};
Per=200 OneYear;result=NDSolve[{
 x0''[t]==Ax[Pc0,Pc1]+Ax[Pc0,Pc2]+Ax[Pc0,Pc3]+Ax[Pc0,Pc4]+Ax[Pc0,Pc5]+Ax[Pc0,Pc6]+Ax[Pc0,Pc7]+Ax[Pc0,Pc8], y0''[t]==Ay[Pc0,Pc1]+Ay[Pc0,Pc2]+Ay[Pc0,Pc3]+Ay[Pc0,Pc4]+Ay[Pc0,Pc5]+Ay[Pc0,Pc6]+Ay[Pc0,Pc7]+Ay[Pc0,Pc8], z0''[t]==Az[Pc0,Pc1]+Az[Pc0,Pc2]+Az[Pc0,Pc3]+Az[Pc0,Pc4]+Az[Pc0,Pc5]+Az[Pc0,Pc6]+Az[Pc0,Pc7]+Az[Pc0,Pc8], x1''[t]==Ax[Pc1,Pc0]+Ax[Pc1,Pc2]+Ax[Pc1,Pc3]+Ax[Pc1,Pc4]+Ax[Pc1,Pc5]+Ax[Pc1,Pc6]+Ax[Pc1,Pc7]+Ax[Pc1,Pc8], y1''[t]==Ay[Pc1,Pc0]+Ay[Pc1,Pc2]+Ay[Pc1,Pc3]+Ay[Pc1,Pc4]+Ay[Pc1,Pc5]+Ay[Pc1,Pc6]+Ay[Pc1,Pc7]+Ay[Pc1,Pc8], z1''[t]==Az[Pc1,Pc0]+Az[Pc1,Pc2]+Az[Pc1,Pc3]+Az[Pc1,Pc4]+Az[Pc1,Pc5]+Az[Pc1,Pc6]+Az[Pc1,Pc7]+Az[Pc1,Pc8], x2''[t]==Ax[Pc2,Pc0]+Ax[Pc2,Pc1]+Ax[Pc2,Pc3]+Ax[Pc2,Pc4]+Ax[Pc2,Pc5]+Ax[Pc2,Pc6]+Ax[Pc2,Pc7]+Ax[Pc2,Pc8], y2''[t]==Ay[Pc2,Pc0]+Ay[Pc2,Pc1]+Ay[Pc2,Pc3]+Ay[Pc2,Pc4]+Ay[Pc2,Pc5]+Ay[Pc2,Pc6]+Ay[Pc2,Pc7]+Ay[Pc2,Pc8],
z2''[t]==Az[Pc2,Pc0]+Az[Pc2,Pc1]+Az[Pc2,Pc3]+Az[Pc2,Pc4]+Az[Pc2,Pc5]+Az[Pc2,Pc6]+Az[Pc2,Pc7]+Az[Pc2,Pc8],
x3''[t]==Ax[Pc3,Pc0]+Ax[Pc3,Pc1]+Ax[Pc3,Pc2]+Ax[Pc3,Pc4]+Ax[Pc3,Pc5]+Ax[Pc3,Pc6]+Ax[Pc3,Pc7]+Ax[Pc3,Pc8], y3''[t]==Ay[Pc3,Pc0]+Ay[Pc3,Pc1]+Ay[Pc3,Pc2]+Ay[Pc3,Pc4]+Ay[Pc3,Pc5]+Ay[Pc3,Pc6]+Ay[Pc3,Pc7]+Ay[Pc3,Pc8], z3''[t]==Az[Pc3,Pc0]+Az[Pc3,Pc1]+Az[Pc3,Pc2]+Az[Pc3,Pc4]+Az[Pc3,Pc5]+Az[Pc3,Pc6]+Az[Pc3,Pc7]+Az[Pc3,Pc8], x4''[t]==Ax[Pc4,Pc0]+Ax[Pc4,Pc1]+Ax[Pc4,Pc2]+Ax[Pc4,Pc3]+Ax[Pc4,Pc5]+Ax[Pc4,Pc6]+Ax[Pc4,Pc7]+Ax[Pc4,Pc8], y4''[t]==Ay[Pc4,Pc0]+Ay[Pc4,Pc1]+Ay[Pc4,Pc2]+Ay[Pc4,Pc3]+Ay[Pc4,Pc5]+Ay[Pc4,Pc6]+Ay[Pc4,Pc7]+Ay[Pc4,Pc8], z4''[t]==Az[Pc4,Pc0]+Az[Pc4,Pc1]+Az[Pc4,Pc2]+Az[Pc4,Pc3]+Az[Pc4,Pc5]+Az[Pc4,Pc6]+Az[Pc4,Pc7]+Az[Pc4,Pc8], x5''[t]==Ax[Pc5,Pc0]+Ax[Pc5,Pc1]+Ax[Pc5,Pc2]+Ax[Pc5,Pc3]+Ax[Pc5,Pc4]+Ax[Pc5,Pc6]+Ax[Pc5,Pc7]+Ax[Pc5,Pc8], y5''[t]==Ay[Pc5,Pc0]+Ay[Pc5,Pc1]+Ay[Pc5,Pc2]+Ay[Pc5,Pc3]+Ay[Pc5,Pc4]+Ay[Pc5,Pc6]+Ay[Pc5,Pc7]+Ay[Pc5,Pc8],
z5''[t]==Az[Pc5,Pc0]+Az[Pc5,Pc1]+Az[Pc5,Pc2]+Az[Pc5,Pc3]+Az[Pc5,Pc4]+Az[Pc5,Pc6]+Az[Pc5,Pc7]+Az[Pc5,Pc8],
x6''[t]==Ax[Pc6,Pc0]+Ax[Pc6,Pc1]+Ax[Pc6,Pc2]+Ax[Pc6,Pc3]+Ax[Pc6,Pc4]+Ax[Pc6,Pc5]+Ax[Pc6,Pc7]+Ax[Pc6,Pc8], 
y6''[t]==Ay[Pc6,Pc0]+Ay[Pc6,Pc1]+Ay[Pc6,Pc2]+Ay[Pc6,Pc3]+Ay[Pc6,Pc4]+Ay[Pc6,Pc5]+Ay[Pc6,Pc7]+Ay[Pc6,Pc8], 
z6''[t]==Az[Pc6,Pc0]+Az[Pc6,Pc1]+Az[Pc6,Pc2]+Az[Pc6,Pc3]+Az[Pc6,Pc4]+Az[Pc6,Pc5]+Az[Pc6,Pc7]+Az[Pc6,Pc8], 
x7''[t]==Ax[Pc7,Pc0]+Ax[Pc7,Pc1]+Ax[Pc7,Pc2]+Ax[Pc7,Pc3]+Ax[Pc7,Pc4]+Ax[Pc7,Pc5]+Ax[Pc7,Pc6]+Ax[Pc7,Pc8], 
y7''[t]==Ay[Pc7,Pc0]+Ay[Pc7,Pc1]+Ay[Pc7,Pc2]+Ay[Pc7,Pc3]+Ay[Pc7,Pc4]+Ay[Pc7,Pc5]+Ay[Pc7,Pc6]+Ay[Pc7,Pc8], 
z7''[t]==Az[Pc7,Pc0]+Az[Pc7,Pc1]+Az[Pc7,Pc2]+Az[Pc7,Pc3]+Az[Pc7,Pc4]+Az[Pc7,Pc5]+Az[Pc7,Pc6]+Az[Pc7,Pc8], 
x8''[t]==Ax[Pc8,Pc0]+Ax[Pc8,Pc1]+Ax[Pc8,Pc2]+Ax[Pc8,Pc3]+Ax[Pc8,Pc4]+Ax[Pc8,Pc5]+Ax[Pc8,Pc6]+Ax[Pc8,Pc7], 
y8''[t]==Ay[Pc8,Pc0]+Ay[Pc8,Pc1]+Ay[Pc8,Pc2]+Ay[Pc8,Pc3]+Ay[Pc8,Pc4]+Ay[Pc8,Pc5]+Ay[Pc8,Pc6]+Ay[Pc8,Pc7], 
z8''[t]==Az[Pc8,Pc0]+Az[Pc8,Pc1]+Az[Pc8,Pc2]+Az[Pc8,Pc3]+Az[Pc8,Pc4]+Az[Pc8,Pc5]+Az[Pc8,Pc6]+Az[Pc8,Pc7], 
xQ''[t]==QAx[PcQ,Pc0],yQ''[t]==QAy[PcQ,Pc0],zQ''[t]==QAz[PcQ,Pc0],
x0[0]==0,y0[0]==0,z0[0]==0,x0'[0]==0,y0'[0]==0,z0'[0]==0,
x1[0]==Px1,y1[0]==Py1,z1[0]==Pz1,x1'[0]==Vx1,y1'[0]==Vy1,z1'[0]==Vz1,
x2[0]==Px2,y2[0]==Py2,z2[0]==Pz2,x2'[0]==Vx2,y2'[0]==Vy2,z2'[0]==Vz2,
x3[0]==Px3,y3[0]==Py3,z3[0]==Pz3,x3'[0]==Vx3,y3'[0]==Vy3,z3'[0]==Vz3,
x4[0]==Px4,y4[0]==Py4,z4[0]==Pz4,x4'[0]==Vx4,y4'[0]==Vy4,z4'[0]==Vz4,
x5[0]==Px5,y5[0]==Py5,z5[0]==Pz5,x5'[0]==Vx5,y5'[0]==Vy5,z5'[0]==Vz5,
x6[0]==Px6,y6[0]==Py6,z6[0]==Pz6,x6'[0]==Vx6,y6'[0]==Vy6,z6'[0]==Vz6,
x7[0]==Px7,y7[0]==Py7,z7[0]==Pz7,x7'[0]==Vx7,y7'[0]==Vy7,z7'[0]==Vz7,
x8[0]==Px8,y8[0]==Py8,z8[0]==Pz8,x8'[0]==Vx8,y8'[0]==Vy8,z8'[0]==Vz8,xQ[0]==PxQ,yQ[0]==PyQ,zQ[0]==PzQ,xQ'[0]==VxQ,yQ'[0]==VyQ,zQ'[0]==VzQ},
{x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8,xQ,yQ,zQ},{t,0,Per},MaxSteps->2000000,AccuracyGoal->10];

(*After running the above simulation, do the following: 
1. Run the function VELANG5
2. Run the Perihelion Finder Script 2
3. Plot the result using ListPlot *)

ListPlot[PeriAng,Joined->True]

(*Get line of best fit - in arcseconds per 100 years*)Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Simulate solar system from Mercury to Neptune (and Sun) using HORIZONS data.  Data positions and velocities are relative to SSB at 15:14 on 2001-Feb-01*)
{M0,M1,M2,M3,M4,M5,M6,M7,M8}={Msun,Mmercury,Mvenus,Mearth,Mmars,Mjupiter,Msaturn,Muranus,Mneptune};Pc0={M0,x0[t],y0[t],z0[t]};
Pc1={M1,x1[t],y1[t],z1[t]};
Pc2={M2,x2[t],y2[t],z2[t]};
Pc3={M3,x3[t],y3[t],z3[t]};
Pc4={M4,x4[t],y4[t],z4[t]};
Pc5={M5,x5[t],y5[t],z5[t]};
Pc6={M6,x6[t],y6[t],z6[t]};
Pc7={M7,x7[t],y7[t],z7[t]};
Pc8={M8,x8[t],y8[t],z8[t]};
PcQ={M1,xQ[t],yQ[t],zQ[t]};
G=Gg;
{Px0,Py0,Pz0}={-6.544700998788778*^8,-7.542720159167981*^8,2.232295004897844*^7};
{Vx0,Vy0,Vz0}={1.430709291641341*^1,-6.129689956131971*^0,-3.231469219383437*^-1};
{Px1,Py1,Pz1}={9.459704017632110*^9,4.403743225392255*^10,2.753061478911016*^9};
{Vx1,Vy1,Vz1}={-5.726823818034963*^4,1.254517593507960*^4,6.282561277101745*^3};
{Px2,Py2,Pz2}={-1.521787736620838*^10,1.058313140886273*^11,2.319715498325974*^9};
{Vx2,Vy2,Vz2}={-3.480265957291748*^4,-4.938326945425076*^3,1.941907289679378*^3};
{Px3,Py3,Pz3}={-1.008728284577093*^11,1.073649125740838*^11,2.245984644959867*^7};
{Vx3,Vy3,Vz3}={-2.230964007501381*^4,-2.037499756826389*^4,-8.625517806555294*^-1};
{Px4,Py4,Pz4}={-2.341553529718963*^11,-6.893016826997817*^10,4.333407438725121*^9};
{Vx4,Vy4,Vz4}={7.714702908531727*^3,-2.119711899205197*^4,-6.334684794491539*^2};
{Px5,Py5,Pz5}={2.339841047987117*^11,7.181908392479482*^11,-8.215593369233102*^9};
{Vx5,Vy5,Vz5}={-1.258371434241999*^4,4.669884216770746*^3,2.622933336930586*^2};
{Px6,Py6,Pz6}={6.761216366578872*^11,1.180910647882344*^12,-4.745229184242111*^10};
{Vx6,Vy6,Vz6}={-8.889132514137824*^3,4.778912397853901*^3,2.702901855284778*^2};
{Px7,Py7,Pz7}={2.310818987293303*^12,-1.890798390154246*^12,-3.696547163527751*^10};
{Vx7,Vy7,Vz7}={4.262398872465041*^3,4.953701195921087*^3,-3.660014040687543*^1};
{Px8,Py8,Pz8}={2.665340229040296*^12,-3.630837682718371*^12,1.334489177487350*^10};
{Vx8,Vy8,Vz8}={4.346847466026690*^3,3.248046287448822*^3,-1.671845424812615*^2};
RQ=RPmercury*0.45;VQ=(Sqrt[((G M0)/RQ)]);
{PxQ,PyQ,PzQ}={RQ,0,0};{VxQ,VyQ,VzQ}={0,VQ,0};
Per=100 OneYear;result=NDSolve[{
 x0''[t]==Ax[Pc0,Pc1]+Ax[Pc0,Pc2]+Ax[Pc0,Pc3]+Ax[Pc0,Pc4]+Ax[Pc0,Pc5]+Ax[Pc0,Pc6]+Ax[Pc0,Pc7]+Ax[Pc0,Pc8], y0''[t]==Ay[Pc0,Pc1]+Ay[Pc0,Pc2]+Ay[Pc0,Pc3]+Ay[Pc0,Pc4]+Ay[Pc0,Pc5]+Ay[Pc0,Pc6]+Ay[Pc0,Pc7]+Ay[Pc0,Pc8], z0''[t]==Az[Pc0,Pc1]+Az[Pc0,Pc2]+Az[Pc0,Pc3]+Az[Pc0,Pc4]+Az[Pc0,Pc5]+Az[Pc0,Pc6]+Az[Pc0,Pc7]+Az[Pc0,Pc8], x1''[t]==Ax[Pc1,Pc0]+Ax[Pc1,Pc2]+Ax[Pc1,Pc3]+Ax[Pc1,Pc4]+Ax[Pc1,Pc5]+Ax[Pc1,Pc6]+Ax[Pc1,Pc7]+Ax[Pc1,Pc8], y1''[t]==Ay[Pc1,Pc0]+Ay[Pc1,Pc2]+Ay[Pc1,Pc3]+Ay[Pc1,Pc4]+Ay[Pc1,Pc5]+Ay[Pc1,Pc6]+Ay[Pc1,Pc7]+Ay[Pc1,Pc8], z1''[t]==Az[Pc1,Pc0]+Az[Pc1,Pc2]+Az[Pc1,Pc3]+Az[Pc1,Pc4]+Az[Pc1,Pc5]+Az[Pc1,Pc6]+Az[Pc1,Pc7]+Az[Pc1,Pc8], x2''[t]==Ax[Pc2,Pc0]+Ax[Pc2,Pc1]+Ax[Pc2,Pc3]+Ax[Pc2,Pc4]+Ax[Pc2,Pc5]+Ax[Pc2,Pc6]+Ax[Pc2,Pc7]+Ax[Pc2,Pc8], y2''[t]==Ay[Pc2,Pc0]+Ay[Pc2,Pc1]+Ay[Pc2,Pc3]+Ay[Pc2,Pc4]+Ay[Pc2,Pc5]+Ay[Pc2,Pc6]+Ay[Pc2,Pc7]+Ay[Pc2,Pc8],
z2''[t]==Az[Pc2,Pc0]+Az[Pc2,Pc1]+Az[Pc2,Pc3]+Az[Pc2,Pc4]+Az[Pc2,Pc5]+Az[Pc2,Pc6]+Az[Pc2,Pc7]+Az[Pc2,Pc8],
x3''[t]==Ax[Pc3,Pc0]+Ax[Pc3,Pc1]+Ax[Pc3,Pc2]+Ax[Pc3,Pc4]+Ax[Pc3,Pc5]+Ax[Pc3,Pc6]+Ax[Pc3,Pc7]+Ax[Pc3,Pc8], y3''[t]==Ay[Pc3,Pc0]+Ay[Pc3,Pc1]+Ay[Pc3,Pc2]+Ay[Pc3,Pc4]+Ay[Pc3,Pc5]+Ay[Pc3,Pc6]+Ay[Pc3,Pc7]+Ay[Pc3,Pc8], z3''[t]==Az[Pc3,Pc0]+Az[Pc3,Pc1]+Az[Pc3,Pc2]+Az[Pc3,Pc4]+Az[Pc3,Pc5]+Az[Pc3,Pc6]+Az[Pc3,Pc7]+Az[Pc3,Pc8], x4''[t]==Ax[Pc4,Pc0]+Ax[Pc4,Pc1]+Ax[Pc4,Pc2]+Ax[Pc4,Pc3]+Ax[Pc4,Pc5]+Ax[Pc4,Pc6]+Ax[Pc4,Pc7]+Ax[Pc4,Pc8], y4''[t]==Ay[Pc4,Pc0]+Ay[Pc4,Pc1]+Ay[Pc4,Pc2]+Ay[Pc4,Pc3]+Ay[Pc4,Pc5]+Ay[Pc4,Pc6]+Ay[Pc4,Pc7]+Ay[Pc4,Pc8], z4''[t]==Az[Pc4,Pc0]+Az[Pc4,Pc1]+Az[Pc4,Pc2]+Az[Pc4,Pc3]+Az[Pc4,Pc5]+Az[Pc4,Pc6]+Az[Pc4,Pc7]+Az[Pc4,Pc8], x5''[t]==Ax[Pc5,Pc0]+Ax[Pc5,Pc1]+Ax[Pc5,Pc2]+Ax[Pc5,Pc3]+Ax[Pc5,Pc4]+Ax[Pc5,Pc6]+Ax[Pc5,Pc7]+Ax[Pc5,Pc8], y5''[t]==Ay[Pc5,Pc0]+Ay[Pc5,Pc1]+Ay[Pc5,Pc2]+Ay[Pc5,Pc3]+Ay[Pc5,Pc4]+Ay[Pc5,Pc6]+Ay[Pc5,Pc7]+Ay[Pc5,Pc8],
z5''[t]==Az[Pc5,Pc0]+Az[Pc5,Pc1]+Az[Pc5,Pc2]+Az[Pc5,Pc3]+Az[Pc5,Pc4]+Az[Pc5,Pc6]+Az[Pc5,Pc7]+Az[Pc5,Pc8],
x6''[t]==Ax[Pc6,Pc0]+Ax[Pc6,Pc1]+Ax[Pc6,Pc2]+Ax[Pc6,Pc3]+Ax[Pc6,Pc4]+Ax[Pc6,Pc5]+Ax[Pc6,Pc7]+Ax[Pc6,Pc8], 
y6''[t]==Ay[Pc6,Pc0]+Ay[Pc6,Pc1]+Ay[Pc6,Pc2]+Ay[Pc6,Pc3]+Ay[Pc6,Pc4]+Ay[Pc6,Pc5]+Ay[Pc6,Pc7]+Ay[Pc6,Pc8], 
z6''[t]==Az[Pc6,Pc0]+Az[Pc6,Pc1]+Az[Pc6,Pc2]+Az[Pc6,Pc3]+Az[Pc6,Pc4]+Az[Pc6,Pc5]+Az[Pc6,Pc7]+Az[Pc6,Pc8], 
x7''[t]==Ax[Pc7,Pc0]+Ax[Pc7,Pc1]+Ax[Pc7,Pc2]+Ax[Pc7,Pc3]+Ax[Pc7,Pc4]+Ax[Pc7,Pc5]+Ax[Pc7,Pc6]+Ax[Pc7,Pc8], 
y7''[t]==Ay[Pc7,Pc0]+Ay[Pc7,Pc1]+Ay[Pc7,Pc2]+Ay[Pc7,Pc3]+Ay[Pc7,Pc4]+Ay[Pc7,Pc5]+Ay[Pc7,Pc6]+Ay[Pc7,Pc8], 
z7''[t]==Az[Pc7,Pc0]+Az[Pc7,Pc1]+Az[Pc7,Pc2]+Az[Pc7,Pc3]+Az[Pc7,Pc4]+Az[Pc7,Pc5]+Az[Pc7,Pc6]+Az[Pc7,Pc8], 
x8''[t]==Ax[Pc8,Pc0]+Ax[Pc8,Pc1]+Ax[Pc8,Pc2]+Ax[Pc8,Pc3]+Ax[Pc8,Pc4]+Ax[Pc8,Pc5]+Ax[Pc8,Pc6]+Ax[Pc8,Pc7], 
y8''[t]==Ay[Pc8,Pc0]+Ay[Pc8,Pc1]+Ay[Pc8,Pc2]+Ay[Pc8,Pc3]+Ay[Pc8,Pc4]+Ay[Pc8,Pc5]+Ay[Pc8,Pc6]+Ay[Pc8,Pc7], 
z8''[t]==Az[Pc8,Pc0]+Az[Pc8,Pc1]+Az[Pc8,Pc2]+Az[Pc8,Pc3]+Az[Pc8,Pc4]+Az[Pc8,Pc5]+Az[Pc8,Pc6]+Az[Pc8,Pc7], 
xQ''[t]==QAx[PcQ,Pc0],yQ''[t]==QAy[PcQ,Pc0],zQ''[t]==QAz[PcQ,Pc0],
x0[0]==Px0,y0[0]==Py0,z0[0]==Pz0,x0'[0]==Vx0,y0'[0]==Vy0,z0'[0]==Vz0,
x1[0]==Px1,y1[0]==Py1,z1[0]==Pz1,x1'[0]==Vx1,y1'[0]==Vy1,z1'[0]==Vz1,
x2[0]==Px2,y2[0]==Py2,z2[0]==Pz2,x2'[0]==Vx2,y2'[0]==Vy2,z2'[0]==Vz2,
x3[0]==Px3,y3[0]==Py3,z3[0]==Pz3,x3'[0]==Vx3,y3'[0]==Vy3,z3'[0]==Vz3,
x4[0]==Px4,y4[0]==Py4,z4[0]==Pz4,x4'[0]==Vx4,y4'[0]==Vy4,z4'[0]==Vz4,
x5[0]==Px5,y5[0]==Py5,z5[0]==Pz5,x5'[0]==Vx5,y5'[0]==Vy5,z5'[0]==Vz5,
x6[0]==Px6,y6[0]==Py6,z6[0]==Pz6,x6'[0]==Vx6,y6'[0]==Vy6,z6'[0]==Vz6,
x7[0]==Px7,y7[0]==Py7,z7[0]==Pz7,x7'[0]==Vx7,y7'[0]==Vy7,z7'[0]==Vz7,
x8[0]==Px8,y8[0]==Py8,z8[0]==Pz8,x8'[0]==Vx8,y8'[0]==Vy8,z8'[0]==Vz8,xQ[0]==PxQ,yQ[0]==PyQ,zQ[0]==PzQ,xQ'[0]==VxQ,yQ'[0]==VyQ,zQ'[0]==VzQ},
{x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8,xQ,yQ,zQ},{t,0,Per},MaxSteps->2000000,AccuracyGoal->10];

(*After running the above simulation, do the following: 
1. Run the function VELANG6
2. Run the Perihelion Finder Script 2
3. Plot the result using ListPlot *)

ListPlot[PeriAng,Joined->True]

(*Get line of best fit - in arcseconds per 100 years*)Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*
Super Accurate Simulations
*)

(* This next section will be a repeat of the above simulations and done with very high accuracy.  In order to be sure of the results it was considered necessary to simulate Mercury's position to within 10cm.  The perihelion is defined as the nearest point to the Sun and the difference between the perihelion and another point 1 arcsecond along the orbit represents the Sun-Mercury distance increasing by 10cm.
For a comparison of accuracy, Mathematica simulates a single orbit of Mercury with a position error of 15km.  Whereas this new requirement represents an accuracy of 3 microns per orbit, which translates to 8cm over 100 years or 13cm over 200 years. 
Doing this with Mathematica is going to be impossible because it requires a timestep smaller than what is allowed.  Instead the simulation will be done 'manually' via scripts.  Two methods are used to get such extreme accuracies:
1. Use a large number of derivative terms.  In this case 12 terms are used.
2.  It's also difficult to even define the derivatives in algebraic form because they get exponentially complex with each term.  
To overcome the second point, only the 'jerk' term (3rd order derivative or rate of change of acceleration) is calculated.  A short burst of jerk calculations is done very close together (time-wise) and this is called a 'pre run'.  Differences between those prerun steps are used to determine higher-order derivatives, from 4th to 12th order.

(*Simulate Mercury orbiting a fixed Sun and an outer planet in a fixed circular orbit.  The outer planet can be selected by changing the variable pno, e.g. pno=2 is Venus and pno=8 is Neptune.  
The outer planet has a circular orbit with radius equal to its semimajor axis.  This is a high accuracy version of the above simulation.  It has a runtime of about 12 mins*)
Timing[Ff2[mas_,wp_,rp3_]:=-G mas wp/rp3;Ff3[mas_,wp_,rp3_,xp_,yp_,wv_,xv_,yv_]:=-G mas(wv/rp3-3 wp (xp xv+yp yv)/(rp3^(5/3)));Posd[bp_,wp_,tp_]:=derivs[[bp,wp,1]]+derivs[[bp,wp,2]] tp+derivs[[bp,wp,3]](tp^2)/2+derivs[[bp,wp,4]](tp^3)/6+derivs[[bp,wp,5]](tp^4)/24+derivs[[bp,wp,6]](tp^5)/120+derivs[[bp,wp,7]](tp^6)/720+derivs[[bp,wp,8]](tp^7)/5040+derivs[[bp,wp,9]](tp^8)/40320+derivs[[bp,wp,10]](tp^9)/362880+derivs[[bp,wp,11]](tp^10)/3628800+derivs[[bp,wp,12]](tp^11)/39916800+derivs[[bp,wp,13]](tp^12)/479001600;Veld[bp_,wp_,tp_]:=derivs[[bp,wp,2]] +derivs[[bp,wp,3]]tp+derivs[[bp,wp,4]](tp^2)/2+derivs[[bp,wp,5]](tp^3)/6+derivs[[bp,wp,6]](tp^4)/24+derivs[[bp,wp,7]](tp^5)/120+derivs[[bp,wp,8]](tp^6)/720+derivs[[bp,wp,9]](tp^7)/5040+derivs[[bp,wp,10]](tp^8)/40320+derivs[[bp,wp,11]](tp^9)/362880+derivs[[bp,wp,12]](tp^10)/3628800+derivs[[bp,wp,13]](tp^11)/39916800;
VelMerc[tm_]:=((Posd[2,1,tm]-Posd[1,1,tm])(Veld[2,1,tm]-Veld[1,1,tm])+(Posd[2,2,tm]-Posd[1,2,tm])(Veld[2,2,tm]-Veld[1,2,tm]))/√((((Posd[2,1,tm]-Posd[1,1,tm]))^2)+(((Posd[2,2,tm]-Posd[1,2,tm]))^2));
AngArcsec[tm_]:=ArcTan[(Posd[2,1,tm]-Posd[1,1,tm]),(Posd[2,2,tm]-Posd[1,2,tm])]/π 180 3600;
pno=2;G=SetPrecision[Gg,180];M0=SetPrecision[Msun,180];M1=SetPrecision[Mmercury,180];M2=SetPrecision[Pmass[[pno]],180];R1=SetPrecision[RPmercury,180];R2=SetPrecision[Psemimaj[[pno]],180];e1=SetPrecision[Emercury,180];
nbods=3;frbod=2;tobod=2;ndr=12;PeriAng={};firstperi=True;
Px1=R1;Py1=0;Vx1=0;Vy1=(Sqrt[(1+e1)G M0/R1]);Px2=R2;Py2=0;Vx2=0;Vy2=(Sqrt[(G M0)/R2]);
permerc=2π(Sqrt[(((R1/(1-e1)))^3)/(G M0)]);
Px2=SetPrecision[Px2,180];Vy2=SetPrecision[Vy2,180];posn={{0,0},{Px1,Py1},{ Px2, Py2}};velo={{0,0},{Vx1,Vy1},{Vx2,Vy2}};masses=SetPrecision[{M0,M1,M2},180];prerun=Array[0&,{nbods,2,ndr-2,4}];derivs=Array[0&,{nbods,2,ndr+1}];
count=0;ndiv=132 ;dt=permerc/ndiv;ht=SetPrecision[(10^(-17))*permerc/(2π),180];tt=0;
While[tt<100OneYear,
posn=SetPrecision[posn,180];velo=SetPrecision[velo,180];
(*determine "pre-run" - a series of jerk calculations done very close together*)
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
prerun[[b,w,1,1]]=posn[[b,w]];
prerun[[b,w,1,2]]=velo[[b,w]];
prerun[[b,w,1,3]]=0;
prerun[[b,w,1,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((posn[[b,1]]-posn[[bs,1]]))^2)+(((posn[[b,2]]-posn[[bs,2]]))^2)))^(3/2));
prerun[[b,w,1,3]]+=Ff2[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3];
prerun[[b,w,1,4]]+=Ff3[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3,posn[[b,1]]-posn[[bs,1]],posn[[b,2]]-posn[[bs,2]],velo[[b,w]]-velo[[bs,w]],velo[[b,1]]-velo[[bs,1]],velo[[b,2]]-velo[[bs,2]]];
];
];
];
];
prerun[[3,1,1,1]]=posn[[3,1]];
prerun[[3,2,1,1]]=posn[[3,2]];
prerun[[3,1,1,2]]=velo[[3,1]];
prerun[[3,2,1,2]]=velo[[3,2]];
For[p=2,p≤ndr-2,p++,
tt+=ht;
ps=p-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
Ff2C=0;
Ff3C=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]]))^2)+(((prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]]))^2)))^(3/2));
Ff2C+=Ff2[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3];
Ff3C+=Ff3[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3,prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]],prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]],prerun[[b,w,ps,2]]-prerun[[bs,w,ps,2]],prerun[[b,1,ps,2]]-prerun[[bs,1,ps,2]],prerun[[b,2,ps,2]]-prerun[[bs,2,ps,2]]];
];
];
prerun[[b,w,p,1]]=prerun[[b,w,ps,1]]+prerun[[b,w,ps,2]]ht+Ff2C (ht^2)/2+Ff3C (ht^3)/6;
prerun[[b,w,p,2]]=prerun[[b,w,ps,2]]+Ff2C ht+Ff3C (ht^2)/2;
];
];
prerun[[3,1,p,1]]=Px2 Cos[tt Vy2/Px2];
prerun[[3,2,p,1]]=Px2 Sin[tt Vy2/Px2];
prerun[[3,1,p,2]]=-Vy2 Sin[tt Vy2/Px2];
prerun[[3,2,p,2]]=Vy2 Cos[tt Vy2/Px2];
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
prerun[[b,w,p,3]]=0;
prerun[[b,w,p,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,p,1]]-prerun[[bs,1,p,1]]))^2)+(((prerun[[b,2,p,1]]-prerun[[bs,2,p,1]]))^2)))^(3/2));
prerun[[b,w,p,3]]+=Ff2[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3];prerun[[b,w,p,4]]+=Ff3[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3,prerun[[b,1,p,1]]-prerun[[bs,1,p,1]],prerun[[b,2,p,1]]-prerun[[bs,2,p,1]],prerun[[b,w,p,2]]-prerun[[bs,w,p,2]],prerun[[b,1,p,2]]-prerun[[bs,1,p,2]],prerun[[b,2,p,2]]-prerun[[bs,2,p,2]]];
];
];
];
];
];
(*determine derivitives from prerun calcs*)
ndm1=ndr-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
derivs[[b,w,1]]=prerun[[b,w,ndm1-1,1]];derivs[[b,w,2]]=prerun[[b,w,ndm1-1,2]];derivs[[b,w,3]]=prerun[[b,w,ndm1-1,3]];derivs[[b,w,4]]=prerun[[b,w,ndm1-1,4]];
derivs[[b,w,5]]=(prerun[[b,w,ndm1-1,4]]-prerun[[b,w,ndm1-2,4]])/ht;
derivs[[b,w,6]]=(prerun[[b,w,ndm1-1,4]]-2prerun[[b,w,ndm1-2,4]]+prerun[[b,w,ndm1-3,4]])/(ht^2);
derivs[[b,w,7]]=(prerun[[b,w,ndm1-1,4]]-3prerun[[b,w,ndm1-2,4]]+3prerun[[b,w,ndm1-3,4]]-prerun[[b,w,ndm1-4,4]])/(ht^3);
derivs[[b,w,8]]=(prerun[[b,w,ndm1-1,4]]-4prerun[[b,w,ndm1-2,4]]+6prerun[[b,w,ndm1-3,4]]-4prerun[[b,w,ndm1-4,4]]+prerun[[b,w,ndm1-5,4]])/(ht^4);derivs[[b,w,9]]=(prerun[[b,w,ndm1-1,4]]-5prerun[[b,w,ndm1-2,4]]+10prerun[[b,w,ndm1-3,4]]-10prerun[[b,w,ndm1-4,4]]+5prerun[[b,w,ndm1-5,4]]-prerun[[b,w,ndm1-6,4]])/(ht^5);
derivs[[b,w,10]]=(prerun[[b,w,ndm1-1,4]]-6prerun[[b,w,ndm1-2,4]]+15prerun[[b,w,ndm1-3,4]]-20prerun[[b,w,ndm1-4,4]]+15prerun[[b,w,ndm1-5,4]]-6prerun[[b,w,ndm1-6,4]]+prerun[[b,w,ndm1-7,4]])/(ht^6);
derivs[[b,w,11]]=(prerun[[b,w,ndm1-1,4]]-7prerun[[b,w,ndm1-2,4]]+21prerun[[b,w,ndm1-3,4]]-35prerun[[b,w,ndm1-4,4]]+35prerun[[b,w,ndm1-5,4]]-21prerun[[b,w,ndm1-6,4]]+7prerun[[b,w,ndm1-7,4]]-prerun[[b,w,ndm1-8,4]])/(ht^7);
derivs[[b,w,12]]=(prerun[[b,w,ndm1-1,4]]-8prerun[[b,w,ndm1-2,4]]+28prerun[[b,w,ndm1-3,4]]-56prerun[[b,w,ndm1-4,4]]+70prerun[[b,w,ndm1-5,4]]-56prerun[[b,w,ndm1-6,4]]+28prerun[[b,w,ndm1-7,4]]-8prerun[[b,w,ndm1-8,4]]+prerun[[b,w,ndm1-9,4]])/(ht^8);
derivs[[b,w,13]]=(prerun[[b,w,ndm1-1,4]]-9prerun[[b,w,ndm1-2,4]]+36prerun[[b,w,ndm1-3,4]]-84prerun[[b,w,ndm1-4,4]]+126prerun[[b,w,ndm1-5,4]]-126prerun[[b,w,ndm1-6,4]]+84prerun[[b,w,ndm1-7,4]]-36prerun[[b,w,ndm1-8,4]]+9prerun[[b,w,ndm1-9,4]]-prerun[[b,w,ndm1-10,4]])/(ht^9);
];
];
(*determine new position/velocity simulation points*)
derivs=SetPrecision[derivs,31];
dt=SetPrecision[permerc/ndiv-(ndr-3)ht,31];
tt+=dt;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
posn[[b,w]]=Posd[b,w,dt];velo[[b,w]]=Veld[b,w,dt];
];
];
tt=SetPrecision[tt,180];
posn[[3,1]]=Px2 Cos[tt Vy2/Px2];
posn[[3,2]]=Px2 Sin[tt Vy2/Px2];
velo[[3,1]]=-Vy2 Sin[tt Vy2/Px2];
velo[[3,2]]=Vy2 Cos[tt Vy2/Px2];
(*determine perihelion location of Mercury*)
If[posn[[2,1]]>0.98Px1&&posn[[2,2]]>Py1-0.1Px1&&posn[[2,2]]<Py1+0.1Px1,
awyvel0=VelMerc[0];awyvel1=VelMerc[dt];
If[firstperi==True||(awyvel0<0&&awyvel1>0),
crosdt2=-awyvel0 dt/(awyvel1-awyvel0);
awyvel2=VelMerc[crosdt2];
If[awyvel2<0,
crosdt3=crosdt2-awyvel2( dt-crosdt2)/(awyvel1-awyvel2);
,
crosdt3=-awyvel0 crosdt2/(awyvel2-awyvel0);
];
perang=AngArcsec[crosdt3];
crostim=tt-dt+crosdt3;
AppendTo[PeriAng,{SetPrecision[crostim/OneYear,5],SetPrecision[perang,5]}];
];
firstperi=False;
];
count+=1;
];]

(*Venus*)ListPlot[PeriAng,Joined->True]

Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Earth*)ListPlot[PeriAng,Joined->True]

Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Mars*)ListPlot[PeriAng,Joined->True]

Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Jupiter*)ListPlot[PeriAng,Joined->True]

Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Saturn*)ListPlot[PeriAng[[1;;429]],Joined->True](*103years*)

Normal[LinearModelFit[PeriAng[[1;;429]],x,x]][[2]]100/x

(*Uranus*)ListPlot[PeriAng,Joined->True]

Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Neptune*)ListPlot[PeriAng,Joined->True]

Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Total = 293.0+96.0+2.4+155.3+7.3+0.13+0.03 = 554.1 *)

(*Simulate Mercury orbiting a fixed Sun and surrounded by Jupiter and Saturn.  The outer planets will have a circular orbit with a radius equal to their semimajor axis.  Runtime approx 20mins.*)Timing[Ff2[mas_,wp_,rp3_]:=-G mas wp/rp3;Ff3[mas_,wp_,rp3_,xp_,yp_,wv_,xv_,yv_]:=-G mas(wv/rp3-3 wp (xp xv+yp yv)/(rp3^(5/3)));Posd[bp_,wp_,tp_]:=derivs[[bp,wp,1]]+derivs[[bp,wp,2]] tp+derivs[[bp,wp,3]](tp^2)/2+derivs[[bp,wp,4]](tp^3)/6+derivs[[bp,wp,5]](tp^4)/24+derivs[[bp,wp,6]](tp^5)/120+derivs[[bp,wp,7]](tp^6)/720+derivs[[bp,wp,8]](tp^7)/5040+derivs[[bp,wp,9]](tp^8)/40320+derivs[[bp,wp,10]](tp^9)/362880+derivs[[bp,wp,11]](tp^10)/3628800+derivs[[bp,wp,12]](tp^11)/39916800+derivs[[bp,wp,13]](tp^12)/479001600;Veld[bp_,wp_,tp_]:=derivs[[bp,wp,2]] +derivs[[bp,wp,3]]tp+derivs[[bp,wp,4]](tp^2)/2+derivs[[bp,wp,5]](tp^3)/6+derivs[[bp,wp,6]](tp^4)/24+derivs[[bp,wp,7]](tp^5)/120+derivs[[bp,wp,8]](tp^6)/720+derivs[[bp,wp,9]](tp^7)/5040+derivs[[bp,wp,10]](tp^8)/40320+derivs[[bp,wp,11]](tp^9)/362880+derivs[[bp,wp,12]](tp^10)/3628800+derivs[[bp,wp,13]](tp^11)/39916800;
VelMerc[tm_]:=((Posd[2,1,tm]-Posd[1,1,tm])(Veld[2,1,tm]-Veld[1,1,tm])+(Posd[2,2,tm]-Posd[1,2,tm])(Veld[2,2,tm]-Veld[1,2,tm]))/√((((Posd[2,1,tm]-Posd[1,1,tm]))^2)+(((Posd[2,2,tm]-Posd[1,2,tm]))^2));
AngArcsec[tm_]:=ArcTan[(Posd[2,1,tm]-Posd[1,1,tm]),(Posd[2,2,tm]-Posd[1,2,tm])]/π 180 3600;
G=SetPrecision[Gg,180];{M0,M1,M2,M3}=SetPrecision[{Msun,Mmercury,Mjupiter,Msaturn},180];
{R1,R2,R3}=SetPrecision[{RPmercury,RAjupiter,RAsaturn},180];e1=SetPrecision[Emercury,180];
nbods=4;frbod=2;tobod=2;ndr=12;PeriAng={};firstperi=True;
{Px1,Px2,Px3}={R1,R2,R3};{Py1,Py2,Py3}={0,0,0};
{Vx1,Vx2,Vx3}={0,0,0};{Vy1,Vy2,Vy3}={(Sqrt[(1+e1)G M0/R1]),(Sqrt[(G M0)/R2]),(Sqrt[(G M0)/R3])};
permerc=2π(Sqrt[(((R1/(1-e1)))^3)/(G M0)]);
Px2=SetPrecision[Px2,180];Vy2=SetPrecision[Vy2,180];posn={{0,0},{Px1,Py1},{ Px2, Py2},{ Px3, Py3}};velo={{0,0},{Vx1,Vy1},{Vx2,Vy2},{Vx3,Vy3}};masses=SetPrecision[{M0,M1,M2,M3},180];
strXpos={0,Px1,Px2,Px3};strYvel={0,Vy1,Vy2,Vy3};prerun=Array[0&,{nbods,2,ndr-2,4}];derivs=Array[0&,{nbods,2,ndr+1}];
count=0;ndiv=132;dt=permerc/ndiv;ht=SetPrecision[(10^(-17))*permerc/(2π),180];tt=0;
While[tt<115OneYear,
posn=SetPrecision[posn,180];velo=SetPrecision[velo,180];
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
prerun[[b,w,1,1]]=posn[[b,w]];
prerun[[b,w,1,2]]=velo[[b,w]];
prerun[[b,w,1,3]]=0;
prerun[[b,w,1,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((posn[[b,1]]-posn[[bs,1]]))^2)+(((posn[[b,2]]-posn[[bs,2]]))^2)))^(3/2));
prerun[[b,w,1,3]]+=Ff2[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3];
prerun[[b,w,1,4]]+=Ff3[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3,posn[[b,1]]-posn[[bs,1]],posn[[b,2]]-posn[[bs,2]],velo[[b,w]]-velo[[bs,w]],velo[[b,1]]-velo[[bs,1]],velo[[b,2]]-velo[[bs,2]]];
];
];
];
];
For[b=tobod+1,b≤nbods,b++,
For[w=1,w≤2,w++,
prerun[[b,w,1,1]]=posn[[b,w]];
prerun[[b,w,1,2]]=velo[[b,w]];
];
];
For[p=2,p≤ndr-2,p++,
tt+=ht;
ps=p-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
Ff2C=0;
Ff3C=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]]))^2)+(((prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]]))^2)))^(3/2));
Ff2C+=Ff2[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3];
Ff3C+=Ff3[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3,prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]],prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]],prerun[[b,w,ps,2]]-prerun[[bs,w,ps,2]],prerun[[b,1,ps,2]]-prerun[[bs,1,ps,2]],prerun[[b,2,ps,2]]-prerun[[bs,2,ps,2]]];
];
];
prerun[[b,w,p,1]]=prerun[[b,w,ps,1]]+prerun[[b,w,ps,2]]ht+Ff2C (ht^2)/2+Ff3C (ht^3)/6;
prerun[[b,w,p,2]]=prerun[[b,w,ps,2]]+Ff2C ht+Ff3C (ht^2)/2;
];
];
For[b=tobod+1,b≤nbods,b++,
velDdis=strYvel[[b]]/strXpos[[b]];
prerun[[b,1,p,1]]=strXpos[[b]] Cos[tt velDdis];
prerun[[b,2,p,1]]=strXpos[[b]] Sin[tt velDdis];
prerun[[b,1,p,2]]=-strYvel[[b]] Sin[tt velDdis];
prerun[[b,2,p,2]]=strYvel[[b]] Cos[tt velDdis];
];
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
prerun[[b,w,p,3]]=0;
prerun[[b,w,p,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,p,1]]-prerun[[bs,1,p,1]]))^2)+(((prerun[[b,2,p,1]]-prerun[[bs,2,p,1]]))^2)))^(3/2));
prerun[[b,w,p,3]]+=Ff2[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3];prerun[[b,w,p,4]]+=Ff3[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3,prerun[[b,1,p,1]]-prerun[[bs,1,p,1]],prerun[[b,2,p,1]]-prerun[[bs,2,p,1]],prerun[[b,w,p,2]]-prerun[[bs,w,p,2]],prerun[[b,1,p,2]]-prerun[[bs,1,p,2]],prerun[[b,2,p,2]]-prerun[[bs,2,p,2]]];
];
];
];
];
];
ndm1=ndr-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
derivs[[b,w,1]]=prerun[[b,w,ndm1-1,1]];derivs[[b,w,2]]=prerun[[b,w,ndm1-1,2]];derivs[[b,w,3]]=prerun[[b,w,ndm1-1,3]];derivs[[b,w,4]]=prerun[[b,w,ndm1-1,4]];
derivs[[b,w,5]]=(prerun[[b,w,ndm1-1,4]]-prerun[[b,w,ndm1-2,4]])/ht;
derivs[[b,w,6]]=(prerun[[b,w,ndm1-1,4]]-2prerun[[b,w,ndm1-2,4]]+prerun[[b,w,ndm1-3,4]])/(ht^2);
derivs[[b,w,7]]=(prerun[[b,w,ndm1-1,4]]-3prerun[[b,w,ndm1-2,4]]+3prerun[[b,w,ndm1-3,4]]-prerun[[b,w,ndm1-4,4]])/(ht^3);
derivs[[b,w,8]]=(prerun[[b,w,ndm1-1,4]]-4prerun[[b,w,ndm1-2,4]]+6prerun[[b,w,ndm1-3,4]]-4prerun[[b,w,ndm1-4,4]]+prerun[[b,w,ndm1-5,4]])/(ht^4);derivs[[b,w,9]]=(prerun[[b,w,ndm1-1,4]]-5prerun[[b,w,ndm1-2,4]]+10prerun[[b,w,ndm1-3,4]]-10prerun[[b,w,ndm1-4,4]]+5prerun[[b,w,ndm1-5,4]]-prerun[[b,w,ndm1-6,4]])/(ht^5);
derivs[[b,w,10]]=(prerun[[b,w,ndm1-1,4]]-6prerun[[b,w,ndm1-2,4]]+15prerun[[b,w,ndm1-3,4]]-20prerun[[b,w,ndm1-4,4]]+15prerun[[b,w,ndm1-5,4]]-6prerun[[b,w,ndm1-6,4]]+prerun[[b,w,ndm1-7,4]])/(ht^6);
derivs[[b,w,11]]=(prerun[[b,w,ndm1-1,4]]-7prerun[[b,w,ndm1-2,4]]+21prerun[[b,w,ndm1-3,4]]-35prerun[[b,w,ndm1-4,4]]+35prerun[[b,w,ndm1-5,4]]-21prerun[[b,w,ndm1-6,4]]+7prerun[[b,w,ndm1-7,4]]-prerun[[b,w,ndm1-8,4]])/(ht^7);
derivs[[b,w,12]]=(prerun[[b,w,ndm1-1,4]]-8prerun[[b,w,ndm1-2,4]]+28prerun[[b,w,ndm1-3,4]]-56prerun[[b,w,ndm1-4,4]]+70prerun[[b,w,ndm1-5,4]]-56prerun[[b,w,ndm1-6,4]]+28prerun[[b,w,ndm1-7,4]]-8prerun[[b,w,ndm1-8,4]]+prerun[[b,w,ndm1-9,4]])/(ht^8);
derivs[[b,w,13]]=(prerun[[b,w,ndm1-1,4]]-9prerun[[b,w,ndm1-2,4]]+36prerun[[b,w,ndm1-3,4]]-84prerun[[b,w,ndm1-4,4]]+126prerun[[b,w,ndm1-5,4]]-126prerun[[b,w,ndm1-6,4]]+84prerun[[b,w,ndm1-7,4]]-36prerun[[b,w,ndm1-8,4]]+9prerun[[b,w,ndm1-9,4]]-prerun[[b,w,ndm1-10,4]])/(ht^9);
];
];
derivs=SetPrecision[derivs,31];
dt=SetPrecision[permerc/ndiv-(ndr-3)ht,31];
tt+=dt;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
posn[[b,w]]=Posd[b,w,dt];velo[[b,w]]=Veld[b,w,dt];
];
];
tt=SetPrecision[tt,180];
For[b=tobod+1,b≤nbods,b++,
velDdis=strYvel[[b]]/strXpos[[b]];
posn[[b,1]]=strXpos[[b]] Cos[tt velDdis];
posn[[b,2]]=strXpos[[b]] Sin[tt velDdis];
velo[[b,1]]=-strYvel[[b]] Sin[tt velDdis];
velo[[b,2]]=strYvel[[b]] Cos[tt velDdis];
];
If[posn[[2,1]]>0.98Px1&&posn[[2,2]]>Py1-0.1Px1&&posn[[2,2]]<Py1+0.1Px1,
awyvel0=VelMerc[0];awyvel1=VelMerc[dt];
If[firstperi==True||(awyvel0<0&&awyvel1>0),
crosdt2=-awyvel0 dt/(awyvel1-awyvel0);
awyvel2=VelMerc[crosdt2];
If[awyvel2<0,
crosdt3=crosdt2-awyvel2( dt-crosdt2)/(awyvel1-awyvel2);
,
crosdt3=-awyvel0 crosdt2/(awyvel2-awyvel0);
];
perang=AngArcsec[crosdt3];
crostim=tt-dt+crosdt3;
AppendTo[PeriAng,{SetPrecision[crostim/OneYear,5],SetPrecision[perang,5]}];
];
firstperi=False;
];
count+=1;
];]

(*Jupiter+Saturn*)ListPlot[PeriAng[[1;;469]],Joined->True](*112.7years*)

Normal[LinearModelFit[PeriAng[[1;;469]],x,x]][[2]]100/x

(*Simulate Mercury orbiting a fixed Sun and surrounded by all the outer planets in fixed orbits.  The outer planets will have a circular orbit with a radius equal to their semimajor axis.  Runtime approx 74mins.*)Timing[Ff2[mas_,wp_,rp3_]:=-G mas wp/rp3;Ff3[mas_,wp_,rp3_,xp_,yp_,wv_,xv_,yv_]:=-G mas(wv/rp3-3 wp (xp xv+yp yv)/(rp3^(5/3)));Posd[bp_,wp_,tp_]:=derivs[[bp,wp,1]]+derivs[[bp,wp,2]] tp+derivs[[bp,wp,3]](tp^2)/2+derivs[[bp,wp,4]](tp^3)/6+derivs[[bp,wp,5]](tp^4)/24+derivs[[bp,wp,6]](tp^5)/120+derivs[[bp,wp,7]](tp^6)/720+derivs[[bp,wp,8]](tp^7)/5040+derivs[[bp,wp,9]](tp^8)/40320+derivs[[bp,wp,10]](tp^9)/362880+derivs[[bp,wp,11]](tp^10)/3628800+derivs[[bp,wp,12]](tp^11)/39916800+derivs[[bp,wp,13]](tp^12)/479001600;Veld[bp_,wp_,tp_]:=derivs[[bp,wp,2]] +derivs[[bp,wp,3]]tp+derivs[[bp,wp,4]](tp^2)/2+derivs[[bp,wp,5]](tp^3)/6+derivs[[bp,wp,6]](tp^4)/24+derivs[[bp,wp,7]](tp^5)/120+derivs[[bp,wp,8]](tp^6)/720+derivs[[bp,wp,9]](tp^7)/5040+derivs[[bp,wp,10]](tp^8)/40320+derivs[[bp,wp,11]](tp^9)/362880+derivs[[bp,wp,12]](tp^10)/3628800+derivs[[bp,wp,13]](tp^11)/39916800;
VelMerc[tm_]:=((Posd[2,1,tm]-Posd[1,1,tm])(Veld[2,1,tm]-Veld[1,1,tm])+(Posd[2,2,tm]-Posd[1,2,tm])(Veld[2,2,tm]-Veld[1,2,tm]))/√((((Posd[2,1,tm]-Posd[1,1,tm]))^2)+(((Posd[2,2,tm]-Posd[1,2,tm]))^2));
AngArcsec[tm_]:=ArcTan[(Posd[2,1,tm]-Posd[1,1,tm]),(Posd[2,2,tm]-Posd[1,2,tm])]/π 180 3600;
G=SetPrecision[Gg,180];{M0,M1,M2,M3,M4,M5,M6,M7,M8}=SetPrecision[{Msun,Mmercury,Mvenus,Mearth,Mmars,Mjupiter,Msaturn,Muranus,Mneptune},180];
{R1,R2,R3,R4,R5,R6,R7,R8}=SetPrecision[{RPmercury,RAvenus,RAearth,RAmars,RAjupiter,RAsaturn,RAuranus,RAneptune},180];e1=SetPrecision[Emercury,180];
nbods=9;frbod=2;tobod=2;ndr=12;PeriAng={};firstperi=True;
{Px1,Px2,Px3,Px4,Px5,Px6,Px7,Px8}={R1,R2,R3,R4,R5,R6,R7,R8};{Py1,Py2,Py3,Py4,Py5,Py6,Py7,Py8}={0,0,0,0,0,0,0,0};
{Vx1,Vx2,Vx3,Vx4,Vx5,Vx6,Vx7,Vx8}={0,0,0,0,0,0,0,0};{Vy1,Vy2,Vy3,Vy4,Vy5,Vy6,Vy7,Vy8}={(Sqrt[(1+e1)G M0/R1]),(Sqrt[(G M0)/R2]),(Sqrt[(G M0)/R3]),(Sqrt[(G M0)/R4]),(Sqrt[(G M0)/R5]),(Sqrt[(G M0)/R6]),(Sqrt[(G M0)/R7]),(Sqrt[(G M0)/R8])};
permerc=2π(Sqrt[(((R1/(1-e1)))^3)/(G M0)]);posn2={{Px1,Py1}};
posn={{0,0},{Px1,Py1},{ Px2,Py2},{Px3,Py3},{Px4,Py4},{Px5,Py5},{Px6,Py6},{Px7,Py7},{Px8,Py8}};velo={{0,0},{Vx1,Vy1},{Vx2,Vy2},{Vx3,Vy3},{Vx4,Vy4},{Vx5,Vy5},{Vx6,Vy6},{Vx7,Vy7},{Vx8,Vy8}};masses=SetPrecision[{M0,M1,M2,M3,M4,M5,M6,M7,M8},180];
strXpos={0,Px1,Px2,Px3,Px4,Px5,Px6,Px7,Px8};strYvel={0,Vy1,Vy2,Vy3,Vy4,Vy5,Vy6,Vy7,Vy8};prerun=Array[0&,{nbods,2,ndr-2,4}];derivs=Array[0&,{nbods,2,ndr+1}];
count=0;ndiv=132;dt=permerc/ndiv;ht=SetPrecision[(10^(-17))*permerc/(2π),180];tt=0;
While[tt<150OneYear,
posn=SetPrecision[posn,180];velo=SetPrecision[velo,180];
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
prerun[[b,w,1,1]]=posn[[b,w]];
prerun[[b,w,1,2]]=velo[[b,w]];
prerun[[b,w,1,3]]=0;
prerun[[b,w,1,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((posn[[b,1]]-posn[[bs,1]]))^2)+(((posn[[b,2]]-posn[[bs,2]]))^2)))^(3/2));
prerun[[b,w,1,3]]+=Ff2[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3];
prerun[[b,w,1,4]]+=Ff3[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3,posn[[b,1]]-posn[[bs,1]],posn[[b,2]]-posn[[bs,2]],velo[[b,w]]-velo[[bs,w]],velo[[b,1]]-velo[[bs,1]],velo[[b,2]]-velo[[bs,2]]];
];
];
];
];
For[b=tobod+1,b≤nbods,b++,
For[w=1,w≤2,w++,
prerun[[b,w,1,1]]=posn[[b,w]];
prerun[[b,w,1,2]]=velo[[b,w]];
];
];
For[p=2,p≤ndr-2,p++,
tt+=ht;
ps=p-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
Ff2C=0;
Ff3C=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]]))^2)+(((prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]]))^2)))^(3/2));
Ff2C+=Ff2[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3];
Ff3C+=Ff3[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3,prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]],prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]],prerun[[b,w,ps,2]]-prerun[[bs,w,ps,2]],prerun[[b,1,ps,2]]-prerun[[bs,1,ps,2]],prerun[[b,2,ps,2]]-prerun[[bs,2,ps,2]]];
];
];
prerun[[b,w,p,1]]=prerun[[b,w,ps,1]]+prerun[[b,w,ps,2]]ht+Ff2C (ht^2)/2+Ff3C (ht^3)/6;
prerun[[b,w,p,2]]=prerun[[b,w,ps,2]]+Ff2C ht+Ff3C (ht^2)/2;
];
];
For[b=tobod+1,b≤nbods,b++,
velDdis=strYvel[[b]]/strXpos[[b]];
prerun[[b,1,p,1]]=strXpos[[b]] Cos[tt velDdis];
prerun[[b,2,p,1]]=strXpos[[b]] Sin[tt velDdis];
prerun[[b,1,p,2]]=-strYvel[[b]] Sin[tt velDdis];
prerun[[b,2,p,2]]=strYvel[[b]] Cos[tt velDdis];
];
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
prerun[[b,w,p,3]]=0;
prerun[[b,w,p,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,p,1]]-prerun[[bs,1,p,1]]))^2)+(((prerun[[b,2,p,1]]-prerun[[bs,2,p,1]]))^2)))^(3/2));
prerun[[b,w,p,3]]+=Ff2[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3];prerun[[b,w,p,4]]+=Ff3[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3,prerun[[b,1,p,1]]-prerun[[bs,1,p,1]],prerun[[b,2,p,1]]-prerun[[bs,2,p,1]],prerun[[b,w,p,2]]-prerun[[bs,w,p,2]],prerun[[b,1,p,2]]-prerun[[bs,1,p,2]],prerun[[b,2,p,2]]-prerun[[bs,2,p,2]]];
];
];
];
];
];
ndm1=ndr-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
derivs[[b,w,1]]=prerun[[b,w,ndm1-1,1]];derivs[[b,w,2]]=prerun[[b,w,ndm1-1,2]];derivs[[b,w,3]]=prerun[[b,w,ndm1-1,3]];derivs[[b,w,4]]=prerun[[b,w,ndm1-1,4]];
derivs[[b,w,5]]=(prerun[[b,w,ndm1-1,4]]-prerun[[b,w,ndm1-2,4]])/ht;
derivs[[b,w,6]]=(prerun[[b,w,ndm1-1,4]]-2prerun[[b,w,ndm1-2,4]]+prerun[[b,w,ndm1-3,4]])/(ht^2);
derivs[[b,w,7]]=(prerun[[b,w,ndm1-1,4]]-3prerun[[b,w,ndm1-2,4]]+3prerun[[b,w,ndm1-3,4]]-prerun[[b,w,ndm1-4,4]])/(ht^3);
derivs[[b,w,8]]=(prerun[[b,w,ndm1-1,4]]-4prerun[[b,w,ndm1-2,4]]+6prerun[[b,w,ndm1-3,4]]-4prerun[[b,w,ndm1-4,4]]+prerun[[b,w,ndm1-5,4]])/(ht^4);derivs[[b,w,9]]=(prerun[[b,w,ndm1-1,4]]-5prerun[[b,w,ndm1-2,4]]+10prerun[[b,w,ndm1-3,4]]-10prerun[[b,w,ndm1-4,4]]+5prerun[[b,w,ndm1-5,4]]-prerun[[b,w,ndm1-6,4]])/(ht^5);
derivs[[b,w,10]]=(prerun[[b,w,ndm1-1,4]]-6prerun[[b,w,ndm1-2,4]]+15prerun[[b,w,ndm1-3,4]]-20prerun[[b,w,ndm1-4,4]]+15prerun[[b,w,ndm1-5,4]]-6prerun[[b,w,ndm1-6,4]]+prerun[[b,w,ndm1-7,4]])/(ht^6);
derivs[[b,w,11]]=(prerun[[b,w,ndm1-1,4]]-7prerun[[b,w,ndm1-2,4]]+21prerun[[b,w,ndm1-3,4]]-35prerun[[b,w,ndm1-4,4]]+35prerun[[b,w,ndm1-5,4]]-21prerun[[b,w,ndm1-6,4]]+7prerun[[b,w,ndm1-7,4]]-prerun[[b,w,ndm1-8,4]])/(ht^7);
derivs[[b,w,12]]=(prerun[[b,w,ndm1-1,4]]-8prerun[[b,w,ndm1-2,4]]+28prerun[[b,w,ndm1-3,4]]-56prerun[[b,w,ndm1-4,4]]+70prerun[[b,w,ndm1-5,4]]-56prerun[[b,w,ndm1-6,4]]+28prerun[[b,w,ndm1-7,4]]-8prerun[[b,w,ndm1-8,4]]+prerun[[b,w,ndm1-9,4]])/(ht^8);
derivs[[b,w,13]]=(prerun[[b,w,ndm1-1,4]]-9prerun[[b,w,ndm1-2,4]]+36prerun[[b,w,ndm1-3,4]]-84prerun[[b,w,ndm1-4,4]]+126prerun[[b,w,ndm1-5,4]]-126prerun[[b,w,ndm1-6,4]]+84prerun[[b,w,ndm1-7,4]]-36prerun[[b,w,ndm1-8,4]]+9prerun[[b,w,ndm1-9,4]]-prerun[[b,w,ndm1-10,4]])/(ht^9);
];
];
derivs=SetPrecision[derivs,31];
dt=SetPrecision[permerc/ndiv-(ndr-3)ht,31];
tt+=dt;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
posn[[b,w]]=Posd[b,w,dt];velo[[b,w]]=Veld[b,w,dt];
];
];
tt=SetPrecision[tt,180];
For[b=tobod+1,b≤nbods,b++,
velDdis=strYvel[[b]]/strXpos[[b]];
posn[[b,1]]=strXpos[[b]] Cos[tt velDdis];
posn[[b,2]]=strXpos[[b]] Sin[tt velDdis];
velo[[b,1]]=-strYvel[[b]] Sin[tt velDdis];
velo[[b,2]]=strYvel[[b]] Cos[tt velDdis];
];
If[posn[[2,1]]>0.98Px1&&posn[[2,2]]>Py1-0.1Px1&&posn[[2,2]]<Py1+0.1Px1,
awyvel0=VelMerc[0];awyvel1=VelMerc[dt];
If[firstperi==True||(awyvel0<0&&awyvel1>0),
crosdt2=-awyvel0 dt/(awyvel1-awyvel0);
awyvel2=VelMerc[crosdt2];
If[awyvel2<0,
crosdt3=crosdt2-awyvel2( dt-crosdt2)/(awyvel1-awyvel2);
,
crosdt3=-awyvel0 crosdt2/(awyvel2-awyvel0);
];
perang=AngArcsec[crosdt3];
crostim=tt-dt+crosdt3;
AppendTo[PeriAng,{SetPrecision[crostim/OneYear,5],SetPrecision[perang,5]}];
];
firstperi=False;
];
count+=1;
];]

ListPlot[PeriAng[[1;;619]],Joined->True](*148.8years*)

Normal[LinearModelFit[PeriAng[[1;;619]],x,x]][[2]]100/x

(*Simulate Mercury orbiting a floating Sun and an outer planet in a circular orbit.  The outer planet can be selected by changing the variable pno, e.g. pno=2 is Venus and pno=8 is Neptune.  
The outer planet has a circular orbit with radius equal to its semimajor axis.  It's initial velocity will be boosted to ensure it is circular.  Runtime approx 30min.*)
Timing[Ff2[mas_,wp_,rp3_]:=-G mas wp/rp3;Ff3[mas_,wp_,rp3_,xp_,yp_,wv_,xv_,yv_]:=-G mas(wv/rp3-3 wp (xp xv+yp yv)/(rp3^(5/3)));Posd[bp_,wp_,tp_]:=derivs[[bp,wp,1]]+derivs[[bp,wp,2]] tp+derivs[[bp,wp,3]](tp^2)/2+derivs[[bp,wp,4]](tp^3)/6+derivs[[bp,wp,5]](tp^4)/24+derivs[[bp,wp,6]](tp^5)/120+derivs[[bp,wp,7]](tp^6)/720+derivs[[bp,wp,8]](tp^7)/5040+derivs[[bp,wp,9]](tp^8)/40320+derivs[[bp,wp,10]](tp^9)/362880+derivs[[bp,wp,11]](tp^10)/3628800+derivs[[bp,wp,12]](tp^11)/39916800+derivs[[bp,wp,13]](tp^12)/479001600;Veld[bp_,wp_,tp_]:=derivs[[bp,wp,2]] +derivs[[bp,wp,3]]tp+derivs[[bp,wp,4]](tp^2)/2+derivs[[bp,wp,5]](tp^3)/6+derivs[[bp,wp,6]](tp^4)/24+derivs[[bp,wp,7]](tp^5)/120+derivs[[bp,wp,8]](tp^6)/720+derivs[[bp,wp,9]](tp^7)/5040+derivs[[bp,wp,10]](tp^8)/40320+derivs[[bp,wp,11]](tp^9)/362880+derivs[[bp,wp,12]](tp^10)/3628800+derivs[[bp,wp,13]](tp^11)/39916800;
VelMerc[tm_]:=((Posd[2,1,tm]-Posd[1,1,tm])(Veld[2,1,tm]-Veld[1,1,tm])+(Posd[2,2,tm]-Posd[1,2,tm])(Veld[2,2,tm]-Veld[1,2,tm]))/√((((Posd[2,1,tm]-Posd[1,1,tm]))^2)+(((Posd[2,2,tm]-Posd[1,2,tm]))^2));
AngArcsec[tm_]:=ArcTan[(Posd[2,1,tm]-Posd[1,1,tm]),(Posd[2,2,tm]-Posd[1,2,tm])]/π 180 3600;
pno=5;G=SetPrecision[Gg,180];{M0,M1,M2}=SetPrecision[{Msun,Mmercury,Pmass[[pno]]},180];
{R1,R2}=SetPrecision[{RPmercury,Psemimaj[[pno]]},180];e1=SetPrecision[Emercury,180];
nbods=3;frbod=1;tobod=nbods;ndr=12;PeriAng={};firstperi=True;
{Px1,Px2}={R1,R2};{Py1,Py2}={0,0};
{Vx1,Vx2}={0,0};{Vy1,Vy2}={(Sqrt[(1+e1)G M0/R1]),SetPrecision[PveladjS[[pno]],180](Sqrt[(G M0)/R2])};
permerc=2π(Sqrt[(((R1/(1-e1)))^3)/(G M0)]);
posn={{0,0},{Px1,Py1},{ Px2,Py2}};velo={{0,0},{Vx1,Vy1},{Vx2,Vy2}};masses={M0,M1,M2};
prerun=Array[0&,{nbods,2,ndr-2,4}];derivs=Array[0&,{nbods,2,ndr+1}];
count=0;ndiv=132;dt=permerc/ndiv;ht=SetPrecision[(10^(-17))*permerc/(2π),180];tt=0;
While[tt<101OneYear,
posn=SetPrecision[posn,180];velo=SetPrecision[velo,180];
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
prerun[[b,w,1,1]]=posn[[b,w]];
prerun[[b,w,1,2]]=velo[[b,w]];
prerun[[b,w,1,3]]=0;
prerun[[b,w,1,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((posn[[b,1]]-posn[[bs,1]]))^2)+(((posn[[b,2]]-posn[[bs,2]]))^2)))^(3/2));
prerun[[b,w,1,3]]+=Ff2[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3];
prerun[[b,w,1,4]]+=Ff3[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3,posn[[b,1]]-posn[[bs,1]],posn[[b,2]]-posn[[bs,2]],velo[[b,w]]-velo[[bs,w]],velo[[b,1]]-velo[[bs,1]],velo[[b,2]]-velo[[bs,2]]];
];
];
];
];
For[p=2,p≤ndr-2,p++,
tt+=ht;
ps=p-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
Ff2C=0;
Ff3C=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]]))^2)+(((prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]]))^2)))^(3/2));
Ff2C+=Ff2[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3];
Ff3C+=Ff3[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3,prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]],prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]],prerun[[b,w,ps,2]]-prerun[[bs,w,ps,2]],prerun[[b,1,ps,2]]-prerun[[bs,1,ps,2]],prerun[[b,2,ps,2]]-prerun[[bs,2,ps,2]]];
];
];
prerun[[b,w,p,1]]=prerun[[b,w,ps,1]]+prerun[[b,w,ps,2]]ht+Ff2C (ht^2)/2+Ff3C (ht^3)/6;
prerun[[b,w,p,2]]=prerun[[b,w,ps,2]]+Ff2C ht+Ff3C (ht^2)/2;
];
];
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
prerun[[b,w,p,3]]=0;
prerun[[b,w,p,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,p,1]]-prerun[[bs,1,p,1]]))^2)+(((prerun[[b,2,p,1]]-prerun[[bs,2,p,1]]))^2)))^(3/2));
prerun[[b,w,p,3]]+=Ff2[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3];prerun[[b,w,p,4]]+=Ff3[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3,prerun[[b,1,p,1]]-prerun[[bs,1,p,1]],prerun[[b,2,p,1]]-prerun[[bs,2,p,1]],prerun[[b,w,p,2]]-prerun[[bs,w,p,2]],prerun[[b,1,p,2]]-prerun[[bs,1,p,2]],prerun[[b,2,p,2]]-prerun[[bs,2,p,2]]];
];
];
];
];
];
ndm1=ndr-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
derivs[[b,w,1]]=prerun[[b,w,ndm1-1,1]];derivs[[b,w,2]]=prerun[[b,w,ndm1-1,2]];derivs[[b,w,3]]=prerun[[b,w,ndm1-1,3]];derivs[[b,w,4]]=prerun[[b,w,ndm1-1,4]];
derivs[[b,w,5]]=(prerun[[b,w,ndm1-1,4]]-prerun[[b,w,ndm1-2,4]])/ht;
derivs[[b,w,6]]=(prerun[[b,w,ndm1-1,4]]-2prerun[[b,w,ndm1-2,4]]+prerun[[b,w,ndm1-3,4]])/(ht^2);
derivs[[b,w,7]]=(prerun[[b,w,ndm1-1,4]]-3prerun[[b,w,ndm1-2,4]]+3prerun[[b,w,ndm1-3,4]]-prerun[[b,w,ndm1-4,4]])/(ht^3);
derivs[[b,w,8]]=(prerun[[b,w,ndm1-1,4]]-4prerun[[b,w,ndm1-2,4]]+6prerun[[b,w,ndm1-3,4]]-4prerun[[b,w,ndm1-4,4]]+prerun[[b,w,ndm1-5,4]])/(ht^4);derivs[[b,w,9]]=(prerun[[b,w,ndm1-1,4]]-5prerun[[b,w,ndm1-2,4]]+10prerun[[b,w,ndm1-3,4]]-10prerun[[b,w,ndm1-4,4]]+5prerun[[b,w,ndm1-5,4]]-prerun[[b,w,ndm1-6,4]])/(ht^5);
derivs[[b,w,10]]=(prerun[[b,w,ndm1-1,4]]-6prerun[[b,w,ndm1-2,4]]+15prerun[[b,w,ndm1-3,4]]-20prerun[[b,w,ndm1-4,4]]+15prerun[[b,w,ndm1-5,4]]-6prerun[[b,w,ndm1-6,4]]+prerun[[b,w,ndm1-7,4]])/(ht^6);
derivs[[b,w,11]]=(prerun[[b,w,ndm1-1,4]]-7prerun[[b,w,ndm1-2,4]]+21prerun[[b,w,ndm1-3,4]]-35prerun[[b,w,ndm1-4,4]]+35prerun[[b,w,ndm1-5,4]]-21prerun[[b,w,ndm1-6,4]]+7prerun[[b,w,ndm1-7,4]]-prerun[[b,w,ndm1-8,4]])/(ht^7);
derivs[[b,w,12]]=(prerun[[b,w,ndm1-1,4]]-8prerun[[b,w,ndm1-2,4]]+28prerun[[b,w,ndm1-3,4]]-56prerun[[b,w,ndm1-4,4]]+70prerun[[b,w,ndm1-5,4]]-56prerun[[b,w,ndm1-6,4]]+28prerun[[b,w,ndm1-7,4]]-8prerun[[b,w,ndm1-8,4]]+prerun[[b,w,ndm1-9,4]])/(ht^8);
derivs[[b,w,13]]=(prerun[[b,w,ndm1-1,4]]-9prerun[[b,w,ndm1-2,4]]+36prerun[[b,w,ndm1-3,4]]-84prerun[[b,w,ndm1-4,4]]+126prerun[[b,w,ndm1-5,4]]-126prerun[[b,w,ndm1-6,4]]+84prerun[[b,w,ndm1-7,4]]-36prerun[[b,w,ndm1-8,4]]+9prerun[[b,w,ndm1-9,4]]-prerun[[b,w,ndm1-10,4]])/(ht^9);
];
];
derivs=SetPrecision[derivs,31];
dt=SetPrecision[permerc/ndiv-(ndr-3)ht,31];
tt+=dt;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
posn[[b,w]]=Posd[b,w,dt];velo[[b,w]]=Veld[b,w,dt];
];
];
If[posn[[2,1]]-posn[[1,1]]>0.98Px1&&posn[[2,2]]-posn[[1,2]]>Py1-0.1Px1&&posn[[2,2]]-posn[[1,2]]<Py1+0.1Px1,
awyvel0=VelMerc[0];awyvel1=VelMerc[dt];
If[firstperi==True||(awyvel0<0&&awyvel1>0),
crosdt2=-awyvel0 dt/(awyvel1-awyvel0);
awyvel2=VelMerc[crosdt2];
If[awyvel2<0,
crosdt3=crosdt2-awyvel2( dt-crosdt2)/(awyvel1-awyvel2);
,
crosdt3=-awyvel0 crosdt2/(awyvel2-awyvel0);
];
perang=AngArcsec[crosdt3];
crostim=tt-dt+crosdt3;
AppendTo[PeriAng,{SetPrecision[crostim/OneYear,5],SetPrecision[perang,5]}];
];
firstperi=False;
];
count+=1;
];]

(*Venus*)ListPlot[PeriAng,Joined->True]

Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Earth*)ListPlot[PeriAng,Joined->True]

Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Mars*)ListPlot[PeriAng,Joined->True]

Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Jupiter*)ListPlot[PeriAng[[1;;434]],Joined->True](*104.3y*)

Normal[LinearModelFit[PeriAng[[1;;434]],x,x]][[2]]100/x

(*Saturn*)ListPlot[PeriAng[[1;;463]],Joined->True](*111.25y*)

Normal[LinearModelFit[PeriAng[[1;;463]],x,x]][[2]](*x100=7.6*)

(*Uranus*)ListPlot[PeriAng,Joined->True](*230y*)

Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Neptune*)ListPlot[PeriAng[[1;;855]],Joined->True](*205.7y*)

Normal[LinearModelFit[PeriAng[[1;;855]],x,x]][[2]]100/x

(*Total = 292.9+95.9+2.4+157.0+7.6+0.14+0.04 = 556.0 *)

(*Simulate Mercury orbiting a floating Sun and an outer planet in an eccentric orbit.  The outer planet can be selected by changing the variable pno, e.g. pno=2 is Venus and pno=8 is Neptune.  
The outer planet starts at its perihelion point.  It's initial velocity will be boosted to ensure the period is correct.  Runtime approx 30min.*)
Timing[Ff2[mas_,wp_,rp3_]:=-G mas wp/rp3;Ff3[mas_,wp_,rp3_,xp_,yp_,wv_,xv_,yv_]:=-G mas(wv/rp3-3 wp (xp xv+yp yv)/(rp3^(5/3)));Posd[bp_,wp_,tp_]:=derivs[[bp,wp,1]]+derivs[[bp,wp,2]] tp+derivs[[bp,wp,3]](tp^2)/2+derivs[[bp,wp,4]](tp^3)/6+derivs[[bp,wp,5]](tp^4)/24+derivs[[bp,wp,6]](tp^5)/120+derivs[[bp,wp,7]](tp^6)/720+derivs[[bp,wp,8]](tp^7)/5040+derivs[[bp,wp,9]](tp^8)/40320+derivs[[bp,wp,10]](tp^9)/362880+derivs[[bp,wp,11]](tp^10)/3628800+derivs[[bp,wp,12]](tp^11)/39916800+derivs[[bp,wp,13]](tp^12)/479001600;Veld[bp_,wp_,tp_]:=derivs[[bp,wp,2]] +derivs[[bp,wp,3]]tp+derivs[[bp,wp,4]](tp^2)/2+derivs[[bp,wp,5]](tp^3)/6+derivs[[bp,wp,6]](tp^4)/24+derivs[[bp,wp,7]](tp^5)/120+derivs[[bp,wp,8]](tp^6)/720+derivs[[bp,wp,9]](tp^7)/5040+derivs[[bp,wp,10]](tp^8)/40320+derivs[[bp,wp,11]](tp^9)/362880+derivs[[bp,wp,12]](tp^10)/3628800+derivs[[bp,wp,13]](tp^11)/39916800;
VelMerc[tm_]:=((Posd[2,1,tm]-Posd[1,1,tm])(Veld[2,1,tm]-Veld[1,1,tm])+(Posd[2,2,tm]-Posd[1,2,tm])(Veld[2,2,tm]-Veld[1,2,tm]))/√((((Posd[2,1,tm]-Posd[1,1,tm]))^2)+(((Posd[2,2,tm]-Posd[1,2,tm]))^2));
AngArcsec[tm_]:=ArcTan[(Posd[2,1,tm]-Posd[1,1,tm]),(Posd[2,2,tm]-Posd[1,2,tm])]/π 180 3600;
pno=2;G=SetPrecision[Gg,180];{M0,M1,M2}=SetPrecision[{Msun,Mmercury,Pmass[[pno]]},180];
{R1,R2}=SetPrecision[{RPmercury,Pperihel[[pno]]},180];e1=SetPrecision[Emercury,180];e2=SetPrecision[Peccen[[pno]],180];
nbods=3;frbod=1;tobod=nbods;ndr=12;PeriAng={};firstperi=True;
{Px1,Px2}={R1,R2};{Py1,Py2}={0,0};
{Vx1,Vx2}={0,0};{Vy1,Vy2}={SetPrecision[PveladjP[[1]],180](Sqrt[(1+e1)G M0/R1]),SetPrecision[PveladjP[[pno]],180](Sqrt[(1+e2)(G M0)/R2])};
permerc=2π(Sqrt[(((R1/(1-e1)))^3)/(G M0)]);
posn={{0,0},{Px1,Py1},{ Px2,Py2}};velo={{0,0},{Vx1,Vy1},{Vx2,Vy2}};masses={M0,M1,M2};
prerun=Array[0&,{nbods,2,ndr-2,4}];derivs=Array[0&,{nbods,2,ndr+1}];
count=0;ndiv=132;dt=permerc/ndiv;ht=SetPrecision[(10^(-17))*permerc/(2π),180];tt=0;
While[tt<100OneYear,
posn=SetPrecision[posn,180];velo=SetPrecision[velo,180];
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
prerun[[b,w,1,1]]=posn[[b,w]];
prerun[[b,w,1,2]]=velo[[b,w]];
prerun[[b,w,1,3]]=0;
prerun[[b,w,1,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((posn[[b,1]]-posn[[bs,1]]))^2)+(((posn[[b,2]]-posn[[bs,2]]))^2)))^(3/2));
prerun[[b,w,1,3]]+=Ff2[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3];
prerun[[b,w,1,4]]+=Ff3[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3,posn[[b,1]]-posn[[bs,1]],posn[[b,2]]-posn[[bs,2]],velo[[b,w]]-velo[[bs,w]],velo[[b,1]]-velo[[bs,1]],velo[[b,2]]-velo[[bs,2]]];
];
];
];
];
For[p=2,p≤ndr-2,p++,
tt+=ht;
ps=p-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
Ff2C=0;
Ff3C=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]]))^2)+(((prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]]))^2)))^(3/2));
Ff2C+=Ff2[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3];
Ff3C+=Ff3[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3,prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]],prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]],prerun[[b,w,ps,2]]-prerun[[bs,w,ps,2]],prerun[[b,1,ps,2]]-prerun[[bs,1,ps,2]],prerun[[b,2,ps,2]]-prerun[[bs,2,ps,2]]];
];
];
prerun[[b,w,p,1]]=prerun[[b,w,ps,1]]+prerun[[b,w,ps,2]]ht+Ff2C (ht^2)/2+Ff3C (ht^3)/6;
prerun[[b,w,p,2]]=prerun[[b,w,ps,2]]+Ff2C ht+Ff3C (ht^2)/2;
];
];
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
prerun[[b,w,p,3]]=0;
prerun[[b,w,p,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,p,1]]-prerun[[bs,1,p,1]]))^2)+(((prerun[[b,2,p,1]]-prerun[[bs,2,p,1]]))^2)))^(3/2));
prerun[[b,w,p,3]]+=Ff2[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3];prerun[[b,w,p,4]]+=Ff3[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3,prerun[[b,1,p,1]]-prerun[[bs,1,p,1]],prerun[[b,2,p,1]]-prerun[[bs,2,p,1]],prerun[[b,w,p,2]]-prerun[[bs,w,p,2]],prerun[[b,1,p,2]]-prerun[[bs,1,p,2]],prerun[[b,2,p,2]]-prerun[[bs,2,p,2]]];
];
];
];
];
];
ndm1=ndr-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
derivs[[b,w,1]]=prerun[[b,w,ndm1-1,1]];derivs[[b,w,2]]=prerun[[b,w,ndm1-1,2]];derivs[[b,w,3]]=prerun[[b,w,ndm1-1,3]];derivs[[b,w,4]]=prerun[[b,w,ndm1-1,4]];
derivs[[b,w,5]]=(prerun[[b,w,ndm1-1,4]]-prerun[[b,w,ndm1-2,4]])/ht;
derivs[[b,w,6]]=(prerun[[b,w,ndm1-1,4]]-2prerun[[b,w,ndm1-2,4]]+prerun[[b,w,ndm1-3,4]])/(ht^2);
derivs[[b,w,7]]=(prerun[[b,w,ndm1-1,4]]-3prerun[[b,w,ndm1-2,4]]+3prerun[[b,w,ndm1-3,4]]-prerun[[b,w,ndm1-4,4]])/(ht^3);
derivs[[b,w,8]]=(prerun[[b,w,ndm1-1,4]]-4prerun[[b,w,ndm1-2,4]]+6prerun[[b,w,ndm1-3,4]]-4prerun[[b,w,ndm1-4,4]]+prerun[[b,w,ndm1-5,4]])/(ht^4);derivs[[b,w,9]]=(prerun[[b,w,ndm1-1,4]]-5prerun[[b,w,ndm1-2,4]]+10prerun[[b,w,ndm1-3,4]]-10prerun[[b,w,ndm1-4,4]]+5prerun[[b,w,ndm1-5,4]]-prerun[[b,w,ndm1-6,4]])/(ht^5);
derivs[[b,w,10]]=(prerun[[b,w,ndm1-1,4]]-6prerun[[b,w,ndm1-2,4]]+15prerun[[b,w,ndm1-3,4]]-20prerun[[b,w,ndm1-4,4]]+15prerun[[b,w,ndm1-5,4]]-6prerun[[b,w,ndm1-6,4]]+prerun[[b,w,ndm1-7,4]])/(ht^6);
derivs[[b,w,11]]=(prerun[[b,w,ndm1-1,4]]-7prerun[[b,w,ndm1-2,4]]+21prerun[[b,w,ndm1-3,4]]-35prerun[[b,w,ndm1-4,4]]+35prerun[[b,w,ndm1-5,4]]-21prerun[[b,w,ndm1-6,4]]+7prerun[[b,w,ndm1-7,4]]-prerun[[b,w,ndm1-8,4]])/(ht^7);
derivs[[b,w,12]]=(prerun[[b,w,ndm1-1,4]]-8prerun[[b,w,ndm1-2,4]]+28prerun[[b,w,ndm1-3,4]]-56prerun[[b,w,ndm1-4,4]]+70prerun[[b,w,ndm1-5,4]]-56prerun[[b,w,ndm1-6,4]]+28prerun[[b,w,ndm1-7,4]]-8prerun[[b,w,ndm1-8,4]]+prerun[[b,w,ndm1-9,4]])/(ht^8);
derivs[[b,w,13]]=(prerun[[b,w,ndm1-1,4]]-9prerun[[b,w,ndm1-2,4]]+36prerun[[b,w,ndm1-3,4]]-84prerun[[b,w,ndm1-4,4]]+126prerun[[b,w,ndm1-5,4]]-126prerun[[b,w,ndm1-6,4]]+84prerun[[b,w,ndm1-7,4]]-36prerun[[b,w,ndm1-8,4]]+9prerun[[b,w,ndm1-9,4]]-prerun[[b,w,ndm1-10,4]])/(ht^9);
];
];
derivs=SetPrecision[derivs,31];
dt=SetPrecision[permerc/ndiv-(ndr-3)ht,31];
tt+=dt;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
posn[[b,w]]=Posd[b,w,dt];velo[[b,w]]=Veld[b,w,dt];
];
];
If[posn[[2,1]]-posn[[1,1]]>0.98Px1&&posn[[2,2]]-posn[[1,2]]>Py1-0.1Px1&&posn[[2,2]]-posn[[1,2]]<Py1+0.1Px1,
awyvel0=VelMerc[0];awyvel1=VelMerc[dt];
If[firstperi==True||(awyvel0<0&&awyvel1>0),
crosdt2=-awyvel0 dt/(awyvel1-awyvel0);
awyvel2=VelMerc[crosdt2];
If[awyvel2<0,
crosdt3=crosdt2-awyvel2( dt-crosdt2)/(awyvel1-awyvel2);
,
crosdt3=-awyvel0 crosdt2/(awyvel2-awyvel0);
];
perang=AngArcsec[crosdt3];
crostim=tt-dt+crosdt3;
AppendTo[PeriAng,{SetPrecision[crostim/OneYear,5],SetPrecision[perang,5]}];
];
firstperi=False;
];
count+=1;
];]

(*Venus*)ListPlot[PeriAng,Joined->True]

Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Earth*)ListPlot[PeriAng,Joined->True]

Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Mars*)ListPlot[PeriAng,Joined->True]

Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Jupiter*)ListPlot[PeriAng[[1;;434]],Joined->True](*104.3y*)

Normal[LinearModelFit[PeriAng[[1;;434]],x,x]][[2]]100/x

(*Saturn*)ListPlot[PeriAng[[1;;463]],Joined->True](*111.3y*)

Normal[LinearModelFit[PeriAng[[1;;463]],x,x]][[2]]100/x

(*Uranus*)ListPlot[PeriAng[[1;;432]],Joined->True](*103.8y*)

Normal[LinearModelFit[PeriAng[[1;;432]],x,x]][[2]]100/x

(*Neptune*)ListPlot[PeriAng[[1;;516]],Joined->True](*124.0y*)

Normal[LinearModelFit[PeriAng[[1;;516]],x,x]][[2]]100/x

(*Total = 285.1+91.7+2.0+153.8+7.5+0.14+0.04 = 540.3 *)

(*Simulate Mercury orbiting a floating Sun and all outer planets from Venus to Neptune with eccentric orbits.  The orbits will be rotated to match the actual planet's orbits but won't be inclined.  Runtime approx 5.6hours*)
Timing[Ff2[mas_,wp_,rp3_]:=-G mas wp/rp3;Ff3[mas_,wp_,rp3_,xp_,yp_,wv_,xv_,yv_]:=-G mas(wv/rp3-3 wp (xp xv+yp yv)/(rp3^(5/3)));Posd[bp_,wp_,tp_]:=derivs[[bp,wp,1]]+derivs[[bp,wp,2]] tp+derivs[[bp,wp,3]](tp^2)/2+derivs[[bp,wp,4]](tp^3)/6+derivs[[bp,wp,5]](tp^4)/24+derivs[[bp,wp,6]](tp^5)/120+derivs[[bp,wp,7]](tp^6)/720+derivs[[bp,wp,8]](tp^7)/5040+derivs[[bp,wp,9]](tp^8)/40320+derivs[[bp,wp,10]](tp^9)/362880+derivs[[bp,wp,11]](tp^10)/3628800+derivs[[bp,wp,12]](tp^11)/39916800+derivs[[bp,wp,13]](tp^12)/479001600;Veld[bp_,wp_,tp_]:=derivs[[bp,wp,2]] +derivs[[bp,wp,3]]tp+derivs[[bp,wp,4]](tp^2)/2+derivs[[bp,wp,5]](tp^3)/6+derivs[[bp,wp,6]](tp^4)/24+derivs[[bp,wp,7]](tp^5)/120+derivs[[bp,wp,8]](tp^6)/720+derivs[[bp,wp,9]](tp^7)/5040+derivs[[bp,wp,10]](tp^8)/40320+derivs[[bp,wp,11]](tp^9)/362880+derivs[[bp,wp,12]](tp^10)/3628800+derivs[[bp,wp,13]](tp^11)/39916800;
VelMerc[tm_]:=((Posd[2,1,tm]-Posd[1,1,tm])(Veld[2,1,tm]-Veld[1,1,tm])+(Posd[2,2,tm]-Posd[1,2,tm])(Veld[2,2,tm]-Veld[1,2,tm]))/√((((Posd[2,1,tm]-Posd[1,1,tm]))^2)+(((Posd[2,2,tm]-Posd[1,2,tm]))^2));
OrbUntran[x_,y_,z_,ι_,Ω_,ω_]:={x(Cos[ι] Cos[ω] Cos[Ω]-Sin[ω] Sin[Ω])+y(Cos[Ω] Sin[ω]+Cos[ι] Cos[ω] Sin[Ω])+z(Cos[ω] Sin[ι]),x(-Cos[ι] Cos[Ω] Sin[ω]-Cos[ω] Sin[Ω])+y(Cos[ω] Cos[Ω]-Cos[ι] Sin[ω] Sin[Ω])+z(-Sin[ι] Sin[ω]),x(-Cos[Ω] Sin[ι])+y(-Sin[ι] Sin[Ω])+z(Cos[ι])};OrbUntranDeg[x_,y_,z_,ι_,Ω_,ω_]:=OrbUntran[x,y,z,ι π/180,Ω π/180,ω π/180];AngArcsec[tm_]:=ArcTan[OrbUntranDeg[(Posd[2,1,tm]-Posd[1,1,tm]),(Posd[2,2,tm]-Posd[1,2,tm]),0,0ιmercury,Ωmercury,ωmercury][[1]]
,OrbUntranDeg[(Posd[2,1,tm]-Posd[1,1,tm]),(Posd[2,2,tm]-Posd[1,2,tm]),0,0ιmercury,Ωmercury,ωmercury][[2]]]/π 180 3600;
G=SetPrecision[Gg,180];nbods=9;frbod=1;tobod=nbods;ndr=12;PeriAng={};firstperi=True;
{M0,M1,M2,M3,M4,M5,M6,M7,M8}=SetPrecision[{Msun,Mmercury,Mvenus,Mearth,Mmars,Mjupiter,Msaturn,Muranus,Mneptune},180];{R1,R2,R3,R4,R5,R6,R7,R8}=SetPrecision[{RPmercury,RPvenus,RPearth,RPmars,RPjupiter,RPsaturn,RPuranus,RPneptune},180];velint=Array[0&,8];
For[i=1,i≤8,i++,velint[[i]]=PveladjP[[i]](Sqrt[(1+ Peccen[[i]])((G Msun)/(Pperihel[[i]]))])];{V1,V2,V3,V4,V5,V6,V7,V8}=velint;
{Px1,Py1,Pz1}=OrbTransXYZdeg[R1,0,0,0ιmercury,Ωmercury,ωmercury];{Vx1,Vy1,Vz1}=OrbTransXYZdeg[0,V1,0,0ιmercury,Ωmercury,ωmercury];
{Px2,Py2,Pz2}=OrbTransXYZdeg[R2,0,0,0ιvenus,Ωvenus,ωvenus];{Vx2,Vy2,Vz2}=OrbTransXYZdeg[0,V2,0,0ιvenus,Ωvenus,ωvenus];
{Px3,Py3,Pz3}=OrbTransXYZdeg[R3,0,0,0ιearth,Ωearth,ωearth];{Vx3,Vy3,Vz3}=OrbTransXYZdeg[0,V3,0,0ιearth,Ωearth,ωearth];
{Px4,Py4,Pz4}=OrbTransXYZdeg[R4,0,0,0ιmars,Ωmars,ωmars];
{Vx4,Vy4,Vz4}=OrbTransXYZdeg[0,V4,0,0ιmars,Ωmars,ωmars];
{Px5,Py5,Pz5}=OrbTransXYZdeg[R5,0,0,0ιjupiter,Ωjupiter,ωjupiter];{Vx5,Vy5,Vz5}=OrbTransXYZdeg[0,V5,0,0ιjupiter,Ωjupiter,ωjupiter];
{Px6,Py6,Pz6}=OrbTransXYZdeg[R6,0,0,0ιsaturn,Ωsaturn,ωsaturn];{Vx6,Vy6,Vz6}=OrbTransXYZdeg[0,V6,0,0ιsaturn,Ωsaturn,ωsaturn];
{Px7,Py7,Pz7}=OrbTransXYZdeg[R7,0,0,0ιuranus,Ωuranus,ωuranus];{Vx7,Vy7,Vz7}=OrbTransXYZdeg[0,V7,0,0ιuranus,Ωuranus,ωuranus];
{Px8,Py8,Pz8}=OrbTransXYZdeg[R8,0,0,0ιneptune,Ωneptune,ωneptune];{Vx8,Vy8,Vz8}=OrbTransXYZdeg[0,V8,0,0ιneptune,Ωneptune,ωneptune];
permerc=2π(Sqrt[(((R1/(1-Emercury)))^3)/(G M0)]);
posn={{0,0},{Px1,Py1},{ Px2,Py2},{ Px3,Py3},{ Px4,Py4},{ Px5,Py5},{ Px6,Py6},{ Px7,Py7},{ Px8,Py8}};velo={{0,0},{Vx1,Vy1},{Vx2,Vy2},{Vx3,Vy3},{Vx4,Vy4},{Vx5,Vy5},{Vx6,Vy6},{Vx7,Vy7},{Vx8,Vy8}};masses={M0,M1,M2,M3,M4,M5,M6,M7,M8};
prerun=Array[0&,{nbods,2,ndr-2,4}];derivs=Array[0&,{nbods,2,ndr+1}];
count=0;ndiv=132;dt=permerc/ndiv;ht=SetPrecision[(10^(-17))*permerc/(2π),180];tt=0;
While[tt<100OneYear,
posn=SetPrecision[posn,180];velo=SetPrecision[velo,180];
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
prerun[[b,w,1,1]]=posn[[b,w]];
prerun[[b,w,1,2]]=velo[[b,w]];
prerun[[b,w,1,3]]=0;
prerun[[b,w,1,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((posn[[b,1]]-posn[[bs,1]]))^2)+(((posn[[b,2]]-posn[[bs,2]]))^2)))^(3/2));
prerun[[b,w,1,3]]+=Ff2[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3];
prerun[[b,w,1,4]]+=Ff3[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3,posn[[b,1]]-posn[[bs,1]],posn[[b,2]]-posn[[bs,2]],velo[[b,w]]-velo[[bs,w]],velo[[b,1]]-velo[[bs,1]],velo[[b,2]]-velo[[bs,2]]];
];
];
];
];
For[p=2,p≤ndr-2,p++,
tt+=ht;
ps=p-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
Ff2C=0;
Ff3C=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]]))^2)+(((prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]]))^2)))^(3/2));
Ff2C+=Ff2[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3];
Ff3C+=Ff3[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3,prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]],prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]],prerun[[b,w,ps,2]]-prerun[[bs,w,ps,2]],prerun[[b,1,ps,2]]-prerun[[bs,1,ps,2]],prerun[[b,2,ps,2]]-prerun[[bs,2,ps,2]]];
];
];
prerun[[b,w,p,1]]=prerun[[b,w,ps,1]]+prerun[[b,w,ps,2]]ht+Ff2C (ht^2)/2+Ff3C (ht^3)/6;
prerun[[b,w,p,2]]=prerun[[b,w,ps,2]]+Ff2C ht+Ff3C (ht^2)/2;
];
];
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
prerun[[b,w,p,3]]=0;
prerun[[b,w,p,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,p,1]]-prerun[[bs,1,p,1]]))^2)+(((prerun[[b,2,p,1]]-prerun[[bs,2,p,1]]))^2)))^(3/2));
prerun[[b,w,p,3]]+=Ff2[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3];prerun[[b,w,p,4]]+=Ff3[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3,prerun[[b,1,p,1]]-prerun[[bs,1,p,1]],prerun[[b,2,p,1]]-prerun[[bs,2,p,1]],prerun[[b,w,p,2]]-prerun[[bs,w,p,2]],prerun[[b,1,p,2]]-prerun[[bs,1,p,2]],prerun[[b,2,p,2]]-prerun[[bs,2,p,2]]];
];
];
];
];
];
ndm1=ndr-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
derivs[[b,w,1]]=prerun[[b,w,ndm1-1,1]];derivs[[b,w,2]]=prerun[[b,w,ndm1-1,2]];derivs[[b,w,3]]=prerun[[b,w,ndm1-1,3]];derivs[[b,w,4]]=prerun[[b,w,ndm1-1,4]];
derivs[[b,w,5]]=(prerun[[b,w,ndm1-1,4]]-prerun[[b,w,ndm1-2,4]])/ht;
derivs[[b,w,6]]=(prerun[[b,w,ndm1-1,4]]-2prerun[[b,w,ndm1-2,4]]+prerun[[b,w,ndm1-3,4]])/(ht^2);
derivs[[b,w,7]]=(prerun[[b,w,ndm1-1,4]]-3prerun[[b,w,ndm1-2,4]]+3prerun[[b,w,ndm1-3,4]]-prerun[[b,w,ndm1-4,4]])/(ht^3);
derivs[[b,w,8]]=(prerun[[b,w,ndm1-1,4]]-4prerun[[b,w,ndm1-2,4]]+6prerun[[b,w,ndm1-3,4]]-4prerun[[b,w,ndm1-4,4]]+prerun[[b,w,ndm1-5,4]])/(ht^4);derivs[[b,w,9]]=(prerun[[b,w,ndm1-1,4]]-5prerun[[b,w,ndm1-2,4]]+10prerun[[b,w,ndm1-3,4]]-10prerun[[b,w,ndm1-4,4]]+5prerun[[b,w,ndm1-5,4]]-prerun[[b,w,ndm1-6,4]])/(ht^5);
derivs[[b,w,10]]=(prerun[[b,w,ndm1-1,4]]-6prerun[[b,w,ndm1-2,4]]+15prerun[[b,w,ndm1-3,4]]-20prerun[[b,w,ndm1-4,4]]+15prerun[[b,w,ndm1-5,4]]-6prerun[[b,w,ndm1-6,4]]+prerun[[b,w,ndm1-7,4]])/(ht^6);
derivs[[b,w,11]]=(prerun[[b,w,ndm1-1,4]]-7prerun[[b,w,ndm1-2,4]]+21prerun[[b,w,ndm1-3,4]]-35prerun[[b,w,ndm1-4,4]]+35prerun[[b,w,ndm1-5,4]]-21prerun[[b,w,ndm1-6,4]]+7prerun[[b,w,ndm1-7,4]]-prerun[[b,w,ndm1-8,4]])/(ht^7);
derivs[[b,w,12]]=(prerun[[b,w,ndm1-1,4]]-8prerun[[b,w,ndm1-2,4]]+28prerun[[b,w,ndm1-3,4]]-56prerun[[b,w,ndm1-4,4]]+70prerun[[b,w,ndm1-5,4]]-56prerun[[b,w,ndm1-6,4]]+28prerun[[b,w,ndm1-7,4]]-8prerun[[b,w,ndm1-8,4]]+prerun[[b,w,ndm1-9,4]])/(ht^8);
derivs[[b,w,13]]=(prerun[[b,w,ndm1-1,4]]-9prerun[[b,w,ndm1-2,4]]+36prerun[[b,w,ndm1-3,4]]-84prerun[[b,w,ndm1-4,4]]+126prerun[[b,w,ndm1-5,4]]-126prerun[[b,w,ndm1-6,4]]+84prerun[[b,w,ndm1-7,4]]-36prerun[[b,w,ndm1-8,4]]+9prerun[[b,w,ndm1-9,4]]-prerun[[b,w,ndm1-10,4]])/(ht^9);
];
];
derivs=SetPrecision[derivs,31];
dt=SetPrecision[permerc/ndiv-(ndr-3)ht,31];
tt+=dt;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤2,w++,
posn[[b,w]]=Posd[b,w,dt];velo[[b,w]]=Veld[b,w,dt];
];
];
If[posn[[2,2]]-posn[[1,2]]>0.98(Py1-0)&&posn[[2,1]]-posn[[1,1]]>(Px1-0)-0.1(Py1-0)&&posn[[2,1]]-posn[[1,1]]<(Px1-0)+0.1(Py1-0),
awyvel0=VelMerc[0];awyvel1=VelMerc[dt];
If[firstperi==True||(awyvel0<0&&awyvel1>0),
crosdt2=-awyvel0 dt/(awyvel1-awyvel0);
awyvel2=VelMerc[crosdt2];
If[awyvel2<0,
crosdt3=crosdt2-awyvel2( dt-crosdt2)/(awyvel1-awyvel2);
,
crosdt3=-awyvel0 crosdt2/(awyvel2-awyvel0);
];
perang=AngArcsec[crosdt3];
crostim=tt-dt+crosdt3;
AppendTo[PeriAng,{SetPrecision[crostim/OneYear,5],SetPrecision[perang,5]}];
];
firstperi=False;
];
count+=1;(*AppendTo[posn2,{posn[[2,1]],posn[[2,2]]}];*)
];]

ListPlot[PeriAng,Joined->True]

Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x

(*Simulate Mercury orbiting a floating Sun and all outer planets from Venus to Neptune with eccentric orbits.  The orbits will be rotated and inclined to match the actual planets' orbits.  This is a 3D situation.  Runtime approx 17.5 hours*)Timing[Ff2[mas_,wp_,rp3_]:=-G mas wp/rp3;Ff3[mas_,wp_,rp3_,xp_,yp_,zp_,wv_,xv_,yv_,zv_]:=-G mas(wv/rp3-3 wp (xp xv+yp yv+zp zv)/(rp3^(5/3)));
Posd[bp_,wp_,tp_]:=derivs[[bp,wp,1]]+derivs[[bp,wp,2]] tp+derivs[[bp,wp,3]](tp^2)/2+derivs[[bp,wp,4]](tp^3)/6+derivs[[bp,wp,5]](tp^4)/24+derivs[[bp,wp,6]](tp^5)/120+derivs[[bp,wp,7]](tp^6)/720+derivs[[bp,wp,8]](tp^7)/5040+derivs[[bp,wp,9]](tp^8)/40320+derivs[[bp,wp,10]](tp^9)/362880+derivs[[bp,wp,11]](tp^10)/3628800+derivs[[bp,wp,12]](tp^11)/39916800+derivs[[bp,wp,13]](tp^12)/479001600;Veld[bp_,wp_,tp_]:=derivs[[bp,wp,2]] +derivs[[bp,wp,3]]tp+derivs[[bp,wp,4]](tp^2)/2+derivs[[bp,wp,5]](tp^3)/6+derivs[[bp,wp,6]](tp^4)/24+derivs[[bp,wp,7]](tp^5)/120+derivs[[bp,wp,8]](tp^6)/720+derivs[[bp,wp,9]](tp^7)/5040+derivs[[bp,wp,10]](tp^8)/40320+derivs[[bp,wp,11]](tp^9)/362880+derivs[[bp,wp,12]](tp^10)/3628800+derivs[[bp,wp,13]](tp^11)/39916800;
VelMerc[tm_]:=((Posd[2,1,tm]-Posd[1,1,tm])(Veld[2,1,tm]-Veld[1,1,tm])+(Posd[2,2,tm]-Posd[1,2,tm])(Veld[2,2,tm]-Veld[1,2,tm])+(Posd[2,3,tm]-Posd[1,3,tm])(Veld[2,3,tm]-Veld[1,3,tm]))/√((((Posd[2,1,tm]-Posd[1,1,tm]))^2)+(((Posd[2,2,tm]-Posd[1,2,tm]))^2)+(((Posd[2,3,tm]-Posd[1,3,tm]))^2));
OrbUntran[x_,y_,z_,ι_,Ω_,ω_]:={x(Cos[ι] Cos[ω] Cos[Ω]-Sin[ω] Sin[Ω])+y(Cos[Ω] Sin[ω]+Cos[ι] Cos[ω] Sin[Ω])+z(Cos[ω] Sin[ι]),x(-Cos[ι] Cos[Ω] Sin[ω]-Cos[ω] Sin[Ω])+y(Cos[ω] Cos[Ω]-Cos[ι] Sin[ω] Sin[Ω])+z(-Sin[ι] Sin[ω]),x(-Cos[Ω] Sin[ι])+y(-Sin[ι] Sin[Ω])+z(Cos[ι])};OrbUntranDeg[x_,y_,z_,ι_,Ω_,ω_]:=OrbUntran[x,y,z,ι π/180,Ω π/180,ω π/180];AngArcsec[tm_]:=ArcTan[OrbUntranDeg[(Posd[2,1,tm]-Posd[1,1,tm]),(Posd[2,2,tm]-Posd[1,2,tm]),(Posd[2,3,tm]-Posd[1,3,tm]),ιmercury,Ωmercury,ωmercury][[1]]
,OrbUntranDeg[(Posd[2,1,tm]-Posd[1,1,tm]),(Posd[2,2,tm]-Posd[1,2,tm]),(Posd[2,3,tm]-Posd[1,3,tm]),ιmercury,Ωmercury,ωmercury][[2]]]/π 180 3600;G=SetPrecision[Gg,180];nbods=9;frbod=1;tobod=nbods;ndr=12;PeriAng={};firstperi=True;
{M0,M1,M2,M3,M4,M5,M6,M7,M8}=SetPrecision[{Msun,Mmercury,Mvenus,Mearth,Mmars,Mjupiter,Msaturn,Muranus,Mneptune},180];{R1,R2,R3,R4,R5,R6,R7,R8}=SetPrecision[{RPmercury,RPvenus,RPearth,RPmars,RPjupiter,RPsaturn,RPuranus,RPneptune},180];velint=Array[0&,8];
For[i=1,i≤8,i++,velint[[i]]=PveladjP[[i]](Sqrt[(1+ Peccen[[i]])((G Msun)/(Pperihel[[i]]))])];{V1,V2,V3,V4,V5,V6,V7,V8}=velint;
{Px1,Py1,Pz1}=OrbTransXYZdeg[R1,0,0,ιmercury,Ωmercury,ωmercury];{Vx1,Vy1,Vz1}=OrbTransXYZdeg[0,V1,0,ιmercury,Ωmercury,ωmercury];
{Px2,Py2,Pz2}=OrbTransXYZdeg[R2,0,0,ιvenus,Ωvenus,ωvenus];{Vx2,Vy2,Vz2}=OrbTransXYZdeg[0,V2,0,ιvenus,Ωvenus,ωvenus];
{Px3,Py3,Pz3}=OrbTransXYZdeg[R3,0,0,ιearth,Ωearth,ωearth];{Vx3,Vy3,Vz3}=OrbTransXYZdeg[0,V3,0,ιearth,Ωearth,ωearth];
{Px4,Py4,Pz4}=OrbTransXYZdeg[R4,0,0,ιmars,Ωmars,ωmars];
{Vx4,Vy4,Vz4}=OrbTransXYZdeg[0,V4,0,ιmars,Ωmars,ωmars];
{Px5,Py5,Pz5}=OrbTransXYZdeg[R5,0,0,ιjupiter,Ωjupiter,ωjupiter];{Vx5,Vy5,Vz5}=OrbTransXYZdeg[0,V5,0,ιjupiter,Ωjupiter,ωjupiter];
{Px6,Py6,Pz6}=OrbTransXYZdeg[R6,0,0,ιsaturn,Ωsaturn,ωsaturn];{Vx6,Vy6,Vz6}=OrbTransXYZdeg[0,V6,0,ιsaturn,Ωsaturn,ωsaturn];
{Px7,Py7,Pz7}=OrbTransXYZdeg[R7,0,0,ιuranus,Ωuranus,ωuranus];{Vx7,Vy7,Vz7}=OrbTransXYZdeg[0,V7,0,ιuranus,Ωuranus,ωuranus];
{Px8,Py8,Pz8}=OrbTransXYZdeg[R8,0,0,ιneptune,Ωneptune,ωneptune];{Vx8,Vy8,Vz8}=OrbTransXYZdeg[0,V8,0,ιneptune,Ωneptune,ωneptune];
posn={{ 0,0,0},{Px1,Py1,Pz1},{ Px2,Py2,Pz2},{ Px3,Py3,Pz3},{ Px4,Py4,Pz4},{ Px5,Py5,Pz5},{ Px6,Py6,Pz6},{ Px7,Py7,Pz7},{ Px8,Py8,Pz8}};velo={{0,0,0},{Vx1,Vy1,Vz1},{Vx2,Vy2,Vz2},{Vx3,Vy3,Vz3},{Vx4,Vy4,Vz4},{Vx5,Vy5,Vz5},{Vx6,Vy6,Vz6},{Vx7,Vy7,Vz7},{Vx8,Vy8,Vz8}};masses={M0,M1,M2,M3,M4,M5,M6,M7,M8};
prerun=Array[0&,{nbods,3,ndr-2,4}];derivs=Array[0&,{nbods,3,ndr+1}];
permerc=SetPrecision[2π(Sqrt[(((RPmercury/(1-Emercury)))^3)/(G Msun)]),180];count=0;ndiv=132;dt=permerc/ndiv;ht=SetPrecision[(10^(-17))*permerc/(2π),180];tt=0;
While[tt<200OneYear,
posn=SetPrecision[posn,180];velo=SetPrecision[velo,180];
For[b=frbod,b≤tobod,b++,
For[w=1,w≤3,w++,
prerun[[b,w,1,1]]=posn[[b,w]];
prerun[[b,w,1,2]]=velo[[b,w]];
prerun[[b,w,1,3]]=0;
prerun[[b,w,1,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((posn[[b,1]]-posn[[bs,1]]))^2)+(((posn[[b,2]]-posn[[bs,2]]))^2)+(((posn[[b,3]]-posn[[bs,3]]))^2)))^(3/2));
prerun[[b,w,1,3]]+=Ff2[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3];
prerun[[b,w,1,4]]+=Ff3[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3,posn[[b,1]]-posn[[bs,1]],posn[[b,2]]-posn[[bs,2]],posn[[b,3]]-posn[[bs,3]],velo[[b,w]]-velo[[bs,w]],velo[[b,1]]-velo[[bs,1]],velo[[b,2]]-velo[[bs,2]],velo[[b,3]]-velo[[bs,3]]];
];
];
];
];
For[p=2,p≤ndr-2,p++,
tt+=ht;
ps=p-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤3,w++,
Ff2C=0;
Ff3C=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]]))^2)+(((prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]]))^2)+(((prerun[[b,3,ps,1]]-prerun[[bs,3,ps,1]]))^2)))^(3/2));
Ff2C+=Ff2[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3];
Ff3C+=Ff3[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3,prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]],prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]],prerun[[b,3,ps,1]]-prerun[[bs,3,ps,1]],prerun[[b,w,ps,2]]-prerun[[bs,w,ps,2]],prerun[[b,1,ps,2]]-prerun[[bs,1,ps,2]],prerun[[b,2,ps,2]]-prerun[[bs,2,ps,2]],prerun[[b,3,ps,2]]-prerun[[bs,3,ps,2]]];
];
];
prerun[[b,w,p,1]]=prerun[[b,w,ps,1]]+prerun[[b,w,ps,2]]ht+Ff2C (ht^2)/2+Ff3C (ht^3)/6;
prerun[[b,w,p,2]]=prerun[[b,w,ps,2]]+Ff2C ht+Ff3C (ht^2)/2;
];
];
For[b=frbod,b≤tobod,b++,
For[w=1,w≤3,w++,
prerun[[b,w,p,3]]=0;
prerun[[b,w,p,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,p,1]]-prerun[[bs,1,p,1]]))^2)+(((prerun[[b,2,p,1]]-prerun[[bs,2,p,1]]))^2)+(((prerun[[b,3,p,1]]-prerun[[bs,3,p,1]]))^2)))^(3/2));
prerun[[b,w,p,3]]+=Ff2[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3];prerun[[b,w,p,4]]+=Ff3[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3,prerun[[b,1,p,1]]-prerun[[bs,1,p,1]],prerun[[b,2,p,1]]-prerun[[bs,2,p,1]],prerun[[b,3,p,1]]-prerun[[bs,3,p,1]],prerun[[b,w,p,2]]-prerun[[bs,w,p,2]],prerun[[b,1,p,2]]-prerun[[bs,1,p,2]],prerun[[b,2,p,2]]-prerun[[bs,2,p,2]],prerun[[b,3,p,2]]-prerun[[bs,3,p,2]]];
];
];
];
];
];
ndm1=ndr-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤3,w++,
derivs[[b,w,1]]=prerun[[b,w,ndm1-1,1]];derivs[[b,w,2]]=prerun[[b,w,ndm1-1,2]];derivs[[b,w,3]]=prerun[[b,w,ndm1-1,3]];derivs[[b,w,4]]=prerun[[b,w,ndm1-1,4]];
derivs[[b,w,5]]=(prerun[[b,w,ndm1-1,4]]-prerun[[b,w,ndm1-2,4]])/ht;
derivs[[b,w,6]]=(prerun[[b,w,ndm1-1,4]]-2prerun[[b,w,ndm1-2,4]]+prerun[[b,w,ndm1-3,4]])/(ht^2);
derivs[[b,w,7]]=(prerun[[b,w,ndm1-1,4]]-3prerun[[b,w,ndm1-2,4]]+3prerun[[b,w,ndm1-3,4]]-prerun[[b,w,ndm1-4,4]])/(ht^3);
derivs[[b,w,8]]=(prerun[[b,w,ndm1-1,4]]-4prerun[[b,w,ndm1-2,4]]+6prerun[[b,w,ndm1-3,4]]-4prerun[[b,w,ndm1-4,4]]+prerun[[b,w,ndm1-5,4]])/(ht^4);derivs[[b,w,9]]=(prerun[[b,w,ndm1-1,4]]-5prerun[[b,w,ndm1-2,4]]+10prerun[[b,w,ndm1-3,4]]-10prerun[[b,w,ndm1-4,4]]+5prerun[[b,w,ndm1-5,4]]-prerun[[b,w,ndm1-6,4]])/(ht^5);
derivs[[b,w,10]]=(prerun[[b,w,ndm1-1,4]]-6prerun[[b,w,ndm1-2,4]]+15prerun[[b,w,ndm1-3,4]]-20prerun[[b,w,ndm1-4,4]]+15prerun[[b,w,ndm1-5,4]]-6prerun[[b,w,ndm1-6,4]]+prerun[[b,w,ndm1-7,4]])/(ht^6);
derivs[[b,w,11]]=(prerun[[b,w,ndm1-1,4]]-7prerun[[b,w,ndm1-2,4]]+21prerun[[b,w,ndm1-3,4]]-35prerun[[b,w,ndm1-4,4]]+35prerun[[b,w,ndm1-5,4]]-21prerun[[b,w,ndm1-6,4]]+7prerun[[b,w,ndm1-7,4]]-prerun[[b,w,ndm1-8,4]])/(ht^7);
derivs[[b,w,12]]=(prerun[[b,w,ndm1-1,4]]-8prerun[[b,w,ndm1-2,4]]+28prerun[[b,w,ndm1-3,4]]-56prerun[[b,w,ndm1-4,4]]+70prerun[[b,w,ndm1-5,4]]-56prerun[[b,w,ndm1-6,4]]+28prerun[[b,w,ndm1-7,4]]-8prerun[[b,w,ndm1-8,4]]+prerun[[b,w,ndm1-9,4]])/(ht^8);
derivs[[b,w,13]]=(prerun[[b,w,ndm1-1,4]]-9prerun[[b,w,ndm1-2,4]]+36prerun[[b,w,ndm1-3,4]]-84prerun[[b,w,ndm1-4,4]]+126prerun[[b,w,ndm1-5,4]]-126prerun[[b,w,ndm1-6,4]]+84prerun[[b,w,ndm1-7,4]]-36prerun[[b,w,ndm1-8,4]]+9prerun[[b,w,ndm1-9,4]]-prerun[[b,w,ndm1-10,4]])/(ht^9);
];
];
derivs=SetPrecision[derivs,31];
dt=SetPrecision[permerc/ndiv-(ndr-3)ht,31];
tt+=dt;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤3,w++,
posn[[b,w]]=Posd[b,w,dt];velo[[b,w]]=Veld[b,w,dt];
];
];
If[posn[[2,2]]-posn[[1,2]]>0.98(Py1-0)&&posn[[2,1]]-posn[[1,1]]>(Px1-0)-0.1(Py1-0)&&posn[[2,1]]-posn[[1,1]]<(Px1-0)+0.1(Py1-0),
awyvel0=VelMerc[0];awyvel1=VelMerc[dt];
If[firstperi==True||(awyvel0<0&&awyvel1>0),
crosdt2=-awyvel0 dt/(awyvel1-awyvel0);
awyvel2=VelMerc[crosdt2];
If[awyvel2<0,
crosdt3=crosdt2-awyvel2( dt-crosdt2)/(awyvel1-awyvel2);
,
crosdt3=-awyvel0 crosdt2/(awyvel2-awyvel0);
];
perang=AngArcsec[crosdt3];
crostim=tt-dt+crosdt3;
AppendTo[PeriAng,{SetPrecision[crostim/OneYear,5],SetPrecision[perang,5]}];
];
firstperi=False;
];
count+=1;
(*print year counter to see progress*)If[Mod[count,548]==0,Print[SetPrecision[tt/OneYear*54816/54800,4]]];
];]

ListPlot[PeriAng[[1;;829]],Joined->True]

Normal[LinearModelFit[PeriAng[[1;;829]],x,x]][[2]]100/x

(*Simulate Mercury orbiting a floating Sun and all outer planets from Venus to Neptune with eccentric orbits.  The orbits will be rotated and inclined to match the actual planets' orbits.  This is the same as above but this time the planets start at their aphelion position.*)Timing[Ff2[mas_,wp_,rp3_]:=-G mas wp/rp3;Ff3[mas_,wp_,rp3_,xp_,yp_,zp_,wv_,xv_,yv_,zv_]:=-G mas(wv/rp3-3 wp (xp xv+yp yv+zp zv)/(rp3^(5/3)));
Posd[bp_,wp_,tp_]:=derivs[[bp,wp,1]]+derivs[[bp,wp,2]] tp+derivs[[bp,wp,3]](tp^2)/2+derivs[[bp,wp,4]](tp^3)/6+derivs[[bp,wp,5]](tp^4)/24+derivs[[bp,wp,6]](tp^5)/120+derivs[[bp,wp,7]](tp^6)/720+derivs[[bp,wp,8]](tp^7)/5040+derivs[[bp,wp,9]](tp^8)/40320+derivs[[bp,wp,10]](tp^9)/362880+derivs[[bp,wp,11]](tp^10)/3628800+derivs[[bp,wp,12]](tp^11)/39916800+derivs[[bp,wp,13]](tp^12)/479001600;Veld[bp_,wp_,tp_]:=derivs[[bp,wp,2]] +derivs[[bp,wp,3]]tp+derivs[[bp,wp,4]](tp^2)/2+derivs[[bp,wp,5]](tp^3)/6+derivs[[bp,wp,6]](tp^4)/24+derivs[[bp,wp,7]](tp^5)/120+derivs[[bp,wp,8]](tp^6)/720+derivs[[bp,wp,9]](tp^7)/5040+derivs[[bp,wp,10]](tp^8)/40320+derivs[[bp,wp,11]](tp^9)/362880+derivs[[bp,wp,12]](tp^10)/3628800+derivs[[bp,wp,13]](tp^11)/39916800;
VelMerc[tm_]:=((Posd[2,1,tm]-Posd[1,1,tm])(Veld[2,1,tm]-Veld[1,1,tm])+(Posd[2,2,tm]-Posd[1,2,tm])(Veld[2,2,tm]-Veld[1,2,tm])+(Posd[2,3,tm]-Posd[1,3,tm])(Veld[2,3,tm]-Veld[1,3,tm]))/√((((Posd[2,1,tm]-Posd[1,1,tm]))^2)+(((Posd[2,2,tm]-Posd[1,2,tm]))^2)+(((Posd[2,3,tm]-Posd[1,3,tm]))^2));
OrbUntran[x_,y_,z_,ι_,Ω_,ω_]:={x(Cos[ι] Cos[ω] Cos[Ω]-Sin[ω] Sin[Ω])+y(Cos[Ω] Sin[ω]+Cos[ι] Cos[ω] Sin[Ω])+z(Cos[ω] Sin[ι]),x(-Cos[ι] Cos[Ω] Sin[ω]-Cos[ω] Sin[Ω])+y(Cos[ω] Cos[Ω]-Cos[ι] Sin[ω] Sin[Ω])+z(-Sin[ι] Sin[ω]),x(-Cos[Ω] Sin[ι])+y(-Sin[ι] Sin[Ω])+z(Cos[ι])};OrbUntranDeg[x_,y_,z_,ι_,Ω_,ω_]:=OrbUntran[x,y,z,ι π/180,Ω π/180,ω π/180];AngArcsec[tm_]:=ArcTan[OrbUntranDeg[(Posd[2,1,tm]-Posd[1,1,tm]),(Posd[2,2,tm]-Posd[1,2,tm]),(Posd[2,3,tm]-Posd[1,3,tm]),ιmercury,Ωmercury,ωmercury][[1]]
,OrbUntranDeg[(Posd[2,1,tm]-Posd[1,1,tm]),(Posd[2,2,tm]-Posd[1,2,tm]),(Posd[2,3,tm]-Posd[1,3,tm]),ιmercury,Ωmercury,ωmercury][[2]]]/π 180 3600;G=SetPrecision[Gg,180];nbods=9;frbod=1;tobod=nbods;ndr=12;PeriAng={};firstperi=True;
{M0,M1,M2,M3,M4,M5,M6,M7,M8}=SetPrecision[{Msun,Mmercury,Mvenus,Mearth,Mmars,Mjupiter,Msaturn,Muranus,Mneptune},180];posint=Array[0&,8];
For[i=1,i≤8,i++,posint[[i]]=((1+Peccen[[i]])/(1-Peccen[[i]]))Pperihel[[i]]];{R1,R2,R3,R4,R5,R6,R7,R8}=posint;For[i=1,i≤8,i++,velint[[i]]=PveladjA[[i]](Sqrt[(1- Peccen[[i]])((G Msun)/(posint[[i]]))])];{V1,V2,V3,V4,V5,V6,V7,V8}=velint;R1=RPmercury;V1=PveladjP[[1]](Sqrt[(1+ Emercury)((G Msun)/R1)]);
{Px1,Py1,Pz1}=OrbTransXYZdeg[R1,0,0,ιmercury,Ωmercury,ωmercury];{Vx1,Vy1,Vz1}=OrbTransXYZdeg[0,V1,0,ιmercury,Ωmercury,ωmercury];
{Px2,Py2,Pz2}=OrbTransXYZdeg[R2,0,0,ιvenus,Ωvenus,ωvenus+180];{Vx2,Vy2,Vz2}=OrbTransXYZdeg[0,V2,0,ιvenus,Ωvenus,ωvenus+180];
{Px3,Py3,Pz3}=OrbTransXYZdeg[R3,0,0,ιearth,Ωearth,ωearth+180];{Vx3,Vy3,Vz3}=OrbTransXYZdeg[0,V3,0,ιearth,Ωearth,ωearth+180];
{Px4,Py4,Pz4}=OrbTransXYZdeg[R4,0,0,ιmars,Ωmars,ωmars+180];
{Vx4,Vy4,Vz4}=OrbTransXYZdeg[0,V4,0,ιmars,Ωmars,ωmars+180];
{Px5,Py5,Pz5}=OrbTransXYZdeg[R5,0,0,ιjupiter,Ωjupiter,ωjupiter+180];{Vx5,Vy5,Vz5}=OrbTransXYZdeg[0,V5,0,ιjupiter,Ωjupiter,ωjupiter+180];
{Px6,Py6,Pz6}=OrbTransXYZdeg[R6,0,0,ιsaturn,Ωsaturn,ωsaturn+180];{Vx6,Vy6,Vz6}=OrbTransXYZdeg[0,V6,0,ιsaturn,Ωsaturn,ωsaturn+180];
{Px7,Py7,Pz7}=OrbTransXYZdeg[R7,0,0,ιuranus,Ωuranus,ωuranus+180];{Vx7,Vy7,Vz7}=OrbTransXYZdeg[0,V7,0,ιuranus,Ωuranus,ωuranus+180];
{Px8,Py8,Pz8}=OrbTransXYZdeg[R8,0,0,ιneptune,Ωneptune,ωneptune+180];{Vx8,Vy8,Vz8}=OrbTransXYZdeg[0,V8,0,ιneptune,Ωneptune,ωneptune+180];
posn={{ 0,0,0},{Px1,Py1,Pz1},{ Px2,Py2,Pz2},{ Px3,Py3,Pz3},{ Px4,Py4,Pz4},{ Px5,Py5,Pz5},{ Px6,Py6,Pz6},{ Px7,Py7,Pz7},{ Px8,Py8,Pz8}};velo={{0,0,0},{Vx1,Vy1,Vz1},{Vx2,Vy2,Vz2},{Vx3,Vy3,Vz3},{Vx4,Vy4,Vz4},{Vx5,Vy5,Vz5},{Vx6,Vy6,Vz6},{Vx7,Vy7,Vz7},{Vx8,Vy8,Vz8}};masses={M0,M1,M2,M3,M4,M5,M6,M7,M8};
prerun=Array[0&,{nbods,3,ndr-2,4}];derivs=Array[0&,{nbods,3,ndr+1}];
permerc=SetPrecision[2π(Sqrt[(((RPmercury/(1-Emercury)))^3)/(G Msun)]),180];count=0;ndiv=132;dt=permerc/ndiv;ht=SetPrecision[(10^(-17))*permerc/(2π),180];tt=0;
While[tt<200OneYear,
posn=SetPrecision[posn,180];velo=SetPrecision[velo,180];
For[b=frbod,b≤tobod,b++,
For[w=1,w≤3,w++,
prerun[[b,w,1,1]]=posn[[b,w]];
prerun[[b,w,1,2]]=velo[[b,w]];
prerun[[b,w,1,3]]=0;
prerun[[b,w,1,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((posn[[b,1]]-posn[[bs,1]]))^2)+(((posn[[b,2]]-posn[[bs,2]]))^2)+(((posn[[b,3]]-posn[[bs,3]]))^2)))^(3/2));
prerun[[b,w,1,3]]+=Ff2[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3];
prerun[[b,w,1,4]]+=Ff3[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3,posn[[b,1]]-posn[[bs,1]],posn[[b,2]]-posn[[bs,2]],posn[[b,3]]-posn[[bs,3]],velo[[b,w]]-velo[[bs,w]],velo[[b,1]]-velo[[bs,1]],velo[[b,2]]-velo[[bs,2]],velo[[b,3]]-velo[[bs,3]]];
];
];
];
];
For[p=2,p≤ndr-2,p++,
tt+=ht;
ps=p-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤3,w++,
Ff2C=0;
Ff3C=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]]))^2)+(((prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]]))^2)+(((prerun[[b,3,ps,1]]-prerun[[bs,3,ps,1]]))^2)))^(3/2));
Ff2C+=Ff2[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3];
Ff3C+=Ff3[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3,prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]],prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]],prerun[[b,3,ps,1]]-prerun[[bs,3,ps,1]],prerun[[b,w,ps,2]]-prerun[[bs,w,ps,2]],prerun[[b,1,ps,2]]-prerun[[bs,1,ps,2]],prerun[[b,2,ps,2]]-prerun[[bs,2,ps,2]],prerun[[b,3,ps,2]]-prerun[[bs,3,ps,2]]];
];
];
prerun[[b,w,p,1]]=prerun[[b,w,ps,1]]+prerun[[b,w,ps,2]]ht+Ff2C (ht^2)/2+Ff3C (ht^3)/6;
prerun[[b,w,p,2]]=prerun[[b,w,ps,2]]+Ff2C ht+Ff3C (ht^2)/2;
];
];
For[b=frbod,b≤tobod,b++,
For[w=1,w≤3,w++,
prerun[[b,w,p,3]]=0;
prerun[[b,w,p,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,p,1]]-prerun[[bs,1,p,1]]))^2)+(((prerun[[b,2,p,1]]-prerun[[bs,2,p,1]]))^2)+(((prerun[[b,3,p,1]]-prerun[[bs,3,p,1]]))^2)))^(3/2));
prerun[[b,w,p,3]]+=Ff2[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3];prerun[[b,w,p,4]]+=Ff3[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3,prerun[[b,1,p,1]]-prerun[[bs,1,p,1]],prerun[[b,2,p,1]]-prerun[[bs,2,p,1]],prerun[[b,3,p,1]]-prerun[[bs,3,p,1]],prerun[[b,w,p,2]]-prerun[[bs,w,p,2]],prerun[[b,1,p,2]]-prerun[[bs,1,p,2]],prerun[[b,2,p,2]]-prerun[[bs,2,p,2]],prerun[[b,3,p,2]]-prerun[[bs,3,p,2]]];
];
];
];
];
];
ndm1=ndr-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤3,w++,
derivs[[b,w,1]]=prerun[[b,w,ndm1-1,1]];derivs[[b,w,2]]=prerun[[b,w,ndm1-1,2]];derivs[[b,w,3]]=prerun[[b,w,ndm1-1,3]];derivs[[b,w,4]]=prerun[[b,w,ndm1-1,4]];
derivs[[b,w,5]]=(prerun[[b,w,ndm1-1,4]]-prerun[[b,w,ndm1-2,4]])/ht;
derivs[[b,w,6]]=(prerun[[b,w,ndm1-1,4]]-2prerun[[b,w,ndm1-2,4]]+prerun[[b,w,ndm1-3,4]])/(ht^2);
derivs[[b,w,7]]=(prerun[[b,w,ndm1-1,4]]-3prerun[[b,w,ndm1-2,4]]+3prerun[[b,w,ndm1-3,4]]-prerun[[b,w,ndm1-4,4]])/(ht^3);
derivs[[b,w,8]]=(prerun[[b,w,ndm1-1,4]]-4prerun[[b,w,ndm1-2,4]]+6prerun[[b,w,ndm1-3,4]]-4prerun[[b,w,ndm1-4,4]]+prerun[[b,w,ndm1-5,4]])/(ht^4);derivs[[b,w,9]]=(prerun[[b,w,ndm1-1,4]]-5prerun[[b,w,ndm1-2,4]]+10prerun[[b,w,ndm1-3,4]]-10prerun[[b,w,ndm1-4,4]]+5prerun[[b,w,ndm1-5,4]]-prerun[[b,w,ndm1-6,4]])/(ht^5);
derivs[[b,w,10]]=(prerun[[b,w,ndm1-1,4]]-6prerun[[b,w,ndm1-2,4]]+15prerun[[b,w,ndm1-3,4]]-20prerun[[b,w,ndm1-4,4]]+15prerun[[b,w,ndm1-5,4]]-6prerun[[b,w,ndm1-6,4]]+prerun[[b,w,ndm1-7,4]])/(ht^6);
derivs[[b,w,11]]=(prerun[[b,w,ndm1-1,4]]-7prerun[[b,w,ndm1-2,4]]+21prerun[[b,w,ndm1-3,4]]-35prerun[[b,w,ndm1-4,4]]+35prerun[[b,w,ndm1-5,4]]-21prerun[[b,w,ndm1-6,4]]+7prerun[[b,w,ndm1-7,4]]-prerun[[b,w,ndm1-8,4]])/(ht^7);
derivs[[b,w,12]]=(prerun[[b,w,ndm1-1,4]]-8prerun[[b,w,ndm1-2,4]]+28prerun[[b,w,ndm1-3,4]]-56prerun[[b,w,ndm1-4,4]]+70prerun[[b,w,ndm1-5,4]]-56prerun[[b,w,ndm1-6,4]]+28prerun[[b,w,ndm1-7,4]]-8prerun[[b,w,ndm1-8,4]]+prerun[[b,w,ndm1-9,4]])/(ht^8);
derivs[[b,w,13]]=(prerun[[b,w,ndm1-1,4]]-9prerun[[b,w,ndm1-2,4]]+36prerun[[b,w,ndm1-3,4]]-84prerun[[b,w,ndm1-4,4]]+126prerun[[b,w,ndm1-5,4]]-126prerun[[b,w,ndm1-6,4]]+84prerun[[b,w,ndm1-7,4]]-36prerun[[b,w,ndm1-8,4]]+9prerun[[b,w,ndm1-9,4]]-prerun[[b,w,ndm1-10,4]])/(ht^9);
];
];
derivs=SetPrecision[derivs,31];
dt=SetPrecision[permerc/ndiv-(ndr-3)ht,31];
tt+=dt;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤3,w++,
posn[[b,w]]=Posd[b,w,dt];velo[[b,w]]=Veld[b,w,dt];
];
];
If[posn[[2,2]]-posn[[1,2]]>0.98(Py1-0)&&posn[[2,1]]-posn[[1,1]]>(Px1-0)-0.1(Py1-0)&&posn[[2,1]]-posn[[1,1]]<(Px1-0)+0.1(Py1-0),
awyvel0=VelMerc[0];awyvel1=VelMerc[dt];
If[firstperi==True||(awyvel0<0&&awyvel1>0),
crosdt2=-awyvel0 dt/(awyvel1-awyvel0);
awyvel2=VelMerc[crosdt2];
If[awyvel2<0,
crosdt3=crosdt2-awyvel2( dt-crosdt2)/(awyvel1-awyvel2);
,
crosdt3=-awyvel0 crosdt2/(awyvel2-awyvel0);
];
perang=AngArcsec[crosdt3];
crostim=tt-dt+crosdt3;
AppendTo[PeriAng,{SetPrecision[crostim/OneYear,5],SetPrecision[perang,5]}];
];
firstperi=False;
];
count+=1;
(*print year counter to see progress*)If[Mod[count,548]==0,Print[SetPrecision[tt/OneYear*54816/54800,4]]];
];]

ListPlot[PeriAng[[1;;820]],Joined->True]

Normal[LinearModelFit[PeriAng[[1;;820]],x,x]][[2]]100/x

(*Simulate solar system from Mercury to Neptune (and Sun) using HORIZONS data.  Data positions and velocities are relative to SSB at 15:14 on 2001-Feb-01.  Runtime approx 8.9hours*)Timing[Ff2[mas_,wp_,rp3_]:=-G mas wp/rp3;Ff3[mas_,wp_,rp3_,xp_,yp_,zp_,wv_,xv_,yv_,zv_]:=-G mas(wv/rp3-3 wp (xp xv+yp yv+zp zv)/(rp3^(5/3)));
Posd[bp_,wp_,tp_]:=derivs[[bp,wp,1]]+derivs[[bp,wp,2]] tp+derivs[[bp,wp,3]](tp^2)/2+derivs[[bp,wp,4]](tp^3)/6+derivs[[bp,wp,5]](tp^4)/24+derivs[[bp,wp,6]](tp^5)/120+derivs[[bp,wp,7]](tp^6)/720+derivs[[bp,wp,8]](tp^7)/5040+derivs[[bp,wp,9]](tp^8)/40320+derivs[[bp,wp,10]](tp^9)/362880+derivs[[bp,wp,11]](tp^10)/3628800+derivs[[bp,wp,12]](tp^11)/39916800+derivs[[bp,wp,13]](tp^12)/479001600;Veld[bp_,wp_,tp_]:=derivs[[bp,wp,2]] +derivs[[bp,wp,3]]tp+derivs[[bp,wp,4]](tp^2)/2+derivs[[bp,wp,5]](tp^3)/6+derivs[[bp,wp,6]](tp^4)/24+derivs[[bp,wp,7]](tp^5)/120+derivs[[bp,wp,8]](tp^6)/720+derivs[[bp,wp,9]](tp^7)/5040+derivs[[bp,wp,10]](tp^8)/40320+derivs[[bp,wp,11]](tp^9)/362880+derivs[[bp,wp,12]](tp^10)/3628800+derivs[[bp,wp,13]](tp^11)/39916800;
VelMerc[tm_]:=((Posd[2,1,tm]-Posd[1,1,tm])(Veld[2,1,tm]-Veld[1,1,tm])+(Posd[2,2,tm]-Posd[1,2,tm])(Veld[2,2,tm]-Veld[1,2,tm])+(Posd[2,3,tm]-Posd[1,3,tm])(Veld[2,3,tm]-Veld[1,3,tm]))/√((((Posd[2,1,tm]-Posd[1,1,tm]))^2)+(((Posd[2,2,tm]-Posd[1,2,tm]))^2)+(((Posd[2,3,tm]-Posd[1,3,tm]))^2));
OrbUntran[x_,y_,z_,ι_,Ω_,ω_]:={x(Cos[ι] Cos[ω] Cos[Ω]-Sin[ω] Sin[Ω])+y(Cos[Ω] Sin[ω]+Cos[ι] Cos[ω] Sin[Ω])+z(Cos[ω] Sin[ι]),x(-Cos[ι] Cos[Ω] Sin[ω]-Cos[ω] Sin[Ω])+y(Cos[ω] Cos[Ω]-Cos[ι] Sin[ω] Sin[Ω])+z(-Sin[ι] Sin[ω]),x(-Cos[Ω] Sin[ι])+y(-Sin[ι] Sin[Ω])+z(Cos[ι])};OrbUntranDeg[x_,y_,z_,ι_,Ω_,ω_]:=OrbUntran[x,y,z,ι π/180,Ω π/180,ω π/180];AngArcsec[tm_]:=ArcTan[OrbUntranDeg[(Posd[2,2,tm]-Posd[1,2,tm]),-(Posd[2,1,tm]-Posd[1,1,tm]),(Posd[2,3,tm]-Posd[1,3,tm]),ιmercury,Ωmercury,ωmercury-90][[1]]
,OrbUntranDeg[(Posd[2,2,tm]-Posd[1,2,tm]),-(Posd[2,1,tm]-Posd[1,1,tm]),(Posd[2,3,tm]-Posd[1,3,tm]),ιmercury,Ωmercury,ωmercury-90][[2]]]/π 180 3600;
G=SetPrecision[Gg,180];{M0,M1,M2,M3,M4,M5,M6,M7,M8}=SetPrecision[{Msun,Mmercury,Mvenus,Mearth,Mmars,Mjupiter,Msaturn,Muranus,Mneptune},180];
nbods=9;frbod=1;tobod=nbods;ndr=12;PeriAng={};firstperi=True;
{Px0,Py0,Pz0}={-6.544700998788778*^8,-7.542720159167981*^8,2.232295004897844*^7};
{Vx0,Vy0,Vz0}={1.430709291641341*^1,-6.129689956131971*^0,-3.231469219383437*^-1};
{Px1,Py1,Pz1}={9.459704017632110*^9,4.403743225392255*^10,2.753061478911016*^9};
{Vx1,Vy1,Vz1}={-5.726823818034963*^4,1.254517593507960*^4,6.282561277101745*^3};
{Px2,Py2,Pz2}={-1.521787736620838*^10,1.058313140886273*^11,2.319715498325974*^9};
{Vx2,Vy2,Vz2}={-3.480265957291748*^4,-4.938326945425076*^3,1.941907289679378*^3};
{Px3,Py3,Pz3}={-1.008728284577093*^11,1.073649125740838*^11,2.245984644959867*^7};
{Vx3,Vy3,Vz3}={-2.230964007501381*^4,-2.037499756826389*^4,-8.625517806555294*^-1};
{Px4,Py4,Pz4}={-2.341553529718963*^11,-6.893016826997817*^10,4.333407438725121*^9};
{Vx4,Vy4,Vz4}={7.714702908531727*^3,-2.119711899205197*^4,-6.334684794491539*^2};
{Px5,Py5,Pz5}={2.339841047987117*^11,7.181908392479482*^11,-8.215593369233102*^9};
{Vx5,Vy5,Vz5}={-1.258371434241999*^4,4.669884216770746*^3,2.622933336930586*^2};
{Px6,Py6,Pz6}={6.761216366578872*^11,1.180910647882344*^12,-4.745229184242111*^10};
{Vx6,Vy6,Vz6}={-8.889132514137824*^3,4.778912397853901*^3,2.702901855284778*^2};
{Px7,Py7,Pz7}={2.310818987293303*^12,-1.890798390154246*^12,-3.696547163527751*^10};
{Vx7,Vy7,Vz7}={4.262398872465041*^3,4.953701195921087*^3,-3.660014040687543*^1};
{Px8,Py8,Pz8}={2.665340229040296*^12,-3.630837682718371*^12,1.334489177487350*^10};
{Vx8,Vy8,Vz8}={4.346847466026690*^3,3.248046287448822*^3,-1.671845424812615*^2};
posn={{ Px0,Py0,Pz0},{Px1,Py1,Pz1},{ Px2,Py2,Pz2},{ Px3,Py3,Pz3},{ Px4,Py4,Pz4},{ Px5,Py5,Pz5},{ Px6,Py6,Pz6},{ Px7,Py7,Pz7},{ Px8,Py8,Pz8}};velo={{Vx0,Vy0,Vz0},{Vx1,Vy1,Vz1},{Vx2,Vy2,Vz2},{Vx3,Vy3,Vz3},{Vx4,Vy4,Vz4},{Vx5,Vy5,Vz5},{Vx6,Vy6,Vz6},{Vx7,Vy7,Vz7},{Vx8,Vy8,Vz8}};masses={M0,M1,M2,M3,M4,M5,M6,M7,M8};
prerun=Array[0&,{nbods,3,ndr-2,4}];derivs=Array[0&,{nbods,3,ndr+1}];
permerc=SetPrecision[2π(Sqrt[(((RPmercury/(1-Emercury)))^3)/(G Msun)]),180];count=0;ndiv=132;dt=permerc/ndiv;ht=SetPrecision[(10^(-17))*permerc/(2π),180];tt=0;
While[tt<100OneYear,
posn=SetPrecision[posn,180];velo=SetPrecision[velo,180];
For[b=frbod,b≤tobod,b++,
For[w=1,w≤3,w++,
prerun[[b,w,1,1]]=posn[[b,w]];
prerun[[b,w,1,2]]=velo[[b,w]];
prerun[[b,w,1,3]]=0;
prerun[[b,w,1,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((posn[[b,1]]-posn[[bs,1]]))^2)+(((posn[[b,2]]-posn[[bs,2]]))^2)+(((posn[[b,3]]-posn[[bs,3]]))^2)))^(3/2));
prerun[[b,w,1,3]]+=Ff2[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3];
prerun[[b,w,1,4]]+=Ff3[masses[[bs]],posn[[b,w]]-posn[[bs,w]],rdist3,posn[[b,1]]-posn[[bs,1]],posn[[b,2]]-posn[[bs,2]],posn[[b,3]]-posn[[bs,3]],velo[[b,w]]-velo[[bs,w]],velo[[b,1]]-velo[[bs,1]],velo[[b,2]]-velo[[bs,2]],velo[[b,3]]-velo[[bs,3]]];
];
];
];
];
For[p=2,p≤ndr-2,p++,
tt+=ht;
ps=p-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤3,w++,
Ff2C=0;
Ff3C=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]]))^2)+(((prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]]))^2)+(((prerun[[b,3,ps,1]]-prerun[[bs,3,ps,1]]))^2)))^(3/2));
Ff2C+=Ff2[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3];
Ff3C+=Ff3[masses[[bs]],prerun[[b,w,ps,1]]-prerun[[bs,w,ps,1]],rdist3,prerun[[b,1,ps,1]]-prerun[[bs,1,ps,1]],prerun[[b,2,ps,1]]-prerun[[bs,2,ps,1]],prerun[[b,3,ps,1]]-prerun[[bs,3,ps,1]],prerun[[b,w,ps,2]]-prerun[[bs,w,ps,2]],prerun[[b,1,ps,2]]-prerun[[bs,1,ps,2]],prerun[[b,2,ps,2]]-prerun[[bs,2,ps,2]],prerun[[b,3,ps,2]]-prerun[[bs,3,ps,2]]];
];
];
prerun[[b,w,p,1]]=prerun[[b,w,ps,1]]+prerun[[b,w,ps,2]]ht+Ff2C (ht^2)/2+Ff3C (ht^3)/6;
prerun[[b,w,p,2]]=prerun[[b,w,ps,2]]+Ff2C ht+Ff3C (ht^2)/2;
];
];
For[b=frbod,b≤tobod,b++,
For[w=1,w≤3,w++,
prerun[[b,w,p,3]]=0;
prerun[[b,w,p,4]]=0;
For[bs=1,bs≤nbods,bs++,
If[bs≠b,
rdist3=((((((prerun[[b,1,p,1]]-prerun[[bs,1,p,1]]))^2)+(((prerun[[b,2,p,1]]-prerun[[bs,2,p,1]]))^2)+(((prerun[[b,3,p,1]]-prerun[[bs,3,p,1]]))^2)))^(3/2));
prerun[[b,w,p,3]]+=Ff2[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3];prerun[[b,w,p,4]]+=Ff3[masses[[bs]],prerun[[b,w,p,1]]-prerun[[bs,w,p,1]],rdist3,prerun[[b,1,p,1]]-prerun[[bs,1,p,1]],prerun[[b,2,p,1]]-prerun[[bs,2,p,1]],prerun[[b,3,p,1]]-prerun[[bs,3,p,1]],prerun[[b,w,p,2]]-prerun[[bs,w,p,2]],prerun[[b,1,p,2]]-prerun[[bs,1,p,2]],prerun[[b,2,p,2]]-prerun[[bs,2,p,2]],prerun[[b,3,p,2]]-prerun[[bs,3,p,2]]];
];
];
];
];
];
ndm1=ndr-1;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤3,w++,
derivs[[b,w,1]]=prerun[[b,w,ndm1-1,1]];derivs[[b,w,2]]=prerun[[b,w,ndm1-1,2]];derivs[[b,w,3]]=prerun[[b,w,ndm1-1,3]];derivs[[b,w,4]]=prerun[[b,w,ndm1-1,4]];
derivs[[b,w,5]]=(prerun[[b,w,ndm1-1,4]]-prerun[[b,w,ndm1-2,4]])/ht;
derivs[[b,w,6]]=(prerun[[b,w,ndm1-1,4]]-2prerun[[b,w,ndm1-2,4]]+prerun[[b,w,ndm1-3,4]])/(ht^2);
derivs[[b,w,7]]=(prerun[[b,w,ndm1-1,4]]-3prerun[[b,w,ndm1-2,4]]+3prerun[[b,w,ndm1-3,4]]-prerun[[b,w,ndm1-4,4]])/(ht^3);
derivs[[b,w,8]]=(prerun[[b,w,ndm1-1,4]]-4prerun[[b,w,ndm1-2,4]]+6prerun[[b,w,ndm1-3,4]]-4prerun[[b,w,ndm1-4,4]]+prerun[[b,w,ndm1-5,4]])/(ht^4);derivs[[b,w,9]]=(prerun[[b,w,ndm1-1,4]]-5prerun[[b,w,ndm1-2,4]]+10prerun[[b,w,ndm1-3,4]]-10prerun[[b,w,ndm1-4,4]]+5prerun[[b,w,ndm1-5,4]]-prerun[[b,w,ndm1-6,4]])/(ht^5);
derivs[[b,w,10]]=(prerun[[b,w,ndm1-1,4]]-6prerun[[b,w,ndm1-2,4]]+15prerun[[b,w,ndm1-3,4]]-20prerun[[b,w,ndm1-4,4]]+15prerun[[b,w,ndm1-5,4]]-6prerun[[b,w,ndm1-6,4]]+prerun[[b,w,ndm1-7,4]])/(ht^6);
derivs[[b,w,11]]=(prerun[[b,w,ndm1-1,4]]-7prerun[[b,w,ndm1-2,4]]+21prerun[[b,w,ndm1-3,4]]-35prerun[[b,w,ndm1-4,4]]+35prerun[[b,w,ndm1-5,4]]-21prerun[[b,w,ndm1-6,4]]+7prerun[[b,w,ndm1-7,4]]-prerun[[b,w,ndm1-8,4]])/(ht^7);
derivs[[b,w,12]]=(prerun[[b,w,ndm1-1,4]]-8prerun[[b,w,ndm1-2,4]]+28prerun[[b,w,ndm1-3,4]]-56prerun[[b,w,ndm1-4,4]]+70prerun[[b,w,ndm1-5,4]]-56prerun[[b,w,ndm1-6,4]]+28prerun[[b,w,ndm1-7,4]]-8prerun[[b,w,ndm1-8,4]]+prerun[[b,w,ndm1-9,4]])/(ht^8);
derivs[[b,w,13]]=(prerun[[b,w,ndm1-1,4]]-9prerun[[b,w,ndm1-2,4]]+36prerun[[b,w,ndm1-3,4]]-84prerun[[b,w,ndm1-4,4]]+126prerun[[b,w,ndm1-5,4]]-126prerun[[b,w,ndm1-6,4]]+84prerun[[b,w,ndm1-7,4]]-36prerun[[b,w,ndm1-8,4]]+9prerun[[b,w,ndm1-9,4]]-prerun[[b,w,ndm1-10,4]])/(ht^9);
];
];
derivs=SetPrecision[derivs,31];
dt=SetPrecision[permerc/ndiv-(ndr-3)ht,31];
tt+=dt;
For[b=frbod,b≤tobod,b++,
For[w=1,w≤3,w++,
posn[[b,w]]=Posd[b,w,dt];velo[[b,w]]=Veld[b,w,dt];
];
];
If[posn[[2,2]]-posn[[1,2]]>0.98(Py1-Py0)&&posn[[2,1]]-posn[[1,1]]>(Px1-Px0)-0.1(Py1-Py0)&&posn[[2,1]]-posn[[1,1]]<(Px1-Px0)+0.1(Py1-Py0),
awyvel0=VelMerc[0];awyvel1=VelMerc[dt];
If[firstperi==True||(awyvel0<0&&awyvel1>0),
crosdt2=-awyvel0 dt/(awyvel1-awyvel0);
awyvel2=VelMerc[crosdt2];
If[awyvel2<0,
crosdt3=crosdt2-awyvel2( dt-crosdt2)/(awyvel1-awyvel2);
,
crosdt3=-awyvel0 crosdt2/(awyvel2-awyvel0);
];
perang=AngArcsec[crosdt3];
crostim=tt-dt+crosdt3;
AppendTo[PeriAng,{SetPrecision[crostim/OneYear,5],SetPrecision[perang,5]}];
];
firstperi=False;
];
count+=1;
(*print year counter to see progress*)If[Mod[count,548]==0,Print[SetPrecision[tt/OneYear*54816/54800,4]]];
];]

PeriAng(*HORIZONS*)

ListPlot[PeriAng,Joined->True]

Normal[LinearModelFit[PeriAng,x,x]][[2]]100/x