[n,x,knotyp,t]=pchkt(n,x,knotyp,t);

function [n,x,knotyp,t]=pchkt(n,x,knotyp,t);
%***BEGIN PROLOGUE  PCHKT
%***SUBSIDIARY
%***PURPOSE  Compute B-spline knot sequence for PCHBS.
%***LIBRARY   SLATEC (PCHIP)
%***CATEGORY  E3
%***TYPE      SINGLE PRECISION (PCHKT-S, DPCHKT-D)
%***AUTHOR  Fritsch, F. N., (LLNL)
%***DESCRIPTION
%
%     Set a knot sequence for the B-spline representation of a PCH
%     function with breakpoints X.  All knots will be at least double.
%     Endknots are set as:
%        (1) quadruple knots at endpoints if KNOTYP=0;
%        (2) extrapolate the length of end interval if KNOTYP=1;
%        (3) periodic if KNOTYP=2.
%
%  Input arguments:  N, X, KNOTYP.
%  Output arguments:  T.
%
%  Restrictions/assumptions:
%     1. N.GE.2 .  (not checked)
%     2. X(i).LT.X(i+1), i=1,...,N .  (not checked)
%     3. 0.LE.KNOTYP.LE.2 .  (Acts like KNOTYP=0 for any other value.)
%
%***SEE ALSO  PCHBS
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   870701  DATE WRITTEN
%   900405  Converted Fortran to upper case.
%   900410  Converted prologue to SLATEC 4.0 format.
%   900410  Minor cosmetic changes.
%   930514  Changed NKNOTS from an output to an input variable.  (FNF)
%   930604  Removed unused variable NKNOTS from argument list.  (FNF)
%***end PROLOGUE  PCHKT
%
%*Internal Notes:
%
%  Since this is subsidiary to PCHBS, which validates its input before
%  calling, it is unnecessary for such validation to be done here.
%
%**end
%
%  Declare arguments.
%
persistent hbeg hend j k ndim ; 

x_shape=size(x);x=reshape(x,1,[]);
t_shape=size(t);t=reshape(t,1,[]);
%
%  Declare local variables.
%
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(ndim), ndim=0; end;
if isempty(hbeg), hbeg=0; end;
if isempty(hend), hend=0; end;
%***FIRST EXECUTABLE STATEMENT  PCHKT
%
%  Initialize.
%
ndim = fix(2.*n);
%
%  Set interior knots.
%
j = 1;
for k = 1 : n;
j = fix(j + 2);
t(j) = x(k);
t(j+1) = t(j);
end; k = fix(n+1);
%     Assertion:  At this point T(3),...,T(NDIM+2) have been set and
%                 J=NDIM+1.
%
%  Set end knots according to KNOTYP.
%
hbeg = x(2) - x(1);
hend = x(n) - x(n-1);
if( knotyp==1 )
%          Extrapolate.
t(2) = x(1) - hbeg;
t(ndim+3) = x(n) + hend;
elseif( knotyp==2 ) ;
%          Periodic.
t(2) = x(1) - hend;
t(ndim+3) = x(n) + hbeg;
else;
%          Quadruple end knots.
t(2) = x(1);
t(ndim+3) = x(n);
end;
t(1) = t(2);
t(ndim+4) = t(ndim+3);
%
%  Terminate.
%
%------------- LAST LINE OF PCHKT FOLLOWS ------------------------------
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
end
%DECK PCHNGS