function y=pccp1(x0,y0,x1,y1,x2,y2,x3,y3,xi)
% calculates interpolated value yi at location xi (x1.le.xi.lex2)
% by piecewise continuous cubic polynomial method.  References:
% Akima, H J, J. Assoc. Comp. Mach., v17, 589-602, 1970.
% Wiggins, R A, BSSA, v66, 2077-2081, 1976.
%
eps=0.0001;
% determine slopes
sx1=(y1-y0)/(x1-x0);
sx2=(y2-y1)/(x2-x1);
sx3=(y3-y2)/(x3-x2);
% weight slopes
w1=1/max([abs(sx1) eps]);
w2=1/max([abs(sx2) eps]);
w3=1/max([abs(sx3) eps]);
s1=(w1*sx1 + w2*sx2)/(w1 + w2);
s2=(w2*sx2 + w3*sx3)/(w2 + w3);
% evaluate polynomial at x=xi
p0=y1;
p1=s1;
p2=(3*(y2-y1)/(x2-x1) - 2*s1 - s2)/(x2-x1);
p3=(s1 + s2 - 2*(y2-y1)/(x2-x1))/(x2-x1)^2;
xmd=xi-x1;
y=p0 + p1*xmd + p2*xmd^2 + p3*xmd^3;