% Active walker model (Chapter 7, Box B)
% Written by David Sumpter, corrections by Qi Ma

clear

worldsize=60;
distancebetween=26;
distanceup=40;


%Resistence initially 1.
R=ones(worldsize,worldsize);



%Number of walkers
numwalkers1=50;
numwalkers2=50;

%Positions of targets and source
tx=[ones(numwalkers1,1)*(worldsize-distancebetween)/2 ; ones(numwalkers2,1)*((worldsize+distancebetween)/2)];
ty=ones(numwalkers1+numwalkers2,1)*(worldsize+distanceup)/2;
sx=ones(numwalkers1+numwalkers2,1)*worldsize/2;
sy=ones(numwalkers1+numwalkers2,1)*(worldsize-distanceup)/2;

wx=ceil(rand(numwalkers1+numwalkers2,1)*55);
wy=ceil(rand(numwalkers1+numwalkers2,1)*55);
px=ones(numwalkers1+numwalkers2,1);
py=ones(numwalkers1+numwalkers2,1);


k=0.5;
lambda=0.99;
kappa=0.8;
mu=2;

alpha=400;
    
simtime=150000

%Walker parameters
for t=1:simtime
   
    %Pick a random walker
    i=ceil(rand*(numwalkers1+numwalkers2));
    
    %Reflective boundary conditions
    if wx(i)<=2 
        wx(i)=2;
    end
    if wx(i)>=worldsize-1; 
        wx(i)=worldsize-1;
    end
    if wy(i)<=2 
        wy(i)=2;
    end
    if wy(i)>=worldsize-1
        wy(i)=worldsize-1;
    end

    
    %Got to the target so change direction
    if (wx(i)==tx(i)) & (wy(i)==ty(i))
        tempx=tx(i);
        tempy=ty(i);
        tx(i)=sx(i);
        ty(i)=sy(i);
        sx(i)=tempx;
        sy(i)=tempy;
    end
        
    
    %Find direction to target
    pleft=max(wx(i)-tx(i),0)/(abs(wy(i)-ty(i))+abs(wx(i)-tx(i)));
    pright=max(tx(i)-wx(i),0)/(abs(wy(i)-ty(i))+abs(wx(i)-tx(i)));
    pdown=max(wy(i)-ty(i),0)/(abs(wy(i)-ty(i))+abs(wx(i)-tx(i)));
    pup=max(ty(i)-wy(i),0)/(abs(wy(i)-ty(i))+abs(wx(i)-tx(i)));

    %Find the resistence in the directions
    rtot=(R(wx(i)-1,wy(i)))+(R(wx(i)+1,wy(i)))+(R(wx(i),wy(i)+1))+(R(wx(i),wy(i)-1));
    rleft=(R(wx(i)-1,wy(i)));
    rright=(R(wx(i)+1,wy(i)));
    rup=(R(wx(i),wy(i)+1));
    rdown=(R(wx(i),wy(i)-1));
    
    probtot=(pleft*(1+alpha*rleft)^mu+k)+(pright*(1+alpha*rright)^mu+k)+(pdown*(1+alpha*rdown)^mu+k)+(pup*(1+alpha*rup)^mu+k);
    probleft=(pleft*(1+alpha*rleft)^mu+k)/probtot;
    probright=(pright*(1+alpha*rright)^mu+k)/probtot;
    probdown=(pdown*(1+alpha*rdown)^mu+k)/probtot;
    probup=(pup*(1+alpha*rup)^mu+k)/probtot;

    r=rand;

    
    if r<probleft
        wx(i)=wx(i)-1;
    elseif r<probleft+probright
        wx(i)=wx(i)+1;
    elseif r<probleft+probright+probdown
        wy(i)=wy(i)-1;
    else
        wy(i)=wy(i)+1;
    end
    
    R(wx(i),wy(i))=R(wx(i),wy(i))+kappa*(1-R(wx(i),wy(i))); 

    R=lambda*R; 
    
 
    if mod(t,50)==0
        figure(2)
        plot((worldsize-distancebetween)/2,(worldsize+distanceup)/2,'ro')
        hold on
        plot((worldsize+distancebetween)/2,(worldsize+distanceup)/2,'ro')
        hold on
        plot((worldsize)/2,(worldsize-distanceup)/2,'ro')
        hold on
        plot(wx,wy,'x')
        hold off
        
        axis([5 55 5 55])
        M(t/50)=getframe;
    end

    
end

movie2avi(M, 'trail.avi', 'compression', 'Cinepak');