% simulating crashing of asteroids and resulting rotation
clear all
clf

% assume uniform density of pi
dimp = 400;
D = zeros(dimp,dimp,6);
d = pi; % density
ra = 30; % radius of initial asteroid
rb = 12; %radius of incoming asteroid
cai = dimp/2; %center of asteroid a x coordinate
caj = dimp/2; %center of asteroid a y coordinate
cbj = 140; cbrj = round(cbj);%x coordinate
cbi = 20; cbri = round(cbi);% y coordinate
vbi = 2; % y velocity
vbj = 1; %x velocity
dt = 1;
f=0

while ge(norm([(cai-cbri),(caj-cbrj)]),(ra+rb))
    for i = 1:dimp
        for j = 1:dimp
            if le(((i-cai)^2+(j-caj)^2),ra^2) % plot initial circle
                D(i,j,1) = 1; %fill binary
            end
            if le(((i-cbri)^2+(j-cbrj)^2),rb^2)
                D(i,j,1)=1;
            end
        end
    end
    
    imagesc(D(:,:,1));
    set(gca,'YDir','normal');
    colormap(1-gray);
    axis equal;
    axis tight;
    saveas(gcf, ['Asteroids___',sprintf('%03d',f),'.png'])
    pause(.005)
    f=f+1
    
    D(:,:,1)=0;
    cbi= cbi+(vbi*dt); %x displacement of center of b
    cbj= cbj+(vbj*dt); %y displacement of center of b
    cbri = round(cbi);
    cbrj = round(cbj);
end

cmi_new = cai-((cai-cbi)*rb^2)/(ra^2+rb^2)
cmj_new = caj-((caj-cbj)*rb^2)/(ra^2+rb^2)
dcai = cai-cmi_new
dcaj = caj-cmj_new
dcbi = cbi-cmi_new
dcbj = cbj-cmj_new

for i = 1:dimp
    for j = 1:dimp
        if le(((i-cai)^2+(j-caj)^2),ra^2) % plot initial circle
            D(i,j,1) = 1; %fill binary
        end
        if le(((i-cbri)^2+(j-cbrj)^2),rb^2)
            D(i,j,1)=1;
        end
        if D(i,j,1) == 1
            D(i,j,5) = ((i-cmi_new)^2+(j-cmj_new)^2);
            D(i,j,6) = D(i,j,5)*d; %gives MR^2 of particles)
        end
    end
end

Inew = sum(sum((D(:,:,6))));
D(:,:,1)=0;

v=norm([vbi;vbj]);
P= rb^2*v;
L = (P/(norm([dcbi;dcbj])));

w= 100000*L/Inew
atheta= atan(((cai-cmi_new)/(caj-cmj_new)));
btheta= atan(((cai-cmi_new)/(caj-cmj_new)))+pi;
for t=1:150
    atheta=atheta+w*dt;
    btheta=btheta+w*dt;
    cai = cmi_new+dcai*(-cos(atheta));
    caj = cmj_new+dcaj*(-sin(atheta));
    cbi = cmi_new+dcbi*(cos(btheta));
    cbj = cmj_new+dcbj*(sin(btheta));
    cbri = round(cbi);
    cbrj = round(cbj);
    cari = round(cai);
    carj = round(caj);
    
    for i = 1:dimp
        for j = 1:dimp
            if le(((i-cari).^2+(j-carj).^2),ra.^2) % plot initial circle
                D(i,j,1) = 1; %fill binary
            end
            if le(((i-cbri).^2+(j-cbrj).^2),rb.^2)
                D(i,j,1)=D(i,j,1)+1;
            end
        end
    end
    
    imagesc(D(:,:,1));
    set(gca,'YDir','normal');
    colormap(1-gray);
    axis equal;
    axis tight;
    saveas(gcf, ['Asteroids___',sprintf('%03d',f),'.png'])
    pause(.005)
    f=f+1
    D(:,:,1)=0;
end

