% Mussel Movement IBM

% Supporting Online Material for
% Lévy walks evolve through interaction between movement and 
% environmental complexity

% Monique de Jager*, Franz J. Weissing, Peter M.J. Herman, Bart A. Nolet, 
% Johan van de Koppel
% *To whom correspondence should be addressed. 
% Email: M.deJager@nioo.knaw.nl;

% Last updated: May 8th 2011.

% Model description:

% Individual-based model of mussel movement with 
% density-dependent probability of moving. 

% np individuals are placed randomly on a 2-dimensional grid (of size 
% gridsize by gridsize);
% the locations of these individuals are stored in two vectors: X and Y. 

% Parameters needed for the simulation are:
% - the Levy exponent of the movement strategy (mu);
% - the number of iterations (steps);
% - the radius of a mussel (msize);
% - the radius in which local mussel density is calculated (D1);
% - the radius in which long-range mussel density is calculated (D2);
% - the minimal walking distance (xmin);
% - size of the grid (gridsize);
% - number of individuals (np);

% Variables:
% - vector of X-coordinates (X);
% - vector of Y-coordinates (Y);
% - the degree of patterning (pattern); 
% - matrices to calculate distances with (AX, BX, AY, BY);
% - matrix of distances between mussels (D); 
% - vector of local densities (ssd);
% - vector of long-range densities (lsd);
% - vector of probability of moving (pmove);
% - vector of step lengths (L);
% - vector of turning angles (a);
% - vector of goal X and Y coordinates (Xgoal, Ygoal);
% - vectors of new X and Y coordinates (Xnew, Ynew);
% - matrices to calculate collisions with (AXgoal, BXgoal, AYgoal, 
% BYgoal);
% - matrix of relative intersection points of line A-Agoal with point 
% Bgoal(r);
    % r is the intersection point of line A-Agoal and Bgoal,
    % which is used for finding the perpendicular projection (P)
    % of Bgoal on line A-Agoal;
    % r has the following meaning:
        % r = 0 -> P = A;
        % r = 1 -> P = Agoal;
        % r < 0 -> P is on backward extension of A-Agoal;
        % r > 1 -> P is on forward extension of A-Agoal;
        % 0 < r < 1 -> P is interior to A-Agoal;
% - matrix to help calculate r (r1);
% - matrices with perpendicular projection points of Bnew on line A-Agoal 
% (PX, PY);
% - matrix with distances from line A-Agoal to point Bgoal (Dq);
% - matrix with fractions of step lengths that are collision-free 
% (FreeL);
% - vector with minimum fraction of step lengths that is collision-free 
% (Lmin);

% 1. According to a probability of moving (which is
% dependent on mussel density at two spatial scales), mussels decide
% whether they will move or not. 

% 2. For each mussel, a random turning angle is drawn from a uniform 
% distribution, and a random step length is drawn 
% from a Pareto distribution with parameter mu, to resemble a Levy walk 
% with a scaling exponent mu.

% 3. To account for collisions, the fraction of the step length that can 
% be traversed without bumping into others is calculated.

% 4. The mussels move according to their probability of moving, the 
% turning angle, the step length, and the fraction of this step length 
% that can be moved without collisions. 

% 5. When a mussel moves outside the grid, the X and Y coordinates are 
% reset to random values. 

% 6. The locations of the mussels are plotted on the 2-D grid
% steps 1 to 6 are repeated for each iteration. 

clear;

% parameter values:
msize = 1; 
D1 = 1.1;
D2 =  7.5;
% (D1 and D2 have been derived from experimental data)
xmin = 0.1;
gridsize = 30; 
np = 1500;
ToPlot = 1;                 % plots the distribution of mussels when 1
WriteDataToFile = 0;        % writes the data to a specified .csv file when 1
FileName = 'filename.csv';    % specify the file name here

mu = 2;
 
% state variables:
X = rand(1,np).*gridsize; 
Y = rand(1,np).*gridsize;

