% GEOS 309 demo code
% Alice Bradley
% Williams College



% set up initial random map of penguins
initial_penguins = double(rand(180,360) >0.5);


% you will need to turn your land mask sideways
landseamask = ncread('lsmask.nc', 'mask')';

% make the land NaN
initial_penguins(landseamask == 0) = NaN;

%start your penguins!
penguins = initial_penguins;
figure(1)

% looping through 2000 days
for n = 1:2000

    next_penguins = penguins; % copy this to get your edge penguins to stay still
    
for r = 2:179 % we're only looping through the non-edge penguins
    for c = 2:359 % rows, then columns
        
        %gather the set of local penguins
        peer_penguins = [penguins(r-1,c), penguins(r+1,c),...
            penguins(r,c-1),penguins(r,c+1)];
        
        %check if they're coastal
        if isnan(sum(peer_penguins))
            next_penguins(r,c) = penguins(r,c);
        %if not, check if they should flip up
        elseif sum(peer_penguins) > 2
            next_penguins(r,c) = 1;
        % if they should flop down
        elseif sum(peer_penguins) < 2
            next_penguins(r,c) = 0;
        % or the random flop    
        else 
            next_penguins(r,c) = rand(1) > 0.51; % this controls the ...
            %likelihood of penguins to flip
        end
        
    end
end

%calculate what fraction of non-land penguins are up
penguin_frac(n) = nanmean(penguins(:));

% this makes the 20 subplots in 4 rows and 5 columns, numbering across
% subplot(4,5,n)
% imagesc(penguins)

penguins = next_penguins;

end


figure(10)
plot(penguin_frac)
xlabel('day')
ylabel('Fraction belly-up')

