% tcell_flow_png.m
%
%  Discussion:
%
%    This MATLAB script file reads the T-Cell flow data:
%
%      geometry (XY values at nodes, assumed to be in 'xy.txt') 
%      flow (UV values at nodes, in a sequence of files starting with 'up001.txt')
%
%    and plots the velocity vectors (U,V)(X,Y), and saves each plot as a
%    PNG file, presumably so the PNG files can be gathered into an
%    animation.
%
%    The file plots either the velocity vector field, or the velocity
%    direction field, depending on the value of the internal logical
%    parameter "normalized".
%
%    The MATLAB routine quiver internally scales the vectors, but this can be
%    adjusted by using a value of SCALE that is not 1.
%
%    For the unnormalized case, the velocity field is computed to a fixed
%    scale determined by finding the largest velocity vector over time.
%
%    This routine requires some auxiliary routines in order to "increment"
%    the name of the current velocity data file to get the name of the next
%    one, with the main such routine being called FILE_NAME_INC.
%
%  Modified:
%
%    05 July 2003
%
%  Author:
%
%    John Burkardt
%
  FALSE = 0;
  TRUE = 1;
  scale = 1.0;
  normalized = TRUE;

  nframes = 500;
%
%  Get the XY coordinates of the nodes.
%
  load xy.txt;

  x = xy(:,1);
  y = xy(:,2);
%
%  Thin the data.
%
  thin_factor = 4;
  thin_dex = thin_index ( x, y, thin_factor );
  thin_num = length ( thin_dex );
  x = x(thin_dex);
  y = y(thin_dex);
%
%  For unnormalized plots, you need to make sure that a fixed scale is preserved from 
%  step to step.  The only way I can see to do this requires that I determine the
%  maximum velocity over all time, and then append one extract node and velocity
%  to the data structure.
%
  if ( normalized == FALSE )

    upfile = 'up000.txt';
    vnorm_max = 0.0;

    for i = 1 : nframes
      fprintf ( 1, '  Checking file %d.\n', i );
      upfile = file_name_inc ( upfile );
      uv = load ( upfile );
      u = uv(thin_dex,1);
      v = uv(thin_dex,2);
      norm = sqrt ( u.^2 + v.^2 );
      vnorm_max = max ( vnorm_max, max ( norm ) );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Maximum visible velocity magnitude is %f\n', vnorm_max );

  end
%
%  Set the coordinates of the boundary lines, and an extra set of 
%  lines that will be invisible, but which make some space for the 
%  file name to be displayed within the plot area.
%
  bx1 = [ 0.00, 1.00, 1.00, 0.00, 0.00, 0.00 ];
  by1 = [ 1.00, 1.00, 0.99, 0.99, 0.99, 1.00 ];

  bx2 = [ 0.00, 0.25, 0.25, 0.75, 0.75, 1.00, 1.00, 0.74, 0.74, 0.26, 0.26, 0.00, 0.00];
  by2 = [ 0.50, 0.50, 0.00, 0.00, 0.50, 0.50, 0.51, 0.51, 0.01, 0.01, 0.51, 0.51, 0.50];

  bx3 = [ -0.10,  1.10, 1.10, -0.10, -0.10 ];
  by3 = [ -0.10, -0.10, 1.10,  1.10, -0.10 ];
%
%  Set the name of the 0th (nonexistent) velocity file.
%  All the velocity files with have this format, with the
%  numeric part of the name incremented to get the next one.
%  In particular, the first velocity file is called UP001.TXT.
%
  upfile = 'up000.txt';
  pngfile = 'up000.png';

  for i = 1 : nframes
    
    fprintf ( 1, '  Converting file %d.\n', i );
  
    upfile = file_name_inc ( upfile );
    pngfile = file_name_inc ( pngfile );
    uv = load ( upfile );
    u = uv(thin_dex,1);
    v = uv(thin_dex,2);

    if ( normalized == TRUE )
      norm = sqrt ( u.^2 + v.^2 );
      nonzero = find ( norm ~= 0.0 );
      u(nonzero) = u(nonzero) ./ norm(nonzero);
      v(nonzero) = v(nonzero) ./ norm(nonzero);
      vnorm_max = 1.0;
    end

    x(thin_num+1) = 0.00;
    y(thin_num+1) = 1.05;
    u(thin_num+1) = vnorm_max;
    v(thin_num+1) = 0.0;

    quiver ( x, y, u, v, scale )
%
%  Draw the boundary lines, and the invisible bounding line.
%
    line ( bx1, by1, 'color', 'r' )
    line ( bx2, by2, 'color', 'r' )
    line ( bx3, by3, 'color', 'w' )
  
    axis equal
    if ( normalized == TRUE )
      title ( 'T-Cell Direction Field' )
    else
      title ( 'T-Cell Flow Field' )
    end
    text ( 0.420, 1.03, upfile )
%
%  Here's where we take a snapshot of the current image, and save it to 
%  a PNG file.
%
    F = getframe;
    [X,map] = frame2im ( F );
    imwrite ( X, pngfile, 'PNG' );

  end