(* To be appended to a Mathematica code after the FHN equations have been solved and before the output of number of excited cells *)

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}];

(*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]]};

(*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]]]]]]];

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]]};

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}];

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

(* Duplicate Connections *)

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];

(* Matrix Plot *)

FinalConnectionMatrix = MatrixPlot[ConnectionsMatrix,ColorRules->{-3->Green,-2->Black,-1->Gray,0->Gray},ColorFunction->(Blend[{White,Red},#]&),AspectRatio->1];


Export["Network Topology/GenerationAnalysis_Iso_Prob_"<>ToString[prob]<>"_Pert_"<>ToString[kick]<>"_Iteration_"<>ToString[JJ]<>".xls",Gaverage];

Export["Network Topology/GenerationAnalysis_Iso_Image_Prob_"<>ToString[prob]<>"_Pert_"<>ToString[kick]<>"_Iteration_"<>ToString[JJ]<>".jpg",FinalConnectionMatrix];