(* Contains code for calculating correlation function, and creating a network diagram *)

Monitor[Monitor[Monitor[Monitor[

(* Cycling through connection probabilities *)

For[prob=1,prob≥0.6, prob=prob-0.1,

(* parameter values *)

A = 3; 
ϵ = 0.2;
α = 3; (* Changes height of cubic nullcline peak *)
γ = 0.05;   (* 1/γ is the gradient of the straight-line nullcline *)
v0 = 0.4; (* -1/γ v0 is the y intercept of the straight-line nullcline *)
w0 = 0.4;
T = 25;
n = 25;
m = 25;

(* initialising parameters *)

kick = 1;
Corr=Table[i,{i,1,n-1}];
K = Table[{i,j},{i,1,n},{j,1,m}];
CorrelationFunc=Table[i,{i,1,100}];

(* iterations *)

For[ee=1,ee≤100,ee++,


initial     = Solve[wi == A*vi*(vi-α)*(1-vi)-w0 && wi == (vi-v0)/γ, {vi,wi}, Reals];
initialline = Solve[wj == A*vj*(vj-α)*(1-vj)-w0 && wj == wi/.initial[[1]],{vj,wj},Reals];
localmax    = FindMaxValue[A*vi*(vi-α)*(1-vi)-w0, {vi,vj/.initialline[[3]]}]-wi/.initial[[1]];

kickleftN = Floor[(n-1)/2];
kickrightN = Floor[(n+3)/2];	
kickcellN = Floor[(n+1)/2];

kickleftM = Floor[(m-1)/2];
kickrightM = Floor[(m+3)/2];	
kickcellM = Floor[(m+1)/2];

maxw1 = 0;

minscale=-0.25;
maxscale=1.2;


(* Set up connections *)

RRh = Table[RandomChoice[{1-prob,prob}->{0,RandomReal[{0.5,10}]}],{n},{m+1}];
RRv = Table[RandomChoice[{1-prob,prob}->{0,RandomReal[{0.5,10}]}],{n+1},{m}];

Kmean =Mean[Join[Flatten[RRh],Flatten[RRv]]];

RRh = ReplacePart[RRh,{_,1|m+1}->0];
RRv = ReplacePart[RRv,{1|n+1,_}->0];

(* Connections around central cell are not zero *)

While[Mean[{RRh[[kickcellN,kickcellM]],RRh[[kickcellN,kickcellM+1]],RRv[[kickcellN,kickcellM]],RRv[[kickcellN+1,kickcellM]]}]==0,

RRh = Table[RandomChoice[{1-prob,prob}->{0,RandomReal[{0.5,10}]}],{n},{m+1}];
RRv = Table[RandomChoice[{1-prob,prob}->{0,RandomReal[{0.5,10}]}],{n+1},{m}];

Kmean =Mean[Join[Flatten[RRh],Flatten[RRv]]];

RRh = ReplacePart[RRh,{_,1|m+1}->0];
RRv = ReplacePart[RRv,{1|n+1,_}->0]];

(* First iteration - are connections appropriate *)

fitz = NDSolve[Join[
{v[1,1]'[t]==((A*v[1,1][t]*(v[1,1][t]-α)*(1-v[1,1][t])-w[1,1][t]-w0)/ϵ)
+RRh[[1,2]]*(v[1,2][t]-v[1,1][t])
+RRv[[2,1]]*(v[2,1][t]-v[1,1][t])},

			Table[v[1,i]'[t]==((A*v[1,i][t]*(v[1,i][t]-α)*(1-v[1,i][t])-w[1,i][t]-w0)/ϵ)
+RRh[[1,i]]*(v[1,i-1][t]-v[1,i][t])
+RRv[[2,i]]*(v[2,i][t]-v[1,i][t])
+RRh[[1,i+1]]*(v[1,i+1][t]-v[1,i][t]),{i,2,m-1}],
			
			{v[1,m]'[t]==((A*v[1,m][t]*(v[1,m][t]-α)*(1-v[1,m][t])-w[1,m][t]-w0)/ϵ)
+RRh[[1,m]]*(v[1,m-1][t]-v[1,m][t])
+RRv[[2,m]]*(v[2,m][t]-v[1,m][t])},

			Table[v[i,1]'[t]==((A*v[i,1][t]*(v[i,1][t]-α)*(1-v[i,1][t])-w[i,1][t]-w0)/ϵ)
+RRv[[i,1]]*(v[i-1,1][t]-v[i,1][t])
+RRh[[i,2]]*(v[i,2][t]-v[i,1][t])
+RRv[[i+1,1]]*(v[i+1,1][t]-v[i,1][t]),{i,2,n-1}],

			Table[v[i,m]'[t]==((A*v[i,m][t]*(v[i,m][t]-α)*(1-v[i,m][t])-w[i,m][t]-w0)/ϵ)
+RRv[[i,m]]*(v[i-1,m][t]-v[i,m][t])
+RRh[[i,m]]*(v[i,m-1][t]-v[i,m][t])
+RRv[[i+1,m]]*(v[i+1,m][t]-v[i,m][t]),{i,2,n-1}],

			{v[n,1]'[t]==((A*v[n,1][t]*(v[n,1][t]-α)*(1-v[n,1][t])-w[n,1][t]-w0)/ϵ)
+RRh[[n,2]]*(v[n,2][t]-v[n,1][t])
+RRv[[n,1]]*(v[n-1,1][t]-v[n,1][t])},

			Table[v[n,i]'[t]==((A*v[n,i][t]*(v[n,i][t]-α)*(1-v[n,i][t])-w[n,i][t]-w0)/ϵ)
+RRh[[n,i]]*(v[n,i-1][t]-v[n,i][t])
+RRv[[n,i]]*(v[n-1,i][t]-v[n,i][t])
+RRh[[n,i+1]]*(v[n,i+1][t]-v[n,i][t]),{i,2,m-1}],

			{v[n,m]'[t]==((A*v[n,m][t]*(v[n,m][t]-α)*(1-v[n,m][t])-w[n,m][t]-w0)/ϵ)
+RRh[[n,m]]*(v[n,m-1][t]-v[n,m][t])
+RRv[[n,m]]*(v[n-1,m][t]-v[n,m][t])},

Flatten[Table[v[i,j]'[t]== ((A*v[i,j][t]*(v[i,j][t]-α)*(1-v[i,j][t])-w[i,j][t]-w0)/ϵ)
+RRv[[i,j]]*(v[i-1,j][t]-v[i,j][t])
+RRh[[i,j]]*(v[i,j-1][t]-v[i,j][t])
+RRv[[i+1,j]]*(v[i+1,j][t]-v[i,j][t])
+RRh[[i,j+1]]*(v[i,j+1][t]-v[i,j][t]),{i,2,n-1},{j,2,m-1}]],

			Flatten[Table[w[i,j]'[t]==v[i,j][t]-γ*w[i,j][t]-v0,{i,1,n},{j,1,m}]],

			Flatten[Table[v[i,j][0]==Evaluate[initial[[1,1,2]]],{i,1,kickleftN},{j,1,m}]],
			Flatten[Table[v[i,j][0]==Evaluate[initial[[1,1,2]]],{i,kickrightN,n},{j,1,m}]],
			Flatten[Table[v[kickcellN,j][0]==Evaluate[initial[[1,1,2]]],{j,1,kickleftM}]],
			Flatten[Table[v[kickcellN,j][0]==Evaluate[initial[[1,1,2]]],{j,kickrightM,m}]],
			{v[kickcellN,kickcellM][0]==kick+Evaluate[initial[[1,1,2]]]},
			Flatten[Table[w[i,j][0]==Evaluate[initial[[1,2,2]]],{i,1,n},{j,1,m}]]],

			Flatten[Join[Table[v[i,j],{i,1,n},{j,1,m}],Table[w[i,j],{i,1,n},{j,1,m}]]],

			{t,0,T}, MaxSteps->100000];

CC=Table[Max[Last/@Level[Cases[ParametricPlot[Evaluate[{v[i,j][t],w[i,j][t]}/.fitz[[1]]],{t,0,T},PlotStyle->{{Purple,Thickness[Large]}},Background->White,PlotRange->{{-1,3.5},{-3,10}}]
,_Line,Infinity],{-2}]],{i,1,n},{j,1,m}];

Numberexcited =Count[Positive[Flatten[CC]],True];

(* No excitation *)

While[Numberexcited>(10/100)*n n,

RRh = Table[RandomChoice[{1-prob,prob}->{0,RandomReal[{0.5,10}]}],{n},{m+1}];
RRv = Table[RandomChoice[{1-prob,prob}->{0,RandomReal[{0.5,10}]}],{n+1},{m}];

Kmean =Mean[Join[Flatten[RRh],Flatten[RRv]]];

RRh = ReplacePart[RRh,{_,1|m+1}->0];
RRv = ReplacePart[RRv,{1|n+1,_}->0];


K[[kickcellN,kickcellM]]=Mean[{RRh[[kickcellN,kickcellM]],RRh[[kickcellN,kickcellM+1]],RRv[[kickcellN,kickcellM]],RRv[[kickcellN+1,kickcellM]]}];

If[K[[kickcellN,kickcellM]]==0,

RRh = Table[RandomChoice[{1-prob,prob}->{0,RandomReal[{0.5,10}]}],{n},{m+1}];
RRv = Table[RandomChoice[{1-prob,prob}->{0,RandomReal[{0.5,10}]}],{n+1},{m}];

Kmean =Mean[Join[Flatten[RRh],Flatten[RRv]]];

RRh = ReplacePart[RRh,{_,1|m+1}->0];
RRv = ReplacePart[RRv,{1|n+1,_}->0];

K[[kickcellN,kickcellM]]=Mean[{RRh[[kickcellN,kickcellM]],RRh[[kickcellN,kickcellM+1]],RRv[[kickcellN,kickcellM]],RRv[[kickcellN+1,kickcellM]]}]];

(* Closed Boundary Conditions *)



fitz = NDSolve[Join[
{v[1,1]'[t]==((A*v[1,1][t]*(v[1,1][t]-α)*(1-v[1,1][t])-w[1,1][t]-w0)/ϵ)
+RRh[[1,2]]*(v[1,2][t]-v[1,1][t])
+RRv[[2,1]]*(v[2,1][t]-v[1,1][t])},

			Table[v[1,i]'[t]==((A*v[1,i][t]*(v[1,i][t]-α)*(1-v[1,i][t])-w[1,i][t]-w0)/ϵ)
+RRh[[1,i]]*(v[1,i-1][t]-v[1,i][t])
+RRv[[2,i]]*(v[2,i][t]-v[1,i][t])
+RRh[[1,i+1]]*(v[1,i+1][t]-v[1,i][t]),{i,2,m-1}],
			
			{v[1,m]'[t]==((A*v[1,m][t]*(v[1,m][t]-α)*(1-v[1,m][t])-w[1,m][t]-w0)/ϵ)
+RRh[[1,m]]*(v[1,m-1][t]-v[1,m][t])
+RRv[[2,m]]*(v[2,m][t]-v[1,m][t])},

			Table[v[i,1]'[t]==((A*v[i,1][t]*(v[i,1][t]-α)*(1-v[i,1][t])-w[i,1][t]-w0)/ϵ)
+RRv[[i,1]]*(v[i-1,1][t]-v[i,1][t])
+RRh[[i,2]]*(v[i,2][t]-v[i,1][t])
+RRv[[i+1,1]]*(v[i+1,1][t]-v[i,1][t]),{i,2,n-1}],

			Table[v[i,m]'[t]==((A*v[i,m][t]*(v[i,m][t]-α)*(1-v[i,m][t])-w[i,m][t]-w0)/ϵ)
+RRv[[i,m]]*(v[i-1,m][t]-v[i,m][t])
+RRh[[i,m]]*(v[i,m-1][t]-v[i,m][t])
+RRv[[i+1,m]]*(v[i+1,m][t]-v[i,m][t]),{i,2,n-1}],

			{v[n,1]'[t]==((A*v[n,1][t]*(v[n,1][t]-α)*(1-v[n,1][t])-w[n,1][t]-w0)/ϵ)
+RRh[[n,2]]*(v[n,2][t]-v[n,1][t])
+RRv[[n,1]]*(v[n-1,1][t]-v[n,1][t])},

			Table[v[n,i]'[t]==((A*v[n,i][t]*(v[n,i][t]-α)*(1-v[n,i][t])-w[n,i][t]-w0)/ϵ)
+RRh[[n,i]]*(v[n,i-1][t]-v[n,i][t])
+RRv[[n,i]]*(v[n-1,i][t]-v[n,i][t])
+RRh[[n,i+1]]*(v[n,i+1][t]-v[n,i][t]),{i,2,m-1}],

			{v[n,m]'[t]==((A*v[n,m][t]*(v[n,m][t]-α)*(1-v[n,m][t])-w[n,m][t]-w0)/ϵ)
+RRh[[n,m]]*(v[n,m-1][t]-v[n,m][t])
+RRv[[n,m]]*(v[n-1,m][t]-v[n,m][t])},

Flatten[Table[v[i,j]'[t]== ((A*v[i,j][t]*(v[i,j][t]-α)*(1-v[i,j][t])-w[i,j][t]-w0)/ϵ)
+RRv[[i,j]]*(v[i-1,j][t]-v[i,j][t])
+RRh[[i,j]]*(v[i,j-1][t]-v[i,j][t])
+RRv[[i+1,j]]*(v[i+1,j][t]-v[i,j][t])
+RRh[[i,j+1]]*(v[i,j+1][t]-v[i,j][t]),{i,2,n-1},{j,2,m-1}]],

			Flatten[Table[w[i,j]'[t]==v[i,j][t]-γ*w[i,j][t]-v0,{i,1,n},{j,1,m}]],

			Flatten[Table[v[i,j][0]==Evaluate[initial[[1,1,2]]],{i,1,kickleftN},{j,1,m}]],
			Flatten[Table[v[i,j][0]==Evaluate[initial[[1,1,2]]],{i,kickrightN,n},{j,1,m}]],
			Flatten[Table[v[kickcellN,j][0]==Evaluate[initial[[1,1,2]]],{j,1,kickleftM}]],
			Flatten[Table[v[kickcellN,j][0]==Evaluate[initial[[1,1,2]]],{j,kickrightM,m}]],
			{v[kickcellN,kickcellM][0]==kick+Evaluate[initial[[1,1,2]]]},
			Flatten[Table[w[i,j][0]==Evaluate[initial[[1,2,2]]],{i,1,n},{j,1,m}]]],

			Flatten[Join[Table[v[i,j],{i,1,n},{j,1,m}],Table[w[i,j],{i,1,n},{j,1,m}]]],

			{t,0,T}, MaxSteps->100000];

CC=Table[Max[Last/@Level[Cases[ParametricPlot[Evaluate[{v[i,j][t],w[i,j][t]}/.fitz[[1]]],{t,0,T},PlotStyle->{{Purple,Thickness[Large]}},Background->White,PlotRange->{{-1,3.5},{-3,10}}]
,_Line,Infinity],{-2}]],{i,1,n},{j,1,m}];

Numberexcited =Count[Positive[Flatten[CC]],True]];

(* end of while loop *)

(* calculated the connections per generation *)

Locationmatrix = Table[{i,j},{i,1,n},{j,1,m}];

Generations = Table[Total[Abs[Locationmatrix[[kickcellN,kickcellM]]-Locationmatrix[[i,j]]]],{i,1,n},{j,1,m}];

Gaverage=Table[i,{i,1,n-1}];

NumberExcited=Table[i,{i,1,n-1}];

ExciteRatio=Table[i,{i,1,n-1}];

K[[kickcellN,kickcellM]] = Mean[{RRh[[kickcellN,kickcellM]],RRh[[kickcellN,kickcellM+1]],RRv[[kickcellN,kickcellM]],RRv[[kickcellN+1,kickcellM]]}];

(*Surrounding kicked cell*)

GenerationCoordinates=Position[Generations,1];

DD1={1};

For[s=1,s≤Length[GenerationCoordinates],s++,

i=GenerationCoordinates[[s,1]];
j=GenerationCoordinates[[s,2]];

DD1=Append[DD1,Max[Last/@Level[Cases[ParametricPlot[Evaluate[{v[i,j][t],w[i,j][t]}/.fitz[[1]]],{t,0,T},PlotStyle->{{Purple,Thickness[Large]}},Background->White,PlotRange->{{-1,3.5},{-3,10}}]
,_Line,Infinity],{-2}]]]];

NumberExcited[[1]]=Count[Positive[Flatten[Delete[DD1,1]]],True];

ExciteRatio[[1]]=NumberExcited[[1]]/Length[GenerationCoordinates];

Gaverage[[1]] = {1,Mean[{RRh[[kickcellN,kickcellM]],RRh[[kickcellN,kickcellM+1]],RRv[[kickcellN,kickcellM]],RRv[[kickcellN+1,kickcellM]]}],ExciteRatio[[1]]};

(* Surrounding kicked cell - Correlation Function*)

K[[kickcellN-1,kickcellM]] = RRv[[kickcellN,kickcellM]];
K[[kickcellN,kickcellM-1]] = RRh[[kickcellN,kickcellM]];
K[[kickcellN,kickcellM+1]] = RRh[[kickcellN,kickcellM+1]];
K[[kickcellN+1,kickcellM]] = RRv[[kickcellN+1,kickcellM]];

Corr[[1]]=Mean[{(K[[kickcellN,kickcellM]]-Kmean)(K[[kickcellN-1,kickcellM]]-Kmean),(K[[kickcellN,kickcellM]]-Kmean)(K[[kickcellN,kickcellM-1]]-Kmean),(K[[kickcellN,kickcellM]]-Kmean)(K[[kickcellN,kickcellM+1]]-Kmean),(K[[kickcellN,kickcellM]]-Kmean)(K[[kickcellN+1,kickcellM]]-Kmean)}];


(*Middle Generations *)

For[ii=2,ii≤n-2,ii++,

GenerationCoordinates=Position[Generations,ii];

valuestore={1};

For[GG=1,GG≤Length[GenerationCoordinates],GG++,

If[GenerationCoordinates[[GG,1]]-kickcellN<0,
valuestore=Append[valuestore,RRv[[GenerationCoordinates[[GG,1]]+1,GenerationCoordinates[[GG,2]]]]]];

If[GenerationCoordinates[[GG,1]]-kickcellN>0,
valuestore=Append[valuestore,RRv[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]]]]];

If[GenerationCoordinates[[GG,2]]-kickcellN<0,
valuestore = Append[valuestore,RRh[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]+1]]]];

If[GenerationCoordinates[[GG,2]]-kickcellN>0,
valuestore=Append[valuestore,RRh[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]]]]];

If[(GenerationCoordinates[[GG,1]]-kickcellN<0) &&(GenerationCoordinates[[GG,2]]-kickcellM==0),
K[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]]]=RRv[[GenerationCoordinates[[GG,1]]+1,GenerationCoordinates[[GG,2]]]]];

If[(GenerationCoordinates[[GG,1]]-kickcellN>0) && (GenerationCoordinates[[GG,2]]-kickcellM== 0),
K[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]]]=RRv[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]]]];

