%    Deformation during motion
%    (1) Plane stagnation flow: Aris Exercise 4.42.1
%    (2) Pure shear, Couette flow
     clear all
     clf reset
    for Case=1:3
     figure(Case)
     cases=['stagf';'shear';'EX413'];
     titles=['Stagnation Flow';'Couette Flow   ';'Section 4.13   '];
     axes=[[0 16 0 4]; [0 10  0 4]; [0 8  0 25]];
     nty=4;
     xin=[[1 2]' [2 2]' [2 3]' [1 3]'];  % Initial element
     M=moviein(20);
     tmovie=[.1: .1: 2]; % times of observation
     t0=0;               % time particle originated
     for m=1:20      % frame of movie
      tf=tmovie(m);  % time of observation
      for n=1:nty    % particle trajectories
         x0=xin(:,n);
         tspan=[t0 tf];
       [t,x]=ode23(cases(Case,:),tspan,x0);
       lt=length(t);
       plot(x(:,1),x(:,2),'r',x(lt,1),x(lt,2),'g*')  % particle path
       axis(axes(Case,:))
       xlabel('X1')
       ylabel('X2')
       title(titles(Case,:))
       hold on
       if n == 1        % Shape of element
        xstreak=x(lt,1);
        ystreak=x(lt,2);
       else
        xstreak=[xstreak x(lt,1)];
        ystreak=[ystreak x(lt,2)];
       end
       if n == nty
        xstreak=[xstreak xstreak(1)];
        ystreak=[ystreak ystreak(1)];
       end
       clear t,x ;
      end
     plot(xstreak,ystreak,'g')
     plot([xin(1,:) xin(1,1)],[xin(2,:) xin(2,1)], 'b')
     clear xstreak,ystreak ;
     M(:,m)=getframe;
     hold off
     end
     movie(M,2,5)
     'Press any key to continue', pause
     clear M ;
    end
