function visualize(u,T,fVis)
% Visualizes the trajectories generated from a system of ODEs for the
% n-body problem in 2D.
%
%   Input:
%     u     Matrix. It is expected that the i-th column u(:,i) represents
%           the value of the trajectory at time t_(i-1) and that for
%           0 <= k < n the values of the column describe:
%             u(1+k*4,i)  : x-position of body k+1
%             u(2+k*4,i)  : x-velocity of body k+1
%             u(3+k*4,i)  : y-position of body k+1
%             u(4+k*4,i)  : y-velocity of body k+1
%           
%     T     Scalar. The time up to which the trajectory was simulated.
%
%     fVis  True/False value (optional). Sets if plots should be visible
%           or not. Default is FALSE.
%
%   Output:
%     * No variable output.
%     * The following files are generated:
%         - trajectory_2D.png: Trajectories as seen on the Euclidean plane.
%         - trajectory_3D.png: Trajectories including the time component.
%         - trajectory.avi: Animation of the trajectories.
%
%   Remark:
%     * So far, this function does not use the velocity information at all.
%     * WARNING: The computation of the animation may be relatively time
%         consuming. Feel free to increase the frame skip or to remove the
%         corresponding part entirely. Also feel free to adapt the
%         frameRate and speedUp to your needs.
%
  fprintf('%s visualize.m was invoked.\n', datestr(now, '[HH:MM:SS]'));

  % Set skip. (Is used later to skip some frames during video generation.)
  % skip = 0 means no skip at all.
  skip = 50;
  
  % Set speedUp for the movie; total movie time is roughly T/speedUp seconds.
  speedUp = 100;
  
  % Set visibility.
  if(nargin < 3)
    fVis = false;
  end
  if(~fVis)
    sVis  = 'off';
  else
    sVis  = 'on';
  end
  
  % Extract general information.
  N   = size(u,2);      % Number of iterates.
  n   = size(u,1)/4;    % Number of bodies.

  % Extract position coordinates at all times.
  X   = u(1:4:end,:);
  Y   = u(3:4:end,:);



  %%%
  % Plot 1: 2D trajectories.
  %%%
  
  % Actual plotting.
  h2D = figure('Visible',sVis);
  for i = 1:n
    x = X(i,:);
    y = Y(i,:);
    z = 0*x;
    col = linspace(0,T,N);
    surface([x;x],[y;y],[z;z],[col;col],...
          'facecol','no',...
          'edgecol','interp',...
          'linew',2);
    hold on;
  end
  xlabel('x');
  ylabel('y');
  hc = colorbar;
  set(get(hc,'title'),'string','time');
  title(sprintf('2D plot of trajectories up to time T=%.2g',T));
  
  % Save plot and close if desired.
  print(h2D, '-dpng', 'trajectory_2D.png');
  if(~fVis)
    close(h2D);
  end



  %%%
  % Plot 2: 3D trajectories.
  %%%
  
  % Actual plotting.
  h3D = figure('Visible',sVis);
  for i = 1:n
    plot3(X(i,:),Y(i,:),linspace(0,T,N));
    hold on;
  end
  title(sprintf('3D plot of trajectories up to time T=%.2g',T));
  xlabel('x');
  ylabel('y');
  zlabel('time');
  
  % Save plot and close if desired.
  print(h3D, '-dpng', 'trajectory_3D.png');
  if(~fVis)
    close(h3D);
  end


  %%%
  % Plot 3: Animated 2D trajectories.
  %%%
  h = figure('Visible', 'off');

  lim_x = [-1+min(X(:)), 1+max(X(:))];
  lim_y = [-1+min(Y(:)), 1+max(Y(:))];

  % Compute and save video frames.
  v = VideoWriter('trajectory.avi');
  v.FrameRate = ceil(length(1:(skip+1):N)*speedUp/T); % Roughly try to match the frame rate with time.
  fprintf('%s Set movie frame rate to %d.\n', datestr(now, '[HH:MM:SS]'), v.FrameRate);
  open(v);
  for k = 1:(skip+1):N % Skip some frames to speed up loop a little.
    scatter(X(:,k),Y(:,k));
    title(sprintf('Frame %d/%d', k, N));
    xlim(lim_x);
    ylim(lim_y);
    writeVideo(v,getframe);
  end
  close(v);
  close(h);
  fprintf('%s Done computing frames.\n', datestr(now, '[HH:MM:SS]'));
end
