%simulation for Tom Seeley
%
%--initial burst of scouts returns slowly
%--bee becomes scout randomly at low rate
%--bee becomse scont at higher rate if within sight of dancer
%--bee dance decriments in time
%Quality max is 300 seconds
%Sswarm size is 4000


%wipe the slate
clear all;
%make sure each run is random
rand('state',sum(100*clock));

%====================================================
%show swarm?
MakeGraphics = 0;
%====================================================
%make a movie
MakeMovie = 0; %change to zero to supress the movie
MovieFolder = 'C:\Documents and Settings\bruce land\My Documents\Matlab\Seeley\';
MovieName = 'test1';
if MakeMovie==1
    mov = avifile([MovieFolder,MovieName,'.avi'],...
        'fps',10,'compression','cinepak');
end
%====================================================
%define time and space
dt = 1/10; %minute
endTime = 30*60; %default end time in minutes
timeMax = endTime %modified if terminated early
t=0:dt:endTime;

%define the swarm space
xSwarm = .3; %meter
ySwarm = .3; %meter

for run=1:100

    %====================================================
    nBee = 200;
    %Initialize Bee state:
    for i=1:nBee
        %alive?
        bee(i).IsAlive = 1;
        %Scout?
        bee(i).IsScout = 0;
        %is the scout out of the nest
        bee(i).IsScoutOut = 0;
        %is the scout out of the nest
        bee(i).IsResting = 0;
        %the next time at which the bee will have an event to handle
        bee(i).NextEventTime = inf;
        %the event to handle-- see list below
        bee(i).NextEvent = 0;
        %position if in swarm and not dancing
        bee(i).xPosition = rand*xSwarm;
        bee(i).yPosition = rand*ySwarm;
        %Probability of becoming scout (per minute) spontaneously
        %units are about 1 scout goes out every 1 hours (60 mins)
        %1/1000 bees * 1/60 min
        bee(i).ProbBecomeScoutSpont = 1/200 * 1/60 ; %per minute
        %Probability of becoming scout (per minute) with dance
        %high qual site dance 3.3+2.8 +...=12.6 min
        %want gain of >1 \
        %For quality=0.37; 0.015 will not sustain scouts, but 0.025 will
        bee(i).ProbBecomeScoutDance = 0.1; %per minute
        %distance to nearest scout
        bee(i).DistToScout = inf;
        %Dance duration is current strength of waggle dance
        bee(i).DanceDuration = 0; %units minutes of dance
        %DanceDecriment is change in dance duration per return
        bee(i).DanceDecrement = 0.5; %units are min of dance/return
        %range of effect of dance
        bee(i).DanceRange = .02; %meters
        %time after dance when bee does nothing
        bee(i).RestDuration = 30; %min
        %Flying speed
        %bee(i).FlyRate = 250; %meters/minute
        %Nest destination
        bee(i).NestDestination = 0; %integer
        %time spent at a potential nest
        bee(i).NestEvalTime = max(20,40+5*randn) ; %min mean=40, sd=5
    end %bee init

    %====================================================
    %possible events list
    ReturnToSwarm = 1;
    ScoutAgain = 2;
    EndDance = 3;
    Rest = 4;

    %====================================================
    %initialize nest state:
    nNest = 4;
    QualityMax = 300/60; %minutes
    %the time of the first dance for a very good nest site (200 sec)

    nest(1).distance = 100; %meters
    nest(1).direction = 90; %degrees
    nest(1).quality = 1.0 ; % varies from 0 to 1
    nest(1).DanceColor = 'red'; %color of the bee when dancing

    %set to 0.37 quality to match paper
    nest(2).distance = 100; %meters
    nest(2).direction = 270; %degrees
    nest(2).quality = 0.75; % varies from 0 to 1
    nest(2).DanceColor = 'blue';

    nest(3).distance = 100; %meters
    nest(3).direction = 180; %degrees
    nest(3).quality = 0.75; % varies from 0 to 1
    nest(3).DanceColor = 'green';

    nest(4).distance = 100; %meters
    nest(4).direction = 0; %degrees
    nest(4).quality = 0.75; % varies from 0 to 1
    nest(4).DanceColor = 'cyan';

    %====================================================
    %repeats fo the same parameters


    %====================================================
    %variables to capture for later analysis
    nScouts(1:length(t),1:nNest) = 0;
    nDancers(1:length(t),1:nNest) = 0;
    nAtNest(1:length(t),1:nNest) = 0;
    %counter is for output data
    count = 1;

    %====================================================
    %construct the graphics objects
    %filled circle for each bee -- default color yellow
    %square region for swarm
    if MakeGraphics
        figure(1)
        clf; set(gcf,'doublebuffer','on')
        beeSize=10;
        beeColor = [.5,.5,.1];
        ax = axes('position',[.05 .05 .9 .9]);
        set(ax,'xlim',[-.05 xSwarm+.05], 'ylim', [-.05 ySwarm+.05]);
        axis equal
        axis tight
        axis off
        hold on
        rectangle('position',[0 0 xSwarm ySwarm],'facecolor',[1 1 1]*.7);

        for i=1:nBee
            bee(i).IconHandle = plot(bee(i).xPosition, bee(i).yPosition, ...
                'markerfacecolor', beeColor, ...
                'markeredgecolor', 'k', ...
                'marker', 'o', ...
                'markersize',beeSize);
        end
    end



    %====================================================
    %Initial scout burst at t=0
    nInitialScout = 6 ;
    time = 0;
    for i=1:nInitialScout
        bee(i).IsScout = 1;
        bee(i).IsScoutOut =1 ;
        bee(i).NestDestination = randint(1,1,nNest)+1;
        %bees arrivals are separated in time by 30 mins
        bee(i).NextEventTime = time + bee(i).NestEvalTime + (i-1)*30 ;
        bee(i).NextEvent = ReturnToSwarm;
        bee(i).DanceDuration = QualityMax * nest(bee(i).NestDestination).quality;
    end

    %====================================================

    %Start time loop
    for time=t

        if ~MakeGraphics & (time-fix(time)<.01)
            time
        end
        %find the current scouts
        scoutList = find([bee.IsScout] & ~[bee.IsScoutOut] & ~[bee.IsResting]);

        %for each bee
        for i=1:nBee

            %====================
            %check to see if bee spontaneously converts to a scout
            if (bee(i).IsScout==0) & ...
                    (rand < bee(i).ProbBecomeScoutSpont*dt)
                bee(i).IsScout = 1;
                bee(i).IsScoutOut =1 ;
                bee(i).NestDestination = randint(1,1,nNest)+1;
                bee(i).NextEventTime = time + bee(i).NestEvalTime;
                bee(i).NextEvent = ReturnToSwarm;
                bee(i).DanceDuration = QualityMax * nest(bee(i).NestDestination).quality;
                if MakeGraphics
                    set(bee(i).IconHandle, 'visible', 'off')
                end
            end %bee spont converts to scout

            %====================
            %check to see if bee sees dance and converts to a scout
            %get distance to each scout
            if (bee(i).IsScout==0) & ...
                    (~isempty(scoutList)) %is at leat one scout in the swarm?
                [bee(i).DistToScout,ss] =  min(max ...
                    (abs(bee(i).xPosition-[bee(scoutList).xPosition]), ...
                    abs(bee(i).yPosition-[bee(scoutList).yPosition])));
                %now do the check for distance and roll the dice
                if (rand < bee(i).ProbBecomeScoutDance*dt) & ...
                        (bee(i).DistToScout <= bee(i).DanceRange)
                    bee(i).NestDestination=bee(scoutList(ss)).NestDestination ;
                    bee(i).DistToScout = inf; %reset distance
                    bee(i).IsScout = 1;
                    bee(i).IsScoutOut =1 ;
                    bee(i).NextEventTime = time + bee(i).NestEvalTime;
                    bee(i).NextEvent = ReturnToSwarm;
                    bee(i).DanceDuration = QualityMax * nest(bee(i).NestDestination).quality;
                    if MakeGraphics
                        set(bee(i).IconHandle, 'visible', 'off')
                    end
                end %bee converts to scout
            end %check to make sure there is at least one scout
            %====================
            %--is there a pending event?
            if (bee(i).IsScout==1) & ...
                    (time >= bee(i).NextEventTime)
                bee(i).NextEventTime = inf; %reset event waiting time

                %===event is Return=====
                if (bee(i).NextEvent==ReturnToSwarm)
                    bee(i).IsScoutOut =0 ;
                    bee(i).xPosition = rand*xSwarm;
                    bee(i).yPosition = rand*ySwarm;
                    if MakeGraphics
                        set(bee(i).IconHandle, 'visible', 'on',...
                            'marker','o', ...
                            'markerfacecolor',nest(bee(i).NestDestination).DanceColor, ...
                            'xdata', bee(i).xPosition, 'ydata', bee(i).yPosition);
                    end
                    %see if dance is done
                    if bee(i).DanceDuration<=0
                        bee(i).NextEvent = EndDance;
                        bee(i).NextEventTime = time ;
                    else
                        bee(i).NextEventTime = time + bee(i).DanceDuration;
                        bee(i).NextEvent = Rest;
                    end
                    %update dance time for the next dance
                    bee(i).DanceDuration = bee(i).DanceDuration - ....
                        bee(i).DanceDecrement;
                    %end return event

                    %===event is rest=====
                elseif (bee(i).NextEvent==Rest)
                    bee(i).IsScoutOut = 0 ;
                    bee(i).IsResting = 1;
                    bee(i).NextEventTime = time + bee(i).RestDuration;
                    bee(i).NextEvent = ScoutAgain;
                    if MakeGraphics
                        set(bee(i).IconHandle, 'visible', 'on',...
                            'markerfacecolor',nest(bee(i).NestDestination).DanceColor, ...
                            'marker','v');
                    end
                    %end rest event

                    %===event is go back to nest site again=====
                elseif (bee(i).NextEvent==ScoutAgain)
                    bee(i).IsScoutOut = 1 ;
                    bee(i).IsResting = 0;
                    bee(i).NextEventTime = time + bee(i).NestEvalTime;
                    bee(i).NextEvent = ReturnToSwarm;
                    if MakeGraphics
                        set(bee(i).IconHandle, 'visible', 'off')
                    end
                    %end scout again event

                    %===event is dance end=====
                elseif (bee(i).NextEvent==EndDance)
                    bee(i).xPosition = rand*xSwarm;
                    bee(i).yPosition = rand*ySwarm;
                    if MakeGraphics
                        set(bee(i).IconHandle, 'visible', 'on',...
                            'markerfacecolor', beeColor, ...
                            'marker','o', ...
                            'xdata', bee(i).xPosition, 'ydata', bee(i).yPosition);
                    end
                    bee(i).NextEvent = 0;
                    bee(i).IsScout=0;
                    %dance end event

                end %event selection ELSEIF sequence

            end %pending event check
            %====================
        end %bee update

        %update the statistics
        scoutList = find([bee.IsScout]);
        for jj=1:nNest
            nScouts(count,jj) = ...
                length(find([bee(scoutList).NestDestination]==jj));
            nDancers(count,jj) = ...
                length(find([bee(scoutList).NestDestination]==jj &...
                ([bee(scoutList).IsScoutOut]==0) & ...
                ([bee(scoutList).IsResting]==0) ));
            nResting(count,jj) = ...
                length(find([bee(scoutList).NestDestination]==jj &...
                [bee(scoutList).IsResting]==1 ));
            nAtNest(count,jj) = length(find([bee(scoutList).NestDestination]==jj &...
                [bee(scoutList).IsScoutOut]==1 ));
        end
        if max(nAtNest(count,:))>100
            timeMax = time;
            break %get out of the loop
        end
        count = count +1;

        %manage the figure window
        if MakeGraphics
            title(['Time=',num2str(time,'%6.2f')]);
            drawnow
        end
        if MakeMovie==1
            F = getframe(gcf);
            mov = addframe(mov,F);
        end
    end %time loop
    %===end the time loop===============================================

    %Firm up the movie file
    if MakeMovie==1
        mov = close(mov);
    end

    %===save some data=========================================
    %save all of the arrays generated by the sim
    
    figfile = ['S',num2str(nBee*bee(1).ProbBecomeScoutSpont),...
        'D', num2str(bee(1).ProbBecomeScoutDance), ...
        'N',num2str(nInitialScout), ...
        'R',num2str(run)]
    figdir = ...
    'C:\Documents and Settings\Bruce\My Documents\m-files\Seeley\Version14\R1\';
    save([figdir,figfile,'.mat'] ,'nScouts','nDancers','nAtNest');
end   
