(* Planet orbiting the Sun  *)
(* Assuming perfect circular orbit*)
(*Initial condition *)
R=aE=1;iini=0;GM=4Pi^2;delta=0.01;
T=1.0;
Deltat=delta*T;ilast=IntegerPart[T/Deltat];
x[0]=aE;y[0]=0;vx[0]=0;vy[0]=2Pi;

For[i=iini,i≤ilast,i++,
r[i]=Sqrt[(x[i]^2+y[i]^2)];
vx[i+1]=vx[i]-Deltat*(4Pi^2/r[i]^3)*x[i];
x[i+1]=x[i]+vx[i+1]*Deltat;
vy[i+1]=vy[i]-Deltat*(4Pi^2/r[i]^3)*y[i];
y[i+1]=y[i]+vy[i+1]*Deltat;
];(*end for i*)
txy=Table[{x[i],y[i]},{i,iini,ilast}];
ListPlot[txy,AspectRatio->1,PlotLabel->"locus of Earth around Sun ",AxesLabel->{"x (in A.U.)","y (in A.U) "}];


lp[i_]:=ListPlot[{{x[i],y[i]}},
PlotRange->{{-1.1aE,1.1aE},{-1.1aE,1.1aE}},
AspectRatio->1,PlotLabel->"locus of Earth around Sun ",AxesLabel->{"x (in A.U.)","y (in A.U) "}];

Print[Manipulate[lp[i],{i,iini,ilast,1}]];(*end of manipulate*)

(* Planet orbiting the Sun  *)
(* Assuming non-perfect circular orbit, i.e., elliptical *)
(*Initial condition *)
beta=1.0;gamma=120Degree;R=aE=1;iini=0;GM=4Pi^2;delta=0.001;T=1.0;
Deltat=delta*T;ilast=1400(*IntegerPart[T/Deltat]*);
x[0]=aE;y[0]=0;vx[0]=2Pi*Sin[gamma]*beta;
vy[0]=2Pi*Cos[gamma]*beta;
For[i=iini,i≤ilast,i++,
r[i]=Sqrt[(x[i]^2+y[i]^2)];vx[i+1]=vx[i]-Deltat*(4Pi^2/r[i]^3)*x[i];
x[i+1]=x[i]+vx[i+1]*Deltat;vy[i+1]=vy[i]-Deltat*(4Pi^2/r[i]^3)*y[i];
y[i+1]=y[i]+vy[i+1]*Deltat;];(*end for i*)
txy=Table[{x[i],y[i]},{i,iini,ilast}];
ListPlot[txy,AspectRatio->1,PlotLabel->"locus of Earth around Sun ",AxesLabel->{"x (in A.U.)","y (in A.U) "}]
lp[i_]:=ListPlot[{{x[i],y[i]}},PlotRange->{{-1.1aE,1.1aE},{-1.1aE,1.1aE}},
AspectRatio->1,PlotLabel->"locus of Earth around Sun ",AxesLabel->{"x (in A.U.)","y (in A.U) "}];
Print[Manipulate[lp[i],{i,iini,ilast,1}]];

(* Planet orbiting the Sun  *)
(* Assuming non-perfect circular orbit, i.e., elliptical *)
(* Improved initial conditions *)
(*Initial condition *)
beta=1.0;gamma=80Degree;
aE=1;iini=0;GM=4Pi^2;delta=0.0001;

(*Planet dependent parameters *)
a0=1.58aE;
(*Planet dependent parameters *)

T=a0^(1.5);Deltat=delta*T;
ilast=5*IntegerPart[T/Deltat];

x[0]=a0;y[0]=0;
v0=2Pi/a0^(1/2);
vx[0]=v0*Sin[gamma]*beta;
vy[0]=v0*Cos[gamma]*beta;

For[i=iini,i≤ilast,i++,
r[i]=Sqrt[(x[i]^2+y[i]^2)];
vx[i+1]=vx[i]-Deltat*(4Pi^2/r[i]^3)*x[i];
x[i+1]=x[i]+vx[i+1]*Deltat;vy[i+1]=vy[i]-Deltat*(4Pi^2/r[i]^3)*y[i];
y[i+1]=y[i]+vy[i+1]*Deltat;];(*end for i*)
txy=Table[{x[i],y[i]},{i,iini,ilast}];
ListPlot[txy,AspectRatio->1,PlotLabel->"locus of Earth around Sun ",AxesLabel->{"x (in A.U.)","y (in A.U) "}]
lp[i_]:=ListPlot[{{x[i],y[i]}},PlotRange->{{-1.1aE,1.1aE},{-1.1aE,1.1aE}},
AspectRatio->1,PlotLabel->"locus of Earth around Sun ",AxesLabel->{"x (in A.U.)","y (in A.U) "}];
Print[Manipulate[lp[i],{i,iini,ilast,1}]];