%circular_orbits
% Plots, and animates, circular orbits of planets about a star.
% LAST UPDATED by Dr French. Sept 2023.
function circular_orbits
M1 = 1; M2 = 1; %Mass of stars in solar masses
r1 = 5.2; r2 = 13.1; %Obital radius of planet in AU
T1 = sqrt( r1^3 / M1 ); %Calculate orbital period in Yr
T2 = sqrt( r2^3 / M2 ); %Calculate orbital period in Yr
Tmax = max([T1,T2]);
[x1,y1,t] = orbit_calculator(r1,M1,T1,Tmax);  %Calculate trajectories
[x2,y2,t] = orbit_calculator(r2,M2,T2,Tmax); N = length(t);

%Set up axes and title
axes('nextplot','add');
tt = title( {['M1=',num2str(M1),', r1=',num2str(r1),'AU, T1=',num2str(T2),'Yr'],...
   ['M2=',num2str(M2),', r2=',num2str(r2),'AU, T2=',num2str(T2),'Yr'] } );

%Plot lines between planet positions
for n=1:N
    plot( [x1(n),x2(n)],[y1(n),y2(n)],'k-');
end

%Plot trajectories
plot(x1,y1,'b-'); axis equal; grid on; box on; xlabel('x/AU'); ylabel('y/AU')
plot(x2,y2,'r-'); star = plot(0,0,'r*');  %Add star 

%Plot planet initial positions
p1 = plot(x1(1),y1(1),'b*'); p2 = plot(x2(1),y2(1),'r*'); %Plot planets
set(p1,'markersize',15); set(p1,'markersize',15); set(star,'markersize',15);

%ANIMATE!
stop = 0; n = 1; 
while stop~=1
    n = n+1; 
    if n>N
        n = 1;
    end
    set( p1, 'xdata', x1(n), 'ydata', y1(n) );
    set( p2, 'xdata', x2(n), 'ydata', y2(n) );
    drawnow; pause(0.01);
end

%%

%Orbit calculator
function [x,y,t] = orbit_calculator(r,M,T,Tmax)
t = linspace( 0, Tmax, 200 ); %Array or vector times /Yr
theta = 2*pi*t/T;  %Polar angle (rad) from the x axis
x = r*cos(theta); y = r*sin(theta); %Vectors for x,y coordinates

%End of code