%close all;
clear all;

%n=2;
n=input('Enter number of vortices (2<=n<=10) >> ');
nts=input('Enter number of timesteps (~200) >> ');

%g=[-1 -1];
for i=1:n;
   l1=['Enter strength of vortex ' num2str(i) ' >> '];
   g(i)=input(l1);
   l2=['Enter x position of vortex ' num2str(i) ' >> '];
   l3=['Enter y position of vortex ' num2str(i) ' >> '];
   x(i)=input(l2);
   y(i)=input(l3);
end;

z=max(max(abs(x)),max(abs(y)));
mc=['r*';'b*';'g*';'m*';'y*';'c*';'k*';'r^';'b^'];
lc=['r- ';'b- ';'g- ';'m- ';'y- ';'c- ';'k- ';'r--';'b--'];

figure;hold on;
for i=1:n;
   plot(x(i),y(i),mc(i,:));
end;

for i=1:n;
   u(i)=0;
   v(i)=0;
   for j=1:n;
      if i~=j;
         s=1/((x(i)-x(j))^2+(y(i)-y(j))^2);
         u(i)=u(i)-g(j)*(y(i)-y(j))*s;
         v(i)=v(i)+g(j)*(x(i)-x(j))*s;
      end;
   end;
end;
w=max(max(abs(u)),max(abs(v)));

dt=.1*z/w;

for i=1:n;
   xn(i)=x(i)+u(i)*dt;
   yn(i)=y(i)+v(i)*dt;
end;

for i=1:n;
   un(i)=0;
   vn(i)=0;
   for j=1:n;
      if i~=j;
         s=1/((xn(i)-xn(j))^2+(yn(i)-yn(j))^2);
         un(i)=un(i)-g(j)*(yn(i)-yn(j))*s;
         vn(i)=vn(i)+g(j)*(xn(i)-xn(j))*s;
      end;
   end;
end;

for i=1:n;
   xn(i)=x(i)+0.5*(u(i)+un(i))*dt;
   yn(i)=y(i)+0.5*(v(i)+vn(i))*dt;
end;

xw(:,1)=x';
xw(:,2)=xn';
yw(:,1)=y';
yw(:,2)=yn';

for k=1:nts;
x=xn;
y=yn;

for i=1:n;
   u(i)=0;
   v(i)=0;
   for j=1:n;
      if i~=j;
         s=1/((x(i)-x(j))^2+(y(i)-y(j))^2);
         u(i)=u(i)-g(j)*(y(i)-y(j))*s;
         v(i)=v(i)+g(j)*(x(i)-x(j))*s;
      end;
   end;
end;
w=max(max(abs(u)),max(abs(v)));

dt=.1*z/w;

for i=1:n;
   xn(i)=x(i)+u(i)*dt;
   yn(i)=y(i)+v(i)*dt;
end;

for i=1:n;
   un(i)=0;
   vn(i)=0;
   for j=1:n;
      if i~=j;
         s=1/((xn(i)-xn(j))^2+(yn(i)-yn(j))^2);
         un(i)=un(i)-g(j)*(yn(i)-yn(j))*s;
         vn(i)=vn(i)+g(j)*(xn(i)-xn(j))*s;
      end;
   end;
end;

for i=1:n;
   xn(i)=x(i)+0.5*(u(i)+un(i))*dt;
   yn(i)=y(i)+0.5*(v(i)+vn(i))*dt;
end;

xw(:,k+2)=xn';
yw(:,k+2)=yn';
end;

for i=1:n;
  plot(xw(i,:),yw(i,:),lc(i,:)); 
end;

%set(gca, 'DataAspectRatio', [1 1 1]);
%axis off;
box on;grid on;

