function [NextLattice,Ctot,Dtot] = MixDetNowakLiao(InitialLattice,NeighborFilterfft,AutoNeighborFilterfft,Vt,Vr,T,R,P,S,base)
%   ===============================================
%   David Liao
%   2008 January 27
%   Ising Trio: Two strains influence each other's growth rates on a 3-d
%   lattice
%   ===============================================
%
%   NextLattice = NowakLiao(InitialLattice,Vt,Vr,T,R,P,S)
%
%   InitialLattice is a 3-d array containing any of three values
%       +1  Cooperator indicated by matrix symbol C
%        0  Vacancy indicated by matrix symbol V
%       -1  Defector indicated by matrix symbol D

%   Enumerate three lattices each containing a specific site type
L = size(InitialLattice);

C = max(InitialLattice,0);
D = max(-InitialLattice,0);
V = ones(L) - C - D;

%   ===============
%   COUNT NEIGHBORS
%   ===============

%   Count the nearest neighbors of a lattice site using fast convolution.
%   Build a filter.  The filter can be easily changed by specifying the
%   locations of non-zero elements.  
% NeighborFilter = zeros(L);
% NeighborFilter(2,1,1) = 1; NeighborFilter(1,2,1) = 1; NeighborFilter(1,1,2) = 1;
% NeighborFilter(L(1),1,1) = 1; NeighborFilter(1,L(2),1) = 1; NeighborFilter(1,1,L(3)) = 1;

%   Convolve the filter and the individual C, V, and D lattices in the
%   frequency domain.  The third filter equation is unnecessary if the
%   filters for each type of neighbor are the same, as they are now.  These
%   calculations have a computational cost O(N log N).  They're probably
%   the bottleneck of the calculation.  
NeighborsC = ifftn(NeighborFilterfft.*fftn(C));
NeighborsD = ifftn(NeighborFilterfft.*fftn(D));
NeighborsV = ifftn(NeighborFilterfft.*fftn(V));

%   ================================================
%   CONTESTANTS FEED BY PLAYING GAMES WITH NEIGHBORS
%   ================================================

%   Calculate the nutrition each site obtains.  Cooperators and defectors
%   can acquire nutrition.  Vacancies do not feed.  
UCLocal = C.*(R*NeighborsC + S*NeighborsD + Vr*NeighborsV);
UDLocal = D.*(T*NeighborsC + P*NeighborsD + Vt*NeighborsV);
%   This is the logarithm of the reproduction rate of each bacterium

%   ========================
%   NEW CODE: EXPONENTIATION
%   ========================
%   Exponentiating the previous equations gives the number of progeny each
%   site can produce (counting a sustained bacterium as one progeny)
UCLocal = C.*(base.^UCLocal);
UDLocal = D.*(base.^UDLocal);

%   Now add the number of progeny that a lattice site can receive from
%   itself and its neighbors.  

%   Add the nutrition that the cooperators and defectors acquire in the
%   neighborhood.  The self-site is the only addition necessary to the
%   neighbor filter already declared.  
% NeighborFilter(1,1,1) = 1;

UC = ifftn(AutoNeighborFilterfft.*fftn(UCLocal))/7;  % Added division by 7
UD = ifftn(AutoNeighborFilterfft.*fftn(UDLocal))/7;  % Added division by 7

%   =====================================================
%   CONTESTANTS SPEND NEIGHBORHOOD NUTRITION TO REPRODUCE
%   =====================================================

%UCD = max(UC+UD,1);
UCD = UC + UD;
UC = UC./UCD;
UD = UD./UCD;

%   Make a random lattice
% 
% FinalRand = rand(L);
% NextC = (FinalRand < UC);
% NextCD = (FinalRand < (UC + UD));

NextC = (UC > UD);
NextCD = ones(L);
NextLattice = 2*NextC - NextCD;

%   Redistribute bacteria at random

Ctot = sum(sum(sum(NextC)));
Dtot = sum(sum(sum(NextCD - NextC)));
CDtot = Ctot + Dtot;
pc = Ctot/CDtot;

RandDistribute = rand(L);
NextCRand = (RandDistribute < pc);
NextDRand = ones(L) - NextCRand;
NextLattice = NextCRand - NextDRand;