%------------------------------------------------------------
% Name: Laurent series
% Description:
%   Evaluates the Laurent series coefficients via the Cauchy
%   integral, and plots the resulting power series.
% Dependencies:
%   cif.m: Cauchy integral formula evaluator
% Parameters (in order):
%   fhandle: Handle to a function of single complex variable
%   N:       Number of terms to evaluate. (2*N+1 terms will be
%             evaluated, centered at the constant term)
%   z0:      Expansion point
%   r:       Radius of loop
% Returns:
%   pos_coeff:   The positive coefficients + constant term
%   neg_coeff:   Negative coefficients
%------------------------------------------------------------
function [pos_coeff neg_coeff] = laurent(fhandle,N,z0,r)

pos_coeff = zeros(1,N+1);
neg_coeff = zeros(1,N);

% Constant term:
pos_coeff(1) = cif(fhandle,z0,0,z0,r);

for j = 1:N
    pos_coeff(j+1) = cif(fhandle,z0, j,z0,r);
    neg_coeff(j)   = cif(fhandle,z0,-j,z0,r);
end

% View results
%------------------------------------------------------------
L   = 100;
x   = linspace(z0-r,z0+r,L);
y_l = zeros(1,L); % Laurent
y_a = zeros(1,L); % Actual

for k = 1:L
    sum = pos_coeff(1);
    for j = 1:N
        sum = sum + pos_coeff(j+1) * (x(k)-z0)^j;
        sum = sum +   neg_coeff(j) * (x(k)-z0)^(-j);
    end
    y_l(k) = sum;
    y_a(k) = fhandle(x(k));
end

plot(x,y_a,'b');
hold on;
plot(x,y_l,'r');
M = max(y_a);
m = min(y_a);
axis([z0-r z0+r m M]);
grid on;
xlabel('x');
ylabel('y');
title('Original function (blue) vs. Laurent expansion (red)');