%kepler
% A visualization of Kepler's Laws or orbital motion. For simplicity, we'll
% consider a model of a planet about a star, whose mass is much greater
% than the mass of the planet. i.e. we can ignore the motion of the star
% about the barycentre of the system (which is effectively the centre of
% the star).
%
% LAST UPDATED by Andy French June 2021

function kepler

AU = 1.495978707e11;    %Astronomical unit /m
G = 6.674e-11;          %Universal gravitational constant /Nm^2kg^-2
Ms = 1.988e30;          %Solar mass /kg
Yr = 365.25*24*3600;    %Solar year /s

%Semi-major axis
a = 5*AU;

%Star mass /solar masses
M = 2*Ms;

%Planet mass /kg
m = M/1e6;

%Orbital eccentricity
ecc = 0.6;

%Orbit time interval fractions for Kepler II visualization
dt = 0.1;
t1 = 0.1; t2 = t1 + dt;
t3 = 0.45; t4 = t3 + dt;

%Fontsize for plots
fsize = 18;   

%Number of angle evaluation points
N = 10000;

%

%Determine orbital period /s
P = sqrt( (a^3)*(4*pi^2)/(G*(m+M)) );

%Determine total energy of system /J
E = -G*M*m/(2*a);

%Determine rate of area swepts (m^2 /s)
dA_by_dt = 0.5*sqrt( G*(m+M)*( 1-ecc^2)*a );

%Determine angular momentum J (kgm^2s^-1)
J = sqrt( G*(m^2)*(M^2)*(1-ecc^2)*a/(m+M) );

%Determine polar angle for elliptical orbit /radians
theta = linspace(0,2*pi,N);

%Determine radius of elliptical orbit
r = a*(1-ecc^2)./( 1 - ecc*cos(theta) );

%Determine orbit x,y coordinates (/m)
x = r.*cos(theta); y = r.*sin(theta);

%Determine orbit time /s
t = zeros(1,N);
dtheta = 2*pi/(N-1);
for n=2:N
    t(n) = t(n-1) + sqrt( ( ( a*(1-ecc^2) )^3 )/(G*(m+M)) )*dtheta*...
        ( 1 - ecc*cos( theta(n) ) )^(-2) ;
end

%Find angles which correspond to Kepler II visualization times
n1 = floor( interp1( t, 1:N, t1*P ) );
n2 = floor( interp1( t, 1:N, t2*P ) );
n3 = floor( interp1( t, 1:N, t3*P ) );
n4 = floor( interp1( t, 1:N, t4*P ) );

%Find radial and tangential orbital velocities (ms^-1)
rdot = -sqrt( G*(m+M)/(a*(1-ecc^2)) )*ecc*sin(theta);
thetadot = sqrt( G*(m+M)/(a*(1-ecc^2)) )*( 1 - ecc*cos(theta) )./r;

%Find x and y velocities
vx = rdot.*cos(theta) - r.*thetadot.*sin(theta);
vy = rdot.*sin(theta) + r.*thetadot.*cos(theta);
v = sqrt( vx.^2 + vy.^2 );

%

%Plot orbit
plot(x/AU,y/AU,'b-','linewidth',2); axis equal; hold on;
set(gca,'fontsize',fsize); grid on;
xlabel('x /AU'); ylabel('y /AU');
title(['a=',num2str(a/AU),'AU, M=',num2str(M/Ms),...
    ', \epsilon=',num2str(ecc),', P=',num2str(P/Yr,3),'years.']);

%Overlay Kepler II visualization
patch( [0,r(n1:n2).*cos(theta(n1:n2)),0]/AU,...
    [0,r(n1:n2).*sin(theta(n1:n2)),0]/AU, 'g','facealpha',0.5,'edgecolor','none' );
patch( [0,r(n3:n4).*cos(theta(n3:n4)),0]/AU,...
    [0,r(n3:n4).*sin(theta(n3:n4)),0]/AU, 'g','facealpha',0.5,'edgecolor','none'  );

%Overlay star
set(gcf,'renderer','zbuffer')
plot( 0,0,'r*','markersize',80);
plot( 0,0,'ro','markersize',50,'markerfacecolor','r' );

%Print PNG and close
print( gcf, 'kepler.png','-dpng','-r600');
clf;

%Plot orbital speed (km/s) vs orbit angle
plot( theta*180/pi, v/1e3, 'r-','linewidth',1 );
set(gca,'fontsize',fsize); grid on;
xlabel('\theta^o'); ylabel('Orbital speed /kms^{-1}'); xlim([0,360]);
title(['a=',num2str(a/AU),'AU, M=',num2str(M/Ms),...
    ', \epsilon=',num2str(ecc),', P=',num2str(P/Yr,3),'years.']);
print( gcf, 'kepler v vs theta.png','-dpng','-r600');
clf;

%Plot orbit time vs angle
plot( theta*180/pi, t/Yr, 'b-','linewidth',1 );
set(gca,'fontsize',fsize); grid on;
xlabel('\theta^o'); ylabel('Time /Yrs'); xlim([0,360]);
title(['a=',num2str(a/AU),'AU, M=',num2str(M/Ms),...
    ', \epsilon=',num2str(ecc),', P=',num2str(P/Yr,3),'years.']);
print( gcf, 'kepler t vs theta.png','-dpng','-r600');
close(gcf);

%End of code
