%
%  Monte-Carlo simulation of chemotactic insect locomotion
%  Ramdom chemotactic movement towards a patch located at (xPatch, yPatch)
%
%  Victor Matveev, June 22, 2010
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mapSize = 600;
figure('Position',[1 200 mapSize mapSize]);

xPatch = 100;  % Chemical source patch center, x coordinate 
yPatch = 100;  % Chemical source patch center, y coordinate 
r0 = 50;       % Patch size (radius of highest concentration)
R = 200;       % Patch maximal spread

chem = zeros(mapSize);  % Array storing chemical concentration
                        
for x = 1:mapSize       % Calculate chemical concentration at each map pixel
    for y = 1:mapSize
        r = sqrt((x-xPatch)^2+(y-yPatch)^2);  % Distance from patch center
        if r < r0
            chem(x,y) = 1;  % Chemical concenration = 1 when r < r0
        elseif r < R
            chem(x,y) = log(r/R)/log(r0/R);  % Log dependence on distance
        else
            chem(x,y) = 0;  % Chemical concenration = 0 when r > R
        end
    end
end
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N = 100;  % Number of insects

ang = rand(N,1)*2*pi;   % Initial polar angles
rad = rand(N,1)*100;    % Initial distance from the enclosure center

xAnt = mapSize/2 + rad .* cos(ang);  % Calculate cartesian positions
yAnt = mapSize/2 + rad .* sin(ang);   

dirAnt = rand(N,1)*2*pi; % Initial insect direction

for i = 1 : N
    chemOld(i) = chem(floor(xAnt(i))+1, floor(yAnt(i))+1);
end;


for j = 1:4000
 
    xAnt = xAnt + cos(dirAnt) .* rand(N,1) * 4; 
    yAnt = yAnt + sin(dirAnt) .* rand(N,1) * 4;

    xAnt(xAnt > mapSize) = 1;      % Don't let anyone escape!
    xAnt(xAnt < 0) = mapSize - 1;  % Don't let anyone escape!
    yAnt(yAnt > mapSize) = 1;      % Don't let anyone escape!
    yAnt(yAnt < 0) = mapSize - 1;  % Don't let anyone escape!
       
    for i = 1 : N
        chemNew(i) = chem(floor(xAnt(i))+1, floor(yAnt(i))+1);
    
        if chemNew(i) >= chemOld(i)    
            dirAnt(i) = dirAnt(i) + randn * 0.25;  % Happy ants stick to old direction
        else
            dirAnt(i) = rand * 2 * pi; % Unhappy ants try something new
        end;
    end;

    chemOld = chemNew;
    
    plot(xAnt, yAnt, '.');
    title(sprintf('Chemotactic movement towards a patch at (%d,%d)', xPatch, yPatch));
    axis([0 mapSize 0 mapSize]);
    drawnow;
    %pause(0.1);
end
