function [NextLattice,NextFoodLattice,NextGarbageLattice] = DoubleBrothF01(InitialLattice,MaskLattice,InitialFoodLattice,InitialGarbageLattice,fc,fd,g,rc,rd,rg,rho,LysingRate,DFchannel,DGchannel,DFres,FoodGaussianfft,FoodDiffusionSigma,GarbageGaussianfft,GarbageDiffusionSigma,NeighborFilterfft,MHPCoords)
%   ===============================================
%   David Liao
%   2008 February 10
%   DoubleBrothF: Two bacterial strains influence each other's growth
%   rates by competing for nutrients in solution.  Two MHPs interact
%   through a small channel.  
%   ===============================================
%
%   [NextLattice,NextFoodLattice,NextGarbageLattice] =
%   DoubleBrothF01(InitialLattice,MaskLattice,InitialFoodLattice,InitialGar
%   bageLattice,fc,fd,g,rc,rd,rg,rho,LysingRate,DFchannel,DGchannel,DFres,F
%   oodGaussianfft,FoodDiffusionSigma,GarbageGaussianfft,GarbageDiffusionSi
%   gma,NeighborFilterfft,MHPCoords)
%
%   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

%   ===============================================
%   CALCULATE FEEDING ON ENTIRE LATTICE WITHOUT SEPARATING MHPS FROM
%   CHANNELS
%   ===============================================

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

C = max(InitialLattice,0);
D = max(-InitialLattice,0);

%   Each bacterium feeds from its own position according to the "hunger"
%   characteristic of its metabolic strategy
LocalFoodC = C.*min(fc,InitialFoodLattice);
LocalFoodD = D.*min(fd,InitialFoodLattice);
LocalGarbageD = D.*min(g,InitialGarbageLattice);

%   This filter has precisely 7 elements equal +1: the six nearest
%   neighbors, and the center site they surround
% 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;
% NeighborFilter(1,1,1) = 1;
% NeighborFilterfft = fftn(NeighborFilter);

%   THE NEIGHBOR FILTER MUST BE BUILT TO THE SAME SIZE AS THE FULL LATTICE,
%   WHICH INCLUDES THE SILICON WALLS AND FREE SOLUTION.  CALCULATE THE
%   FILTER OUTSIDE ONCE TO REDUCE COMPUTATIONAL COST.  

%   Determine the average nutrition acquired by cooperators and defectors
%   in the neighborhood defined by the neighborhood filter.  Multiplication
%   in frequency space corresponds to convolution in position space.  
ZeroPaddedFoodC = [zeros(1,L(2)+2,L(3)); zeros(L(1),1,L(3)) LocalFoodC zeros(L(1),1,L(3)); zeros(1,L(2)+2,L(3))];
ZeroPaddedFoodD = [zeros(1,L(2)+2,L(3)); zeros(L(1),1,L(3)) LocalFoodD zeros(L(1),1,L(3)); zeros(1,L(2)+2,L(3))];
ZeroPaddedGarbageD = [zeros(1,L(2)+2,L(3)); zeros(L(1),1,L(3)) LocalGarbageD zeros(L(1),1,L(3)); zeros(1,L(2)+2,L(3))];

PaddednFoodC = ifftn(NeighborFilterfft.*fftn(ZeroPaddedFoodC))/7;
PaddednFoodD = ifftn(NeighborFilterfft.*fftn(ZeroPaddedFoodD))/7;
PaddednGarbageD = ifftn(NeighborFilterfft.*fftn(ZeroPaddedGarbageD))/7;

nFoodC = PaddednFoodC(2:L(1)+1,2:L(2)+1,:);
nFoodD = PaddednFoodD(2:L(1)+1,2:L(2)+1,:);
nGarbageD = PaddednGarbageD(2:L(1)+1,2:L(2)+1,:);

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

%   Assume that bacteria reproduce with probability proportional to their
%   ingested nutrition
ProbC = nFoodC/(rc*fc);
ProbD = nFoodD/(rd*fd) + nGarbageD/(rg*g);

%   Reproduction probability can fill each lattice site at most once
P = max(1,ProbC + ProbD);
ProbC = ProbC./P;
ProbD = ProbD./P;

%   Progeny distribute spatially with some randomness
ReproductionRand = rand(L);
NextPartC = (ReproductionRand < ProbC);
NextPartCD = (ReproductionRand < (ProbC + ProbD));
NextUnMasked = 2*NextPartC - NextPartCD;