If[(GenerationCoordinates[[GG,1]]-kickcellN== 0) && (GenerationCoordinates[[GG,2]]-kickcellM<0),
K[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]]]=RRh[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]]]];

If[(GenerationCoordinates[[GG,1]]-kickcellN== 0) &&(GenerationCoordinates[[GG,2]]-kickcellM>0),
K[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]]]=RRh[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]+1]]];

If[(GenerationCoordinates[[GG,1]]-kickcellN<0) && (GenerationCoordinates[[GG,2]]-kickcellM<0),
K[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]]]=Mean[{RRv[[GenerationCoordinates[[GG,1]]+1,GenerationCoordinates[[GG,2]]]],RRh[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]+1]]}]];

If[(GenerationCoordinates[[GG,1]]-kickcellN>0) && (GenerationCoordinates[[GG,2]]-kickcellM<0),
K[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]]]=Mean[{RRv[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]]],RRh[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]+1]]}]];

If[(GenerationCoordinates[[GG,1]]-kickcellN<0) &&(GenerationCoordinates[[GG,2]]-kickcellM>0),
K[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]]]=Mean[{RRh[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]]],RRv[[GenerationCoordinates[[GG,1]]+1,GenerationCoordinates[[GG,2]]]]}]];

If[(GenerationCoordinates[[GG,1]]-kickcellN>0 )&& (GenerationCoordinates[[GG,2]]-kickcellM >0),
K[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]]]=Mean[{RRh[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]]],RRv[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]]]}]]];

