Code covered by the BSD License

# Interpolation Utilities

### Joe Henning (view profile)

22 May 2012 (Updated )

A variety of interpolation utilities

cakima(x, y, xi)
```function [yi, ypi, yppi] = cakima(x, y, xi)

% CAKIMA 1-D piecewise cubic Akima interpolation
%    CAKIMA(X,Y,XI) interpolates to find YI, the values of
%    the underlying function Y at the points in the array XI, using
%    piecewise cubic Akima interpolation.  X and Y must be vectors
%    of length N.
%
%    [YI,YPI,YPPI] = CAKIMA() also returns the interpolated
%    quadratic derivative and linear second derivative of the
%    underlying function Y at points XI.

% Joe Henning - Fall 2011

n = length(x);
u = zeros(1,n+3);
yp = zeros(1,n);

% Calculate Akima coefficients and derivatives
for i = 1:n-1
% Shift i to i+2
u(i+2) = (y(i+1)-y(i))/(x(i+1)-x(i));
end

u(n+2) = 2.0*u(n+1) - u(n);
u(n+3) = 2.0*u(n+2) - u(n+1);
u(2) = 2.0*u(3) - u(4);
u(1) = 2.0*u(2) - u(3);

for i = 1:n
a = abs(u(i+3) - u(i+2));
b = abs(u(i+1) - u(i));
if ((a+b) ~= 0)
yp(i) = (a*u(i+1) + b*u(i+2))/(a+b);
else
yp(i) = (u(i+2) + u(i+1))/2.0;
end
end

for i = 1:length(xi)
% Find the right place in the table by means of a bisection.
klo = 1;
khi = n;
while (khi-klo > 1)
k = fix((khi+klo)/2.0);
if (x(k) > xi(i))
khi = k;
else
klo = k;
end
end

b = x(khi) - x(klo);
if (b == 0.0)
fprintf('??? Bad x input to cakima ==> x values must be distinct\n');
yi(i) = NaN;
ypi(i) = NaN;
yppi(i) = NaN;
continue;
end

% Evaluate Akima polynomial
a = xi(i) - x(klo);
yi(i) = y(klo) + yp(klo)*a + (3.0*u(klo+2) - 2.0*yp(klo) - yp(klo+1))*a*a/b + (yp(klo) + yp(klo+1) - 2.0*u(klo+2))*a*a*a/(b*b);

% Differentiate to find the second-order interpolant
ypi(i) = yp(klo) + (3.0*u(klo+2) - 2.0*yp(klo) - yp(klo+1))*2*a/b + (yp(klo) + yp(klo+1) - 2.0*u(klo+2))*3*a*a/(b*b);

% Differentiate to find the first-order interpolant
yppi(i) = (3.0*u(klo+2) - 2.0*yp(klo) - yp(klo+1))*2/b + (yp(klo) + yp(klo+1) - 2.0*u(klo+2))*6*a/(b*b);
end
```