%pluto_and_charon_animation
% Model of two planets in mutually circulart orbits. In this case the dwarf
% planet Pluto and its largest moon Charon.
%
% LAST UPDATED by Andy French June 2021

function pluto_and_charon_animation

%Pluto mass in kg
M = 1.303e22;

%Charon mass in kg
m = 1.586e21; 

%Planet separation a /m
a = 1.9596e7;

%G (Universal Gravitational constant) /Nm^2*kg^-2
G = 6.67e-11;

%Fontsize for graphs
fsize = 18;

%

%Pluto orbit radius /m
R = a*m/(m+M);

%Charon radius /m
r = a*M/(m+M);

%Period from Kepler III /day
P = sqrt( (4*pi^2)* (a^3)/(G*(m+M)) );  %Period in s
Pdays = P/(24*3600); %Period in days

%Angular speed (rad/s)
w = sqrt( G*(m+M)/a^3 );

%

%Define angles 0 to 2*pi
theta = linspace(0,2*pi,100);

%Define times /s
t = theta/w;

%x,y cordinates of Pluto and Charon (/km)
X = R*cos(theta)/1e3; Y = R*sin(theta)/1e3;
x = r*cos(theta + pi)/1e3; y = -r*sin(theta)/1e3;

%Speeds (in m/s)
V = 2*pi*R/P; v = 2*pi*r/P;

%Plot orbits
axes('fontsize',fsize);
plot(X,Y,'r-',x,y,'b'); xlabel('x/km'); ylabel('y/km');
axis equal; hold on; grid on;
pluto = plot( X(1),Y(1),'r*');
charon = plot( x(1),y(1),'b*');
joining_line = plot( [x(1),X(1)], [y(1),Y(1) ],'k--' );
legend({'Pluto','Charon'})
title(['Pluto and Charon orbit. Period = ',num2str(Pdays,3),'days.']);

%Save a PNG file
print( gcf, 'pluto & charon.png','-dpng','-r300' );

%Animation!
stop = 0; n=1;
while stop==0
    n = n+1;
    if n>length(theta)
        n=1;
    end
    set( pluto, 'xdata', X(n), 'ydata', Y(n) );
    set( charon, 'xdata', x(n), 'ydata', y(n) );
    set( joining_line, 'xdata', [x(n),X(n)], 'ydata', [y(n),Y(n)] );
    drawnow
end

%End of code




