function yi=msspline1(x0,y0,x1,y1,x2,y2,x3,y3,xi)
% for x1<=xi<=x2, use mean slope cubic spline method 
% (Ku, C C, 1973) to find interpolated value of y.
slope1=(y1-y0)/(x1-x0);
slope2=(y2-y1)/(x2-x1);
slope3=(y3-y2)/(x3-x2);
h=x2-x1;
u=atan(slope1);
v=atan(slope2);
w=atan(slope3);
c0=y1;
c1=tan(.5*(u+v));
d=tan(.5*(v+w));
a=y2-y1-h*c1;
b=d-c1;
c2=(3*a - h*b)/(h^2);
c3=(a - (h^2)*c2)/(h^3);
xv=xi-x1;
yi=c0 + c1*xv + c2*xv^2 + c3*xv^3;