%  quick discussion of inverse linear interpolation

%  adaptive piecewise linear interpolation
%   here is a version of the book's pwLAdapt that treats points in the plane
%   as one item, rather than as two.

>> type pwladapt

  function xy = pwLAdapt(fname,xyL,xyR,delta,hmin)
%  xy = pwLAdapt(fname,xyL,xyR,delta,hmin)
% Adaptively determines interpolation points for a piecewise linear
% approximation to the function specified.
%
%  fname  is a string that specifies a function of the form  y = f(u).
%  xyL =: [xL;f(xL)] and  xyR  =: [xR;f(xR)] specify the interval [xL .. xR]
%  and also provide the corresponding function values.
%  delta  is a (positive) uppper bound on the error.
%  hmin/2 is a positive lower bound on the interval length.
%
%  xy(1,:) =: x is an vector with the property that
%              xL = x(1) < ... < x(end) = xR,
%  and xy(2,:) is the corresponding vector of values of f. 
%
%  Each subinterval is either <= hmin in length or has the property
%  that at its midpoint m, |f(m) - L(m)| <= delta, where
%  L(x) is the line that connects (xR,fR) with (xL,fL) 
%           
% 
   if (xyR(1)-xyL(1)) <= hmin % Subinterval is acceptable
      xy = [xyL, xyR];

   else
      mid  = (xyL(1)+xyR(1))/2; fmid = feval(fname,mid);

      if (abs(((xyL(2)+xyR(2))/2) - fmid) <= delta ) % Subinterval acceptable
         xy = [xyL, xyR];

      else             % Produce left and right partitions, then synthesize
         xyLeft   = pwLAdapt(fname,xyL,[mid;fmid],delta,hmin);
         xyRight  = pwLAdapt(fname,[mid;fmid],xyR,delta,hmin);
         xy = [xyLeft, xyRight(:,2:end)];
      end
   end

%   Then discussed cubic Hermite interpolation, as cubic interpolation to data
% at four points, as the left two coalesce and the right two coalesce.
% The book's script ShowHermite gives a nice illustration:

>> showhermite

%  Discussed the construction of the cubic Hermite polynomial interpolant using
%  divided difference table and the Newton form, and using the definition
%   f[t_0, ...,t_k] = f^{(k)}(t_0))/k!   in case  t_0 = t_1 = ... t_k .

% Here is the use of the resulting formulas for the construction of the 
% piecewise cubic Hermite interpolant:

>> type pwch

function pc = pwch(x,y,s)
%
% Pre: x = strictly increasing sequence of length n
%      y, s = sequence of the same orientation and length as x
%
% Post: pc = a description, ready for use with ppval, of the piecewise
%      cubic function  q  that satisfies
%          q(x_i) = y_i,  q'(x_i) = s_i,  i=1:n .
 
h = diff(x); dy = diff(y)./h;
dzzh = (dy-s(1:end-1))./h; dzhh = (s(2:end)-dy)./h;
pc = mkpp(x,[(dzhh-dzzh)./h 2*dzzh-dzhh s(1:end-1) y(1:end-1)]);

%% Note how it differs from the book's pwch command, in that it outputs just one
%% item, for direct use in matlab's ppval (rather than four vectors, two which
%% one then has to remember to adjoin also the sequence x of breakpoints
%% when one wants to evaluate the interpolant using the book's evaluator).

%% Finally, started on spline interpolation.
diary off