%% This is the main code 
% This code calculates fitness values of all strategies, S and s, D and d
% S - signalling females, s - non-signalling females
% D - discriminating males (odeying signals), d - non-discriminating males 

function [Ws,WS,Wd,WD,x,y]=mainresults(scenario,T,tau,c,accuracy)
[x,y]=meshgrid(linspace(0.001,0.999,accuracy),linspace(0.001,0.999,accuracy+1)); % Values of x and y - proportions of S females and D males
md=(1-y)/tau; mD=y./(T*x*tau); % Mating probability with d males
maxmatings=100;
Ad=zeros([accuracy+1 accuracy maxmatings+1]); Ed=Ad; Ld=Ad; AD=Ad; F=Ad; 
for k=0:maxmatings
    k1=k+1; % the one that works for indexing
    Ad(:,:,k1)=poisspdf(k,md);
    Ed(:,:,k1)=poisspdf(k,md*(1-T));
    Ld(:,:,k1)=poisspdf(k,md*T);
    AD(:,:,k1)=poisspdf(k,mD*T);
    F(:,:,k1)=fecundity(c,k);
end

disp('Computing Ws and Wds...')
Ws=sum(Ad.*F,3);
Wds=(1-x)./(1-y).*Ws;

disp('Computing WS...')
WS=zeros(size(x));
for i=maxmatings:-1:0
    i1=i+1;
    for j=maxmatings:-1:0
        j1=j+1;
        for k=maxmatings:-1:0
            k1=k+1; WS=WS+Ed(:,:,i1).*Ld(:,:,j1).*AD(:,:,k1).*fecundity(c,i+j+k);
        end
    end
end

disp('Computing WdS1...')
WdS1=x./(1-y).*Ld(1).*AD(1).*sum(Ed.*F,3);

disp('Computing WdS2...')
WdS2=zeros(size(x));
for i=maxmatings:-1:0
    i1=i+1;
    for j=maxmatings:-1:1
        j1=j+1; WdS2=WdS2+Ed(:,:,i1).*Ld(:,:,j1).*fecundity(c,i+j);
    end
end
WdS2=x./(1-y).*AD(1).*WdS2;

disp('Computing WdS3 and WDS3...')
switch scenario
    case 0, % no male is any better than any other male in sperm competition
        WdS3=zeros(size(x)); WDS3=zeros(size(x));
        for i=maxmatings:-1:0
            i1=i+1;
            for j=maxmatings:-1:0
                j1=j+1;
                for k=maxmatings:-1:1
                    k1=k+1; 
                    WdS3=WdS3+Ed(i1).*Ld(j1).*Ad(k1).*fecundity(c,i+j+k).*j./(j+k);
                    WDS3=WDS3+Ed(i1).*Ld(j1).*Ad(k1).*fecundity(c,i+j+k).*k./(j+k);
                end
            end
        end
        WdS3=x./(1-y).*WdS3; WDS3=x./y.*WDS3;
    case 1, % D males beat d males in sperm competition
    WdS3=zeros(size(x)); WDS3=zeros(size(x));
        for i=maxmatings:-1:0
            i1=i+1;
            for j=maxmatings:-1:0
                j1=j+1;
                for k=maxmatings:-1:1
                    k1=k+1; WDS3=WDS3+Ed(:,:,i1).*Ld(:,:,j1).*AD(:,:,k1).*fecundity(c,i+j+k);
                end
            end
        end
        WDS3=x./y.*WDS3;
end

Wd=Wds+WdS1+WdS2+WdS3;
WD=WDS3;