%------------------------------------------------------------
% Name: Cauchy integral formula
% Description:
%   Evaluates the closed loop integral of:
%       1/(2*pi*i) * f(z)/(z-Z)^(j+1)
%   along a circular contour centered at z0 with radius r, 
%   evaluated once in the positive direction.
% Parameters (in order):
%   fhandle: Handle to a function of single complex variable
%   Z:       Appears in the term (z-Z)
%   j:       Power of (z-Z)
%   z0:      Center of loop
%   r:       Radius of loop
% Returns:
%   integ:   The result of integral
%------------------------------------------------------------
function integ = cif(fhandle,Z,j,z0,r)

% Number of segments
%------------------------------------------------------------
N = 50000;

theta   = linspace(0,2*pi,N);
d_theta = theta(2)-theta(1);
zeta    = r*exp(i*theta) + z0;
d_zeta  = r*i*exp(i*theta)*d_theta;

integ   = 0;
for k = 1:N
    term  = fhandle(zeta(k))/(zeta(k)-Z)^(j+1)*d_zeta(k);
    integ = integ + term;
end

integ   = 1/(2*pi*i)*integ;
integ   = clean(integ);

% Truncate numerical trailing
%------------------------------------------------------------
function z_clean = clean(z)
eps = 0.001;
if abs(real(z)) < eps
    z = imag(z)*i;     % Note to self: Recall that imag(z) returns
                       %               a real value!
end
if abs(imag(z)) < eps
    z = real(z);
end
z_clean = z;