clear;
dominance="S"; % Em1>Em2>Ef
clutch_size_max=[100 100 100 100]; % Constant clutch size
larval_weekly_survival=0.85; %Survival before density dependence
T=0.25; % how long one mating takes
L=1; % length of the mating season within a night
dev_dist="dev_unif"; %Development distribution uniform
init_pop=10000;
tmax=44; %duration of a season 
niche_shift=2; %Age at niche shift, but will be ignored if desity_dependence="B"
ddparameter=[1e-8*1 1e-8 1e-8]; % Strength of density dependent mortality when no niche shift,If niche shift:overwritten in OneSeason.m
hybrid_fitness=NaN ; % Only difference among the 3 files producing fig 3 and fig 4
%if hybrid_fitness is NaN then strains don't mix at all, allows looking at coexistence chances 'purely ecologically'

% these are the two that we vary; this'll create 4 x 3 triangles
density_dependence="B"; %When "B": ignore niche_shift in line 11
b=NaN; % mate search efficiency; if you set this to NaN then we assume that everyone mates
% above, b is a vectorbecause this script creates 1 plot for each of them (see loop below)

figure(1); clf;
columns=length(density_dependence);
rows=length(b);
subplotno=0;

for i=1:rows
    for j=1:columns
        subplotno=subplotno+1;
        [i j]
        gridaccuracy=31; %for lines/arrows in figure
        [x,y,z]=computetriangle(gridaccuracy,dominance,clutch_size_max,density_dependence(j),ddparameter,hybrid_fitness,larval_weekly_survival,b(i),T,L,dev_dist,init_pop,tmax,niche_shift);
        % the following 'eval' commands adapt the variable names to the
        % current values of i and j (they create x11 y11 z11, then x12 y12
        % z12, and so on
        eval(['x' int2str(i) int2str(j) '=x;']);
        eval(['y' int2str(i) int2str(j) '=y;']);
        eval(['z' int2str(i) int2str(j) '=z;']);

        %figure(1); subplot(rows,columns*2,subplotno); plottriangle(x,y,z,titletext); drawnow;
            subplotno=subplotno+1;
            
        gridaccuracy=31*3; %For heatmap in figure
        [x,y,z]=computetriangle(gridaccuracy,dominance,clutch_size_max,density_dependence(j),ddparameter,hybrid_fitness,larval_weekly_survival,b(i),T,L,dev_dist,init_pop,tmax,niche_shift);

        eval(['hx' int2str(i) int2str(j) '=x;']);
        eval(['hy' int2str(i) int2str(j) '=y;']);
        eval(['hz' int2str(i) int2str(j) '=z;']);

          

    end
end

save('DataFig2_b_variables.mat')


clear;
dominance="S"; % Em1>Em2>Ef
clutch_size_max= [60 90 120 150]; % Constant clutch size
larval_weekly_survival=0.85; %Survival before density dependence
T=0.25; % how long one mating takes
L=1; % length of the mating season within a night
dev_dist="dev_unif"; %Development distribution uniform
init_pop=10000;
tmax=44; %duration of a season 
niche_shift=2; %Age at niche shift, but will be ignored if desity_dependence="B"
ddparameter=[1e-8*1 1e-8 1e-8]; % Strength of density dependent mortality when no niche shift,If niche shift:overwritten in OneSeason.m
hybrid_fitness=NaN ; % Only difference among the 3 files producing fig 3 and fig 4
%if hybrid_fitness is NaN then strains don't mix at all, allows looking at coexistence chances 'purely ecologically'

% these are the two that we vary; this'll create 4 x 3 triangles
density_dependence="B"; %When "B": ignore niche_shift in line 11
b=NaN; % mate search efficiency; if you set this to NaN then we assume that everyone mates
% above, b is a vectorbecause this script creates 1 plot for each of them (see loop below)

figure(1); clf;
columns=length(density_dependence);
rows=length(b);
subplotno=0;

for i=1:rows
    for j=1:columns
        subplotno=subplotno+1;
        [i j]
        gridaccuracy=31; %for lines/arrows in figure
        [x,y,z]=computetriangle(gridaccuracy,dominance,clutch_size_max,density_dependence(j),ddparameter,hybrid_fitness,larval_weekly_survival,b(i),T,L,dev_dist,init_pop,tmax,niche_shift);
        % the following 'eval' commands adapt the variable names to the
        % current values of i and j (they create x11 y11 z11, then x12 y12
        % z12, and so on
        eval(['x' int2str(i) int2str(j) '=x;']);
        eval(['y' int2str(i) int2str(j) '=y;']);
        eval(['z' int2str(i) int2str(j) '=z;']);

        %figure(1); subplot(rows,columns*2,subplotno); plottriangle(x,y,z,titletext); drawnow;
            subplotno=subplotno+1;
            
        gridaccuracy=31*3; %For heatmap in figure
        [x,y,z]=computetriangle(gridaccuracy,dominance,clutch_size_max,density_dependence(j),ddparameter,hybrid_fitness,larval_weekly_survival,b(i),T,L,dev_dist,init_pop,tmax,niche_shift);

        eval(['hx' int2str(i) int2str(j) '=x;']);
        eval(['hy' int2str(i) int2str(j) '=y;']);
        eval(['hz' int2str(i) int2str(j) '=z;']);

          

    end
end

save('DataFig2_c_variables.mat')


clear;
dominance="S"; % Em1>Em2>Ef
clutch_size_max= [40 80 120 160]; % Constant clutch size
larval_weekly_survival=0.85; %Survival before density dependence
T=0.25; % how long one mating takes
L=1; % length of the mating season within a night
dev_dist="dev_unif"; %Development distribution uniform
init_pop=10000;
tmax=44; %duration of a season 
niche_shift=2; %Age at niche shift, but will be ignored if desity_dependence="B"
ddparameter=[1e-8*1 1e-8 1e-8]; % Strength of density dependent mortality when no niche shift,If niche shift:overwritten in OneSeason.m
hybrid_fitness=NaN ; % Only difference among the 3 files producing fig 3 and fig 4
%if hybrid_fitness is NaN then strains don't mix at all, allows looking at coexistence chances 'purely ecologically'

% these are the two that we vary; this'll create 4 x 3 triangles
density_dependence="B"; %When "B": ignore niche_shift in line 11
b=NaN; % mate search efficiency; if you set this to NaN then we assume that everyone mates
% above, b is a vectorbecause this script creates 1 plot for each of them (see loop below)

figure(1); clf;
columns=length(density_dependence);
rows=length(b);
subplotno=0;

for i=1:rows
    for j=1:columns
        subplotno=subplotno+1;
        [i j]
        gridaccuracy=31; %for lines/arrows in figure
        [x,y,z]=computetriangle(gridaccuracy,dominance,clutch_size_max,density_dependence(j),ddparameter,hybrid_fitness,larval_weekly_survival,b(i),T,L,dev_dist,init_pop,tmax,niche_shift);
        % the following 'eval' commands adapt the variable names to the
        % current values of i and j (they create x11 y11 z11, then x12 y12
        % z12, and so on
        eval(['x' int2str(i) int2str(j) '=x;']);
        eval(['y' int2str(i) int2str(j) '=y;']);
        eval(['z' int2str(i) int2str(j) '=z;']);

        %figure(1); subplot(rows,columns*2,subplotno); plottriangle(x,y,z,titletext); drawnow;
            subplotno=subplotno+1;
            
        gridaccuracy=31*3; %For heatmap in figure
        [x,y,z]=computetriangle(gridaccuracy,dominance,clutch_size_max,density_dependence(j),ddparameter,hybrid_fitness,larval_weekly_survival,b(i),T,L,dev_dist,init_pop,tmax,niche_shift);

        eval(['hx' int2str(i) int2str(j) '=x;']);
        eval(['hy' int2str(i) int2str(j) '=y;']);
        eval(['hz' int2str(i) int2str(j) '=z;']);

          

    end
end

save('DataFig2_d_variables.mat')