%pluto_and_charon
% Example solution of a Kepler two-body problem. Centre of mass of system
% is deemed to be fixed and set at the origon of the coordinate system.
%
%
% LAST UPDATED by Andy French October 2016

function pluto_and_charon

%%% Astronomical data %%%

%Universal Gravitational constant / m^3kg^-1s^-2
G = 6.67384e-11;

%Mass of Pluto /kg
M = 1.303e22;

%Radius of Pluto /m
Rp = 1187e3;

%Mass of Charon /kg
m = 1.586e21;

%Radius of Charon /m
Rc = 606e3;

%Semi-major axis of orbit /m
a = 19596e3;

%Eccentricity or orbit
ecc = 0.0;

%

%Number of data points
N = 1000;

%Fontsize for graphs
fsize = 18;

%Axis length scale /m
scale = a;

%Axis limits /m
xlimits = [-1 - ecc,1];
ylimits = [-1,1];

%

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

%Calculate orbital period in Earth days
Pdays = P/(24*3600);

%Define polar angles /radians
theta = linspace(0,2*pi,N);

%Define separation r /m
r = a*(1-ecc^2)./( 1-ecc*cos(theta) );

%x and y coordinates /m
xp = (m/(M+m))*r.*cos(theta);
xc = -(M/(M+m))*r.*cos(theta);
yp = (m/(M+m))*r.*sin(theta);
yc = -(M/(M+m))*r.*sin(theta);

%Velocity /ms^-1
Vx = (1-ecc*cos(theta))*sqrt( G*(m+M)/((1-ecc^2)*a) ).*...
    ( -sin(theta) - ( ecc*sin(theta)./(1 - ecc*cos(theta) )).*cos(theta) );
Vy = (1-ecc*cos(theta))*sqrt( G*(m+M)/((1-ecc^2)*a) ).*...
    ( cos(theta) - ( ecc*sin(theta)./(1 - ecc*cos(theta) )).*sin(theta) );
vxp = (m/(M+m))*Vx;
vyp = (m/(M+m))*Vy;
vxc = -(M/(M+m))*Vx;
vyc = -(M/(M+m))*Vy;

%Define orbital time /s over three orbits
t = sqrt( ( ( a*(1-ecc^2) )^3 )/( G*(m+M) ) ) * thetafunc( 0, linspace(0,3*2*pi,N), ecc );

%Plot orbits on an x,y graph
fig = figure('color',[1 1 1],'name','kepler');
plot(xc/scale,yc/scale,xp/scale,yp/scale);
legend({'Charon','Pluto'},'fontsize',fsize);
set(gca,'fontsize',fsize)
axis equal
xlim(xlimits);
ylim(ylimits);
grid on;
xlabel('x (10^7m)','fontsize',fsize)
ylabel('y (10^7m)','fontsize',fsize)
title('Orbits of Pluto and its moon Charon','fontsize',fsize)
print( gcf,'Pluto & Charon xy orbit.png','-dpng','-r300');

%Plot orbital angle /radians vs time /days
clf;
plot( t/(24*3600), linspace(0,3*2*pi,N) );
set(gca,'fontsize',fsize)
grid on;
xlabel('time /days','fontsize',fsize)
ylabel('Polar angle /radians','fontsize',fsize)
title(['Charon orbital period is ',num2str(Pdays),' days'],'fontsize',fsize)
hold on;
xlim([0,3*Pdays])
ylimits = get(gca,'ylim');
xlimits = get(gca,'xlim');
plot( [Pdays,Pdays],ylimits,'r');
plot( 2*[Pdays,Pdays],ylimits,'r');
plot( 3*[Pdays,Pdays],ylimits,'r');
plot( xlimits,2*pi*[1,1] ,'r' );
plot( xlimits,2*2*pi*[1,1],'r');
plot( xlimits,3*2*pi*[1,1],'r');
print( gcf,'angle vs time.png','-dpng','-r300');

%Plot orbital velocities /kms^-1 vs time
clf;
t = sqrt( ( ( a*(1-ecc^2) )^3 )/( G*(m+M) ) ) * thetafunc( 0, theta, ecc );
t = t/(24*3600);
plot(t,vxc/1000,t,vyc/1000,t,vxp/1000,t,vyp/1000);
legend({'Charon vx','Charon vy','Pluto vx','Pluto vy'},'fontsize',fsize);
set(gca,'fontsize',fsize)
grid on;
xlabel('time /days','fontsize',fsize)
ylabel('velocity /kms^{-1}','fontsize',fsize)
title('Orbital velocities of Pluto and its moon Charon','fontsize',fsize)
xlim([0,Pdays])
print( gcf,'Pluto & Charon vx vy vs t.png','-dpng','-r300');

%Plot separation vs time
clf;
t = sqrt( ( ( a*(1-ecc^2) )^3 )/( G*(m+M) ) ) * thetafunc( 0, theta, ecc );
t = t/(24*3600);
plot(t,r/scale);
set(gca,'fontsize',fsize)
grid on;
xlabel('time /days','fontsize',fsize)
ylabel('Separation (10^7m) ','fontsize',fsize)
title('Separation of Charon and Pluto vs time','fontsize',fsize)
xlim([0,Pdays])
print( gcf,'Pluto & Charon separation vs t.png','-dpng','-r300');
close(gcf);

%%


%Eliptic integral function
function f = thetafunc( theta0, theta, ecc )
y = 1./( 1 - ecc*cos(theta) ).^2;
f = integrator( theta, y, theta0, theta );

%%

%integrator
% Integration function using cubic splines.
%
% Syntax: I = integrator( xx, yy, x0, x )
%
% xx, yy     Samples from function y(x) which is being integrated. These
%            must be the same dimensions and have at least three elements each.
% x0         Starting point (i.e. lower limit) of the integral
% x          Upper limit of the integral. This can be an array of any
%            dimensions.
%
% I          Output of integral corresponindg to array x
%
% The function is based on the code written on pages 292 and 294
% of "Mastering Matlab 6" by Duane Hanselman & Bruce Littlefield.

function I = integrator( xx, yy, x0, x )

%Check inputs are the right size
if ( numel(xx) == numel(yy) ) && ( numel(xx) >= 3 )
    
    %Sort xx in ascending order
    [xx,i] = sort(xx);
    yy = yy(i);
    
    %Compute piecewise polynomial based upon cubic splines.
    pp=spline(xx,yy);
    
    %Take apart piecewise polynomial
    [br,co,npy,nco] = unmkpp(pp);
    
    %Scale factors for integration
    sf=nco:-1:1;
    
    %Integral coefficients
    ico=[co./sf(ones(npy,1),:) zeros(npy,1)];
    
    %Integral spline has higher order
    nco=nco+1;
    
    %Set integration constant of all splines to be zero
    ico(1,nco)=0;
    
    %Find constant terms in polynomials
    for k=2:npy
        ico(k,nco)=polyval(ico(k-1,:),br(k)-br(k-1));
    end
    
    %Build pp form for integral
    ppi=mkpp(br,ico);
    
    %Evaluate definite integral between limits of x0 to x
    I = ppval(ppi,x) - ppval(ppi,x0);
    
else
    I = NaN(size(x));
end

%End of code