function plot_rd ( T, X, W, ca )

%*****************************************************************************80
%
%% plot_rd() makes two plots of the reaction-diffusion solution.
%
%  Discussion:
%
%    This function makes a "surface plot" of W(T,X), and
%    aÊplot of the spatial L2 norm of the solution over time.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    04 April 2011
%
%  Author:
%
%   Eugene Cliff
%
%  Reference:
%
%    Jeffrey Borggaard, John Burkardt, John Burns, Eugene Cliff,
%    Working Notes on a Reaction Diffusion Model: a Finite Element Formulation.
%
%  Input:
%
%    real T(nt), time values.
%
%    real X(nx), spatial values.
%
%    real W(NT,NX), solution array.
%
%    real CA, coefficient in nonlinear polynomial, used to label plot.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'plot_rd\n' );
%
% plot at most 100 time points
%
  it = ceil(length(T)/100);
  figure; 
  surf ( X, T(1:it:end), W(1:it:end,:) );
  xlabel('\bf axial location'); 
  ylabel('\bf time');
  zlabel('\bf w')
  vv = axis; 
  axis([X([1 end]) T([1 end]) vv([5 6])]);
 
  title_cell = cell(2,1); % Holds strings for two-line title
  title_cell{1} = ['\bf Reaction Diffusion: n = ' int2str(length(X)-1)];
  title_cell{2} = [' \bf Reaction Function: ' ...
                      fix_sign(ca(1),0)  ...
                      fix_sign(ca(2),1)  ...
                      fix_sign(ca(3),2)  ...
                      fix_sign(ca(4),3)];                           
  title(title_cell)
%
%  Compute and plot time history of solution norm 
%
  temp = zeros(ceil(length(T)/it),1);
  wtmp = W(1:it:end,:);
  for ii=1:length(temp)
    temp(ii) = l2_norm(wtmp(ii,:));
  end

  figure; 
  plot ( T(1:it:end), temp, 'r', 'LineWidth', 2 );
  grid ( 'on' );
  xlabel('time'); 
  ylabel('|| w ||^2'); 
  title(title_cell)
 
  return
end
function val = l2_norm(f)

%*****************************************************************************80
%
%% l2_norm() estimates the L2 norm of F using the trapezoidal rule.
%
%  Discussion:
%
%    The vector F is assumed to contain the values of a function at
%    equally spaced points in the interval [0,1].  In that case,
%    the L2 norm of F can be estimated by the trapezoidal rule.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    04 April 2011
%
%  Author:
%
%   Eugene Cliff
%
%  Input:
%
%    real F(N+1), function values.
%
%  Output:
%
%    real VAL, the integral of F.
%
  n = length ( f ) - 1;

  val = ( 0.5 * f(1)^2 + sum ( f(2:end-1).^2 ) + 0.5 * f(end)^2 ) / n;

  return
end
function str_out = fix_sign(c, p)

%*****************************************************************************80
%
%% fix_sign assembles a string consisting of sign(ca) * abs(ca) w^p
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    04 April 2011
%
%  Author:
%
%   Eugene Cliff
%
%  Input:
%
%    real C
%
%    integer P
%
%  Output:
%
%    string STR_OUT, the string.
%
  ps = int2str(p);

  if abs(c) <= eps
    str_out  =  '';
  elseif c > 0 
    str_out = [' + ' num2str(c)  ' w**' ps];
  else
    str_out = [' - ' num2str(-c) ' w**' ps];
  end

  return
end