%   ==========================================================
%   BOXED ANSWER: NEXT POPULATION LATTICE
%   Reproduction cannot occur in silicon walls
NextLattice = NextUnMasked.*MaskLattice;
%   ==========================================================

%   Remove consumed food from food lattices
InitialFoodLattice = InitialFoodLattice - LocalFoodC - LocalFoodD;
InitialGarbageLattice = InitialGarbageLattice - LocalGarbageD + LocalFoodC*LysingRate;
% InitialGarbageLattice = InitialGarbageLattice - LocalGarbageD + C*LysingRate;


%   ==========================================================
%   BREAKDOWN LATTICE INTO MHP AND CHANNEL PIECES
%   ==========================================================

%   LeftMHP = [r1 c1 d1 size(LeftMHPPieces)]; and analogous for Channel and
%   RightMHP;

LMC = MHPCoords(1,:);   %   Left MHP coordinates
CHC = MHPCoords(2,:);   %   Channel coordinates
RMC = MHPCoords(3,:);   %   Right MHP coordinates

FoodLatticeLeft = InitialFoodLattice(LMC(1):LMC(1)+LMC(4)-1,LMC(2):LMC(2)+LMC(5)-1,LMC(3):LMC(3)+LMC(6)-1);
FoodLatticeChannel = InitialFoodLattice(CHC(1):CHC(1)+CHC(4)-1,CHC(2):CHC(2)+CHC(5)-1,CHC(3):CHC(3)+CHC(6)-1);
FoodLatticeRight = InitialFoodLattice(RMC(1):RMC(1)+RMC(4)-1,RMC(2):RMC(2)+RMC(5)-1,RMC(3):RMC(3)+RMC(6)-1);

GarbageLatticeLeft = InitialGarbageLattice(LMC(1):LMC(1)+LMC(4)-1,LMC(2):LMC(2)+LMC(5)-1,LMC(3):LMC(3)+LMC(6)-1);
GarbageLatticeChannel = InitialGarbageLattice(CHC(1):CHC(1)+CHC(4)-1,CHC(2):CHC(2)+CHC(5)-1,CHC(3):CHC(3)+CHC(6)-1);
GarbageLatticeRight = InitialGarbageLattice(RMC(1):RMC(1)+RMC(4)-1,RMC(2):RMC(2)+RMC(5)-1,RMC(3):RMC(3)+RMC(6)-1);

%   ==========================================================
%   DIFFUSION AT INTERFACES BETWEEN MHPs AND CHANNEL
%   ==========================================================

mu = LMC(4:6);
chi = CHC(4:6);

rt = CHC(1) - LMC(1) + 1;
rb = rt + chi(1) - 1;

%   Calculate diffused food and garbage at MHP-channel interface
TransFoodLeft2Right = DFchannel*(mean(mean(mean(FoodLatticeLeft))) - mean(mean(mean(FoodLatticeRight))));

TransGarbageLeft2Right = DGchannel*(mean(mean(mean(GarbageLatticeLeft))) - mean(mean(mean(GarbageLatticeRight))));

%   Calculate surface density and arrange on a 2-d array
SigmaTransFoodLeft2Right = TransFoodLeft2Right*ones(rb-rt+1,1,mu(3))/((rb - rt + 1)*mu(3));
SigmaTransGarbageLeft2Right = TransGarbageLeft2Right*ones(rb-rt+1,1,mu(3))/((rb - rt + 1)*mu(3));

%   Incorporate the surface arrays into 3-d arrays
PrismTransFoodLeft = zeros(mu);
PrismTransFoodRight = zeros(mu);

PrismTransGarbageLeft = zeros(mu);
PrismTransGarbageRight = zeros(mu);

PrismTransFoodLeft(rt:rb,mu(2),:) = -SigmaTransFoodLeft2Right;
PrismTransFoodRight(rt:rb,1,:) = SigmaTransFoodLeft2Right;

PrismTransGarbageLeft(rt:rb,mu(2),:) = -SigmaTransGarbageLeft2Right;
PrismTransGarbageRight(rt:rb,1,:) = SigmaTransGarbageLeft2Right;

%   ==========================================================
%   RESERVOIR DIFFUSION AT NANOCHANNELS
%   ==========================================================

