a=-1;
t=0:0.075:4*pi;
t1=(0:4:360)*pi/180;
x1=a/2*(t*sqrt(2)+sin(t*sqrt(2)));
y1=a/2*(1-cos(t*sqrt(2)));
x2=a/2*(t*sqrt(2)-sin(t*sqrt(2)));
y2=a/2*(1+cos(t*sqrt(2)));
N=length(x1);
figure(1)
axis([min(min(x1,x2))-a,max(max(x1,x2))+a,-1,a+1])
axis equal
hold on
tic
for i=1:N
    xO=(x1(i)+x2(i))/2;
    yO=(y1(i)+y2(i))/2;
    plot(xO+a*0.5*cos(t),yO+a*0.5*sin(t),'k') 
    
    plot(x1(i),y1(i),'ob')
    plot(x2(i),y2(i),'or')
    drawnow
    
    plot(x1(i),y1(i),'ow')
    plot(x2(i),y2(i),'ow')
    plot(xO+a*0.5*cos(t),yO+a*0.5*sin(t),'w') 
end
toc
comet(x1,y1,0.5)
hold on
comet(x2,y2,0.5)