function stackdot_swarm = swarmplot(inputdot,mindist,circdimension,move_bilateral,projection_dir)
% Function for making swarm plots
% Input:
%     inputdot: row/column vector or N x 2 matrix. If the input is a
%               matrix, it will be projected to the direction defined by the variable
%               `projection_dir` and make the swarm plot for that direction
%     mindist? scalar, the minimal distance between two points when they
%               are not considered as overlapped
%     circdimension: boolean (default false), if set to true, make circular swarm plot;
%               when input is N x 2 matrix, it will calculate the angle of
%               each row vector and make 1D swarm plot but plotted on a
%               circle
%     move_bilateral: boolean (true), if true, the two overlapped points will be
%               move in two direction
%     projection_dir: complex scalar (default 1i), complex number whose
%               phase angle defined the projection direction.

if ~exist('circdimension','var')
    circdimension = [];
end

if ~exist('move_bilateral','var')
    move_bilateral = [];
end

if ~exist('projection_dir','var')
    projection_dir = [];
end

if isempty(projection_dir)
    projection_dir =  exp(1i.*pi/2);
else
    projection_dir = projection_dir./abs(projection_dir);
end

if isempty(move_bilateral)
    move_bilateral = true;
end

if isempty(circdimension)
    circdimension = false;
end

if move_bilateral
    altersign = -1;
else
    altersign = 1;
end

if circdimension
    if sum(size(inputdot)>1) == 1
        inputdotComplex = 1i.*inputdot(:);
    else
        inputdotComplex = inputdot(:,1)+1i.*inputdot(:,2);
        inputdotComplex = angle(inputdotComplex).*1i;
    end
    [sortedDot,reverse_order] = sort(imag(inputdotComplex));
    [~,reverse_order] = sort(reverse_order);
    stackdot_swarm = sortedDot;
else
    if sum(size(inputdot)>1) == 1
        inputdotComplex = 1i.*inputdot(:);
        projection_dir =  exp(1i.*pi/2);
    else
        inputdotComplex = inputdot(:,1)+1i.*inputdot(:,2);
    end
    [sortedDot,reverse_order] = sort(real(inputdotComplex.*conj(projection_dir)));
    [~,reverse_order] = sort(reverse_order);
    stackdot_swarm = sortedDot;
end

movesign = 1i;

for i = 2:size(sortedDot,1)
    movesign = movesign*altersign;
    backtraceidx = find((sortedDot(i) - sortedDot(1:i-1))<mindist,1);
    while any(abs(stackdot_swarm(i) - stackdot_swarm(backtraceidx:i-1))<mindist)
        stackdot_swarm(i) = stackdot_swarm(i) + mindist*movesign;
    end
end

if circdimension
    stackdot_swarm = exp(1i.*real(stackdot_swarm)).*(imag(stackdot_swarm)+1); 
else
    stackdot_swarm = stackdot_swarm.*projection_dir;
end
stackdot_swarm = stackdot_swarm(reverse_order);
end