%binary_stars_and_planet
% Simulates the circular orbit of a planet around two massive stars.
% Ignore mass of the planet in the gravitational interaction between the
% stars.
%
% LAST UPDATED by Andy French 21/9/2023

function binary_stars_and_planet

%% INPUTS %%

%Semi-major axis of mutual star orbit in AU
a = 3; 

%Planet (initial) circular orbit radius about star 1
ap = 0.4*a; 

%Initial angle from x axis (anticlockwise) of planet /radians
theta0 = pi/4; 

%Masses of stars in solar masses
M1 = 5; M2 = 2;

%Initial vy velocity multiplier from mutually circular of stars
k = 1;  

%Number of orbital periods
num_periods = 30; 

%Timestep in years
dt = 0.01; 

%Fontsize
fsize = 18; 

%Axes limits
limit = 1.3*a; 

%Starting period for plot
Pstart = 1.23;

%

%% INITIAL CONDITIONS %%

%Determine period in years via Kepler III for binary stars
T = sqrt( (1/(M1+M2))*a^3 );

%Set initial x,y,vx,vy,t parameters

%Star 1
X1 = -M2*a/(M1+M2); Y1 = 0; VX1 = 0;  VY1 = k*2*pi*X1/T;

%Star 2
X2 = M1*a/(M1+M2); Y2 = 0; VX2 = 0;  VY2 = k*2*pi*X2/T;

%Planet
t = 0; x = X1 + ap*cos(theta0); y = ap*sin(theta0);
vx = -2*pi*sqrt( M1/ap)*sin(theta0); vy = 2*pi*sqrt( M1/ap)*cos(theta0) + VY1;

%Determine initial acceleration of star 1
[aX1,aY1] = gravity( X1,Y1, X2, Y2, M2 );

%Determine initial acceleration of star 2
[aX2,aY2] = gravity( X2,Y2, X1, Y1, M1 );

%Determine initial acceleration of planet
[ax1,ay1] = gravity( x,y, X1, Y1, M1 );
[ax2,ay2] = gravity( x,y, X2, Y2, M2 );
ax = ax1 + ax2; ay = ay1 + ay2;

%

%% COMPUTE ORBITS USING VERLET METHOD %%

n=1;
while t(end) < num_periods * T
    t(n+1) = t(n) + dt;
    
    %x,y update for stars and planet
    X1(n+1) = X1(n) + VX1*dt + 0.5*aX1*dt^2;
    Y1(n+1) = Y1(n) + VY1*dt + 0.5*aY1*dt^2;
    X2(n+1) = X2(n) + VX2*dt + 0.5*aX2*dt^2;
    Y2(n+1) = Y2(n) + VY2*dt + 0.5*aY2*dt^2;
    x(n+1) = x(n) + vx*dt + 0.5*ax*dt^2;
    y(n+1) = y(n) + vy*dt + 0.5*ay*dt^2;
    
    %Store old accelerations
    ax_old = ax; ay_old = ay;
    aX1_old = aX1; aY1_old = aY1;
    aX2_old = aX2; aY2_old = aY2;
    
    %Determine acceleration of star 1
    [aX1,aY1] = gravity( X1(end),Y1(end), X2(end), Y2(end), M2 );
    
    %Determine acceleration of star 2
    [aX2,aY2] = gravity( X2(end),Y2(end), X1(end), Y1(end), M1 );
    
    %Determine acceleration of planet
    [ax1,ay1] = gravity( x(end),y(end), X1(end), Y1(end), M1 );
    [ax2,ay2] = gravity( x(end),y(end), X2(end), Y2(end), M2 );
    ax = ax1 + ax2; ay = ay1 + ay2;
    
    %Work out new velocities
    VX1 = VX1 + 0.5*dt*( aX1_old + aX1 );
    VY1 = VY1 + 0.5*dt*( aY1_old + aY1 );
    VX2 = VX2 + 0.5*dt*( aX2_old + aX2 );
    VY2 = VY2 + 0.5*dt*( aY2_old + aY2 );
    vx = vx + 0.5*dt*( ax_old + ax );
    vy = vy + 0.5*dt*( ay_old + ay );
    n = n + 1;
end

%Determine starting n
n = floor( (Pstart/num_periods)*length(x) );

%Plot trajectory
figure('color',[1 1 1],'name','star and planet','units','normalized',...
    'position',[0.2 0.2 0.5 0.5],'renderer','painters');
plot(x,y,'g'); hold on
plot(X1,Y1,'r'); plot(X2,Y2,'b');
s1 = plot( X1(n),Y1(n),'r*','markersize',30 );
s2 = plot( X2(n),Y2(n),'b*','markersize',30 );
ss1 = plot(X1(n),Y1(n),'ro','markersize',15,'markerfacecolor','r');
ss2 = plot(X2(n),Y2(n),'bo','markersize',15,'markerfacecolor','b');
pp = plot(x(1:n),y(1:n),'k-','linewidth',1);
p = plot(x(n),y(n),'ko','markersize',5,'markerfacecolor','k');
axis equal; axis manual; grid on;
xlim([-limit,limit]); ylim([-limit,limit]);
xlabel('x /AU','fontsize',fsize); ylabel('y /AU','fontsize',fsize);
title(['M1=',num2str(M1),', M2=',num2str(M2),...
    ' T=',num2str(T,3),' years, a=',num2str(a,3),'AU, k=',num2str(k),...
    ', a_p=',num2str(ap,3),'AU.'],'fontsize',fsize)
set(gca,'fontsize',fsize);

%Print a PNG file
print( gcf,'binary stars & planet.png','-dpng','-r300' );

%Animation
stop = 0; N = length(x); n=1;
while stop == 0
    try
        set( pp, 'xdata', x(1:n), 'ydata', y(1:n) );
        set( p, 'xdata', x(n), 'ydata', y(n) );
        set( s1, 'xdata', X1(n), 'ydata', Y1(n) );
        set( s2, 'xdata', X2(n), 'ydata', Y2(n) );
        set( ss1, 'xdata', X1(n), 'ydata', Y1(n) );
        set( ss2, 'xdata', X2(n), 'ydata', Y2(n) );
        drawnow; n = n+1; if n>N; n=1; end
    catch
        return
    end
end

%%

%Newton's law of gravitation. x,y accelerations of a mass at location x,y
%resulting from the force of gravity of mass M1 at location x1,y1.
%Units of distance are AU and time are years.
function [ax,ay] = gravity( x,y, x1, y1, M1 )
r = sqrt( (x-x1)^2 + (y-y1)^2 );
if r > 0.01
    ax = -( 4*pi^2 )*M1*(x-x1)/ ( r^3 ) ;
    ay = -( 4*pi^2 )*M1*(y-y1)/ ( r^3 ) ;
else
    %For really close approaches, don't compute acceleration
    ax=0;
    ay=0;
end

%End of code


