function output = spline(x,y,xx)
%MATLAB Code Generation Library Function

%   Limitations:
%   1. Input X must be strictly increasing.
%   2. Y entries with NaN values are not removed, and if Y contains
%      end slopes, no error will be issued if any are Inf.
%   3. For code generation with two inputs, the output struct cannot be
%      used directly by PPVAL in MATLAB. See MKPP for instructions on how
%      to convert it for use outside of a MATLAB Coder project.
%   4. With vector XX inputs, the orientation of vector outputs may not match
%      MATLAB unless Y is a fixed-size or variable-length vector.

%   Copyright 1984-2018 The MathWorks, Inc.
%#codegen

eml_invariant(nargin >= 2,'MATLAB:minrhs');
if eml_is_const(iscolumn(y)) && iscolumn(y)
    % Make y a row vector.
    pp = splinepp(x,y(:).');
else
    pp = splinepp(x,y);
end
if nargin == 2
    output = pp;
else
    eml_invariant(isreal(xx), ...
        'MATLAB:pchip:ComplexInterpPts');
    output = ppval(pp,xx);
end

%--------------------------------------------------------------------------

function pp = splinepp(x,y)
chckxy(x,y);
nx = coder.internal.indexInt(numel(x));
ndy = coder.internal.indexInt(eml_ndims(y));
npages = coder.internal.indexInt(size(y,ndy));
has_endslopes = (npages == nx + 2);
if nx <= 2 || (nx <= 3 && ~has_endslopes)
    pp = degenerateSpline(x,y);
    return
end
pglen = coder.internal.prodsize(y,'below',ndy);
nxm1 = nx - 1;
ZERO = coder.internal.indexInt(0);
type = coder.internal.scalarEg(x,y);
YISROW = coder.const(eml_is_const(isrow(y)) && isrow(y));
nxm2 = nx - 2;
nxm3 = nx - 3;
szdvdf = size(y);
szs = size(y);
if has_endslopes
    szdvdf(ndy) = size(y,ndy) - 3;
    szs(ndy) = size(y,ndy) - 2;
    yoffset = pglen;
else
    szdvdf(ndy) = size(y,ndy) - 1;
    szs(ndy) = size(y,ndy);
    yoffset = ZERO;
end
if YISROW
    s = coder.nullcopy(zeros(1,szs(2),'like',type));
    dvdf = coder.nullcopy(zeros(1,szdvdf(2),'like',type));
else
    s = coder.nullcopy(zeros(szs,'like',type));
    dvdf = coder.nullcopy(zeros(szdvdf,'like',type));
end
dx = coder.nullcopy(x(1:end-1));
for k = 1:nxm1
    dx(k) = x(k + 1) - x(k);
    dpg = (k - 1)*pglen;
    pg = yoffset + dpg;
    pgp1 = yoffset + k*pglen;
    for j = 1:pglen
        dvdf(dpg + j) = (y(pgp1 + j) - y(pg + j))/dx(k);
    end
end
for k = 2:nxm1
    pg = (k - 1)*pglen;
    pgm1 = (k - 2)*pglen;
    d1 = dx(k);
    d2 = dx(k - 1);
    for j = 1:pglen
        s(pg + j)= 3*(d1*dvdf(pgm1 + j) + d2*dvdf(pg + j));
    end
end
if has_endslopes
    d31 = zeros('like',x);
    dnnm2 = zeros('like',x);
    for j = 1:pglen
        s(j)= dx(2)*y(j);
    end
    pgend2 = (nx + 1)*pglen;
    pglast = nxm1*pglen;
    for j = 1:pglen
        s(pglast + j) = dx(nxm2)*y(pgend2 + j);
    end
else
    d31 = x(3) - x(1);
    dnnm2 = x(nx) - x(nx - 2);
    d1 = dx(1);
    d2 = dx(2);
    for j = 1:pglen
        s(j) = ((d1 + 2*d31)*d2*dvdf(j) + d1*d1*dvdf(pglen + j))/d31;
    end
    pg = nxm1*pglen;
    pgm1 = nxm2*pglen;
    pgm2 = nxm3*pglen;
    d1 = dx(nxm1);
    d2 = dx(nxm2);
    for j = 1:pglen
        s(pg + j)= ((d1 + 2*dnnm2)*d2*dvdf(pgm1 + j) + d1*d1*dvdf(pgm2 + j))/dnnm2;
    end
end
md = coder.nullcopy(x);
md(1) = dx(2);
md(nx) = dx(nxm2);
for k = 2:nxm1
    md(k) = 2*(dx(k) + dx(k - 1));
end
% Solve the tridiagonal system.
% The superdiagonal is [d31,dx(1:nx-2)].
% The main diagonal is md.
% The subdiagonal is [dx(2:nx-1),dnnm2].
% Reduce the subdiagonal.
r = dx(2)/md(1);
md(2) = md(2) - r*d31;
for j = 1:pglen
    s(pglen + j) = s(pglen + j) - r*s(j);
end
for k = 3:nxm1
    r = dx(k)/md(k - 1);
    md(k) = md(k) - r*dx(k - 2);
    pg = (k - 1)*pglen;
    pgm1 = (k - 2)*pglen;
    for j = 1:pglen
        s(pg + j) = s(pg + j) - r*s(pgm1 + j);
    end
end
r = dnnm2/md(nxm1);
md(nx) = md(nx) - r*dx(nxm2);
pg = nxm1*pglen;
pgm1 = nxm2*pglen;
for j = 1:pglen
    s(pg + j) = s(pg + j) - r*s(pgm1 + j);
end
% Now do the backsolve.
pg = nxm1*pglen;
for j = 1:pglen
    s(pg + j) = s(pg + j)/md(nx);
end
for k = nxm1:-1:2
    pg = (k - 1)*pglen;
    pgp1 = k*pglen;
    d1 = dx(k - 1);
    for j = 1:pglen
        s(pg + j) = (s(pg + j) - d1*s(pgp1 + j))/md(k);
    end
end
for j = 1:pglen
    s(j) = (s(j) - d31*s(pglen + j))/md(1);
end
pp = pwchcore(x,y,yoffset,s,dx,dvdf);

%--------------------------------------------------------------------------

function pp = degenerateSpline(x,y)
% For numel(x) == 2 or numel(x) == 3 (without endslopes).
nx = coder.internal.indexInt(numel(x));
ndy = coder.internal.indexInt(eml_ndims(y));
npages = coder.internal.indexInt(size(y,ndy));
pglen = coder.internal.prodsize(y,'below',ndy);
has_endslopes = (npages == nx + 2);
type = coder.internal.scalarEg(x,y);
YISROW = coder.const(eml_is_const(isrow(y)) && isrow(y));
if nx <= 2
    % Note that nx < 2 is excluded by chckxy(x,y), so nx == 2.
    if has_endslopes
        % The interpolant is the cubic Hermite polynomial.
        numterms = 4;
    else
        % The interpolant is a straight line.
        numterms = 2;
    end
else
    % The interpolant is a parabola.
    numterms = 3;
end
% Allocate coefficient array.
if YISROW
    c = coder.nullcopy(zeros(1,numterms,'like',type));
else
    szc = [size(y),0];
    szc(ndy) = 1;
    szc(ndy + 1) = numterms;
    c = coder.nullcopy(zeros(szc,'like',type));
end
if nx <= 2
    if has_endslopes
        % The interpolant is the cubic Hermite polynomial.
        if YISROW
            endslopes = coder.nullcopy(zeros(1,2,'like',y));
        else
            szes = size(y);
            szes(ndy) = 2;
            endslopes = coder.nullcopy(zeros(szes,'like',y));
        end
        pgend = (npages - 1)*pglen;
        for j = 1:pglen
            endslopes(j) = y(j);
            endslopes(pglen + j) = y(pgend + j);
        end
        % We can't use the output of pwchcore directly because the
        % dimensionality of pp.coefs might be less than that of other
        % cases, and then the compiler would make pp.coefs an N-D array
        % with a trailing singleton dimension in this case. That would not
        % be correctly interpreted by ppval.
        pp1 = pwchcore(x,y,pglen,endslopes);
        c(:) = pp1.coefs(:);
    else
        % The interpolant is a straight line.
        dx = x(2) - x(1);
        for j = 1:pglen
            c(j) = (y(pglen + j) - y(j))/dx;
            c(pglen + j) = y(j);
        end
    end
    pp = makepp(x,c);
else % nx == 3 && ~has_endslopes
    d31 = x(3) - x(1);
    d21 = x(2) - x(1);
    d32 = x(3) - x(2);
    for j = 1:pglen
        dvdf1 = (y(pglen + j) - y(j))/d21;
        dvdf2 = (y(2*pglen + j) - y(pglen + j))/d32;
        c(j) = (dvdf2 - dvdf1)/d31;
        c(pglen + j) = dvdf1 - c(j)*d21;
        c(2*pglen + j) = y(j);
    end
    pp = makepp([x(1),x(3)],c);
end

%--------------------------------------------------------------------------
