function interpolateSine(max_degree)
% interpolateSine: Plots some polynomial interpolants of the sin() function.
%
%   interpolateSine(max_degree)
%
% This function constructs some polynomial interpolants of the sin(20 pi x)
% function using the monomial basis, and plots the results.  The
% interpolants start with degree zero (a constant) and work their way up to
% degree max_degree.  Interpolation points start at x = 0 and are chosen
% equally spaced so that the max_degree point is 1.
%
% Also plotted is the error in the interpolant.
%
% Input Parameters:
%
%   max_degree - Maximum degree of polynomial interpolant.  Default = 10.
%
% Output Parameters:
%
%   None.
%
% Ian Mitchell for CPSC 303 (January 10, 2006).

  % Assign default inputs.
  if(nargin < 1)
    max_degree = 10;
  end
  
  % What is our function?
  f = inline('sin(20 * pi * x)', 'x');
  
  % Interpolation abscissae.  Not all will be used for every interpolant.
  % Apply transpose to get column vector.
  xs = linspace(0, 1, max_degree + 1)';
  
  % Interpolation data values.
  ys = feval(f, xs);

  % Get a new figure window.
  figure;
  
  % Loop over the polynomial degrees.  
  for n = 0 : max_degree

    % What data points will we use this time?
    x = xs(1:n+1);
    y = ys(1:n+1);
    
    % Basis function matrix.  Note that Matlab's vander() function
    % returns the columns in the opposite order than that used in class.
    A = vander(x);
    
    % Find coefficients.
    c = A \ y;
    
    % Range of interest for the plot.
    plot_range_x = [ min(x) - 0.25, max(x) + 0.25 ];
    plot_x = linspace(plot_range_x(1), plot_range_x(2), 201);
    
    % Clear the figure window.
    clf;
    
    % Plot the function we are trying to match.
    subplot(2,1,1);
    f_y = feval(f, plot_x);
    hf = plot(plot_x, f_y, 'b-');
    hold on;
    
    % Plot the interpolation points.
    hp = plot(x, y, 'r*');
    
    % Also plot the interpolant between interpolation points, to compare
    % results.
    plot_y = polyval(c, plot_x);
    hi = plot(plot_x, plot_y, 'r:');

    % Label the plot.
    ylabel('f(x)');
    title('Monomial interpolation example - Ian Mitchell for CPSC 303');
    legend('sin(20 \pi x)', 'interpolant');
    
    % Easier to see lines.
    set([ hf; hi ], 'LineWidth', 2);
    set(hp, 'MarkerSize', 9);

    % Tighten up the bounds so that we can see the original function.
    axis([ plot_range_x, -1.2, +1.2 ]);
       
    % Plot the error.  We will only look at magnitude of error, and we will use
    % a separate subplot so that we can examine the error with logarithmic
    % scaling.  Note how error gets much smaller near the interpolation
    % points (it goes to exactly zero at the interpolation points).
    subplot(2,1,2);
    plot_error = abs(feval(f, plot_x) - plot_y);
    he = semilogy(plot_x, plot_error);

    % Matlab seems to use weird axis bounds for this plot, so we will
    % tighten them up.  First, get the current bounds so as to keep the
    % range in the y-axis.
    axis_error = axis;
    % Now adjust the range in the x-axis to agree with the range for the
    % plot above.
    axis_error(1:2) = plot_range_x;
    % And finally, reset the axis.
    axis(axis_error);
    grid on;
    
    % Label the plot.  Both axes share the same label in the x-axis.
    % Only one line in this plot, so y-axis label substitutes for a legend.
    ylabel('error');
    xlabel('x');

    set(he, 'LineWidth', 2);
    
    % Wait for the user to digest the current plot.
    pause;
  end