function [out_pop,N1,N2,Ntot,Vg,P,cline,selfpenalty,z,nohap,meanfitness] =  DiploidMainCont(mu,recrate,K,N,theta,inp_pop,TL)
global M
global loci
global alleffect

n = length(TL);

cline = zeros(loci,M+1);
out_pop = cell(1,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
        [gametes,Ngam,r] = DiploidSelection(inp_pop,N,theta,K,recrate);
        %Step 2: Mutations
        gametes = DiploidMutation(gametes,mu,Ngam);
        %Step 3: Fertilisation
        [inp_pop,Np,Pself,Ntot] = DiploidFertilisation(gametes,Ngam);
        %Step 4: Dispersal
        [inp_pop,N] = DiploidMigration(inp_pop,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);
        %selfpenalty2(l,in) = selfpenalty(l,in)/Np(in);
        %Clines in allele frequencies
    	cline(:,in) = sum(sum(inp_pop{1,in},1),3)/(2*alleffect*N(in));
        %Each haplotype
        haplotypes{1,in} = cat(2,reshape(inp_pop{1,in}(1,:,:),size(inp_pop{1,in},2),size(inp_pop{1,in},3)),reshape(inp_pop{1,in}(2,:,:),size(inp_pop{1,in},2),size(inp_pop{1,in},3)));
        %Unique haplotypes
        nohap(l,in) = size(unique(haplotypes{1,in}','rows'),1);
        %The average fitness in local population
        meanfitness(l,in) = mean(2*exp(squeeze(r{1,in})));
    end
    %Local average phenotype
    z(l,:) = mean(cline);
    %Local genetic variance
    Vg(l,:) = 2*alleffect^2*sum(cline.*(1-cline));
    %Fraction of selfing (0 if selfing completely disallowed)
    P(l,:) = Pself;
end
for in = 1:M+1
    %All genes
    out_pop{1,in} = logical(inp_pop{1,in});
end
end