% 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