function [C P] = taylor(XmdX, deriv, pade,out)
% function [C P] = taylor(XmdX, deriv, pade)
% CALCULATES THE COEFFICIENTS FOR TAYLOR EXPANSIONS TO APPROXIMATE
% DERIVATIVES. INCLUDES PADE SCHEME EXPANSIONS
%
% INPUTS
%   XmdX:       Vector containing values of (x-dx)
%   deriv:      Which derivative to approximate: d^deriv/dx^deriv
%   pade:       Optional input. 2D array, each column containing
%               [{where to evaluate derivite}, {which derivative to use}]
%   out:        If true, scheme is output to the screen
% OUTPUTS
%   C:          Coefficients in front of f(x) terms. Ordered from large to
%                   small
%   P:          Pade coefficients in the order inputted
%   STDOUT - the scheme is printed to standard output if out==true
%
% Example usage:
%1. To recover 2nd order central difference for first order derivatives:
% >> C=taylor([-1 1],1)
% 
%  1/(2h^1) *[1 f(1*dx) + -0 f(0*dx) + -1 f(-1*dx)]
%  		 + h^2 * 0.166667 d^3f^(3)/dx^3
% 
% C =
% 
%     0.5000
%          0
%    -0.5000
%
%2. To recover 2nd order forward difference for second order derivatives:
%>>  C=taylor([3 2 1],2)
% 
%  1/(1h^2) *[-1 f(3*dx) + 4 f(2*dx) + -5 f(1*dx) + 2 f(0*dx)]
%  		 + h^2 * -0.916667 d^4f^(4)/dx^4
% 
% C =
% 
%    -1.0000
%     4.0000
%    -5.0000
%     2.0000
%     
%3. To recover 4th order central difference pade scheme for first order derivatives:
% >> [C P] = taylor([-1 1],1,[-1 1;1 1])
% 
%  1/(1.33333h^1) *[1 f(1*dx) + -0 f(0*dx) + -1 f(-1*dx)] + 
% 	(h^0) *-0.25 f^(1)(-1*dx) + (h^0) *-0.25 f^(1)(1*dx)
%  		 + h^4 * -0.00833333 d^5f^(5)/dx^5
% 
% C =
% 
%     0.7500
%          0
%    -0.7500
% 
% 
% P =
% 
%    -0.2500
%    -0.2500
%3. To recover 4th order central difference pade scheme for second order derivatives:
%>> [C P] = taylor([-1 1],2,[-1 2;1 2])
% 
%  1/(0.833333h^2) *[1 f(1*dx) + -2 f(0*dx) + 1 f(-1*dx)] + 
% 	(h^0) *-0.1 f^(2)(-1*dx) + (h^0) *-0.1 f^(2)(1*dx)
%  		 + h^4 * -0.005 d^6f^(6)/dx^6
% 
% C =
% 
%     1.2000
%    -2.4000
%     1.2000
% 
% 
% P =
% 
%    -0.1000
%    -0.1000
%
%Author: Matt Ueckermann for MIT course 2.29

%check inputs
if nargin<2
    error('Not enough inputs')
end
if nargin<3
    pade=[];
end
if nargin<4
    out=false;
end
%%
%Set RHS of system. Derivative to be solved for has value 1, all else is
%zero.
rhs = zeros(length(XmdX)+size(pade,1),1);
rhs(deriv)=1;
%Calculate the coefficients (term) of the Taylor series expansion
for j=1:length(XmdX)
    xmdx=XmdX(j);
    for i=1:length(XmdX)+3+size(pade,1)
        term(i,j)=xmdx^i/factorial(i);
    end
end
%Calculate the coefficients (term) of the Taylor series expansion for
%derivatives (Pade scheme)
for j=length(XmdX)+1:length(XmdX)+size(pade,1)
    xmdx=pade(j-length(XmdX),1);
    for i=0:length(XmdX)+3-pade(j-length(XmdX),2)+size(pade,1)
        %Need to be careful with offsets here
        term(i+pade(j-length(XmdX),2),j)=xmdx^i/factorial(i);
    end
end

%Error checking
minij=min(size(term));
if minij<length(rhs)
    error('Not enough terms included to solve for derivative\n Add more terms to XmdX')
end
coeffs=term(1:minij,1:minij)\rhs;
% find remainder term:
Rem = term*coeffs;
Rem(deriv) = 0;
dR=find(abs(Rem)>10*(eps),1);
%Separate into Pade and non-pade
divide = length(XmdX);
C = coeffs(1:divide);
P = coeffs(divide+1:end);
%Normalize coefficients for output to screen
denom=1/(min(abs(C)));
C=C*denom;
C=[-sum(C);C];
XmdX = [0 XmdX];
[XmdX S] = sort(XmdX,'descend'); %Sort
C=C(S);
%Print coefficients to screen
if out
fprintf('\n 1/(%gh^%g) *[',denom,deriv)
for j=1:length(C)
    fprintf('%g f(%g*dx)',C(j),XmdX(j))
    if j~= length(C)
        fprintf(' + ')
    end
end
end
C=C/denom; %Take out normalization.
if out
fprintf(']')
if ~isempty(pade)
    %Print pade coefficients to screen
    fprintf(' + \n\t')
    for j=1:length(P)
%         fprintf('(h^%g)/(%g) *',pade(j,2)-deriv,denom)
        fprintf('(h^%g) *',pade(j,2)-deriv)
        fprintf('%g f^(%g)(%g*dx)',P(j),pade(j,2),pade(j,1))
        if j~= length(P)
            fprintf(' + ')
        end
    end
end
%Print remainder to screen
fprintf('\n \t\t + h^%g * %g d^%gf^(%g)/dx^%g\n',dR-deriv,Rem(dR),dR,dR,dR)
end