(* Middle generations - correlation function *)

Corr[[ii]]=Mean[Flatten[Table[(K[[kickcellN,kickcellM]]-Kmean)(K[[GenerationCoordinates[[GG,1]],GenerationCoordinates[[GG,2]]]]-Kmean),{GG,1,Length[GenerationCoordinates]}]]]];


DD2={1};

For[s=1,s≤Length[GenerationCoordinates],s++,

i=GenerationCoordinates[[s,1]];
j=GenerationCoordinates[[s,2]];

DD2=Append[DD2,Max[Last/@Level[Cases[ParametricPlot[Evaluate[{v[i,j][t],w[i,j][t]}/.fitz[[1]]],{t,0,T},PlotStyle->{{Purple,Thickness[Large]}},Background->White,PlotRange->{{-1,3.5},{-3,10}}]
,_Line,Infinity],{-2}]]]];

NumberExcited[[ii]]=Count[Positive[Flatten[Delete[DD2,1]]],True];

ExciteRatio[[ii]]=NumberExcited[[ii]]/Length[GenerationCoordinates];

Gaverage[[ii]]={ii,Mean[Flatten[Delete[valuestore,1]]],ExciteRatio[[ii]]};


(*Dealing with corners *)

corneraverage={1};

GenerationCoordinates=Position[Generations,n-1];

