function []=Simulation_Sexual_Selection(N,r,s,MaxLoops)
%
% Mixed_Games(Population_size, decay, sample_size, MaxLoops)
%    where N is the population size.
%    decay is the loss for having a long tail 
%         [survial = exp(-decay*tail_length) ]
%    s is the number of males each female samples before selecting the best
%    MaxLoops is the number of loops
%

if nargin<4
    MaxLoops=10000;
    if nargin<3
        s=5;
        if nargin<2
            r=0.1;
            if nargin<1
                N=1000;
            end
        end
    end
end

r=abs(r); % just in case they give negative value
%randomly asign males & females
Sex=round(rand(N,1));
%randomly asign tail length
Tail=abs(randn(N,1));

for L=1:MaxLoops
    
    Female=find(Sex==0); Male=find(Sex==1);
    
    %Males die due to tail length
    
    Surv_Males=Male(find(rand(length(Male),1)<exp(-r*Tail(Male))));
    
    if length(Surv_Males)==0
        error('All Males have died out');
    end
    if length(Female)==0
        error('All Females have died out');
    end
    
    % Females choose s random males & pick the best
    
    for i=1:length(Female)
        f=Female(i);
        m=Surv_Males(randsample(length(Surv_Males),min([length(Surv_Males) s])));
        [y j]=max(Tail(m));
        best(f)=m(j);
    end
    
    % Now generate next year's offspring
    
    F=Female(ceil(length(Female)*rand(N,1))); % which females produce each of the N offspring
    Sex=round(rand(N,1));
    Tail=abs((Tail(F)+Tail(best(F)))/2 + randn(N,1)); % tail length is average of male and female + mutation
    
    
    % Now plot the data
    MT=max(Tail); x=[0:0.01:1]*MT;
    hM=hist(Tail(Sex==1),x);
    hF=hist(Tail(Sex==0),x);
    
    h=bar(x,[hM; hF]'/N,'stacked'); 
    hold on;
    H=plot(x,exp(-r*x),'-k','LineWidth',2);
    legend([h H],'Males','Females','Survival');
    axis([0 MT 0 1.1*max(hM+hF)/N]);
    hold off;
    drawnow;
    
end
end