function clayton_chebfun_test01 ( )

%*****************************************************************************80
%
%% CLAYTON_CHEBFUN_TEST01 examines the gaussian correlation.
%
%  Discussion:
%
%    Note that I can't use the original value EIGEN_NUM = 20.  The 
%    EIGS function fails to converge.  In fact, it still fails at 12.
%
%    This code is based substantially on a document by Toby Driscoll.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Toby Driscoll,
%    Mercer's theorem and the Karhunen-Loeve expansion,
%    Oxford University Mathematical Institute,
%    http://www2.maths.ox.ac.uk/chebfun/examples/stats/pdf/MercerKarhunenLoeve.pdf
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CLAYTON_CHEBFUN_TEST01:\n' );
  fprintf ( 1, '  Demonstrate Mercer''s theorem and the KL expansion\n' );
  fprintf ( 1, '  for the gaussian kernel.\n' );
%
%  Set the interval.
%
  a = 0.0;
  b = 20.0;
  fprintf ( 1, '  Using interval [%g,%g]\n', a, b );
  s_num = 21;
  s_vec = linspace ( a, b, s_num );
%
%  FRED is a function from the CHEBFUN library.
%
%  It constructs a "chebop" representing the Fredholm integral operator
%  with kernel K for functions in domain [A,B].
%
  F = fred ( @gaussian_correlation, domain ( [ a, b ] ) );
%
%  EIGS has been extended to be able to compute the eigenvalues Lambda
%  and eigenfunctions Psi of the Fredholm integral operator represented by F.
%
%  The "LM" switch requests that we return the eigenvalues of largest magnitude.
%
%  Each Psi is a CHEBFUN, that is, it takes a real number argument.
%
  eigen_num = 20;
  [ Psi, Lambda ] = eigs ( F, eigen_num, 'lm' );

  eigen_found = length ( Lambda );
  fprintf ( 1, '  Requested %d eigenmodes, computed %d\n', eigen_num, eigen_found );
  eigen_num = min ( eigen_num, eigen_found );
%
%  Print the eigenvalues.
%
  lambda_vec = diag ( Lambda );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   I       Lambda(I)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : eigen_num
    fprintf ( 1, '  %2d  %10g\n', i, lambda_vec(i) );
  end
%
%  Plot selected eigenfunctions.
%
  figure ( 1 )
  subplot ( 5, 1, 1 )
  plot ( Psi(:,1), 'Linewidth', 2 );
  title ( 'gaussian 1.0: eigenfunction PSI(1)' )
  grid on
  subplot ( 5, 1, 2 )
  plot ( Psi(:,2), 'Linewidth', 2 );
  title ( 'gaussian 1.0: eigenfunction PSI(2)' )
  grid on
  subplot ( 5, 1, 3 )
  plot ( Psi(:,3), 'Linewidth', 2 );
  title ( 'gaussian 1.0: eigenfunction PSI(3)' )
  grid on
  subplot ( 5, 1, 4 )
  plot ( Psi(:,4), 'Linewidth', 2 );
  title ( 'gaussian 1.0: eigenfunction PSI(4)' )
  grid on
  subplot ( 5, 1, 5 )
  plot ( Psi(:,5), 'Linewidth', 2 );
  title ( 'gaussian 1.0: eigenfunction PSI(5)' )
  grid on
  print -dpng 'clayton_chebfun_gaussian_1.0.png'
%
%-------------------------------------
%  Repeat for correlation length 0.25.
%-------------------------------------
%
  F = fred ( @gaussian_correlation_025, domain ( [ a, b ] ) );

  eigen_num = 20;
  [ Psi, Lambda ] = eigs ( F, eigen_num, 'lm' );

  eigen_found = length ( Lambda );
  fprintf ( 1, '  Requested %d eigenmodes, computed %d\n', eigen_num, eigen_found );
  eigen_num = min ( eigen_num, eigen_found );

  lambda_vec = diag ( Lambda );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   I       Lambda(I)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : eigen_num
    fprintf ( 1, '  %2d  %10g\n', i, lambda_vec(i) );
  end

  figure ( 2 )
  subplot ( 5, 1, 1 )
  plot ( Psi(:,1), 'Linewidth', 2 );
  title ( 'gaussian 0.25: eigenfunction PSI(1)' )
  grid on
  subplot ( 5, 1, 2 )
  plot ( Psi(:,2), 'Linewidth', 2 );
  title ( 'gaussian 0.25: eigenfunction PSI(2)' )
  grid on
  subplot ( 5, 1, 3 )
  plot ( Psi(:,3), 'Linewidth', 2 );
  title ( 'gaussian 0.25: eigenfunction PSI(3)' )
  grid on
  subplot ( 5, 1, 4 )
  plot ( Psi(:,4), 'Linewidth', 2 );
  title ( 'gaussian 0.25: eigenfunction PSI(4)' )
  grid on
  subplot ( 5, 1, 5 )
  plot ( Psi(:,5), 'Linewidth', 2 );
  title ( 'gaussian 0.25: eigenfunction PSI(5)' )
  grid on
  print -dpng 'clayton_chebfun_gaussian_0.25.png'
  return
end
