(* Define the birth functions *)
RickerB[x_, y_]=Exp[-(x+y)];
(* Pond 1 Birth *)
cS1=1;
cI1=1;
BS[x_, y_]=bS*RickerB[cS1*x,cS1*y];
BI[x_, y_]=bI*RickerB[cI1*x,cI1*y];
(* Pond 2 Birth will be alpha times that of pond 1 *)

(* Define the Model *)
(* Pond 1 *)
NewS1[S1_, I1_, S2_, I2_]=BS[S1, I1]*S1 + pS*Exp[-w1*I1]*S1+pI*r*I1+α*mS*Exp[-w2*I2]*S2+α*mI*r*I2;
NewI1[S1_, I1_, S2_, I2_]=BI[S1, I1]*I1 + pS*(1-Exp[-w1*I1])*S1 + pI*(1-r)*I1+α*mS*(1-Exp[-w2*I2])*S2+α*mI*(1-r)*I2;
(* Pond 2 *)
NewS2[S1_, I1_, S2_, I2_]=α*BS[S2, I2]*S2 + α*pS*Exp[-w2*I2]*S2+α*pI*r*I2+mS*Exp[-w1*I1]*S1+mI*r*I1;
NewI2[S1_, I1_, S2_, I2_]=α*BI[S2, I2]*I2 + α*pS*(1-Exp[-w2*I2])*S2 +α* pI*(1-r)*I2+mS*(1-Exp[-w1*I1])*S1+mI*(1-r)*I1;

(* Define the parameters *)
(* Pond 1 *)
(*bS1=7;*)
(*bI1=4;*)
pS=0.4;
pI=0.3;
mS=0.3;
mI=0.2;
r=0.15;
w1=1;
(* Pond 2 *)
w2=1;

(* Create the Bifurcation Simulation *)
condition = 500;
numPoints = 50;

b=5;

beginα = 0;
endα = Min[(1-mS)/pS,(1-mI)/pI,(1-pS)/mS,(1-pI)/mI];  (* Note: α p + m < 1 and p + α m < 1 *)
stepα=endα/200;
index=1;

listS1 = {};
listI1 = {};
listS2 = {};
listI2 = {};

Do [
numS1=2+10*Random[];
numI1=2+10*Random[];
numS2=2+10*Random[];
numI2=2+10*Random[];

bS:=b;
bI:=b;

Do [  (* Conditioning Loop -- Let the system settle. These are Numbers, not Lists *)
newS1=NewS1[numS1, numI1, numS2, numI2];
newI1=NewI1[numS1, numI1,numS2, numI2];
newS2=NewS2[numS1, numI1, numS2, numI2];
newI2=NewI2[numS1,numI1, numS2, numI2];
numS1=newS1;
numI1=newI1;
numS2=newS2;
numI2=newI2;
, {time, 1, condition}
];

listS1 = Append[listS1, {numS1}];
listI1 =  Append[listI1, {numI1}];
listS2 =  Append[listS2, {numS2}];
listI2 =  Append[listI2, {numI2}];

Do [  (* Generate the points to plot. Store in Lists of Lists *)
newS1=NewS1[numS1, numI1, numS2, numI2];
newI1=NewI1[numS1, numI1,numS2, numI2];
newS2=NewS2[numS1, numI1, numS2, numI2];
newI2=NewI2[numS1,numI1, numS2, numI2];
numS1=newS1;
numI1=newI1;
numS2=newS2;
numI2=newI2;
listS1[[index]]=Append[listS1[[index]],numS1];
listI1[[index]]=Append[listI1[[index]],numI1];
listS2[[index]]=Append[listS2[[index]],numS2];
listI2[[index]]=Append[listI2[[index]],numI2];

, {time, 1, numPoints}
];
index++;

, {α,beginα, endα, stepα}
];