For[GG=1,GG≤Length[GenerationCoordinates],GG++,

If[GenerationCoordinates[[GG]]=={1,1},
corneraverage =Append[corneraverage,{RRh[[1,2]],RRv[[2,1]]}]];

If[GenerationCoordinates[[GG]]== {n,1},
corneraverage =Append[corneraverage,{RRv[[n,1]],RRh[[n,2]]}]];

If[GenerationCoordinates[[GG]]== {1,n},
corneraverage =Append[corneraverage,{RRh[[1,n]],RRv[[2,n]]}]];

If[GenerationCoordinates[[GG]]=={n,n},
corneraverage =Append[corneraverage,{RRv[[n,n]],RRh[[n,n]]}]]];

DD3={1};

For[s=1,s≤Length[GenerationCoordinates],s++,

i=GenerationCoordinates[[s,1]];
j=GenerationCoordinates[[s,2]];

DD3=Append[DD3,Max[Last/@Level[Cases[ParametricPlot[Evaluate[{v[i,j][t],w[i,j][t]}/.fitz[[1]]],{t,0,T},PlotStyle->{{Purple,Thickness[Large]}},Background->White,PlotRange->{{-1,3.5},{-3,10}}]
,_Line,Infinity],{-2}]]]];

NumberExcited[[n-1]]=Count[Positive[Flatten[Delete[DD3,1]]],True];

