function [RandDots, datafilename] = randomdot_simfunc(hardcore_distance, total_fly_number, number_of_arenas, arena_diameter, datapath)
%
% to use:
%     hardcore_distance          = 10;  % mm 
%     total_fly_number           = 50; 
%     number_of_arenas           = 20; 
%     arena_diameter             = 90;  % mm
%     datapath                   = 'D:\';
%     [RandomDots, datafilename] = randomdot_simfunc(hardcore_distance, total_fly_number, number_of_arenas, arena_diameter, datapath);
%
% YZ 201902

datafilename = [datapath 'Randomdots' ...
    '_d' num2str(hardcore_distance),  ...
    '_X' num2str(number_of_arenas),   ...
    '_N' num2str(total_fly_number),   ...
    '_D' num2str(arena_diameter)];

fly_halfwidth = 0.5;  %mm.  half-width of the fly, to control fly from touching the arena
center_point  = [arena_diameter/2,arena_diameter/2]; %(x,y)

flies   = [];
notdone =  1; 
goodrun =  0; 
cc      =  0; % seed for random number generator
    
while notdone
    fly_num      = total_fly_number * 50000; % 1000 for hard-core distance 0-9, 50000 for 10
    rng(cc);
    initflies    = rand(fly_num,2)*arena_diameter; %(x,y). a series of flies with random positions
    insidepoints = find(hypot(initflies(:,1)-center_point(1), initflies(:,2)-center_point(2)) < arena_diameter/2-fly_halfwidth);

    flyset = initflies(insidepoints,:); %keep these inside of the arena
    okfly  = 1; 
    for ffii = 2:length(flyset)
        if (hardcore_distance == 0)   % Poisson Point Process
            okfly = [okfly ffii];     % to keep everything
        else
            M  = pdist2(flyset(okfly,:),flyset(ffii  ,:)) < hardcore_distance;
            if ~sum(M) okfly = [okfly ffii  ]; end %keep it (into okfly), if not too close to existing flies 
        end
        if length(okfly) >= total_fly_number, break; end
    end
    cc = cc + 1;
    if mod(cc,2)==0 fprintf('|%d',cc); end
    if length(okfly)<total_fly_number
        notdone = 1;
    else
        goodrun = goodrun + 1;
        pickfly = okfly(1:total_fly_number);
        flies(goodrun,:,:) = flyset(pickfly,:); %data structure: ArenaID x flyID x coordinates 
        if (goodrun == number_of_arenas) notdone = 0; end
    end
end

nop = 200;
r   = arena_diameter/2;
da  = linspace(0,2*pi, nop);
[cxs, cys] = pol2cart([da da(1)], r);
cxs = cxs + arena_diameter/2;
cys = cys + arena_diameter/2;

[circle.xs, circle.ys,circle.r] = deal(cxs, cys, r);

%these data will also be saved to the data file (distance unit -- mm):
RandDots.style                       = 'YZ 2018';
RandDots.flies                       = flies;
RandDots.fly_num                     = total_fly_number;
RandDots.fly_num_started             = fly_num;
RandDots.number_of_arenas            = number_of_arenas;
RandDots.fly_halfwidth               = fly_halfwidth;
RandDots.hardcore_distance           = hardcore_distance;
RandDots.center_point                = center_point;
RandDots.circle                      = circle;  %(x,y)
RandDots.arena_diameter              = arena_diameter;

save(datafilename, 'RandDots'); 



