clear all
close all

%=========
% path_lines
%
% Draw path lines
%=========

Nsteps=400;

Dl = 0.01;  % spatial step

xmin =-1.0;  % plotting window
xmax = 1.0;
ymin =-1.0;
ymax = 1.0;

Xstart = [0.1, 0.2, 0.3]; % starting points
Ystart = [0.1, 0.2, 0.3];

hold on

%---
for n=1:size(Xstart')    % run over initial points
%---

 X0 = Xstart(n);
 Y0 = Ystart(n);

 Xsave=X0; Ysave=Y0;

 clear Xplot Yplot

 Xplot(1)=X0;
 Yplot(1)=Y0;

 X = X0; Y = Y0; t=0;

%--
 for i=1:Nsteps
%--

  [ux, uy] = velocity(X,Y,t);
  um=sqrt(ux*ux+uy*uy);
  Dt = Dl/um;
  X1 = X+ux*Dt;
  Y1 = Y+uy*Dt;
  t1 = t+Dt;

  [ux1, uy1] = velocity(X1,Y1,t1);
  X = X + 0.5*Dt*(ux+ux1);
  Y = Y + 0.5*Dt*(uy+uy1);
  t = t+Dt;

  Xplot(i+1)=X;
  Yplot(i+1)=Y;

%--- stopping checks:

  closed = sqrt((X-Xsave)^2+(Y-Ysave)^2);  % for a closed line
  if(closed < 0.9*Dl | X<xmin |X>xmax |Y<ymin|Y>ymax) 
   break
  end

end

plot(Xplot,Yplot,'.')

end  % of path lines

%---
% finish
%---

xlabel('x','fontsize',15)
ylabel('y','fontsize',15)
axis equal
set(gca,'fontsize',15)
box