clear all
%The code is part of my capstone project in the 2012 Distinguished
%Fulbright Awards in Teaching hosted by the University of Maryland.
%This code simulate the orbit of a planet around the sun according to
%Feynman (Chapter 9: Newton's laws of Dynamics, in Feynman Lectures on
%Physics (Feynman, Leighton and Sands) ,1963. Addison-Wesley Publishing
%Company, Reading Massachusetts.
%algorithm. with taking into account the mutual perturbations between Uranus and Neptune.
%The initial velocity is very important to give the right ellipse and that
%should be found by try and error as to give the right number in A.U
%We assume that the sun is fixed and located in (0,0) and that GM_s =1 as
%was assumed by Feynman to make the calculations simple.
%Note that is quite straight forward to add another planet (like Jupiter)
%but remember to the right terms in the equations for Neptune and
%Uranus.
%For comments and questions please contact:
%Hezi Yizhaq hezi.yizhaq1@gmail.com
%__________________________________________________________________________
%Initial conditions and parameter definitions
delta=6.2648;%delta is a parameter for adjusting the initial velocity.
delta1=6.5;
eps=0.04; %eps is the time step
eps2=eps/2; %eps2 is used for computing the first time step.
rr1=19.22; % rr1 is the initial distance of Uranus from the Sun.
rr2=30.8; % rr2 is the initial distance of Neptune from the Sun.
x1(1)=rr1; %x1 is the x coordinatee of Uranus.
y1(1)=0; % y1 is the y coordinate of Uranus.
x2(1)=rr2; %x2 is x coordinate of Neptune
y2(1)=0; %y2 is the y coordinate of Neptune
m1=4.373e-5; %m1 is the mass of Uranus (the ratio between Uranus mass and the Sun mass.
m2=5.178e-5; %m2 is the mass of Neptune (the ratio between Neptune mass and the Sun mass.
m12=m1*m2; %m12 is defined to make the calculation faster.
m22=1/m2; % m22 is needed in order to compute the ratio of the gravity forces.
tu=84.01;% Uranus revolution time (in years).
tn=164.8;% Neptune revolution time (in years).
%_______________________________________________________
% Calculating the first time step for Uranus
vx1(1)=0; %Initial velocity in the x direction.
vy1(1)=2*pi*rr1/(tu*delta); %Initial velocity in the y direction.
r1(1)=(x1(1)^2+y1(1)^2)^0.5;% calculating the initial distance from Neptune to the Sun.
r13(1)=1/(r1(1)^3);
run(1)=((x2(1)-x1(1))^2+(y2(1)-y1(1))^2)^0.5;%calculating the initial distance between Uranus and Neptune.
run3(1)=1/(run(1)^3);
vsqr1(1)=vx1(1)^2+vy1(1)^2; %calculation of the squared velocity of Uranus.
v1(1)=vsqr1(1)^0.5;%calculating the initial velocity of Uranus
r1(1)=(x1(1)^2+y1(1)^2)^0.5;
ax1(1)=-x1(1)*r13(1)-m2*(-x2(1)+x1(1))*run3(1);% the acceleration of Uranus in the x direction.
ay1(1)=-y1(1)*r13(1)-m2*(-y2(1)+y1(1))*run3(1);% the acceleration of Uranus in the y direction.
vvx1(1)=vx1(1)+ax1(1)*eps2;%these values are used for calculation of the next time step.
vvy1(1)=vy1(1)+ay1(1)*eps2;%these values are used for calculation of the next time step.
energy1(1)=-m1/r1(1)+0.5*m1*vsqr1(1);%the total initial energy of Uranus
j1(1)=m1*v1(1)*r1(1);%the total angular momentum of Uranus
%_______________________________________________________
%Calculating the first time step for Neptune
%the same explanations written above are applied fro Neptune.
vx2(1)=0;
vy2(1)=2*pi*rr2/(tn*delta1);
r2(1)=(x2(1)^2+y2(1)^2)^0.5;
r23(1)=1/(r2(1)^3);
vsqr2(1)=vx2(1)^2+vy2(1)^2;
v2(1)=vsqr2(1)^0.5;
r2(1)=(x2(1)^2+y2(1)^2)^0.5;
ax2(1)=-x2(1)*r23(1)-m1*(x2(1)-x1(1))*run3(1);;
ay2(1)=-y2(1)*r23(1)-m1*(y2(1)-y1(1))*run3(1);;
vvx2(1)=vx2(1)+ax2(1)*eps2;
vvy2(1)=vy2(1)+ay2(1)*eps2;
energy2(1)=-m2/r2(1)+0.5*m2*vsqr2(1);
%energy12 is the potential energy between Neptune and Uranus
energy12(1)=-m12/run(1);
energy_total(1)=energy1(1)+energy2(1)+energy12(1);% the total energy of the system: Uranus and Neptune.
j2(1)=m2*v2(1)*r2(1);%the angular momentum of Neptune.
f(1)=m22*(run(1)/r1(1))^2;% the ratio of the forces applied by the Sun to Uranus to that of Neptune on Uranus.
x1(2)=x1(1)+eps*vvx1(1);
y1(2)=y1(1)+eps*vvy1(1);
x2(2)=x2(1)+eps*vvx2(1);
y2(2)=y2(1)+eps*vvy2(1);
%end of calculating the first step
%Starting the main loop of the program
%n is the number of time steps
%__________________________________________________________________________
%n=53860;
n=40000;
for i=2:n
%calculating the distances necessary for the calculation.
r1(i)=(x1(i)^2+y1(i)^2)^0.5;
r13(i)=1/(r1(i)^3);
r2(i)=(x2(i)^2+y2(i)^2)^0.5;
r23(i)=1/(r2(i)^3);
run(i)=((x2(i)-x1(i))^2+(y2(i)-y1(i))^2)^0.5;
run3(i)=1/(run(i)^3);
%__________________________________________
%calculating acceleration components
ax1(i)=-x1(i)*r13(i)-m2*(-x2(i)+x1(i))*run3(i);
ay1(i)=-y1(i)*r13(i)-m2*(-y2(i)+y1(i))*run3(i);
ax2(i)=-x2(i)*r23(i)-m1*(x2(i)-x1(i))*run3(i);
ay2(i)=-y2(i)*r23(i)-m1*(y2(i)-y1(i))*run3(i);
%_______________________________________
%Note that the vvx1(i) are calculated in the middle of interval at times
% (0.5, 1.5. 2.5,...)*eps/2 whereas the x and y are computed at
% (1,2,3,...)*eps such as a time delay of eps/2 between velocity and x.
% This important when computing the total energy.
vvx1(i)=vvx1(i-1)+ax1(i)*eps;
vvy1(i)=vvy1(i-1)+ay1(i)*eps;
vvx2(i)=vvx2(i-1)+ax2(i)*eps;
vvy2(i)=vvy2(i-1)+ay2(i)*eps;
%_______________________________________
x1(i+1)=x1(i)+eps*vvx1(i);
y1(i+1)=y1(i)+eps*vvy1(i);
x2(i+1)=x2(i)+eps*vvx2(i);
y2(i+1)=y2(i)+eps*vvy2(i);
vsqr1(i)=vvx1(i)^2+vvy1(i)^2;
v1(i)=vsqr1(i)^0.5;
r1(i)=(x1(i)^2+y1(i)^2)^0.5;
vsqr2(i)=vvx2(i)^2+vvy2(i)^2;
v2(i)=vsqr2(i)^0.5;
r2(i)=(x2(i)^2+y2(i)^2)^0.5;
%_______________________________________
%calculating the total energy, total angular momentum and the ration between
%the gravity forces on Uranus.
energy1(i)=-m1/r1(i)+0.5*m1*vsqr1(i);
j1(i)=m1*v1(i)*r1(i);
energy2(i)=-m2/r2(i)+0.5*m2*vsqr2(i);
energy12(i)=-m12/run(i);
j2(i)=m2*v2(i)*r2(i);
energy_total(i)=energy1(i)+energy2(i)+energy12(i);
f(i)=m22*(run(i)/r1(i))^2;
% The lines below which are commented allow to plot the location of the
% planets in intermediate times.
% if i==(n/4)*j
% j=j+1;
% figure
% t=(0:eps:(n-1)*eps);
% plot(x1,y1,'k','LineWidth',2)
% grid on
% hold on
% plot(x2,y2,'b','LineWidth',2)
% hold on
% plot(x1(i),y1(i),'o')
% hold on
% plot(x2(i),y2(i),'o')
% xlabel('x [AU]')
% ylabel('y[AU]')
% end
end
%____________________________________________________________________
%the part of the code below plots the different results of the
%calculations
% The first figure plot the orbits of Uranus and Neptune.
figure
set(gca,'FontSize',16)
t=(0:eps:(n-1)*eps);
plot(x1,y1,'k','LineWidth',2)
grid on
hold on
plot(x2,y2,'b','LineWidth',2)
hold on
plot(x1(i),y1(i),'o')
hold on
plot(x2(i),y2(i),'o')
xlabel('x [AU]')
ylabel('y[AU]')
axis square
%___________________________________________
% This figure plot the total energy of the system and the total angular
% momentum of the system.
figure
set(gca,'FontSize',16)
subplot(1,2,1)
plot(t,energy_total,'r','LineWidth',2)
grid on
% hold on
% plot(t,energy2,'b','LineWidth',2)
xlabel('Time')
ylabel('Total Energy')
axis square
subplot(1,2,2)
set(gca,'FontSize',16)
plot(t,j1+j2,'r','LineWidth',2)
grid on
% hold on
% plot(t,j2,'b','LineWidth',2)
xlabel('Time')
ylabel('Angular Momentum')
axis square
%__________________________________________________________________
%This figure plot the orbital velocity of the planets as a function of the
%orbit radius
figure
set(gca,'FontSize',16)
plot(r1,v1,'k','LineWidth',2')
grid on
hold on
plot(r2,v2,'b','LineWidth',2')
xlabel('r')
ylabel('Orbital Velocity')
%______________________________________________________________________
%This figure plot the ratio between the gravity forces acting on Uranus by
%the Sun and by Neptune.
figure
set(gca,'FontSize',16)
plot(t,f,'b','LineWidth',2')
grid on
xlabel('Time')
ylabel('Ratio between the Gravity forces on Uranus')