%   Calculate food diffused at reservoir interface
ResFoodTop = DFres*(rho - mean2(FoodLatticeRight(1,:,:)));
ResFoodBot = DFres*(rho - mean2(FoodLatticeRight(mu(1),:,:)));

%   Convert into a surface density on a surface array
SigmaResFoodTop = ResFoodTop*ones(1,mu(2),mu(3))/(mu(2)*mu(3));
SigmaResFoodBot = ResFoodBot*ones(1,mu(2),mu(3))/(mu(2)*mu(3));

%   Incorporate surface array in rectangular prism
PrismResFood = [SigmaResFoodTop; zeros(mu(1)-2,mu(2),mu(3)); SigmaResFoodBot];

%   Put together all diffusive corrections
FoodLatticeLeft = FoodLatticeLeft + PrismTransFoodLeft;
FoodLatticeRight = FoodLatticeRight + PrismTransFoodRight + PrismResFood;

GarbageLatticeLeft = GarbageLatticeLeft + PrismTransGarbageLeft;
GarbageLatticeRight = GarbageLatticeRight + PrismTransGarbageRight;

%   ==========================================================
%   DIFFUSION WITHIN MHPs
%   ==========================================================

%   Create empty region surrounding MHPs for diffusion in the horizontal
%   and vertical directions, but ignore depth

GreatestDiffLength = max(FoodDiffusionSigma,GarbageDiffusionSigma);
ZeroPadTB = zeros(GreatestDiffLength,2*GreatestDiffLength + mu(2),mu(3));
ZeroPadLR = zeros(mu(1),GreatestDiffLength,mu(3));

FoodLatticeLeftExpanded = [ZeroPadTB; ZeroPadLR FoodLatticeLeft ZeroPadLR; ZeroPadTB];
FoodLatticeRightExpanded = [ZeroPadTB; ZeroPadLR FoodLatticeRight ZeroPadLR; ZeroPadTB];
GarbageLatticeLeftExpanded = [ZeroPadTB; ZeroPadLR GarbageLatticeLeft ZeroPadLR; ZeroPadTB];
GarbageLatticeRightExpanded = [ZeroPadTB; ZeroPadLR GarbageLatticeRight ZeroPadLR; ZeroPadTB];

%   Write a filter for food diffusion
CAPMU = size(FoodLatticeLeftExpanded);

%   Diffuse food
FoodLatticeLeftExpanded = ifftn(FoodGaussianfft.*fftn(FoodLatticeLeftExpanded));
FoodLatticeRightExpanded = ifftn(FoodGaussianfft.*fftn(FoodLatticeRightExpanded));

%   Diffuse garbage
GarbageLatticeLeftExpanded = ifftn(GarbageGaussianfft.*fftn(GarbageLatticeLeftExpanded));
GarbageLatticeRightExpanded = ifftn(GarbageGaussianfft.*fftn(GarbageLatticeRightExpanded));

%   Fold back on self horizontally
HorizontalZeroPad = zeros(CAPMU(1),mu(2)-GreatestDiffLength,mu(3));
FoodLatticeLeftExpanded = [FoodLatticeLeftExpanded(:,GreatestDiffLength:-1:1,:) HorizontalZeroPad] + FoodLatticeLeftExpanded(:,GreatestDiffLength+1:GreatestDiffLength + mu(2),:) + [HorizontalZeroPad FoodLatticeLeftExpanded(:,CAPMU(2):-1:GreatestDiffLength+mu(2)+1,:)];
FoodLatticeRightExpanded = [FoodLatticeRightExpanded(:,GreatestDiffLength:-1:1,:) HorizontalZeroPad] + FoodLatticeRightExpanded(:,GreatestDiffLength+1:GreatestDiffLength + mu(2),:) + [HorizontalZeroPad FoodLatticeRightExpanded(:,CAPMU(2):-1:GreatestDiffLength+mu(2)+1,:)];

GarbageLatticeLeftExpanded = [GarbageLatticeLeftExpanded(:,GreatestDiffLength:-1:1,:) HorizontalZeroPad] + GarbageLatticeLeftExpanded(:,GreatestDiffLength+1:GreatestDiffLength + mu(2),:) + [HorizontalZeroPad GarbageLatticeLeftExpanded(:,CAPMU(2):-1:GreatestDiffLength+mu(2)+1,:)];
GarbageLatticeRightExpanded = [GarbageLatticeRightExpanded(:,GreatestDiffLength:-1:1,:) HorizontalZeroPad] + GarbageLatticeRightExpanded(:,GreatestDiffLength+1:GreatestDiffLength + mu(2),:) + [HorizontalZeroPad GarbageLatticeRightExpanded(:,CAPMU(2):-1:GreatestDiffLength+mu(2)+1,:)];

