% Interacting galaxies simulation!
% Concentric rings of masslets orbit about two galactic masses in mutually
% circular orbits.
%
% LAST UPDATED by Andy French June 2021

function galaxies
global d
a = 12;              %Semi-major axis of mutual star orbit in galactic radii
d.M1 = 15.0;          %Mass of galaxy 1
d.M2 = 3.0;          %Mass of galaxy 2
d.dM = 0.1;          %Change in mass when you press the m or n keys (for M1)
dt = 0.05;           %Timestep in galactic years*
nrings = 100;        %Number of rings of masslets
rmin = 1;            %Radius of minimum ring
dr = 0.05;           %Difference in ring radii
m0 = 20;             %Masslets in inner ring - this sets spacing for the outer rings
limit = a;           %Plot axis limits
msize = 0.2;         %Masslet marker size
vm = 1.0;            %Velocity multiplier for VY2 from circular orbit condition
fsize = 18;          %Fontsize
tprint = 0;          %Print frame every tprint galactic years

%* For Milky Way https://en.wikipedia.org/wiki/Milky_Way
% M = 0.8 to 1.5 x10^12 solar masses
% Distance of Sun to galactic centre = 26.4 x 1000 light years
% Sun's galactic rotation period = 240 million years

%Determine period in years via Kepler III of mutual galaxy rotation
T = sqrt( (1/(d.M1+d.M2))*a^3 );

%Set initial x,y,vx,vy,t parameters
d.t = 0;

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

%Galaxy 2
X2 = d.M1*a/(d.M1+d.M2); Y2 = 0; VX2 = 0;  VY2 = vm*2*pi*X2/T;

%Rings
[x1,y1,vx1,vy1] = rings( rmin, dr, nrings, m0, d.M1, X1, Y1, VX1, VY1 );
[x2,y2,vx2,vy2] = rings( rmin, dr, nrings, m0, d.M2, X2, Y2, VX2, VY2 );
num_rings = length(x1);

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

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

%Determine initial acceleration of rings

%Rings about galaxy 1
[ax11,ay11] = gravity( x1,y1, X1, Y1, d.M1 );
[ax12,ay12] = gravity( x1,y1, X2, Y2, d.M2 );
ax1 = ax11 + ax12; ay1 = ay11 + ay12;

%Rings about galaxy 2
[ax21,ay21] = gravity( x2,y2, X1, Y1, d.M1 );
[ax22,ay22] = gravity( x2,y2, X2, Y2, d.M2 );
ax2 = ax21 + ax22; ay2 = ay21 + ay22;

%Initialize trajectory
figure('color',[1 1 1],'name','Galaxies','keypressfcn',@keypressfunc,...
    'units','normalized','position',[0.2 0.2 0.5 0.5],...
    'renderer','painters');
S1 = plot(X1,Y1,'r'); hold on
S2 = plot(X2,Y2,'b');
s1 = plot( X1(1),Y1(1),'r*','markersize',20 );
s2 = plot( X2(1),Y2(1),'b*','markersize',20 );
r1 = plot( x1(1,:),y1(1,:),'m*','markersize',msize );
r2 = plot( x2(1,:),y2(1,:),'g*','markersize',msize );
axis equal; axis manual; grid on;
xlim(2*[-limit,limit]);
ylim(2*[-limit,limit]);
xlabel('x','fontsize',fsize); ylabel('y','fontsize',fsize);
tit = title(['M1=',num2str(d.M1),', M2=',num2str(d.M2),', T=',num2str(T,3),', t=0'],...
    'fontsize',fsize);
set(gca,'fontsize',fsize);

%Print first frame
fname = ['gravity t=',num2str(d.t(end)),' M1=',num2str(d.M1),' M2=',num2str(d.M2)];
fname = strrep(fname,'.','p');
print( gcf, [fname,'.png'],'-dpng','-r300');

