function [ xp, yp ] = interp_chebyshev ( f_string, a, b, n )

%*****************************************************************************80
%
%% interp_chebyshev() is an interactive Chebyshev spaced node interpolation program.
%
%  Discussion:
%
%    This function manages an interactive computation which uses
%    n Chebyshev spaced nodes to interpolate a function f(x)
%    in the interval a <= x <= b.
%
%    The user enters a formula for f(x), the values of a and b,
%    and the number of interpolation points n.
%
%    The program returns a set of 101 pairs of values (xp,yp) that can be used
%    to plot the interpolating polynomial.
%
%    The program can be invoked by a function call, in which case the
%    string specifying f(x) must be quoted:
%      [ xp, yp ] = interp_chebyshev ( 'x^2', -1, 3, 2 );
%    or, if called with no arguments, it will request them:
%      [ xp, yp ] = interp_chebyshev ( );
%      Enter function formula, like x^2: x^2
%      Enter left limit, a: -1
%      Enter right limit, b: 3
%      Enter number of interpolation points n: 2
%
%    The function is specified as a string which is either:
%    * a MATLAB expression using the argument 'x', 
%    * the name of an M-file followed by the argument '(x)'.
%
%    The string should not contain any spaces between symbols, except
%    when it is passed as a function argument in quotes.
%
%    It is not necessary to use the "dot" notation for expressions involving
%    '*', '/', or '^', but it doesn't hurt either.
%
%    Examples of function specifications:
%
%      x^2
%      x.^2
%      3/(x^4+5*x-6)
%      sin(7*x)*sqrt(x)/8
%      wiggle(x)     <-- where "wiggle.m" is a user-provided M file.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    29 June 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string f_string: the function to be interpolated.
%
%    real a, b: the interval endpoints.
%
%    integer n: the number of Chebyshev points to use.
%
%  Output:
%
%    real xp(), yp(): samples of the interpolant, which can be plotted.
%

%
%  Get F, A, B, N.
%
  if ( nargin < 1 )
    f_string = input ( 'Enter function formula, like x^2: ', 's' );
  end
  f_element = vector_to_element ( f_string );
  f = str2func ( [ '@(x)', f_element ] );

  if ( nargin < 2 )
    a = input ( 'Enter left limit a: ' );
  elseif ( ischar ( a ) )
    a = str2num ( a );
  end

  if ( nargin < 3 )
    b = input ( 'Enter right limit b: ' );
  elseif ( ischar ( b ) )
    b = str2num ( b );
  end

  if ( nargin < 4 )
    n = input ( 'Enter number of Chebyshev points n: ' );
  elseif ( ischar ( n ) )
    n = str2num ( n );
  end

  xd = chebyspace ( a, b, n );
  yd = f ( xd );
  dd = divdif ( xd, yd );

  np = 101;
  xp = linspace ( a, b, np );
  yp = interp ( xd, dd, xp );

  return
end
function x = chebyspace ( a, b, n )

%*****************************************************************************80
%
%% chebyspace creates a vector of Chebyshev spaced values.
%
%  Discussion:
%
%    A set of points with Chebyshev spacing is generated for the
%    interval [A,B].
%
%    An R8VEC is a vector of R8's.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 June 2011
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N, the number of entries in the vector.
%
%    real A, B, the interval endpoints.
%
%  Output:
%
%    real X(N,1), a vector of Chebyshev spaced data.
%
  x = zeros ( n, 1 );

  if ( n == 1 )

    x(1) = ( a + b ) / 2.0;

  else

    for i = 1 : n

      theta = ( n - i ) * pi / ( n - 1 );

      c = cos ( theta );

      if ( mod ( n, 2 ) == 1 )
        if ( 2 * i - 1 == n )
          c = 0.0;
        end
      end

      x(i) = ( ( 1.0 - c ) * a  ...
             + ( 1.0 + c ) * b ) ...
             /   2.0;

    end

  end

  return
end
function yd = divdif ( x, y )

%*****************************************************************************80
%
%% divdif computes the divided differences of a set of data.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    29 June 2019
%
%  Reference:
%
%    Kendall Atkinson, Weimin Han,
%    Elementary Numerical Analysis,
%    Wiley, 2004,
%    ISBN: 0471433373,
%    LC: QA297.A83.2004.
%
%  Input:
%
%    x, y, two real vectors of length n, containing the x and y
%    coordinates of data.
%
%  Output:
%
%    yd, a real vector of length n, containing the divided difference
%    table for the data.
%
  yd = y;
  n = length ( x );

  for i = 2 : n
    for j = n : -1 : i
      yd(j) = ( yd(j) - yd(j-1) ) / ( x(j) - x(j-i+1) );
    end
  end

  return
end

function yp = interp ( xd, yd, xp )

%*****************************************************************************80
%
%% interp evaluates an interpolating polynomial in divided difference form.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    29 June 2019
%
%  Reference:
%
%    Kendall Atkinson, Weimin Han,
%    Elementary Numerical Analysis,
%    Wiley, 2004,
%    ISBN: 0471433373,
%    LC: QA297.A83.2004.
%
%  Input:
%
%    xd, yd, two real vectors of length nd, containing the x coordinates of
%    the data, and the divided differences of the data.
%
%    xp, a real array of length np, evaluation points for the interpolating
%    polynomial.
%
%  Output:
%
%    yp, a real vector of length np, the value of the interpolating polynomial
%    at the points xp.
%
  nd = length ( xd );

  yp = yd(nd) * ones ( size ( xp ) );
  for i = nd - 1 : -1 : 1
    yp = yd(i) + ( xp - xd(i) ) .* yp;
  end

  return
end
function f_element = vector_to_element ( f_vector )

%*****************************************************************************80
%
%% vector_to_element rewrites a MATLAB vector function elementwise.
%
%  Discussion:
%
%    In MATLAB, by default, the operators '*', '/', and '^' are interpreted
%    as linear algebraic operations.  In many cases, a user simply desired
%    the elementwise operation instead, which MATLAB indicates using a 
%    dotted notation: '.*', './', and '.^'.
%
%    This function converts a string to a form that uses dotted operators.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    29 June 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string f_vector, a string presumably representing a MATLAB expression.
%
%  Output:
%
%    string f_element, a copy of f_vector in which all operators have been
%    replaced by their dotted (elementwise) form.
%
  n_vector = length ( f_vector );
  f_element = [''];

  c_prev = '';
  for iv = 1 : n_vector
    c = f_vector(iv);
    if ( ( c == '*' || c == '/' || c == '^' ) && ( c_prev ~= '.' ) )
      f_element = [ f_element, ' .', c, ' ' ];
    else
      f_element = [ f_element, c ];
    end
    c_prev = c;
  end

  return
end

