function paraheat_functional_plot ( diffusivity_seed )

%*****************************************************************************80
%
%% paraheat_functional_plot plots paraheat_functional sensor data.
%
%  Discussion:
%
%    This program sets up sample parameters sets "vc" as input to the
%    paraheat_functional() program, and saves the resulting sensor
%    measurements "vs".  
%
%    The sets of pairs (vc,vs) are written out to a file.
%
%    The "vc" values specify a piecewise constant diffusivity function,
%    which is used when solving a steady heat equation in a 2D domain with
%    zero Dirichlet boundary conditions and a right hand side function f(x,y).
%    The computed solution is then sampled at 50 sensor locations, which
%    are returned as the values "vs".
%
%    For details, see the paraheat_functional() program.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    17 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer diffusivity_seed, the seed for the random number
%    generator, used to set the diffusivity values.
%
  addpath ( '../rbf_interp_2d' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'paraheat_functional_plot\n' );
  fprintf ( 1, '  For given piecewise constant diffusivity data,\n' );
  fprintf ( 1, '  use paraheat 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' );

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

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Diffusivity random seed is %d\n', diffusivity_seed );
%
%  Set the diffusivity.
%
  rng ( diffusivity_seed );
  nxc = 4;
  nyc = 4;
  vc = 0.5 + 5 * rand ( nxc, nyc );
%
%  Solve the finite element system, and return samples of the solution.
%
  vs = paraheat_functional ( vc );
%
%  A hack to get sensor locations.
%
  sensor_location_seed = 987654321;
  [ ns, xs, ys ] = sensor_location_define ( sensor_location_seed );
%
%  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_functional_plot\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../rbf_interp_2d' );

  return
end