S1Plot=ListPlot[Transpose[listS1],DataRange->{0,endα} ,  AxesLabel->{"α", "(Subscript[S, 1])"},PlotStyle->{Directive[PointSize[Tiny],  Black]}, AxesOrigin->{0,0}];
I1Plot=ListPlot[Transpose[listI1],DataRange->{0,endα}  ,  AxesLabel->{"α", "(Subscript[I, 1])"},PlotStyle->{Directive[PointSize[Tiny],  Black]}, AxesOrigin->{0,0}];
S2Plot=ListPlot[Transpose[listS2],DataRange->{0,endα}  ,  AxesLabel->{"α", "(Subscript[S, 2])"},PlotStyle->{Directive[PointSize[Tiny],  Black]}, AxesOrigin->{0,0}];
I2Plot=ListPlot[Transpose[listI2],DataRange->{0,endα}  ,  AxesLabel->{"α", "(Subscript[I, 2])"},PlotStyle->{Directive[PointSize[Tiny],  Black]}, AxesOrigin->{0,0}];
GraphicsGrid[{{S1Plot, I1Plot},{S2Plot, I2Plot}},Spacings->{Scaled[0.25],Scaled[0.25]}]

Export["T:/white/research/ICTCM/BifAlphaB5S1_1.pdf",S1Plot]
Export["T:/white/research/ICTCM/BifAlphaB5I1_1.pdf",I1Plot]
Export["T:/white/research/ICTCM/BifAlphaB5S2_1.pdf",S2Plot]
Export["T:/white/research/ICTCM/BifAlphaB5I2_1.pdf",I2Plot]

(* Create the Bifurcation Simulation *)
condition = 500;
numPoints = 50;

b=20;

beginα = 0;
endα = Min[(1-mS)/pS,(1-mI)/pI,(1-pS)/mS,(1-pI)/mI];  (* Note: α p + m < 1 *)
stepα=endα/200;
index=1;

listS1 = {};
listI1 = {};
listS2 = {};
listI2 = {};

Do [
numS1=2+10*Random[];
numI1=2+10*Random[];
numS2=2+10*Random[];
numI2=2+10*Random[];

bS:=b;
bI:=b;

Do [  (* Conditioning Loop -- Let the system settle. These are Numbers, not Lists *)
newS1=NewS1[numS1, numI1, numS2, numI2];
newI1=NewI1[numS1, numI1,numS2, numI2];
newS2=NewS2[numS1, numI1, numS2, numI2];
newI2=NewI2[numS1,numI1, numS2, numI2];
numS1=newS1;
numI1=newI1;
numS2=newS2;
numI2=newI2;
, {time, 1, condition}
];

listS1 = Append[listS1, {numS1}];
listI1 =  Append[listI1, {numI1}];
listS2 =  Append[listS2, {numS2}];
listI2 =  Append[listI2, {numI2}];

Do [  (* Generate the points to plot. Store in Lists of Lists *)
newS1=NewS1[numS1, numI1, numS2, numI2];
newI1=NewI1[numS1, numI1,numS2, numI2];
newS2=NewS2[numS1, numI1, numS2, numI2];
newI2=NewI2[numS1,numI1, numS2, numI2];
numS1=newS1;
numI1=newI1;
numS2=newS2;
numI2=newI2;
listS1[[index]]=Append[listS1[[index]],numS1];
listI1[[index]]=Append[listI1[[index]],numI1];
listS2[[index]]=Append[listS2[[index]],numS2];
listI2[[index]]=Append[listI2[[index]],numI2];

, {time, 1, numPoints}
];
index++;

, {α,beginα, endα, stepα}
];

S1Plot=ListPlot[Transpose[listS1],DataRange->{0,endα} ,  AxesLabel->{"α", "(Subscript[S, 1])"},PlotStyle->{Directive[PointSize[Tiny],  Black]}, AxesOrigin->{0,0}];
I1Plot=ListPlot[Transpose[listI1],DataRange->{0,endα}  ,  AxesLabel->{"α", "(Subscript[I, 1])"},PlotStyle->{Directive[PointSize[Tiny],  Black]}, AxesOrigin->{0,0}];
S2Plot=ListPlot[Transpose[listS2],DataRange->{0,endα}  ,  AxesLabel->{"α", "(Subscript[S, 2])"},PlotStyle->{Directive[PointSize[Tiny],  Black]}, AxesOrigin->{0,0}];
I2Plot=ListPlot[Transpose[listI2],DataRange->{0,endα}  ,  AxesLabel->{"α", "(Subscript[I, 2])"},PlotStyle->{Directive[PointSize[Tiny],  Black]}, AxesOrigin->{0,0}];
GraphicsGrid[{{S1Plot, I1Plot},{S2Plot, I2Plot}},Spacings->{Scaled[0.25],Scaled[0.25]}]

