%    Calculate streamlines, particle paths, and streaklines
     clear all
     clf reset
     global tstream
%    Streamlines at t = tstream
     figure(1)
     tstream=0.5; 
     s0=0; sf=4; % Integration limits
     tspan=[s0 sf];
     xic=[[0 1]' [.3 .7]' [.5 .5]' [.7 .3]' [1. 0]']; % Initial cond.
     for n=1:5
      x0=xic(:,n);
      [s,x]=ode23('stream',tspan,x0);
      plot(x(:,1),x(:,2),'r')
      xlabel('X1')
      ylabel('X2')
      title(['Streamlines at t = ', num2str(tstream)])
      hold on
      clear s,x ;
     end
     'Press any key to continue', pause
%    Particle paths
     clear all
     figure(2)
     t0=0; tf=4; % Integration limits
     tspan=[t0 tf];
     xic=[[0 1]' [.3 .7]' [.5 .5]' [.7 .3]' [1. 0]']; % Initial cond.
     for n=1:5
      x0=xic(:,n);
      [t,x]=ode23('path',tspan,x0);
      plot(x(:,1),x(:,2),'r')
      xlabel('X1')
      ylabel('X2')
      title('Particle Paths')
      hold on
      clear t,x ;
     end
     'Press any key to continue', pause
%    Streaklines
     clear all
     figure(3)
     x0=[.5 .5]';    % Fixed source point y
     M=moviein(20);
     tmovie=[.2: .2: 4]; % times of observation
     ty=[0: .2: 4];  % times particle originated at y
     for m=1:20      % frame of movie
      tf=tmovie(m);  % time of observation
      nty=1;
      while ty(nty) < tf, nty=nty+1; end  % determine number of particles
      for n=1:nty    % particle trajectories
         t0=ty(n)-1.e-4;% time particle originated at y
         tspan=[t0 tf];
       [t,x]=ode23('path',tspan,x0);
       lt=length(t);
       plot(x(:,1),x(:,2),'r',x(lt,1),x(lt,2),'g*')  % particle path
       axis([0 3 0 30])
       xlabel('X1')
       ylabel('X2')
       title(['Particle Paths and Streakline, Y1, Y2 = ',num2str(x0(1)),...
       ', ',num2str(x0(2))])
       hold on
       if n == 1        % streakline
        xstreak=x(lt,1);
        ystreak=x(lt,2);
       else
        xstreak=[xstreak x(lt,1)];
        ystreak=[ystreak x(lt,2)];
       end
       clear t,x ;
      end
     plot(xstreak,ystreak,'g')
     text(.5, 25, [' t = ',num2str(tf)])
     clear xstreak,ystreak ;
     M(:,m)=getframe;
     hold off
     end
     movie(M,2,5)
