function paraheat_pwc_plot ( diffusivity_seed )

%*****************************************************************************80
%
%% paraheat_pwc_plot plots paraheat_pwc sensor data.
%
%  Discussion:
%
%    This program uses the random number seed to choose parameters 
%    vc as input to the paraheat_pwc() program, and 
%    saves the resulting sensor measurements ns, xs, ys, vs.
%
%    For now, the program then simply uses the sensor measurements
%    to try to reconstruct the finite element solution and plot it.  
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer diffusivity_seed, the seed for the random number
%    generator, used to set the diffusivity values.
%
  addpath ( '../paraheat_pwc' );

  timestamp ( );
  fprintf ( 1, 'paraheat_pwc_plot\n' );
  fprintf ( 1, '  For given piecewise constant diffusivity data,\n' );
  fprintf ( 1, '  use paraheat_pwc() to compute an FEM solution u(x,y),\n' );
  fprintf ( 1, '  sample solution values at sensor locations,\n' );
  fprintf ( 1, '  use radial basis functions to reconstruct solution u.\n' );
  fprintf ( 1, '\n' );

  if ( nargin < 1 )

    fprintf ( 1, '  Diffusivity seed not specified by user, using default.\n' );
    diffusivity_seed = 123456789;
  end

  fprintf ( 1, '  Diffusivity random seed is %d\n', diffusivity_seed );
%
%  Set the diffusivity.
%
  rng ( diffusivity_seed );
  nxc = 4;
  nyc = 4;
  vc = 0.5 + 5.0 * rand ( nxc, nyc );
%
%  Solve the finite element system, and return samples of the solution.
%
  [ ns, xs, ys, vs ] = paraheat_pwc ( vc );
