[l,c,tc,a]=pcoef(l,c,tc,a);

function [l,c,tc,a]=pcoef(l,c,tc,a);
persistent fac i ll llp1 llp2 new nr savemlv ; 

if isempty(fac), fac=0; end;
if isempty(savemlv), savemlv=0; end;
if isempty(i), i=0; end;
if isempty(ll), ll=0; end;
if isempty(llp1), llp1=0; end;
if isempty(llp2), llp2=0; end;
if isempty(new), new=0; end;
if isempty(nr), nr=0; end;
%***BEGIN PROLOGUE  PCOEF
%***PURPOSE  Convert the POLFIT coefficients to Taylor series form.
%***LIBRARY   SLATEC
%***CATEGORY  K1A1A2
%***TYPE      SINGLE PRECISION (PCOEF-S, DPCOEF-D)
%***KEYWORDS  CURVE FITTING, DATA FITTING, LEAST SQUARES, POLYNOMIAL FIT
%***AUTHOR  Shampine, L. F., (SNLA)
%           Davenport, S. M., (SNLA)
%***DESCRIPTION
%
%     Written BY L. F. Shampine and S. M. Davenport.
%
%     Abstract
%
%     POLFIT  computes the least squares polynomial fit of degree  L  as
%     a sum of orthogonal polynomials.  PCOEF  changes this fit to its
%     Taylor expansion about any point  C , i.e. writes the polynomial
%     as a sum of powers of (X-C).  Taking  C=0.  gives the polynomial
%     in powers of X, but a suitable non-zero  C  often leads to
%     polynomials which are better scaled and more accurately evaluated.
%
%     The parameters for  PCOEF  are
%
%     INPUT --
%         L -      Indicates the degree of polynomial to be changed to
%                  its Taylor expansion.  To obtain the Taylor
%                  coefficients in reverse order, input  L  as the
%                  negative of the degree desired.  The absolute value
%                  of L  must be less than or equal to NDEG, the highest
%                  degree polynomial fitted by  POLFIT .
%         C -      The point about which the Taylor expansion is to be
%                  made.
%         A -      Work and output array containing values from last
%                  call to  POLFIT .
%
%     OUTPUT --
%         TC -     Vector containing the first LL+1 Taylor coefficients
%                  where LL=ABS(L).  If  L.GT.0 , the coefficients are
%                  in the usual Taylor series order, i.e.
%                    P(X) = TC(1) + TC(2)*(X-C) + ... + TC(N+1)*(X-C)**N
%                  If L .LT. 0, the coefficients are in reverse order,
%                  i.e.
%                    P(X) = TC(1)*(X-C)**N + ... + TC(N)*(X-C) + TC(N+1)
%
%***REFERENCES  L. F. Shampine, S. M. Davenport and R. E. Huddleston,
%                 Curve fitting by polynomials in one variable, Report
%                 SLA-74-0270, Sandia Laboratories, June 1974.
%***ROUTINES CALLED  PVALUE
%***REVISION HISTORY  (YYMMDD)
%   740601  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  PCOEF
%
a_shape=size(a);a=reshape(a,1,[]);
tc_shape=size(tc);tc=reshape(tc,1,[]);
%***FIRST EXECUTABLE STATEMENT  PCOEF
ll = fix(abs(l));
llp1 = fix(ll + 1);
ll_orig=ll;    [ll,dumvar2,c,dumvar4,dumvar5,a]=pvalue(ll,ll,c,tc(1),tc(sub2ind(size(tc),max(2,1)):end),a);    ll(dumvar2~=ll_orig)=dumvar2(dumvar2~=ll_orig);   dumvar4i=find((tc(1))~=(dumvar4));dumvar5i=find((tc(sub2ind(size(tc),max(2,1)):end))~=(dumvar5));   tc(1-1+dumvar4i)=dumvar4(dumvar4i); tc(2-1+dumvar5i)=dumvar5(dumvar5i); 
if( ll>=2 )
fac = 1.0;
for i = 3 : llp1;
fac = fac.*(i-1);
tc(i) = tc(i)./fac;
end; i = fix(llp1+1);
end;
if( l<0 )
nr = fix(fix(llp1./2));
llp2 = fix(ll + 2);
for i = 1 : nr;
savemlv = tc(i);
new = fix(llp2 - i);
tc(i) = tc(new);
tc(new) = savemlv;
end; i = fix(nr+1);
end;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
tc_shape=zeros(tc_shape);tc_shape(:)=tc(1:numel(tc_shape));tc=tc_shape;
end
%DECK PFQAD