%Kevin Lin
%Final Project
%Simulate population of organism in a fixed field.

NX = 10;
NY = 10;

itCount = 200;

Z=zeros(NX,NY);

plotX = zeros(1,itCount);
plotY = zeros(1,itCount);

pGrowth = 1;
eRate = 2;

orgCount = 100;


%Initializing org
org = zeros(4,orgCount);  %Row 1 = x, row 2 = y row 3 = Health, row 4 = dead? Columns = org #
org(3,:) = 30; %Starting health
org(4,:) = 1;  %Dead?
for oCount = 1:orgCount
    org(1,oCount) = ceil(rand()*NX);
    org(2,oCount) = ceil(rand()*NY);
end


%Initial condition for grass, which will be 50 cm high
for i = 1:NX
    Z(i,:) = 50;
end


X = linspace(0,NX,NX);
Y = linspace(0,NY,NY);

for it=1:itCount
    
    %Growth of grass
    for xGrow=1:NX
        for yGrow=1:NY
            Z(xGrow,yGrow) = Z(xGrow,yGrow) + rand()*pGrowth;
        end
    end
    
    %org eats
    for orgEat=1:orgCount
        if org(4,orgEat) ~= 0%  %Dead?
            if org(3,orgEat) > 99
                org(3,orgEat) = 50;
                continue;
            end
            if Z(org(1,orgEat),org(2,orgEat)) > 30     %Only eats anything longer than 30 cm, big mouth
                Z(org(1,orgEat),org(2,orgEat)) = Z(org(1,orgEat),org(2,orgEat)) - eRate;
                org(3,orgEat) = org(3,orgEat) + 10;   %Increase health by 2
            else
                org(3,orgEat) = org(3,orgEat) - 50;
                if org(3,orgEat) <= 0
                    org(4,orgEat) = 0;
                end
            end
        end
    end
    birth = 0;
    %Org moves/Birth
    for orgMove=1:orgCount
        if org(4,orgMove) == 0 && birth < 2  %Create at most 2
            org(4,orgMove) = 1;
            org(3,orgMove) = 30;
            org(1,orgMove) = ceil(rand()*NX);
            org(2,orgMove) = ceil(rand()*NY);
            birth = birth + 1;
        end
        if org(4,orgMove) == 0
            continue;
        end
        xr = org(1,orgMove) + 1;
        xl = org(1,orgMove) - 1;
        xu = org(2,orgMove) + 1;
        xd = org(2,orgMove) - 1;
        
        %Boundary Conditions
        if xr > NX
            xr = NX;
        end
        if xl < 1
            xl = 1;
        end
        if xu > NX
            xu = NX;
        end
        if xd < 1
            xd = 1;
        end
        choose = [Z(xr,org(2,orgMove)),Z(xl,org(2,orgMove)),Z(org(1,orgMove),xu),Z(org(1,orgMove),xd)];
        choose = sort(choose);
        
        %The walk
        p = rand();
        if p < .1
            if choose(1) == Z(xr,org(2,orgMove))
                org(1,orgMove) = xr;
            elseif choose(1) == Z(xl,org(2,orgMove))
                org(1,orgMove) = xl;
            elseif choose(1) == Z(org(1,orgMove),xu)
                org(2,orgMove) = xu;
            else
                org(2,orgMove) = xd;
            end
        elseif p <= .25
            if choose(2) == Z(xr,org(2,orgMove))
                org(1,orgMove) = xr;
            elseif choose(2) == Z(xl,org(2,orgMove))
                org(1,orgMove) = xl;
            elseif choose(2) == Z(org(1,orgMove),xu)
                org(2,orgMove) = xu;
            else
                org(2,orgMove) = xd;
            end
        elseif p <= 0.5
            if choose(3) == Z(xr,org(2,orgMove))
                org(1,orgMove) = xr;
            elseif choose(3) == Z(xl,org(2,orgMove))
                org(1,orgMove) = xl;
            elseif choose(3) == Z(org(1,orgMove),xu)
                org(2,orgMove) = xu;
            else
                org(2,orgMove) = xd;
            end
        else
            if choose(4) == Z(xr,org(2,orgMove))
                org(1,orgMove) = xr;
            elseif choose(4) == Z(xl,org(2,orgMove))
                org(1,orgMove) = xl;
            elseif choose(4) == Z(org(1,orgMove),xu)
                org(2,orgMove) = xu;
            else
                org(2,orgMove) = xd;
            end
        end
    end
    %Graph
    clf;
    cmap = zeros(10,10,3);
    cmap(:,:,2) = 1;
    
    hold on;
    %Subplot
    subplot(1,2,1);
    surf(Z,cmap);
    %colormap cmap;
    n = 1;
    for oprint=1:orgCount
        if org(4,oprint) == 1
            plotOrg(n) = oprint;
            n = n+1;
        end
    end
    a = 0;
    d = 0;
    for alive=1:orgCount
        if org(4,alive) == 1
            a = a + 1;
        end
    end
    plotX(it) = it;
    plotY(it) = a;
    hold on;
    subplot(1,2,1);
    plot(org(1,plotOrg(1:n-1)),org(2,plotOrg(1:n-1)),'o');
    axis([0 NX 0 NY 0 100]);
    view([30 10 7]);
    str = sprintf('Day %d    Alive: %d',it,a);
    title(str);
    
    subplot(1,2,2);
    plot(plotX,plotY,'r-');
    xlabel('Day');
    ylabel('Population');
    title('Population vs Day for a eating rate of 2x growth rate');
    axis([1 200 0 120]);
    set(gcf,'position',[161 198 974 420])
    drawnow
    
    saveas(gcf, ['Eat_',sprintf('%03d',it),'png'])
   % print(gcf, ['2xEat',sprintf('%03d',it),'.png'], '-dpng','-r100')
    
end