%
%  Create the RBF reconstruction.
%
  [ xmin, xmax, ymin, ymax ] = region_box ( );

  filename = sprintf ( 'fem_contour_rbf_%d.png', diffusivity_seed );
  fem_contour_rbf ( ns, xs, ys, vs, xmin, xmax, ymin, ymax, filename );

  filename = sprintf ( 'fem_surface_rbf_%d.png', diffusivity_seed );
  fem_surface_rbf ( ns, xs, ys, vs, xmin, xmax, ymin, ymax, filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'paraheat_pwc_plot\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  timestamp ( );

  rmpath ( '../paraheat_pwc' );

  return
end
function fem_contour_rbf ( ns, xs, ys, vs, xmin, xmax, ymin, ymax, filename );

%*****************************************************************************80
%
%% fem_contour_rbf reconstructs the FEM solution from sensor data.
%
%  Discussion:
%
%    A radial basis function approach is used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NS, the number of sensors.
%
%    Input, real XS(NS), YS(NS), the sensor coordinates.
%
%    Input, real VS(NS), the sensor readings.
%
%    Input, real XMIN, XMAX, YMIN, YMAX, the argument ranges.
%
%    Input, string FILENAME, the filename to use.
%
  figure ( 'Name', 'RBF interpolant to FEM sensor samples' );
  clf ( );
%
%  Determine the weights for the RBF interpolant.
%
  m = 2;
  xys = [ xs, ys ];
  xys = xys';
  volume = prod ( max ( xys, [], 2 ) - min ( xys, [], 2 ) );
  r0 = ( volume / ns ) ^ ( 1.0 / m );
  w = rbf_weight ( m, ns, xys, r0, @phi1, vs );
%
%  Evaluate the RBF interpolant on a grid.
%
  nix = 21;
  niy = 21;
  xi = linspace ( xmin, xmax, nix );
  yi = linspace ( ymin, ymax, niy );
  [ XI, YI ] = meshgrid ( xi, yi );
  nivec = nix * niy;
  xivec = reshape ( XI, nivec, 1 );
  yivec = reshape ( YI, nivec, 1 );
  xyvec = [ xivec'; yivec' ];
  zvec = rbf_interp ( 2, ns, xys, r0, @phi1, w, nivec, xyvec );
  ZI = reshape ( zvec, nix, niy );

  z_max = max ( zvec );
  z_min = min ( zvec );
  z_min = min ( z_min, 0.0 );
  z_range = z_max - z_min;
  caxis ( [ z_min, z_max ] );
%
%  Plot the RBF interpolant as a surface.
%
  [ M, c ] = contour ( XI, YI, ZI, 20 );
  c.LineWidth = 3;
  xlabel ( '<--- X --->', 'Fontsize', 16 );
  ylabel ( '<--- Y --->', 'Fontsize', 16 );
  zlabel ( '<---Z(X,Y)--->', 'Fontsize', 16 );
  filename2 = s_escape_tex ( filename );
  datestring = datestr ( datevec ( now ), 0 );
  title ( { filename2; datestring }, 'Fontsize', 16 );
  colorbar ( );
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s".\n', filename );
  pause ( 5 );
  close ( );

  return
end
function fem_surface_rbf ( ns, xs, ys, vs, xmin, xmax, ymin, ymax, filename )

%*****************************************************************************80
%
%% fem_surface_rbf reconstructs the FEM solution from sensor data.
%
%  Discussion:
%
%    A radial basis function approach is used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NS, the number of sensors.
%
%    Input, real XS(NS), YS(NS), the sensor coordinates.
%
%    Input, real VS(NS), the sensor readings.
%
%    Input, real XMIN, XMAX, YMIN, YMAX, the argument ranges.
%
%    Input, string FILENAME, the filename to use.
%
  figure ( 'Name', 'RBF interpolant to FEM sensor samples' );
  clf ( );
%
%  Determine the weights for the RBF interpolant.
%
  m = 2;
  xys = [ xs, ys ];
  xys = xys';
  volume = prod ( max ( xys, [], 2 ) - min ( xys, [], 2 ) );
  r0 = ( volume / ns ) ^ ( 1.0 / m );
  w = rbf_weight ( m, ns, xys, r0, @phi1, vs );
%
%  Evaluate the RBF interpolant on a grid.
%
  nix = 21;
  niy = 21;
  xi = linspace ( xmin, xmax, nix );
  yi = linspace ( ymin, ymax, niy );
  [ XI, YI ] = meshgrid ( xi, yi );
  nivec = nix * niy;
  xivec = reshape ( XI, nivec, 1 );
  yivec = reshape ( YI, nivec, 1 );
  xyvec = [ xivec'; yivec' ];
  zvec = rbf_interp ( 2, ns, xys, r0, @phi1, w, nivec, xyvec );
  ZI = reshape ( zvec, nix, niy );

  z_max = max ( zvec );
  z_min = min ( zvec );
  z_min = min ( z_min, 0.0 );
  z_range = z_max - z_min;
  caxis ( [ z_min, z_max ] );
%
%  Plot the RBF interpolant as a surface.
%
  surf ( XI, YI, ZI, 'EdgeColor', 'k' );
  xlabel ( '<--- X --->', 'Fontsize', 16 );
  ylabel ( '<--- Y --->', 'Fontsize', 16 );
  zlabel ( '<---Z(X,Y)--->', 'Fontsize', 16 );
  filename2 = s_escape_tex ( filename );
  datestring = datestr ( datevec ( now ), 0 );
  title ( { filename2; datestring }, 'Fontsize', 16 );
  daspect ( [ 1.0, 1.0, z_range ] );
  colorbar ( );
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s".\n', filename );
  pause ( 5 );
  close ( );

  return
end
function v = phi1 ( r, r0 )

%*****************************************************************************80
%
%% PHI1 evaluates the multiquadric radial basis function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    William Press, Brian Flannery, Saul Teukolsky, William Vetterling,
%    Numerical Recipes in FORTRAN: The Art of Scientific Computing,
%    Third Edition,
%    Cambridge University Press, 2007,
%    ISBN13: 978-0-521-88068-8,
%    LC: QA297.N866.
%
%  Parameters:
%
%    Input, real R(), the radial separation.
%    0 < R.
%
%    Input, real R0, a scale factor.
%
%    Output, real V(), the value of the radial basis function.
%
  v = sqrt ( r.^2 + r0^2 );

  return
end
function fi = rbf_interp ( m, nd, xd, r0, phi, w, ni, xi )

%*****************************************************************************80
%
%% RBF_INTERP evaluates a radial basis function interpolant.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    William Press, Brian Flannery, Saul Teukolsky, William Vetterling,
%    Numerical Recipes in FORTRAN: The Art of Scientific Computing,
%    Third Edition,
%    Cambridge University Press, 2007,
%    ISBN13: 978-0-521-88068-8,
%    LC: QA297.N866.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer ND, the number of data points.
%
%    Input, real XD(M,ND), the data points.
%
%    Input, real R0, a scale factor.  R0 should be larger than the typical
%    separation between points, but smaller than the maximum separation.
%    The value of R0 has a significant effect on the resulting interpolant.
%
%    Input, function V = PHI ( R, R0 ), a function handle to evaluate the radial
%    basis functions.
%
%    Input, real W(ND), the weights, as computed by RBF_WEIGHTS.
%
%    Input, integer NI, the number of interpolation points.
%
%    Input, real XI(M,NI), the interpolation points.
%
%    Output, real FI(NI), the interpolated values.
% 
  fi = zeros ( ni, 1 );

  if ( true )

    for i = 1 : ni
      if ( m == 1 )
        d = xd - xi(i);
        r = abs ( d );
        v = phi ( r(1:nd), r0 );
        fi(i) = v' * w;
      else
        d = xd - repmat ( xi(:,i), 1, nd );
        r = sqrt ( sum ( d.^2 ) );
        v = phi ( r(1:nd), r0 );
        fi(i) = v * w;
      end
    end

  else
%
%  Alternative, unvectorized code.
%
    for i = 1 : ni
      for j = 1 : nd
        r = norm ( xi(:,i) - xd(:,j) );
        fi(i) = fi(i) + w(j) * phi ( r, r0 );
      end
    end

  end

  return
end
function w = rbf_weight ( m, nd, xd, r0, phi, fd )

%*****************************************************************************80
%
%% RBF_WEIGHT computes weights for radial basis function interpolation.
%
%  Discussion:
%
%    We assume that there are N (nonsingular) equations in N unknowns.
%
%    However, it should be clear that, if we are willing to do some kind
%    of least squares calculation, we could allow for singularity,
%    inconsistency, or underdetermine systems.  This could be associated
%    with data points that are very close or repeated, a smaller number
%    of data points than function values, or some other ill-conditioning
%    of the system arising from a peculiarity in the point spacing.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    William Press, Brian Flannery, Saul Teukolsky, William Vetterling,
%    Numerical Recipes in FORTRAN: The Art of Scientific Computing,
%    Third Edition,
%    Cambridge University Press, 2007,
%    ISBN13: 978-0-521-88068-8,
%    LC: QA297.N866.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer ND, the number of data points.
%
%    Input, real XD(M,ND), the data points.
%
%    Input, real R0, a scale factor.  R0 should be larger than the typical
%    separation between points, but smaller than the maximum separation.
%    The value of R0 has a significant effect on the resulting interpolant.
%
%    Input, function V = PHI ( R, R0 ), a function handle to evaluate the radial
%    basis functions.
%
%    Input, real FD(ND,1), the function values at the data points.
%
%    Output, real W(ND,1), the weights.
%
  a = zeros ( nd, nd );
 
  for i = 1 : nd
    if ( m == 1 )
      d = xd - xd(i);
      r = abs ( d );
    else
      d = xd - repmat ( xd(:,i), 1, nd );
      r = sqrt ( sum ( d.^2 ) );
    end
    a(i,1:nd) = phi ( r(1:nd), r0 );
  end
%
%  Solve for the weights.
%
  w = a \ fd;

  return
end
function [ xmin, xmax, ymin, ymax ] = region_box ( )

%*****************************************************************************80
%
%% region_box returns the minimum and maximum x and y values of the region.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real XMIN, XMAX, YMIN, YMAX, the argument ranges.
%
  xmin = 0.0;
  xmax = 1.0;
  ymin = 0.0;
  ymax = 1.0;

  return
end
function s2 = s_escape_tex ( s1 )

%*****************************************************************************80
%
%% S_ESCAPE_TEX de-escapes TeX escape sequences.
%
%  Discussion:
%
%    In particular, every occurrence of the characters '\', '_',
%    '^', '{' and '}' will be replaced by '\\', '\_', '\^',
%    '\{' and '\}'.  A TeX interpreter, on seeing these character
%    strings, is then likely to return the original characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S1, the string to be de-escaped.
%
%    Output, string S2, a copy of the string, modified to avoid TeX escapes.
%
  s1_length = length ( s1 );

  s1_pos = 0;
  s2_pos = 0;
  s2 = [];

  while ( s1_pos < s1_length )

    s1_pos = s1_pos + 1;

    if ( s1(s1_pos) == '\' || ...
         s1(s1_pos) == '_' || ...
         s1(s1_pos) == '^' || ...
         s1(s1_pos) == '{' || ...
         s1(s1_pos) == '}' )
      s2_pos = s2_pos + 1;
      s2 = strcat ( s2, '\' );
    end

    s2_pos = s2_pos + 1;
    s2 = strcat ( s2, s1(s1_pos) );

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

