function v = pchip(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.
%   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.

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

eml_invariant(nargin >= 2,'MATLAB:minrhs');
chckxy(x,y);
eml_invariant(nargin < 3 || isreal(xx), ...
    'MATLAB:pchip:ComplexInterpPts');
ndy = coder.internal.indexInt(eml_ndims(y));
nx = coder.internal.indexInt(numel(x));
nxm1 = nx - 1;
ZERO = coder.internal.indexInt(0);
% h = diff(x);
h = coder.nullcopy(zeros(1,nxm1,'like',x));
for k = 1:nxm1
    h(k) = x(k + 1) - x(k);
end
% Compute slopes
if eml_is_const(iscolumn(y)) && iscolumn(y)
    % Handle the column vector case.
    m = coder.internal.indexInt(1);
    szdel = coder.internal.indexInt([size(y,2),size(y,1)-1]);
else
    m = coder.internal.prodsize(y,'below',ndy); % y is psychologically m-by-nx.
    szdel = coder.internal.indexInt(size(y));
    szdel(ndy) = szdel(ndy) - 1;
end
type = coder.internal.scalarEg(x,y);
del = coder.nullcopy(zeros(szdel,'like',type));
% del = diff(y,1,2)./repmat(h,m,1);
for k = 1:nxm1
    c1 = (k - 1)*m;
    c2 = k*m;
    for j = 1:m
        del(c1 + j) = (y(c2 + j) - y(c1 + j))/h(k);
    end
end
slopes = coder.nullcopy(zeros(size(y),'like',type));
% Compute pchipslopes for each row, and separately for real and imaginary
% parts of del if del is not real.
if nx == 2
    % Special case n = 2, use linear interpolation.
    for k = 1:nx
        c1 = (k - 1)*m;
        for j = 1:m
            slopes(c1 + j) = del(j);
        end
    end
else
    nxm2 = nx - 2;
    nxm3 = nx - 3;
    % Compute interior point slopes.
    for k = 1:nxm2
        hs = h(k) + h(k + 1);
        hs3 = 3*hs;
        w1 = (h(k) + hs)/hs3;
        w2 = (h(k + 1) + hs)/hs3;
        c1 = (k - 1)*m;
        c2 = k*m;
        for j = 1:m
            slopes(c2 + j) = interiorSlope(del(c1 + j),del(c2 + j),w1,w2);
        end
    end
    % Compute end point slopes via non-centered, shape-preserving
    % three-point formulae.
    c0 = nxm1*m;
    c1 = nxm2*m;
    c2 = nxm3*m;
    for j = 1:m
        slopes(j) = exteriorSlope(del(j),del(j + m),h(1),h(2));
        slopes(j + c0) = exteriorSlope(del(j + c1),del(j + c2),h(nxm1),h(nxm2));
    end
end
% Compute piecewise cubic Hermite interpolant to those values and slopes
pp = pwchcore(x,y,ZERO,slopes,h,del);
if nargin == 3
    v = ppval(pp,xx);
else
    v = pp;
end

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

function s = exteriorSlope(d1,d2,h1,h2)
if isreal(d1) && isreal(d2)
    s = realExteriorSlope(d1,d2,h1,h2);
else
    s = complex( ...
        realExteriorSlope(real(d1),real(d2),h1,h2), ...
        realExteriorSlope(imag(d1),imag(d2),h1,h2));
end

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

function s = realExteriorSlope(d1,d2,h1,h2)
s = ((2*h1 + h2)*d1 - h1*d2)/(h1 + h2);
signd1 = sign(d1);
signs = sign(s);
if signs ~= signd1
    s(1) = 0;
elseif (signd1 ~= sign(d2)) && (abs(s) > abs(3*d1))
    s(1) = 3*d1;
end

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

function s = interiorSlope(d1,d2,w1,w2)
if isreal(d1) && isreal(d2)
    s = realInteriorSlope(d1,d2,w1,w2);
else
    s = complex( ...
        realInteriorSlope(real(d1),real(d2),w1,w2), ...
        realInteriorSlope(imag(d1),imag(d2),w1,w2));
end

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

function s = realInteriorSlope(d1,d2,w1,w2)
s = zeros(class(coder.internal.scalarEg(d1,d2,w1,w2)));
% If d1 and d2 have the same sign,
% s = dmin/conj(w1*(d1/dmax) + w2*(d2/dmax)), otherwise s = 0.
if d1 < 0
    if d2 <= d1
        % dmin = -d1, dmax = -d2
        s = d1/conj(w1*(d1/d2) + w2);
    elseif d2 < 0
        % dmin = -d2, dmax = -d1
        s = d2/conj(w1 + w2*(d2/d1));
    end
elseif d1 > 0
    if d2 >= d1
        % dmin = d1, dmax = d2
        s = d1/conj(w1*(d1/d2) + w2);
    elseif d2 > 0
        % dmin = d2, dmax = d1
        s = d2/conj(w1 + w2*(d2/d1));
    end
end

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