function [D] = extrapcoefs(N, x)
% function [D] = extrapcoefs(N, x)
% Creates the coefficients used to extrapolate a value from a polynomial
% function. 
%
% INPUTS:
%       N:      The degree of the polynomial
%       x:      (optional) the points where the polynomial is evaluated.
%           x(1): x value where extrapolation is desired
%           x(2 : N + 1): x values of existing data
%         Default value x = [-1, 0, 1, ..., N-1]
%
% OUTPUTS:
%       D:      Extrapolation coefficients, such that:
%           F(x(1)) ~ D * F(x(2 : N + 1)) 
%
% Examples: 
% To linearly extrapolate the known pairs (x,y) = (0, 0), (1,2) and find
% the y value of (-1,?) use:
% >>> D = extrapcoefs(1)
% >>> y = D * [0;2]
%
% To linearly extrapolate the known pairs (x,y) = (0, 0), (1,2) and find
% the y value of (3,?) use:
% >>> D = extrapcoefs(1,[3 0 1])
% >>> y = D * [0;2]
%
% To extrapolate from the known points (x,y) = (0,1), (1,2), (2,1) (3,2) to
% the point (-1,?) use:
% >>> D = extrapcoefs(3)
% >>> y = D * [1;2;1;2]
% >>> plot([-1:3],[y,1,2,1,2],'o')
% >>> hold on
% >>> plot([-1:0.01:3],polyval(polyfit(0:3, [1 2 1 2], 3),[-1:0.01:3]),'r')
% Notice the poor answer: CAREFUL WITH EXTRAPOLATION!
%
% Written by: Matt Ueckermann

%Creating the coefficients is pretty simple:
% Eventually we want to know the y-value at x*, which can be evaluated
% using:
% [x*^0 x*^1 x*^2 ... x*^n] * [a1;a2;...;an] = y*
%
% But for that we need to solve for coefficients of polynomial:
% [x1^0 x1^1 x1^2 ... x1^n;
%  x2^0 x2^1 x2^2 ... x2^n;
%  ...   ...  ...  ... ...;
%  xn^0 xn^1 xn^2 ... xn^n;] * [a1;a2;...;an] = [y1;y2;...;yn]
%
% A = X^-1 * Y
% (X*) * A = (X*) * X^-1 * Y;
%          = D * Y; 

%Create the Vandermonde matrix:
if nargin < 2
    x = -1:N;
end
V = fliplr(vander(x(2 : N + 2)));

x = fliplr(vander(repmat(x(1), 1, N + 1)));
D = x(1,:) / V;


