function pp = Schumaker(FF,ff,tt)
% Schumaker shape-preserving quadratic interpolation spline. 
%
% pp = Schumaker(F,f,t) returns the pp-form of the spline.
%
% F, f and t should be vectors of same size:
% t = [t1 ... tn] are the interpolation nodes;
% F = [F(t1) ... F(tn)] are the function values and
% f = [F'(t1) ... F'(tn)] are the derivative values at the extremes.
% f is optional; if an [] is supplied, derivatives are estimated from F
% and t.

% copyright 2001 Leonardo Rezende.
% This program can be freely copied, used and modified, as long as this
% copyright notice is kept in place. Please inform the author of any
% useful modifications at lrezende@stanford.edu.

n = max(size(tt));

% estimates derivatives if there were not supplied

if isempty(ff),
   L = sqrt( (tt(2:n)-tt(1:n-1)).^2 + (FF(2:n)-FF(1:n-1)).^2);
   d = (FF(2:n)-FF(1:n-1))./(tt(2:n)-tt(1:n-1));
   sb =(d(1:n-2).*d(2:n-1)>0).*(L(1:n-2).*d(1:n-2)+L(2:n-1).*d(2:n-1))./(L(1:n-2)+L(2:n-1));
   if size(sb,1) > size(sb,2),
     ff = [(3*d(1)-sb(1))/2; sb; (3*d(n-1)-sb(n-2))/2];
   else
     ff = [(3*d(1)-sb(1))/2  sb  (3*d(n-1)-sb(n-2))/2];
   end
end

SS = []; 
ts = [];

for j = 1:(n-1);
  [C, C1, fake] = Schum([FF(j) FF(j+1)], [ff(j) ff(j+1)], [tt(j) tt(j+1)]);
  SS = [SS; C;C1];
  ts = [ts; tt(j); fake];
end

ts = [ts; tt(n)];

% SSS = zeros(size(SS));
% for j = 1:size(SS,1)
%   SSS(j,:) = SS(j,:)*[  1        0     0;
%                       -2*ts(j)   1     0;
%                       ts(j)^2 -ts(j) 1]'; 
% end


pp = ppmak(ts,SS,1);

return



function [S, S1, tsi] = Schum(F,f,t)

if sum(f)*(t(2)-t(1)) == 2*(F(2)-F(1)),
   tsi = t(2);
else 
   d = (F(2) -F(1))/(t(2)-t(1));
   if (f(1)-d)*(f(2)-d) >= 0,
       tsi = sum(t)/2;
   elseif abs(f(2)-d) < abs(f(1)-d),
       tsi = (t(1) + (t(2)-t(1))*(f(2)-d)/(f(2)-f(1)));
   else 
       tsi = (t(2) + (t(2)-t(1))*(f(1)-d)/(f(2)-f(1)));
   end
end

a = tsi-t(1); b = t(2)-tsi;
fbar = (2*(F(2)-F(1))-(a*f(1)+b*f(2)))/(t(2)-t(1));

S = [ (fbar-f(1))/(2*a) f(1) F(1)];

if b == 0,
   S1 = []; tsi = [];
else
   S1 = [ (f(2)-fbar)/(2*b) fbar S*[a*a; a; 1]];
end;

return