%   Fold back on self vertically
VerticalZeroPad = zeros(mu(1)-GreatestDiffLength,mu(2),mu(3));
FoodLatticeLeftExpanded = [FoodLatticeLeftExpanded(GreatestDiffLength:-1:1,:,:); VerticalZeroPad] + FoodLatticeLeftExpanded(GreatestDiffLength+1:GreatestDiffLength+mu(1),:,:) + [VerticalZeroPad; FoodLatticeLeftExpanded(CAPMU(1):-1:GreatestDiffLength+mu(1)+1,:,:)];
FoodLatticeRightExpanded = [FoodLatticeRightExpanded(GreatestDiffLength:-1:1,:,:); VerticalZeroPad] + FoodLatticeRightExpanded(GreatestDiffLength+1:GreatestDiffLength+mu(1),:,:) + [VerticalZeroPad; FoodLatticeRightExpanded(CAPMU(1):-1:GreatestDiffLength+mu(1)+1,:,:)];

GarbageLatticeLeftExpanded = [GarbageLatticeLeftExpanded(GreatestDiffLength:-1:1,:,:); VerticalZeroPad] + GarbageLatticeLeftExpanded(GreatestDiffLength+1:GreatestDiffLength+mu(1),:,:) + [VerticalZeroPad; GarbageLatticeLeftExpanded(CAPMU(1):-1:GreatestDiffLength+mu(1)+1,:,:)];
GarbageLatticeRightExpanded = [GarbageLatticeRightExpanded(GreatestDiffLength:-1:1,:,:); VerticalZeroPad] + GarbageLatticeRightExpanded(GreatestDiffLength+1:GreatestDiffLength+mu(1),:,:) + [VerticalZeroPad; GarbageLatticeRightExpanded(CAPMU(1):-1:GreatestDiffLength+mu(1)+1,:,:)];

%   ==========================================================
%   DIFFUSION WITHIN CHANNEL
%   ==========================================================

%   Make linear ramp between concentrations at faces of channel
for i = 1:chi(2)
    FoodLatticeChannel(:,i,:) = FoodLatticeLeftExpanded(rt:rb,mu(2),:) + (FoodLatticeRightExpanded(rt:rb,1,:) - FoodLatticeLeftExpanded(rt:rb,mu(2),:))*i/chi(2);
    GarbageLatticeChannel(:,i,:) = GarbageLatticeLeftExpanded(rt:rb,mu(2),:) + (GarbageLatticeRightExpanded(rt:rb,1,:) - GarbageLatticeLeftExpanded(rt:rb,mu(2),:))*i/chi(2);
end

%   Reassemble food lattice of MHP array with padding, report
NextFoodLattice = zeros(size(InitialFoodLattice));
NextGarbageLattice = zeros(size(InitialGarbageLattice));

NextFoodLattice(LMC(1):LMC(1)+LMC(4)-1,LMC(2):LMC(2)+LMC(5)-1,LMC(3):LMC(3)+LMC(6)-1) = FoodLatticeLeftExpanded;
NextFoodLattice(CHC(1):CHC(1)+CHC(4)-1,CHC(2):CHC(2)+CHC(5)-1,CHC(3):CHC(3)+CHC(6)-1) = FoodLatticeChannel;
NextFoodLattice(RMC(1):RMC(1)+RMC(4)-1,RMC(2):RMC(2)+RMC(5)-1,RMC(3):RMC(3)+RMC(6)-1) = FoodLatticeRightExpanded;

NextGarbageLattice(LMC(1):LMC(1)+LMC(4)-1,LMC(2):LMC(2)+LMC(5)-1,LMC(3):LMC(3)+LMC(6)-1) = GarbageLatticeLeftExpanded;
NextGarbageLattice(CHC(1):CHC(1)+CHC(4)-1,CHC(2):CHC(2)+CHC(5)-1,CHC(3):CHC(3)+CHC(6)-1) = GarbageLatticeChannel;
NextGarbageLattice(RMC(1):RMC(1)+RMC(4)-1,RMC(2):RMC(2)+RMC(5)-1,RMC(3):RMC(3)+RMC(6)-1) = GarbageLatticeRightExpanded;