%                    Function mfile  flow1        
%
%   Call is flow1(corners, T) where corners = [a,b,c,d] is the
%   vector of corner coordinates of the rectangle R where the
%   vector field is to be displayed.  T is the time interval over
%   which the flow is to be followed. 
%       The flow is generated by a vector field [u,v].  u(x,y)
%   and v(x,y) are to be provided in mfiles u.m and v.m. 
%   Also requires the mfile wdot.m
%   After the call, the vector field is plotted, and the program asks
%   for the number of starting points of streamlines. After this
%   number M is entered, hit return. The program then waits for the
%   user to click on the figure at the desired starting points. 
%   The streamline starting at this point is computed using ode45
%   and plotted. This is done M times. 


function out = flow1(corners,T)
 a = corners(1); b = corners(2); c = corners(3); d = corners(4);
 x = linspace(a,b,16); y = linspace(c,d,16);
 [X,Y] = meshgrid(x,y);

 U = u(X,Y); V = v(X,Y);

quiver(X,Y,U,V)
axis([a b c d])
hold on

M = input('enter the number of starting points    ')

m = 1;
while m < M+1
  [aa bb] = ginput(1);
  x0 = aa(1); y0 = bb(1);
  w0 = [x0, y0]';

  [t,w] = ode45('wdot', [0, T], w0);
  plot( w(:,1), w(:,2), 'r')
  nsteps = size(w,1)
  m = m+1;
end


hold off