ExciteRatio[[n-1]]=NumberExcited[[n-1]]/Length[GenerationCoordinates];

Gaverage[[n-1]]={n-1,Mean[Flatten[Delete[corneraverage,1]]],ExciteRatio[[n-1]]};

K[[1,1]]=Mean[{RRh[[1,2]],RRv[[2,1]]}];
K[[n,1]]=Mean[{RRv[[n,1]],RRh[[n,2]]}];
K[[1,n]]=Mean[{RRh[[1,n]],RRv[[2,n]]}];
K[[n,n]]=Mean[{RRv[[n,n]],RRh[[n,n]]}];

(* Corners - correlation function *)

Corr[[n-1]] =Mean[{(K[[kickcellN,kickcellM]]-Kmean)(K[[1,1]]-Kmean),(K[[kickcellN,kickcellM]]-Kmean)(K[[n,1]]-Kmean),(K[[kickcellN,kickcellM]]-Kmean)(K[[1,n]]-Kmean),(K[[kickcellN,kickcellM]]-Kmean)(K[[n,n]]-Kmean)}];

Gaverage=Prepend[Gaverage,{"Generation","Average Coupling Size","Relative Cluster Size"}];

(* Duplicate Connections - for network diagram *)

RRh2=RRh;
RRv2=RRv;

