function [v, vp, vpp] = pqhip(x,yin,u)
%PQHIP  piecewise quintic Hermite interpolation.
%  [v, vp, vpp] = pqhip(x,y,u) finds the piecewise quintic
%  interpolant P(x), with P(x(j)) = y(1,j), P'(x(j)) = y(2,j)
%  and returns v(k) = P(u(k)), vp(k) = P'(u(k)), and vpp(k) = P''(u(k)).
%
%x = sort(x);  nodes  must be sorted 
%  RMC December 2010 modified from pchiptx
 
% -------------------------------------------------------

function d = pchipslopes(h,delta)
%  PCHIPSLOPES  Slopes for shape-preserving Hermite cubic
%  pchipslopes(h,delta) computes d(k) = P'(x(k)).

%  Slopes at interior points
%  delta = diff(y)./diff(x).
%  d(k) = 0 if delta(k-1) and delta(k) have opposite
%         signs or either is zero.
%  d(k) = weighted harmonic mean of delta(k-1) and
%         delta(k) if they have the same sign.

   n = length(h)+1;
   d = zeros(size(h));
   k = find(sign(delta(1:n-2)).*sign(delta(2:n-1))>0)+1;
   w1 = 2*h(k)+h(k-1);
   w2 = h(k)+2*h(k-1);
   d(k) = (w1+w2)./(w1./delta(k-1) + w2./delta(k));

%  Slopes at endpoints

   d(1) = pchipend(h(1),h(2),delta(1),delta(2));
   d(n) = pchipend(h(n-1),h(n-2),delta(n-1),delta(n-2));
end
% -------------------------------------------------------

function d = pchipend(h1,h2,del1,del2)
%  Noncentered, shape-preserving, three-point formula.
   d = ((2*h1+h2)*del1 - h1*del2)/(h1+h2);
   if sign(d) ~= sign(del1)
      d = 0;
   elseif (sign(del1)~=sign(del2))&(abs(d)>abs(3*del1))
      d = 3*del1;
   end
end

   %x = sort(x);  nodes  must be sorted 
   h = diff(x);
   delta = diff(yin(2,:))./h;
   y(1:2,:) = yin(1:2,:);
   y(3,:)   = pchipslopes(h,delta);



%  Piecewise polynomial coefficients are just the data themselves,
%  because we use the Hermite interpolational basis.

   n = length(x);

%  Find subinterval indices k so that x(k) <= u < x(k+1)

   k = ones(size(u));
   for j = 2:n-1
      k(x(j) <= u) = j;
   end

%  Evaluate interpolant, using formulas generated by Maple
%  and modified by hand: * --> .*, ^ --> .^, and / --> ./

   s = u - x(k);
   h = x(k+1)-x(k);
   if any( h <= 0),
       error 'This routine requires strictly increasing real nodes.';
   end
   theta = s./h;

%   v = ((-s + h) .^ 3 .* (6 * s .^ 2 + h .^ 2 + 3 * s .* h) ./ h .^ 5 .* y(1, k)) ...
%       + (s .* (-s + h) .^ 3 .* (h + 3 * s) ./ h .^ 4 .* y(2, k)) ...
%       + (s .^ 2 .* (-s + h) .^ 3 ./ h .^ 3 .* y(3, k)) / 0.2e1 ...
%       + (s .^ 3 .* (6 * s .^ 2 - 15 * s .* h + 10 * h .^ 2) ./ h .^ 5 .* y(1, k + 1)) ...
%       - (s .^ 3 .* (-s + h) .* (-3 * s + 4 * h) ./ h .^ 4 .* y(2, k + 1)) ...
%       + (s .^ 3 .* (-s + h) .^ 2 / h .^ 3 .* y(3, k + 1)) / 0.2e1;
   v = -((theta - 1) .^ 3 .* (6 * theta .^ 2 + 1 + 3 * theta) .* y(1, k)) ...
       - (theta .* h .* (theta - 1) .^ 3 .* (1 + 3 * theta) .* y(2, k)) ...
       - (theta .^ 2 .* h .^ 2 .* (theta - 1) .^ 3 .* y(3, k)) / 0.2e1 ...
       + (theta .^ 3 .* (6 * theta .^ 2 - 15 * theta + 10) .* y(1, k + 1)) ...
       - (theta .^ 3 .* h .* (theta - 1) .* (3 * theta - 4) .* y(2, k + 1)) ...
       + (theta .^ 3 .* h .^ 2 .* (theta - 1) .^ 2 .* y(3, k + 1)) / 0.2e1;


%   vp = -(30 * (-s + h) .^ 2 .* s .^ 2 ./ h .^ 5 .* y(1, k)) ...
%       + ((-s + h) .^ 2 .* (h + 5 * s) .* (h - 3 * s) ./ h .^ 4 .* y(2, k)) ...
%       + (s .* (-s + h) .^ 2 .* (-5 * s + 2 * h) ./ h .^ 3 .* y(3, k)) / 0.2e1 ...
%       + (30 * (-s + h) .^ 2 .* s .^ 2 ./ h .^ 5 .* y(1, k + 1)) ...
%       - (s .^ 2 .* (6 * h - 5 * s) .* (2 * h - 3 * s) ./ h .^ 4 .* y(2, k + 1)) ...
%       + (s .^ 2 .* (-s + h) .* (-5 * s + 3 * h) ./ h .^ 3 .* y(3, k + 1)) / 0.2e1;
    vp = -(30 * theta .^ 2 .* (theta - 1) .^ 2 .* y(1, k))./h ...
        - ((theta - 1) .^ 2 .* (5 * theta + 1) .* (3 * theta - 1) .* y(2, k)) ...
        - (theta .* h .* (theta - 1) .^ 2 .* (5 * theta - 2) .* y(3, k)) / 0.2e1 ...
        + (30 .* theta .^ 2 .* (theta - 1) .^ 2 .* y(1, k + 1))./h ...
        - (theta .^ 2 .* (3 * theta - 2) .* (5 * theta - 6) .* y(2, k + 1)) ...
        + (theta .^ 2 .* h .* (theta - 1) .* (5 * theta - 3) .* y(3, k + 1)) / 0.2e1;

%   vpp = -60 * (-s + h) .* s .* (-2 * s + h) ./ h .^ 5 .* y(1, k) ...
%       - 12 * (-s + h) .* s .* (-5 * s + 3 * h) ./ h .^ 4 .* y(2, k) ...
%       + (-s + h) .* (10 * s .^ 2 - 8 * s .* h + h .^ 2) ./ h .^ 3 .* y(3, k) ...
%       + 60 * (-s + h) .* s .* (-2 * s + h) ./ h .^ 5 .* y(1, k + 1) ...
%       - 12 * s .* (-s + h) .* (-5 * s + 2 * h) ./ h .^ 4 .* y(2, k + 1) ...
%       + s .* (10 * s .^ 2 - 12 * s .* h + 3 * h .^ 2) ./ h .^ 3 .* y(3, k + 1);

    vpp = 60  * (theta - 1) .* theta .* (2 * theta - 1) .* (y(1,k+1)-y(1, k))./ h .^ 2 ...
        - 12 * (theta - 1) .* theta .* (5 * theta - 3) .* y(2, k)./h ...
        - (theta - 1) .* (10 * theta .^ 2 - 8 * theta + 1) .* y(3, k) ...
        - 12 * theta .* (5 * theta - 2) .* (theta - 1) .* y(2, k + 1)./h ...
        + theta .* (10 * theta .^ 2 - 12 * theta + 3) .* y(3, k + 1);
end