%A simple probabilistic accrestion model of modular stromatoporoid
%growth based on that described by Swan and Kershaw (Palaeontology,
%1994, v.37, p. 409-423).  Written by T.Olszewski for G404-Geobiology
%9/7/2000.

clear all
clf
colormap(jet(256));
%color codes
dirt = 120;		%120 = light green
water = 98;		%98 = cyan
band1 = 256;	%256 = dark red
band2 = 165;	%165 = yellow

h = 151;   %height of raster matrix
w = 151;   %width of raster matrix
s = input('input a random number seed (mandatory):');
rand('seed',s);
geo = input('input the strength of geotropism (default = 1):');
sed = input('stength of positive effect of sediment (default = 0):');
depslug = input('sedimentation increment (0 for no deposition):');
depint = input('interval between sedimentation increments:');
depstart = input('start-up interval:');	

[xax,yax]=meshgrid(0:1:(w-1),0:-1:-(h-1));
A = 98*ones(h,w);    %raster matrix; filled with water as default
AA = 98*ones(h,w);   %copy of raster matrix
growth = band1;  %keeps track of growth band color

A(h,1:w) = dirt;     %sets base of model as a sediment layer
A(h-1,round(w/2)) = growth;    %sets initial pixel                      

left = zeros(1,h);
left(1:h) = 2;        %this set of variables initializes the 
left(h-2:h) = round(w/2)-1; %boundaries of the area to be analyzed 
right = zeros(1,h);
right(1:h) = w-1;
right(h-2:h) = round(w/2)+1;
top = h-1;              %initializes top of strom
swi = h;              %initializes sed-h2o interface
rside = round(w/2);		%initializes right side of strom
lside = round(w/2);		%initializes left side of strom

ccount = 0;
scount = 0;
quitflag = 0;
while top > 1 | rside > 0 | lside < w   %stop when strom hits either side or the top
    AA = A;
    ccount = ccount+1;          %color counter
    if ccount == 16; ccount=0; end
    if ccount < 8
    	growth=band1;
    else
    	growth=band2;
    end
    
    scount = scount + 1;        %sedimentation counter
    if depslug > 0
        if (scount >= depstart) & (rem(scount,depint) == 0)
            for b = swi:-1:(swi-depslug)
                for a = 1:w
                    if A(b,a) == water; 
                        AA(b,a) = dirt;
                        A(b,a) = dirt; 
                    end                    
                end
            end    
            swi = swi-depslug;
        end
    end
    if swi < top
    	image(A)
    	break
    end

    y = swi;
    t = top-1;
    while y >= t			%this condition allows strom to grow as long as it is not buried

        x = lside-1;
        r = rside+1;
        while x <= r

            if A(y,x) == water
                n = 0;               %number of adjacent strom pixels
                if A(y-1,x) == band1 | A(y-1,x) == band2; n=n+1; end
                if A(y+1,x) == band1 | A(y+1,x) == band2; n=n+1; end
                if A(y,x-1) == band1 | A(y,x-1) == band2; n=n+1; end
                if A(y,x+1) == band1 | A(y,x+1) == band2; n=n+1; end
                
                if n > 0 & sed > 0 & A(y-1,x) == dirt; n=n+sed; end    %this section handles
                if n > 0 & sed > 0 & A(y+1,x) == dirt; n=n+sed; end    %positive effects of
                if n > 0 & sed > 0 & A(y,x-1) == dirt; n=n+sed; end    %sediment
                if n > 0 & sed > 0 & A(y,x+1) == dirt; n=n+sed; end
                 
                if n == 1 & geo ~= 1                   %this section handles geotropism
                    if A(y-1,x) == band1 | A(y-1,x) == band2 | A(y+1,x) == band1 | A(y+1,x) == band2
                        n = geo;
                    end
                    if A(y,x-1) == band1 | A(y,x-1) == band2 | A(y,x+1) == band1 | A(y,x+1) == band2
                        n = 1/geo;
                    end
                end
                
                if n > 0                
                    if rand <= n/(n+1)
                        AA(y,x) = growth;
                        if lside > x; lside = x; end
                        if rside < x; rside = x; end
                        if top > y; top = y; end
                        if lside == 2; quitflag = 1; end
                        if rside == w-1; quitflag = 1; end
                        if top == 2; quitflag = 1; end
                    end
                end
            end
            x=x+1;
            if quitflag == 1; break; end
        end
        if quitflag == 1; break; end
        y=y-1;
    end

	A = AA;
	image(A)
	axis('equal');
    pause(0.05)

    if quitflag == 1; break; end
end 
text(round(w/2),-5,['seed=',num2str(s,5),' geo=',num2str(geo,5),' sedhelp=',num2str(sed,5),' depslug=',num2str(depslug,5),' depint=',num2str(depint,5),' depstart=',num2str(depstart,5)],'HorizontalAlignment','center');