For[i=n,i≥2,i--,
RRh2=Transpose[Insert[Transpose[RRh2],Transpose[RRh][[i]],i]]];

For[i=n,i≥2,i--,
RRv2=Insert[RRv2,RRv[[i]],i]];

(* Horizontal Connections *)

RRhcoord = Partition[Flatten[Table[{i,2j},{i,1,n},{j,1,n}]],2];
RRhmatrix=Insert[RRh2,-1,RRhcoord];
RRhcoord2 = Table[{i},{i,1,n+1}];
ZeroVector = Flatten[Table[{0},{i,1,3n}]];
RRhmatrix2 = Insert[RRhmatrix,ZeroVector,RRhcoord2];

RRhcoord3=Table[{2i+1},{i,1,n-1}];
RRhmatrix3=Insert[RRhmatrix2,ZeroVector,RRhcoord3];

(* Vertical Connections *)

RRvcoord=Table[{2i},{i,1,n}];
NegCoordValue = Table[-1,{i,1,n}];
RRvmatrix=Insert[RRv2,NegCoordValue,RRvcoord];
ZeroVector2=Table[0,{i,1,3n}];
RRvcoord2 = Table[{i},{i,1,n+1}];
RRvmatrix2 = Transpose[Insert[Transpose[RRvmatrix],ZeroVector2,RRvcoord2]];
RRvcoord3=Table[{2i+1},{i,1,n-1}];
RRvmatrix3=Transpose[Insert[Transpose[RRvmatrix2],ZeroVector2,RRvcoord3]];

