close all; clear all;
N = 7;
T = 1000;
sep = 2.5;
mark = '.k';
animate = 0;

t = linspace(0,2*pi,t_step);

circle = zeros(N, T, 2);
for i = 1:N
  circle(i,:,1) = cos(t*i);
  circle(i,:,2) = sin(t*i);
end

curve = zeros(N, N, T, 2);
for i = 1:N
  for j = 1:N
    curve(i,j,:,1) = circle(i,:,1);
    curve(i,j,:,2) = circle(j,:,2);
  end
end

hold on;
axis([-1 N*sep+1 -(N*sep+1) 1])
axis square
set(gca,'visible','off')

if animate == 0
  for i = 1:N
    plot(circle(i,:,1)+sep*i, circle(i,:,2),       mark)
    plot(circle(i,:,1),       circle(i,:,2)-sep*i, mark)

    for j = 1:N
      plot(curve(i,j,:,1)+sep*i, curve(i,j,:,2)-sep*j, mark)
    end

  end
end

if animate == 1
  for t = 1:T
    xbordercirc = zeros(N,1);
    ybordercirc = zeros(N,1);
    curvecirc = zeros(N,N);

    for i = 1:N
      plot(circle(i,t,1)+sep*i, circle(i,t,2),     '.k')
      plot(circle(i,t,1),     circle(i,t,2)-sep*i, '.k')

      xbordercirc(i) = plot(circle(i,t,1)+sep*i, circle(i,t,2),     'ok');
      ybordercirc(i) = plot(circle(i,t,1),     circle(i,t,2)-sep*i, 'ok');

      for j = 1:N
        plot(curve(i,j,t,1)+sep*i, curve(i,j,t,2)-sep*j, '.k')

        curvecirc(i,j) = plot(curve(i,j,t,1)+sep*i, curve(i,j,t,2)-sep*j, 'ok');
      end

    end

    pause(0.01)
    delete(xbordercirc);
    delete(ybordercirc);
    delete(curvecirc);
  end
end