%
%  rotor simulation (load theta.mat file first)
%

r = 0.95;  uc = complex(1,0);  vc = complex(-1,0);
w1 = exp(i*2*pi/3);  w2 = exp(-i*2*pi/3);

th0 = 0:pi/15:2*pi;  th1 = th(:,2)';  th2 = th(:,3)';
u0 = uc + r*   exp(i*th1);  v0 = vc + r*   exp(i*th2);
u1 = uc + r*w1*exp(i*th1);  v1 = vc + r*w1*exp(i*th2);
u2 = uc + r*w2*exp(i*th1);  v2 = vc + r*w2*exp(i*th2);

figure(1)
plot(th(:,1),th1,'b');  hold on;  plot(th(:,1),th2,'g')
title('\bf rotation');  xlabel('\bf time');  
ylabel('\bf \theta_{right} (b) & \theta_{left} (g)')
pause

grey = [1 1 1]*0.6;

figure(3);  clf
plot([-1 -1 1 1 -1]*2,[-1 1 1 -1 -1],'k');
hold on;  axis equal;  axinfo = axis;

for k = 1:size(th,1)
  hold off;  clf;  axis off; hold on; axis(axinfo)
  plot([u2(k) uc u1(k) uc u0(k)],'color',grey,'linewidth',4)
  plot([v2(k) vc v1(k) vc v0(k)],'color',grey,'linewidth',4)
  pause(0.01)
end