(* Sum the two matrices *)

ConnectionsMatrix = RRhmatrix3+RRvmatrix3;

NegVector = Table[-1,{i,1,3n}];
ConnectionsMatrix = ReplacePart[ConnectionsMatrix,{1}->NegVector];
ConnectionsMatrix = ReplacePart[ConnectionsMatrix,{3n}->NegVector];

For[i=2,i≤3n,i=i+1,
ConnectionsMatrix=ReplacePart[ConnectionsMatrix,{i,1}->-1];
ConnectionsMatrix=ReplacePart[ConnectionsMatrix,{i,3n}->-1]];

For[i=3,i≤3n,i=i+3,
ConnectionsMatrix = ReplacePart[ConnectionsMatrix,Table[{i,3j},{j,1,n}]->-1]];

For[i=4,i≤3n+1,i=i+3,
ConnectionsMatrix = ReplacePart[ConnectionsMatrix,Table[{i,3j},{j,1,n}]->-1]];

For[i=3,i≤3n,i=i+3,
ConnectionsMatrix = ReplacePart[ConnectionsMatrix,Table[{i,3j+1},{j,1,n}]->-1]];

For[i=4,i≤3n+1,i=i+3,
ConnectionsMatrix = ReplacePart[ConnectionsMatrix,Table[{i,3j+1},{j,1,n}]->-1]];

ConnectionsMatrix = ReplacePart[ConnectionsMatrix,{ (3n+1)/2,(3n+1)/2}->-3];

CorrelationFunc[[ee]]=Table[Corr[[ii]],{ii,1,n-1}]];

CorrTable=Table[CorrelationFunc[[ee]],{ee,1,50}];

(* Exports the correlation information *)

Export["CorrelationFunction_Isotropic_No_Excitation_Prob_"<>ToString[prob]<>".xls",CorrTable]


],Numberexcited],ee],kick],prob]

(* Connection matrix images, if required *)

(*FinalConnectionMatrixImage= MatrixPlot[ConnectionsMatrix,ColorRules->{-3->White,-2->Darker[Black],-1->Lighter[Gray],0->Lighter[Gray]},ColorFunction->(Blend[{White,Darker[Gray]},#]&),AspectRatio->1]*)

(*TruncatedMatrix = ConnectionsMatrix[[3*(kickcellN-6);;3*(kickcellN+6),3*(kickcellM-6);;3*(kickcellM+6)]];*)

(*TruncatedConnectionMatrixImage= MatrixPlot[TruncatedMatrix,ColorRules->{-3->White,-2->Darker[Black],-1->Lighter[Gray],0->Lighter[Gray]},ColorFunction->(Blend[{White,Darker[Gray]},#]&),AspectRatio->1]*)

(*Export["/Users/Rachel/Documents/MOAC/PhD/Modelling/FHN Model/Global Network Diagrams/UniformDistribution_FullExcitation_Prob_0.6.jpg",TruncatedConnectionMatrixImage]*)