%---
close all
clear all
hold on
%---

%=========
% streamlines in a rectangular cavity
%
% confined in ax<x<bx
%             ay<y<by
%=========

ax=-1.2; bx=1.2;
ay=-1.0; by=1.0;
Nx=16; Ny=20;

Nmax=400;  % maximum steps along a streamline
Ds=0.020;  % travel distance

%-----
% grid generation
%-----

[glx,gly,gx,gy] = grid_2d (ax,bx,ay,by,Nx,Ny);

%---
% specify the grid velocity
%---

 for i=1:Nx+1
   for j=1:Ny+1
     % extensional
     gux(i,j) = gx(i,j);
     guy(i,j) =-gy(i,j);
     % kovasnay
     px = gx(i,j);
     py = gy(i,j);
      wn  = 2*pi/0.5;
      pxs = px+0.5;
      gux(i,j) = 1.0-cos(wn*py)*exp(-1.0*wn*pxs);
      guy(i,j) =    -sin(wn*py)*exp(-1.0*wn*pxs);
     % rotational
     gux(i,j) =-gy(i,j);
     guy(i,j) = gx(i,j);
   end
 end

%---------------------------
% initial streamline points
%---------------------------

%x0=[-0.499, -0.499, -0.499, -0.499, -0.499, -0.499, -0.499, -0.499];
%y0=[ 0.1,    0.2,    0.3,    0.4,   -0.1,   -0.2,   -0.3,   -0.4];
%Nstr=8;

x0=[0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.2];
Nstr=size(x0');
for i=1:Nstr
 y0(i)=0.001;
end

for l=1:Nstr       % loop over streamlines

  xn = x0(l);
  yn = y0(l);
  xstr(1) = xn; % new point
  ystr(1) = yn; % new point

%----
% integrate by the modified Euler method
%---- 

  for i=2:Nmax   % step in time

   [ux,uy] = rec_2d_int ...
    ...
      (ax,bx ... % x end points
      ,ay,by ... % y end points
      ,Nx,Ny ... % grid size
      ,glx,gly ... % grid lines
      ,gux,guy ... % grid velocity
      ,xn,yn ...      % interpolation point
      );

      Umag = sqrt(ux*ux+uy*uy);
      Dt   = Ds/Umag;   % establish time step
                        %  for constant travel distance
      x1 = xn + Dt*ux;
      y1 = yn + Dt*uy;

   [ux1,uy1] = rec_2d_int ...
    ...
      (ax,bx ... % x end points
      ,ay,by ... % y end points
      ,Nx,Ny ... % grid size
      ,glx,gly ... % grid lines
      ,gux,guy ... % grid velocity
      ,x1,y1 ...      % interpolation point
      );

      Dth = 0.5*Dt;
      xn  = xn + Dth*(ux+ux1);
      yn  = yn + Dth*(uy+uy1);

      xstr(i) = xn;
      ystr(i) = yn;

%----
% check for a closed streamline after 5 steps
%---

     if(i>5) 
       Dist = sqrt((xn-xstr(1))^2+(yn-ystr(1))^2);
       if(Dist<Ds)
         xstr(i+1) = xstr(1);
         ystr(i+1) = ystr(1);
         break;
       end
     end

   end

 plot(xstr,ystr,'r-')
 clear xstr ystr

end

%---
% finish
%---
axis([ax-0.1, bx+0.1, ay-0.1, by+0.1])
axis equal
xlabel('x','fontsize',15)
ylabel('y','fontsize',15)
set(gca,'fontsize',15)
box