steps = 0;
pattern = 0;
Dmoved = 0;
while pattern < 2,
    % 1. Accoding to a probability of moving (which is
    % dependent on mussel density at two spatial scales), mussels decide
    % whether they should move or not. 
        % Calculate distances between all individuals:
        [AX BX] = meshgrid(X);
        [AY BY] = meshgrid(Y);
        D = sqrt((AX - BX).^2 + (AY - BY).^2);

        % Calculcate short-range (ssd) and long-range (lsd) densities for 	  
        % each individual
        ssd = (sum((D < D1)')-1)/(pi * D1^2);
        lsd = (sum((D < D2)')-1)/(pi * D2^2);

        pattern = mean(ssd/lsd);

        % Calculate the probability of moving from these short-range and
        % long-range densities. 
        pmove = (0.6 - 1.3.*ssd + 1.1.*lsd) > rand(1,np);
 
        TotalMove = 0;
        while TotalMove < 50,

% 2. For each mussel, a random turning angle is drawn from a 
% uniform distribution, and a random step length is drawn 
% from a Pareto distribution with parameter mu, to resemble a 
% Levy walk with a scaling exponent mu.

          	% Calculating the walking distances for each individual 		 	
            % in both directions X and Y

                L = xmin./(1 - rand(1,np)).^(1/(mu-1));
                TotalMove = TotalMove + mean(L);
                a = rand(1,np)*2*pi;

                % goal X and Y coordinates 
                Xgoal = X + cos(a).*L;
                Ygoal = Y + sin(a).*L;

           	% 3. To account for collisions, the fraction of the step 			
            % length that can be traversed without bumping into others is 		
            % calculated.    

           	% to prevent mussels bumping into each other:
                [AXgoal BXgoal] = meshgrid(Xgoal);
                [AYgoal BYgoal] = meshgrid(Ygoal);

          	% relative intersection points of line A-Agoal with point 		
            % Bgoal:
                r1 = ((BXgoal - AX).*(AXgoal - AX) + (BYgoal - AY).*(AYgoal - AY));
                r = r1./meshgrid(L).^2;

           	% the perpendicular projection of Bgoal (P) can be found:
                PX = AX + r.*(AXgoal - AX);
                PY = AY + r.*(AYgoal - AY);

           	% with PX and PY, we can calculate the distance from the line
           	% AX-AXgoal to the point Bgoal:   
                Dq = sqrt((PX - BXgoal).^2 + (PY - BYgoal).^2);
           	% If Dq is smaller than msize, the mussel is blocking the 		
            % track

           	% Lmin is the fraction of the step length which can be walked
         	% without crossing other individuals
                FreeL= (r <= 0) + (Dq > msize) + (r > 0).*(Dq <= msize).*r;
                Lmin = min((FreeL)');Lmin(Lmin > 1) = 1;

           	% 4. the mussels move according to their probability of 			
            % moving, the turning angle, the step length, and the 			
            % fraction of this step length that can be moved without 			
            % collisions. 
            
           	% new X and Y values; the steps are truncated when other
           	% individuals are on the path. 
                Xnew = X + pmove.*cos(a).*L.*Lmin;
                Ynew = Y + pmove.*sin(a).*L.*Lmin;
                Dmoved = Dmoved + sum(pmove.*L.*Lmin);
                
          	% 5. when a mussel moves outside the grid, the X and Y 			
            % coordinates are reset to random values. 
           	% if X or Y coord is outside the grid, it will be assigned a 		
            % randomly generated new value of X or Y 
                X = (Xnew < gridsize).*(Xnew > 0).*Xnew + (((Xnew < gridsize).*(Xnew > 0)) == 0).* rand(1,np).*gridsize;
                Y = (Ynew < gridsize).*(Ynew > 0).*Y + (((Ynew < gridsize).*(Ynew > 0)) == 0).* rand(1,np).*gridsize;           

           	% 6. plot the population distribution
            if ToPlot == 1,
                plot(X,Y,'o', 'MarkerEdgeColor','b', 'MarkerFaceColor', 'm', 'Markersize',7), axis([0 gridsize 0 gridsize])
               pause(0.1)
            end            
            steps = steps + 1;
        end
end
if WriteDataToFile == 1,
    fid = fopen(FileName,'a');
        fprintf(fid, '%4.4f, %4.4f, %4.4f\r\n', mu, steps, Dmoved);
    fclose(fid);
end