%Update orbit and animate
d.stop = 0; d.quit = 0; n=1; np = 1;
while d.quit == 0
    if d.stop == 0
        
        %Update time
        d.t(n+1) = d.t(n) + dt;
        
        %x,y update
        
        %Galaxies
        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;
        
        %Rings
        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;
        
        %Store old accelerations
        ax1_old = ax1; ay1_old = ay1;
        ax2_old = ax2; ay2_old = ay2;
        aX1_old = aX1; aY1_old = aY1;
        aX2_old = aX2; aY2_old = aY2;
        
        %Determine acceleration of galaxy 1
        [aX1,aY1] = gravity( X1(end),Y1(end), X2(end), Y2(end), d.M2 );
        
        %Determine acceleration of galaxy 2
        [aX2,aY2] = gravity( X2(end),Y2(end), X1(end), Y1(end), d.M1 );
        
        %Determine acceleration of rings
        
        %Rings about galaxy 1
        [ax11,ay11] = gravity( x1(end,:),y1(end,:), X1(end), Y1(end), d.M1 );
        [ax12,ay12] = gravity( x1(end,:),y1(end,:), X2(end), Y2(end), d.M2 );
        ax1 = ax11 + ax12; ay1 = ay11 + ay12;
        
        %Rings about galaxy 2
        [ax21,ay21] = gravity( x2(end,:),y2(end,:), X1(end), Y1(end), d.M1 );
        [ax22,ay22] = gravity( x2(end,:),y2(end,:), X2(end), Y2(end), d.M2 );
        ax2 = ax21 + ax22; ay2 = ay21 + ay22;
        
        %Work out new velocities
        
        %Galaxies
        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 );
        
        %Rings
        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 );
        
        %Update plot
        n = n + 1;
        if d.quit == 0
            try
                set( r1, 'xdata', x1(n,:), 'ydata', y1(n,:) );
                set( r2, 'xdata', x2(n,:), 'ydata', y2(n,:) );
                set( s1, 'xdata', X1(n), 'ydata', Y1(n) );
                set( s2, 'xdata', X2(n), 'ydata', Y2(n) );
                set( S1, 'xdata', X1(1:n), 'ydata', Y1(1:n) );
                set( S2, 'xdata', X2(1:n), 'ydata', Y2(1:n) );
                set( tit, 'string', ['M1=',num2str(d.M1),', M2=',num2str(d.M2),...
                    ', T=',num2str(T,3),', t=',num2str(d.t(end),3)])
                drawnow;
            catch
                return
            end
        end
        
        %Print frame
        if rem( d.t(n), tprint ) < dt
            fname = ['gravity t=',num2str(d.t(end)),' M1=',num2str(d.M1),' M2=',num2str(d.M2)];
            fname = strrep(fname,'.','p');
            print( gcf, [fname,'.png'],'-dpng','-r300');
        end
    else
        pause(1);
    end
end
close(gcf);

%%

%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 )
rmin = 0.1;
r = sqrt( (x-x1).^2 + (y-y1).^2 );
ax = -( 4*pi^2 )*M1*(x-x1)./(r.^3) ;
ay = -( 4*pi^2 )*M1*(y-y1)./(r.^3) ;
ax( r < rmin ) = 0;
ay( r < rmin ) = 0;


%%

%Rings of masslests about mass M1 at position x1,y1
function [x,y,vx,vy] = rings( rmin, dr, nrings, m0, M1, x1, y1, vx1, vy1 )
x = []; y = []; vx = []; vy = [];
for n=1:nrings
    %Ring radius
    r = rmin +(n-1)*dr;
    
    %Calculate masslets per ring
    m = round( m0*r/rmin );
    
    %Determine initial x,y coordinates
    theta = 0 : 2*pi/m : 2*pi*(m-1)/m;
    x = [x, x1 + r*cos(theta)] ;
    y = [y, y1 + r*sin(theta)] ;
    
    %Determine initial x,y velocities in AU/year
    vx = [vx,vx1-2*pi*sqrt( M1/r )*sin(theta)] ;
    vy = [vy,vy1+2*pi*sqrt( M1/r )*cos(theta)] ;
end

%%

% Key press handler
function keypressfunc( fig,evnt )
global d
if strcmp(get(fig,'currentkey'),'s')==1
    d.stop = 1;
elseif strcmp(get(fig,'currentkey'),'c')==1
    d.stop = 0;
elseif strcmp(get(fig,'currentkey'),'q')==1
    d.stop = 1; d.quit = 1;
    close(fig);
    return
elseif strcmp(get(fig,'currentkey'),'m')==1
    d.M1 = d.M1 + d.dM;
elseif strcmp(get(fig,'currentkey'),'n')==1
    d.M1 = d.M1 - d.dM;
elseif strcmp(get(fig,'currentkey'),'p')==1
    fname = ['gravity t=',num2str(d.t(end)),' M1=',num2str(d.M1),' M2=',num2str(d.M2)];
    fname = strrep(fname,'.','p');
    print( gcf, [fname,'.png'],'-dpng','-r300');
end

%End of code


