function [x,y,z]=computetriangle(gridaccuracy,dominance,clutch_size_max,density_dependence,ddparameter,hybrid_fitness,larval_weekly_survival,b,T,L,dev_dist,init_pop,tmax,niche_shift)

%initialization
strain1=linspace(0,1,gridaccuracy);
strain2=linspace(0,1,gridaccuracy);
[strain1freq,strain2freq]=meshgrid(strain1,strain2);
strain1freq=strain1freq(:);
strain2freq=strain2freq(:);
strain3freq=max(0,1-strain1freq-strain2freq); % strain 3's frequency is what's left over from strain 1 and 2.
initialfreq=[strain1freq strain2freq strain3freq];
% Note that the above leads to some where strain 1 and 2 together add up to more than 100%. They're removed in the next step.
f=find(strain1freq+strain2freq>1); initialfreq(f,:)=[];
seasonality_option=0;
plotnumber=0; % if 0, then OneSeason saves time by not plotting anything, if >0, then it plots stuff into a figure with that number.
% (It's useful to be able to specify that number to be different from figure numbers that are used e.g. here in this script, 
% so OneSeason won't try to plot things into the same figure)
for i=1:size(initialfreq,1)
    InitEggs=init_pop*[initialfreq(i,:) 0 0 0]';
    [SABOE]=TrackDynamics(InitEggs,[],[],tmax,dominance,clutch_size_max,density_dependence,ddparameter,hybrid_fitness,larval_weekly_survival,b,T,L,dev_dist,niche_shift,seasonality_option,plotnumber);
    finalfreq(i,1)=sum(SABOE(end-3:end,1,1)+0.5*SABOE(end-3:end,1,4)+0.5*SABOE(end-3:end,1,5))./sum(sum(SABOE(end-3:end,1,:)));
    finalfreq(i,2)=sum(SABOE(end-3:end,1,2)+0.5*SABOE(end-3:end,1,4)+0.5*SABOE(end-3:end,1,6))./sum(sum(SABOE(end-3:end,1,:)));
    finalfreq(i,3)=sum(SABOE(end-3:end,1,3)+0.5*SABOE(end-3:end,1,5)+0.5*SABOE(end-3:end,1,6))./sum(sum(SABOE(end-3:end,1,:)));
    relativesuccess=finalfreq(i,:)./(initialfreq(i,:)+eps); [tmp,best(i)]=max(relativesuccess);
end

best=best(:); % to make it vertical
% the output will be strains 1, 2, 3 denoted as x y z, with initial
% frequency followed by final frequency and finally (in the 3rd column)
% whether this strain was the biggest winner (1 = yes, 0 = no), that last
% info will be used for choosing plot colours in the 'xmas tree'
x=[initialfreq(:,1) finalfreq(:,1) (best==1)];
y=[initialfreq(:,2) finalfreq(:,2) (best==2)];
z=[initialfreq(:,3) finalfreq(:,3) (best==3)];
