(* ::Package:: *)

  function Sx = PositionCVL(a,t,x0,v0)
% a and t are column m-vectors with 0 = t(1) < t(2) <...< t (m).
% a(i) is the acceleration of an object at time t(i)
% x0 and v0 are the  position and velocity of the object time 0. 
% Sx is the pp-representation of a spline that approxiamtes position.
%
% Compute Sa(t), the spline interpolant of (t(1),a(1)),...,(t (m).a(m))

Sa = spline(t,a);
areaVals = IntervalAreas(Sa);
v = cumsum([v0; areaVals]);
Sv = spline(t,v);
areaVals = IntervalAreas(Sv);
x = cumsum([x0; areaVals]);
Sx = spline(t,x);

  function areaVals = IntervalAreas(S)
% S is a spline in pp-representation 
% If the breakpoints are x(1)<x(2)<...<x(n), then areaVals(k) is the
% integral of S from x(k) to x(k+1)
  
[x,c] = unmkpp(S);
n = length(x);
areaVals = zeros(n-1,1);
for k=1:n-1
  % Integrate S from x(k) to x(k+1)..
  % The kth local cubic is 
  %        c(k,1)(x-x(k))^3 + c(k,2)(x-x (k)^2 + c(k,3)(x-x(k) + c(k,4)
  % Antiderivative:
  %    c(k,1)(x-x(k))^4/4 + c(k,2)(x-x(k))^3/3 + c(k,3)(x-x(k))^2/2 + c(k,4)(x-x(k))
  % So
  del = x(k+1)-x(k);
  areaVals(k) = ((((c (k,1)/4)*del + (c (k,2)/3))*del + c (k,3)/2)*del + c(k,4))*del;
end

% -1 for each redundancy in the method.
% -2~4 if the result is incorrect.
