% Concentric rings of masslets orbit about two galactic masses in mutually
% circular orbits.
function interacting_galaxies
global d
a = 20;              %Semi-major axis of mutual star orbit in galactic radii
d.M1 = 5.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.01;           %Timestep in galactic years*
nrings = 120;        %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 = 1;           %Planet marker size
vm = 1.2;            %Velocity multiplier for VY2 from circular orbit condition
fsize = 18;          %Fontsize

%* 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;

%Galaxies
X1 = -d.M2*a/(d.M1+d.M2); Y1 = 0; 
VX1 = 0;  VY1 = 2*pi*X1/T;
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);

%

%Initialize trajectory
figure('color',[1 1 1],'name','Galaxies','keypressfcn',@keypressfunc,...
    'units','normalized','position',[0.05 0.05 0.85 0.85],...
    '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([-limit,limit]); ylim([-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);

%Update orbit and animate
d.stop = 0; d.quit = 0; n=1;
while d.quit == 0
    if d.stop == 0
        
        %Update time
        d.t(n+1) = d.t(n) + dt;
         
        %Run gravity solver
        
        %Form inputs
        x0 = [ X1(n),X2(n),x1(n,:),x2(n,:) ];
        y0 = [ Y1(n),Y2(n),y1(n,:),y2(n,:) ];
        z0 = zeros( size(x0) );
        vx0 = [ VX1(n),VX2(n),vx1(n,:),vx2(n,:) ];
        vy0 = [ VY1(n),VY2(n),vy1(n,:),vy2(n,:) ];
        vz0 = zeros( size(x0) );
        M = [ d.M1,d.M2,zeros(1,2*num_rings) ];
        
        %Run solver
        [t,x,y,z,vx,vy,vz] = gravity( dt,dt,x0,y0,z0,vx0,vy0,vz0,M );
        
        %Get outputs
        X1(n+1) = x(2,1); X2(n+1) = x(2,2);
        x1(n+1,:) = x(2,3:2+num_rings);
        x2(n+1,:) = x(2,3+num_rings:end);   
        
        Y1(n+1) = y(2,1); Y2(n+1) = y(2,2);
        y1(n+1,:) = y(2,3:2+num_rings);
        y2(n+1,:) = y(2,3+num_rings:end); 
        
        VX1(n+1) = vx(2,1); VX2(n+1) = vx(2,2);
        vx1(n+1,:) = vx(2,3:2+num_rings);
        vx2(n+1,:) = vx(2,3+num_rings:end);    
        
        VY1(n+1) = vy(2,1); VY2(n+1) = vy(2,2);
        vy1(n+1,:) = vy(2,3:2+num_rings);
        vy2(n+1,:) = vy(2,3+num_rings:end); 
        
        %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
    else
        pause(1);
    end
end
close(gcf);

%%

%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


