algorithm1.m

function [ vortices ] = algorithm1( wavefunction)
% Returns the position and charge of vortices found in a 2D wavefunction.
%   Uses Plaquette algorithm (1) to detect vortices in a wavefuntion.
%   The format of the result is rows of:
%   [vortex charge, x-position, y-position] for each vortex.
% Note that the position given is the top left point of the square of
% four points that bounds the vortex.

[x_size, y_size] = size(wavefunction);
centre_x = x_size/2;
centre_y = y_size/2;

vortices_found = 0;
vortices = [];

for x = 2:(x_size)
for y = 1:(y_size-1)
    charge = 0;
    % Save the arg of each point in an anticlockwise direction.
    arg(1) = angle(wavefunction(x, y));
    arg(5) = arg(1);
    arg(2) = angle(wavefunction(x-1, y));
    arg(3) = angle(wavefunction(x-1, y+1));
    arg(4) = angle(wavefunction(x, y+1));

    % Counts the jumps needed to 'phase unwrap'.
    for ind = 1:4
        if arg(ind+1) - arg(ind) <= -pi
            charge = charge+1;
        elseif arg(ind+1) - arg(ind) >= pi
            charge = charge-1;
        end
    end

    % Saves the info in an array if there is a vortex detected 
    if charge ~= 0 
       vortices_found = vortices_found+1;
       vortices(vortices_found, 1) = charge;
       vortices(vortices_found, 2) = y;
        vortices(vortices_found, 3) = x; 
    end
    end
end

end

algorithm2.m

function [ vortices_found, vort_info ] = algorithm2( wavefunction , h)
% Returns the position and charge of vortices found in a 2D wavefunction. 
%   Uses algorithm (2) to detect vortices in a wavefuntion. 
%   The format of the result is number of vortices then rows of:
%   {vortex charge, x-position, y-position} for each vortex.
% % Note that the position given is the top left point of the square of
% % four points that bounds the vortex.

[x_size, y_size] = size(wavefunction);
centre_x = x_size/2;
centre_y = y_size/2;

vortices_found = 0;
vortices = [];

arg = angle(wavefunction);

% Note that dimensions are swapped in Matlab.

argx = unwrap(arg, [], 1);
vx = diff(argx, 1, 2)./h;

argy = unwrap(arg, [], 2);
vy = diff(argy, 1, 1)./h;

omega = (diff(vy, 1, 2) - diff(vx, 1, 1))./h;
[Y, X] = find(abs(omega) > 1);
vortices_found = length(X);
charges = sum(sign(omega(Y, X)));
vort_info = [charges(:), X, Y+1];
end