%Main function for the range expansion simulations
function [lattice,N1,N2,Ntot,Vg,P,clinehist,selfpenalty,z,nohap,meanfitness] =  DiploidMain(mu,recrate,K,N,theta,inplattice,TL)
global M
global loci
global alleffect

n = length(TL);
if n > 5
    %Use first row for shorter simulations and second row for longer (>7
    %days) simulations.
    %LD_snaps = [1,max(2,round(n/80)),max(3,round(n/20)),max(4,round(n/5)),max(5,round(n/2)),n];
    LD_snaps = [1,max(2,round(n/40)),max(3,round(n/10)),max(4,round(n/2.5)),n];
else
    LD_snaps = 1:n;
end
LD_count = 1;

cline = zeros(loci,M+1);
clinehist = zeros(loci,M+1,length(LD_snaps));
lattice = cell(length(LD_snaps),M+1);
haplotypes = cell(1,M+1);
meanfitness = zeros(n,M+1);
nohap = zeros(n,M+1);
z = zeros(n,M+1);
Vg = zeros(n,M+1);
N1 = zeros(n,M+1);
N2 = zeros(n,M+1);
P = zeros(n,M+1);
selfpenalty = zeros(n,M+1);

for l = 1:n
    for k = 1:TL(l)
        %Step 1: Selection and creation of gametes
        [gametes,Ngam,r] = DiploidSelection(inplattice,N,theta,K,recrate);
        %Step 2: Mutations
        gametes = DiploidMutation(gametes,mu,Ngam);
        %Step 3: Fertilisation
        [inplattice,Np,Pself,Ntot] = DiploidFertilisation(gametes,Ngam);
        %Step 4: Dispersal
        [inplattice,N] = DiploidMigration(inplattice,Np,Ntot);
    end
    
    %Local population size
    N1(l,:) = N;
    N2(l,:) = Np;
    
    for in = 1:M+1
        %Penalty for disallowing selfing (0 if selfing allowed at no cost)
        selfpenalty(l,in) = (Np(in) - floor(Ngam(in)/2))/Np(in);
        %Clines in allele frequencies
    	cline(:,in) = sum(sum(inplattice{1,in},1),3)/(2*alleffect*N(in));
        %Each haplotype
        haplotypes{1,in} = cat(2,reshape(inplattice{1,in}(1,:,:),size(inplattice{1,in},2),size(inplattice{1,in},3)),reshape(inplattice{1,in}(2,:,:),size(inplattice{1,in},2),size(inplattice{1,in},3)));
        %Unique haplotypes
        nohap(l,in) = size(unique(haplotypes{1,in}','rows'),1);
        %The average growth rate in local population
        meanfitness(l,in) = mean(2*exp(squeeze(r{1,in})));
    end
    %Local average phenotype
    z(l,:) = mean(cline);
    %Local genetic variance (LE component)
    Vg(l,:) = 2*alleffect^2*sum(cline.*(1-cline));
    %Fraction of selfing (0 if selfing completely disallowed)
    P(l,:) = Pself;
    if l == LD_snaps(LD_count)
        %All genes
        clinehist(:,:,LD_count) = cline;
        for in = 1:M+1
            lattice{LD_count,in} = logical(inplattice{1,in});
        end
        LD_count = LD_count + 1;
    end    
end
end