%
%   XYFLOW1.M  Sample showing various parameters for a simple XY flow field
%
%   This example comes directly from several sample exercises from the textbook
%   by Munson (3rd edition)  --  see examples 4.1, 4.2 and 4.5.  This is also the
%   same problem as Ex. 3.5 in the Fluid Mechanics text by Hibbeler (2015)
%
%   For this steady flow (i.e. time independent) example, the velocity velocity is
%      V  =  Vx i +  Vy j  =  2x i - 2y j   with units of m/s
%
%   File prepared by J. R. White, UMass-Lowell (last update:  Feb. 2017)
%
      
      clear all;   close all;   nfig = 0;
%
%   plot the VELOCITY vector field using Matlab's quiver function
      xf = 6;  yf = 6;  r = [eps xf eps yf];
      Nx = 16;   xx = linspace(0,xf,Nx);
      Ny = 16;   yy = linspace(0,yf,Ny);
      [x,y] = meshgrid(xx,yy);
      Vx = 2*x;   Vy = -2*y;   
%
      nfig = nfig+1;   figure(nfig)
      quiver(x,y,Vx,Vy), axis(r)
      title('XYFlow1:  Velocity Vector Field')
      xlabel('x-direction'),ylabel('y-direction')
%
%   plot the VELOCITY direction field (use unit vectors to just show direction)
      xf = 6;  yf = 6;  r = [eps xf eps yf];
      Nx = 25;   xx = linspace(0,xf,Nx);
      Ny = 25;   yy = linspace(0,yf,Ny);
      [x,y] = meshgrid(xx,yy);
      Vx = 2*x;   Vy = -2*y;   Vmag = sqrt(Vx.^2 + Vy.^2);
      Vxn = Vx./Vmag;   Vyn = Vy./Vmag;
%
      nfig = nfig+1;   figure(nfig)
      quiver(x,y,Vxn,Vyn), axis(r)
      title('XYFlow1:  Velocity Direction Field')
      xlabel('x-direction'),ylabel('y-direction')
%
%   do this again but let's add some lines of constant SPEED (magnitude of velocity)
      nfig = nfig+1;   figure(nfig)
      quiver(x,y,Vxn,Vyn); hold on
      [cs,hh] = contour(x,y,Vmag,[2 4 6 8 10 12]); 
      clabel(cs); set(hh,'LineWidth',2); hold off, axis(r)
      title('XYFlow1:  Velocity Direction Field and Magnitude Contours (m/s)')
      xlabel('x-direction'),ylabel('y-direction')
%
%   do this again but let's add some streamlines this time
%***  this version solves the differential eqn for the streamlines using ODE23 
      nfig = nfig+1;   figure(nfig)
      quiver(x,y,Vxn,Vyn); hold on, axis(r)
      title('XYFlow1:  Velocity Direction Field and Typical Streamlines (using ODE23)')
      xlabel('x-direction'),ylabel('y-direction')
%   find some typical streamlines that are everywhere tangent to the velocity vector
      xo = [0.05 0.2 0.6 1.5 3 4.5];  yo = yf;  Ni = length(xo);
      for n = 1:Ni
        fxy = @(x,y) -y./x;   % this is dy/dx = f(x,y) where f(x,y) = -y/x here
        [xvs,yvs] = ode23(fxy,[xo(n) xf],yo);  % numerical solution using ODE23
        plot(xvs,yvs,'r-','LineWidth',2)
      end
      hold off
%
%   do this again but with the analytical solution for the streamlines 
%***  this version uses the ANALYTICAL solution for the streamlines 
      nfig = nfig+1;   figure(nfig)
      quiver(x,y,Vxn,Vyn); hold on, axis(r)
      title('XYFlow1:  Velocity Direction Field and Typical Streamlines (Analytical Soln)')
      xlabel('x-direction'),ylabel('y-direction')
%   find some typical streamlines that are everywhere tangent to the velocity vector
      xo = [0.05 0.2 0.6 1.5 3 4.5];  yo = yf;  C = yo*xo;  Ni = length(xo);
      for n = 1:Ni
        Nx = 50;   xx = linspace(xo(n),xf,Nx);
        yx = @(x) C(n)./x;   % this is the analytical solution for the streamlines
        plot(xx,yx(xx),'r-','LineWidth',2)
      end
      hold off
%
%   let's plot the ACCELERATION direction field with and without magnitude contours
      Ax = 4*x;   Ay = 4*y;   Amag = sqrt(Ax.^2 + Ay.^2);
      Axn = Ax./Vmag;   Ayn = Ay./Amag;
      nfig = nfig+1;   figure(nfig)
      quiver(x,y,Axn,Ayn), axis(r)
      title('XYFlow1:  Acceleration Direction Field')
      xlabel('x-direction'),ylabel('y-direction')
%
      nfig = nfig+1;   figure(nfig)
      quiver(x,y,Axn,Ayn), axis(r), hold on
      [cs,hh] = contour(x,y,Amag,[4 8 12 16 20 24]); 
      clabel(cs); set(hh,'LineWidth',2); hold off
      title('XYFlow1:  Acceleration Direction Field and Magnitude Contours (m/s^2)')
      xlabel('x-direction'),ylabel('y-direction')
%
%   end of demo