Export["T:/white/research/ICTCM/BifAlphaB20S1_1.pdf",S1Plot]
Export["T:/white/research/ICTCM/BifAlphaB20I1_1.pdf",I1Plot]
Export["T:/white/research/ICTCM/BifAlphaB20S2_1.pdf",S2Plot]
Export["T:/white/research/ICTCM/BifAlphaB20I2_1.pdf",I2Plot]

(* Create the Bifurcation Simulation *)
condition = 1000;
numPoints = 100;

b=60;

beginα = 0;
endα = Min[(1-mS)/pS,(1-mI)/pI,(1-pS)/mS,(1-pI)/mI];  (* Note: α p + m < 1 *)
stepα=endα/200;
index=1;

listS1 = {};
listI1 = {};
listS2 = {};
listI2 = {};

Do [
numS1=2+10*Random[];
numI1=2+10*Random[];
numS2=2+10*Random[];
numI2=2+10*Random[];

bS:=b;
bI:=b;

Do [  (* Conditioning Loop -- Let the system settle. These are Numbers, not Lists *)
newS1=NewS1[numS1, numI1, numS2, numI2];
newI1=NewI1[numS1, numI1,numS2, numI2];
newS2=NewS2[numS1, numI1, numS2, numI2];
newI2=NewI2[numS1,numI1, numS2, numI2];
numS1=newS1;
numI1=newI1;
numS2=newS2;
numI2=newI2;
, {time, 1, condition}
];

listS1 = Append[listS1, {numS1}];
listI1 =  Append[listI1, {numI1}];
listS2 =  Append[listS2, {numS2}];
listI2 =  Append[listI2, {numI2}];

Do [  (* Generate the points to plot. Store in Lists of Lists *)
newS1=NewS1[numS1, numI1, numS2, numI2];
newI1=NewI1[numS1, numI1,numS2, numI2];
newS2=NewS2[numS1, numI1, numS2, numI2];
newI2=NewI2[numS1,numI1, numS2, numI2];
numS1=newS1;
numI1=newI1;
numS2=newS2;
numI2=newI2;
listS1[[index]]=Append[listS1[[index]],numS1];
listI1[[index]]=Append[listI1[[index]],numI1];
listS2[[index]]=Append[listS2[[index]],numS2];
listI2[[index]]=Append[listI2[[index]],numI2];

, {time, 1, numPoints}
];
index++;

, {α,beginα, endα, stepα}
];

S1Plot=ListPlot[Transpose[listS1],DataRange->{0,endα} ,  AxesLabel->{"α", "(Subscript[S, 1])"},PlotStyle->{Directive[PointSize[Tiny],  Black]}, AxesOrigin->{0,0}];
I1Plot=ListPlot[Transpose[listI1],DataRange->{0,endα}  ,  AxesLabel->{"α", "(Subscript[I, 1])"},PlotStyle->{Directive[PointSize[Tiny],  Black]}, AxesOrigin->{0,0}];
S2Plot=ListPlot[Transpose[listS2],DataRange->{0,endα}  ,  AxesLabel->{"α", "(Subscript[S, 2])"},PlotStyle->{Directive[PointSize[Tiny],  Black]}, AxesOrigin->{0,0}];
I2Plot=ListPlot[Transpose[listI2],DataRange->{0,endα}  ,  AxesLabel->{"α", "(Subscript[I, 2])"},PlotStyle->{Directive[PointSize[Tiny],  Black]}, AxesOrigin->{0,0}];
GraphicsGrid[{{S1Plot, I1Plot},{S2Plot, I2Plot}},Spacings->{Scaled[0.25],Scaled[0.25]}]

Export["T:/white/research/ICTCM/BifAlphaB60S1_1.pdf",S1Plot]
Export["T:/white/research/ICTCM/BifAlphaB60I1_1.pdf",I1Plot]
Export["T:/white/research/ICTCM/BifAlphaB60S2_1.pdf",S2Plot]
Export["T:/white/research/ICTCM/BifAlphaB60I2_1.pdf",I2Plot]