[xx,n,x,c,d,work]=polcof(xx,n,x,c,d,work);

function [xx,n,x,c,d,work]=polcof(xx,n,x,c,d,work);
persistent i im1 k km1 km1pi km2n km2npi nm1 nmkp1 npkm1 pone ptwo ; 

if isempty(pone), pone=0; end;
if isempty(ptwo), ptwo=0; end;
if isempty(i), i=0; end;
if isempty(im1), im1=0; end;
if isempty(k), k=0; end;
if isempty(km1), km1=0; end;
if isempty(km1pi), km1pi=0; end;
if isempty(km2n), km2n=0; end;
if isempty(km2npi), km2npi=0; end;
if isempty(nm1), nm1=0; end;
if isempty(nmkp1), nmkp1=0; end;
if isempty(npkm1), npkm1=0; end;
%***BEGIN PROLOGUE  POLCOF
%***PURPOSE  Compute the coefficients of the polynomial fit (including
%            Hermite polynomial fits) produced by a previous call to
%            POLINT.
%***LIBRARY   SLATEC
%***CATEGORY  E1B
%***TYPE      SINGLE PRECISION (POLCOF-S, DPOLCF-D)
%***KEYWORDS  COEFFICIENTS, POLYNOMIAL
%***AUTHOR  Huddleston, R. E., (SNLL)
%***DESCRIPTION
%
%     Written by Robert E. Huddleston, Sandia Laboratories, Livermore
%
%     Abstract
%        subroutine POLCOF computes the coefficients of the polynomial
%     fit (including Hermite polynomial fits ) produced by a previous
%     call to POLINT. The coefficients of the polynomial, expanded about
%     XX, are stored in the array D. The expansion is of the form
%     P(Z) = D(1) + D(2)*(Z-XX) +D(3)*((Z-XX)**2) + ... +
%                                                  D(N)*((Z-XX)**(N-1)).
%     Between the call to POLINT and the call to POLCOF the variable N
%     and the arrays X and C must not be altered.
%
%     *****  INPUT PARAMETERS
%
%     XX   - The point about which the Taylor expansion is to be made.
%
%     N    - ****
%            *     N, X, and C must remain unchanged between the
%     X    - *     call to POLINT or the call to POLCOF.
%     C    - ****
%
%     *****  OUTPUT PARAMETER
%
%     D    - The array of coefficients for the Taylor expansion as
%            explained in the abstract
%
%     *****  STORAGE PARAMETER
%
%     WORK - This is an array to provide internal working storage. It
%            must be dimensioned by at least 2*N in the calling program.
%
%
%     **** Note - There are two methods for evaluating the fit produced
%     by POLINT. You may call POLYVL to perform the task, or you may
%     call POLCOF to obtain the coefficients of the Taylor expansion and
%     then write your own evaluation scheme. Due to the inherent errors
%     in the computations of the Taylor expansion from the Newton
%     coefficients produced by POLINT, much more accuracy may be
%     expected by calling POLYVL as opposed to writing your own scheme.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   890213  DATE WRITTEN
%   891024  Corrected KEYWORD section.  (WRB)
%   891024  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%***end PROLOGUE  POLCOF
%
x_shape=size(x);x=reshape(x,1,[]);
c_shape=size(c);c=reshape(c,1,[]);
d_shape=size(d);d=reshape(d,1,[]);
work_shape=size(work);work=reshape(work,1,[]);
%***FIRST EXECUTABLE STATEMENT  POLCOF
for k = 1 : n;
d(k) = c(k);
end; k = fix(n+1);
if( n==1 )
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
work(1) = 1.0;
pone = c(1);
nm1 = fix(n - 1);
for k = 2 : n;
km1 = fix(k - 1);
npkm1 = fix(n + k - 1);
work(npkm1) = xx - x(km1);
work(k) = work(npkm1).*work(km1);
ptwo = pone + work(k).*c(k);
pone = ptwo;
end; k = fix(n+1);
d(1) = ptwo;
if( n==2 )
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
for k = 2 : nm1;
km1 = fix(k - 1);
km2n = fix(k - 2 + n);
nmkp1 = fix(n - k + 1);
for i = 2 : nmkp1;
km2npi = fix(km2n + i);
im1 = fix(i - 1);
km1pi = fix(km1 + i);
work(i) = work(km2npi).*work(im1) + work(i);
d(k) = d(k) + work(i).*d(km1pi);
end; i = fix(nmkp1+1);
end; k = fix(nm1+1);
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
end
%DECK POLFIT