Денис Савицкий

clc; clear all; close all;

% Start points
x1 = 0; y1 = 0; r1 = 2;
x2 = 6; y2 = 1; r2 = 2;

% Absolute Velocity
v1 = 6;
v2 = 3;

% Velocity vectors
L1 = rand(1,2);
L2 = rand(1,2);

% Time
t = 0:0.02:10;

% Vectors with coordinates by steps (time)
vec1(1,:) = x1 + v1*L1(1)*t; vec1(2,:) = y1 + v1*L1(2)*t;
vec2(1,:) = x2 + v2*L2(1)*t; vec2(2,:) = y2 + v2*L2(2)*t;

% Draw circles
ang = 0:0.01:2*pi;
xp1 = r1*cos(ang); yp1 = r1*sin(ang);
xp2 = r2*cos(ang); yp2 = r2*sin(ang);

for i=1:length(vec1)
    plot(vec1(1,i)+xp1,vec1(2,i)+yp1);
    hold on
    plot(vec2(1,i)+xp2,vec2(2,i)+yp2);
    hold off
    axis([min(min(vec1(1,:), vec2(1,:))), max(max(vec1(1,:), vec2(1,:))), min(min(vec1(2,:), vec2(2,:))), max(max(vec1(2,:), vec2(2,:)))])
    pause(0.02) % delay
    
    if sqrt((vec1(1,i)-vec2(1,i)).^2 + (vec1(2,i)-vec2(2,i)).^2) < r1+r2
        
        disp('COLLAPSE!')
        
        L1 = rot90(L1,1);
        L2 = rot90(L2,1);
        tt = 0:0.02:10;
        
        newvec1(1,:) = vec1(1,i) + v1*L1(1)*tt;
        newvec1(2,:) = vec1(2,i) + v1*L1(2)*tt;
        newvec2(1,:) = vec2(1,i) + v2*L2(1)*tt;
        newvec2(2,:) = vec2(2,i) + v2*L2(2)*tt;
        
        for j=i:length(vec1)
            plot(newvec1(1,j-i+1)+xp1,newvec1(2,j-i+1)+yp1);
            hold on;
            
            plot(newvec2(1,j-i+1)+xp2,newvec2(2,j-i+1)+yp2);
            hold off;
            
            axis([min(min(vec1(1,:), vec2(1,:))), max(max(vec1(1,:), vec2(1,:))), min(min(vec1(2,:), vec2(2,:))), max(max(vec1(2,:), vec2(2,:)))])
            pause(0.1) % delay
        end
        
        return
        
    end
end