function p=shell;
%Malcolm Walsh, Chris King
%Follows the algorithm of 
%Kusch and Markus J. theor. Biol. (1996) 178, 333?340
%Replicating several species of mollusc shell.
high=1000;
wide=500;

%r=1; R=2; w1=5; w2=10; m0=0; m1=.3; p=2*10^-3; d=0;
%r=2; R=0; w1=6; w2=35; m0=0; m1=0.05; p=2*10^-3; d=0.1; %KM fig 11
r=2; R=0; w1=10; w2=48; m0=0; m1=0; p=2*10^-3; d=0;  %KM fig 4
%r=3; R=0; w1=5; w2=12; m0=0; m1=0.22; p=4*10^-3; d=0.19; %KM fig 10


A=[];
B=[];
A1=[];
B1=[];
A2=[];
B2=[];
A=zeros(high,wide);         %inhibitor matrix
B=zeros(high,wide);         %activating matrix

A1=zeros(1,wide);
A2=zeros(1,wide);
B1=zeros(1,wide);
B2=zeros(1,wide);

for j=1:wide;
     k=rand;
    if k>0.5;
       B(1,j)=1;
    else
    end
end

for i=1:(high-1);
     for j=(1+max(r,R)):(wide-max(r,R));

         N=0;
         for h=(j-r):(j+r);    %for all neighbours of the cell
             if B(i,h)==1;    %adds up how many neighbours are occupied
                N=N+1;
             else
             end
         end
         if B(i,j)==1;
             N=N-1;
         else
         end

         if A(i,j)>=1;
             if d==0
                A1(j)=A(i,j)-1;
             else
                A1(j)=round(A(i,j)*(1-d));
             end
         else
             A1(j)=0;
         end

         if B(i,j)==0;
             x=rand;
             if x<p;
                 B1(j)=1;
             else
                 B1(j)=0;
             end
         else
             B1(j)=1;
         end

         if B1(j)==1;
             A2(j)=A1(j)+w1;
         else
             A2(j)=A1(j);
         end

         if B1(j)==0;
             if N>round(m0+m1*A2(j));
                 B2(j)=1;
             else
                 B2(j)=B1(j);
             end
         else
            B2(j)=B1(j);
         end

         average=0;
         dummy=0;
         for g=(j-R):(j+R);    %for all neighbours of the cell
             dummy=dummy+A2(g);
         end
         average=(1/((2*R)+1))*(dummy);
         A(i+1,j)=round(average);

         if A(i+1,j)>=w2;
             B(i+1,j)=0;
         else
             B(i+1,j)=B2(j);
         end
     end
end
for v=1:high;
     for y=1:wide;
     Q(v,y)=B((high+1)-v,y);
     end
end
%Q2=Q(1+max(r,R):(high-max(r,R)),5:wide);
imat=zeros(401,401,3);
Q2=zeros(401,401);
for i=1:401
    Q2(i,:)=Q(401-i+1,50:450);
end
Q2=~Q2;

for i=1:3
    imat(:,:,i)=Q2;
end
image(imat);
%cellgraph(400,400,transpose(Q2))
%cellgraph(wide,high,transpose(Q))