% Chebyshev polynomial package
% programer: t. shores		date: 12/5/00
% description: this package contains a discussion of routines for
% implementing a Chebyshev polynomial package suitable for Octave and
% Matlab. 
% Matlab users will have to unpackage functions of this file into 
% separate function files for each of the routines below, while
% Octave users can simply read in this single file. 
% The basic structure is a row vector a=[a_0,...,a_n] of coefficients
% representing a polynomial p(x)=sum(a_k*T_k(x), k=0..n), where
% T_k(x) is the k-th Chebyshev polynomial.
%
% not a function file
1;

%

function retval = chebyeval(a,x)
% usage: y = chebyeval(a,x)
% description:  given input Chebyshev vector a and vector
% x of absicssas, this routine returns the vector of values of 
% corresponding Chebyshev polynomial at the coordinates of x.
 
% local variables: 
% j,k:  index variables
% m: size of coefficient array
% n: size of node array
% z: 2*x(j)
% xpar,xtmp: temporaries

m = length(a);
n = length(x);
retval = zeros(size(x));
for j = 1:n % loop over entries of x 
  % use Clenshaw's method to evaluate the polynomial
  z = 2*x(j);
  xpar = zeros(1,3);
  for k = 1:m
    xtmp = z*xpar(3) - xpar(2) + a(k);
    xpar = [xpar(2:3), xtmp];
  end;
  retval(j) = (xpar(3)-xpar(1))/2;
end;

endfunction % Matlab users must cut above this line for a function file

%

function retval = chebycoef(fn_, range,nn)
% usage: y = chebycoef('fcn', range,nodenum)
% description:  given input function fcn, range of non-negative indices
% range, and positive integer nn, this routine returns an array y
% the size of range with each entry the corresponding Chebychev
% coefficient calculated by using an nn-point trapezoidal method.
% The function fcn must handle vector 
% arguments correctly: evaluate each coordinate.

% local variables: 
% nrange: size of range
% dtheta: 2/nn
% thetas: vector of cosines
% j:  index variables

retval = zeros(size(range));
nrange = length(range);
dtheta = 2/nn;
costhetas = cos((0:nn)*pi/nn);
for j = 1:nrange  % loop over indices of range 
  retval(j) = dtheta*trapz(feval(fn_,costhetas).*cos((0:nn)*range(j)*pi/nn));
end;

endfunction % Matlab users must cut above this line for a function file

%

function retval = chebyinterp(fn_,n,x)
% usage: z = chebyinterp(fcn,x)
% description:  given function fcn, non-negative integer n, and arbitrary
% vector x, this routine returns z = p(x), where p is the polynomial that 
% interpolates fcn with values y at the zeros of T_n(x), listed in 
% increasing order. The function fcn must handle vector 
% arguments correctly: evaluate each coordinate.

% local variables: 
% xnodes: ordered zeros of T_n(x)
% ynodes: values of fcn at xnodes

retval = zeros(size(x));
xnodes = cos(pi/(2*n)*(2*n-1:-2:1));
ynodes = feval(fn_,xnodes);
retval = Interp(xnodes, DividDif(xnodes,ynodes), x);

endfunction % Matlab users must cut above this line for a function file

%

function retval = chebyder(a)
% usage: b = PPder(a)
% description: given input a vector of Chebyshev polynomial 
% coefficients, this function returns a vector of a Chebyshev
% polynomial of the same size representing the derivative 
% of the function represented by a.

% local variables: 
% j,k:  index variables

% just a stub for now

endfunction % Matlab users must cut above this line for a function file

%

function retval = chebyint(fn_)
% usage retval = chebyint('fcn')
% description:  takes a function of one variable fcn  and returns an 
% approximate value for the definite integral of fcn(x)/sqrt(1-x^2) on 
% interval [-1,1], calculated by using a 5-point Chebyshev-Gauss
% quadrature method .  The function fcn must handle vector 
% arguments correctly: evaluate each coordinate.

% local variables:
% xabsc: Chebyshev abscissas 
% y: values of fcn at abscissas

xabsc = [-0.951056516295154,-0.587785252292473,0.0, \
              0.587785252292473, 0.951056516295154];
y = feval(fn_, xabsc);
retval = pi/5*sum(y);

endfunction % Matlab users must cut above this line for a function file.

%

function retval = DividDif(x,y)
% usage: DividDif(x,y)
% description: for input vector x of abscissas and y of ordinates
% of function f(x), returns vector of differences f[x_1], f[x_1,x_2], etc.

% local variables
% j,k: index variables
% nn,mm: size of input vectors x,y

nn = length(x);
mm = length(y);
if (mm ~= nn)
  disp('There is a size mismatch');
  return;
end;
retval = y;
for j = (2:nn)
  for k =(nn:-1:j)
    retval(k)=(retval(k)-retval(k-1))/(x(k)-x(k-j+1));
  end;
end;

endfunction % Matlab users must cut above this line for a function file.

%

function retval = Interp(x,d,t)
% usage: Interp(x,d,t)
% description: for input vectors x of abscissas, d of divided
% differences of function f, and vector t of arguments, 
% returns p(t) where p is the interpolating polynomial.

% local variables
% j,k: index variables
% nn,mm,pp: size of input vectors x,d,t

nn = length(x);
mm = length(d);
if (mm ~= nn)
  disp('There is a size mismatch');
  return;
end;
pp = length(t);
retval = zeros(size(t));
for k = (1:pp)
  retval(k)= d(nn);
  for j = (nn-1:-1:1)
    retval(k) = d(j) + (t(k) - x(j))*retval(k);
  end;
end;

endfunction % Matlab users must cut above this line for a function file.