(* ::Package:: *)

(* Section 1 *)
(* Section 1.1 *)
(* In this section, users need to input the number of bodies, the number of inputs and the number of carriers according to the setup of the transmission. Generally, they assigning rules of all bodies are: the body number of ground is 1; the body number of input 1 is 2, the body number of input 2 is 3, and so on. Then the body number of planet-carriers are numbered right after the last input. At last, the body of output is the last number to be assigned. In some special cases, the body of input and the planet-carrier may share the same body number. If such cases, the overlap box is checked if they share the same number. Otherwise, leave the box unchecked*)
(*Known condiditons; input by users*)
ClearAll[C1,C2,overlap1,overlap2];
DialogInput[
{TextCell["Number of bodies: "],InputField[Dynamic[numberofbodies],Number,FieldSize->5],TextCell["Number of inputs: "],InputField[Dynamic[numberofinputs],Number,FieldSize->5],
TextCell["Number of outputs: "],InputField[Dynamic[numberofoutputs],Number,FieldSize->5],
TextCell["Number of planet-carriers: "],InputField[Dynamic[numberofcarriers],Number,FieldSize->5],
(*TextCell["Output torque (The last body number, N/m): "],InputField[Dynamic[Tout],Number,FieldSize\[Rule]5],*)
(*TextCell["Slide friction coefficient: "],InputField[Dynamic[\[Mu]f],Number,FieldSize\[Rule]5],*)Row[{Button["Quit",FrontEndExecute@FrontEndToken["EvaluatorQuit"]],Button["Next - Go to the overlap setting",DialogReturn[]]}]}];

DialogInput[{TextCell["Input gear and planetary carrier overlap setting"],
Panel[Grid[{
{Manipulate[If[overlap1,bodynumberofpc1=bodynumberofinput1,bodynumberofpc1!=bodynumberofinput1;"Input1 overlap with a planetary carrier"],{overlap1,{True,False}},LocalizeVariables->False]},
{Manipulate[If[overlap2,bodynumberofpc2=bodynumberofinput2,bodynumberofpc2!=bodynumberofinput2;"Input2 overlap with a planetary carrier"],{overlap2,{True,False}},LocalizeVariables->False]}
},Alignment->Left]],Row[{Button["Clear History",FrontEndExecute@FrontEndToken["EvaluatorQuit"]],Button["Next - Go to the dimension matrix",DialogReturn[]]}]}];

(* Section 1.2 *)
(*In Section 1.2, users have to input the dimensions for matrix C1, and connection types for C2 accodring to the dimension matrix. Generally,for matirx C2,users only need to input the bottom left part from the dimension matrix since the connection types matrix is an identity matrix. Here, G stands for gear mesh connection, and B stands for bearing connection. Note that hey should be typed in as capital letters.*)
(*Dimension Matrix C1*)
With[{numberofbodies=numberofbodies,numberofinputs=numberofinputs,numberofoutputs=numberofoutputs,numberofcarriers=numberofcarriers,Tout=Tout,\[Mu]f=\[Mu]f},
mdimension=1+numberofbodies;
DialogInput[{TextCell["Input the dimension matrix"],
Panel[Grid[{
{With[{rows1=mdimension,cols1=mdimension,fieldSize=5},Manipulate[Column[{"Dimension Matrix = "Dynamic[MatrixForm[C1]]}],{C1,None},{initialized1,None},Evaluate@With[{makeRow1=Function[{rowIndex1},Map[Function[{colIndex1},InputField[Dynamic[C1[[rowIndex1,colIndex1]]],FieldSize->fieldSize]],Range[cols1]],HoldAll]},Grid[Map[makeRow1,Range[rows1]]]],Initialization->(If[Not[initialized3===True],C1=Table[Null,{rows1},{cols1}]];
initialized1=True),LocalizeVariables->False]]}
},Alignment->Left]],
B1=ConstantArray[Null,{mdimension,1}];
Do[B1[[i+1,1]]=Subscript[B, i],{i,1,numberofbodies,1}];
B1=Flatten[B1];
C1[[All,1]]=B1;
C1[[1]]=B1;
Do[C1[[i,i]]=" X ",{i,1,numberofbodies+1,1}];,Row[{Button["Cancel",FrontEndExecute@FrontEndToken["EvaluatorQuit"]],Button["Next - Go to the connection type matrix",DialogReturn[]]}]}];

DialogInput[{TextCell["Input the connection type matrix"],
Panel[Grid[{
{With[{rows2=mdimension,cols2=mdimension,fieldSize=3},Manipulate[Column[{"Connection Type Matrix = "Dynamic[MatrixForm[C2]]}],{C2,None},{initialized2,None},Evaluate@With[{makeRow2=Function[{rowIndex2},Map[Function[{colIndex2},InputField[Dynamic[C2[[rowIndex2,colIndex2]]],FieldSize->fieldSize]],Range[cols2]],HoldAll]},Grid[Map[makeRow2,Range[rows2]]]],Initialization->(If[Not[initialized3===True],C2=Table[Null,{rows2},{cols2}]];
initialized2=True),LocalizeVariables->False]]}
},Alignment->Left]],

Xcell=Reap[Do[Sow[Subscript[Xcell, i]=DiagonalMatrix[Table[" X ",{mdimension-i}],i]],{i,mdimension-1}]]//Last;
Xcell=Total[Transpose[Flatten[Xcell,1]]]/.{0->Null};
C2[[2;;mdimension,1;;mdimension]]=Xcell[[1;;mdimension-1,1;;mdimension]];

C2[[All,1]]=B1;
C2[[1]]=B1;
Do[C2[[i,i]]=" X ",{i,1,numberofbodies+1,1}];,Row[{Button["Cancel",FrontEndExecute@FrontEndToken["EvaluatorQuit"]],Button["Next - Start efficiency evaluation",DialogReturn[]]}]}];

(*DialogInput[{TextCell["Input the bearing dimensions"],
Panel[Grid[{
{With[{rows3=mdimension,cols3=mdimension,fieldSize=5},Manipulate[Column[{"Bearing dimensions = "Dynamic[MatrixForm[BD]]}],{BD,None},{initialized3,None},Evaluate@With[{makeRow3=Function[{rowIndex3},Map[Function[{colIndex3},InputField[Dynamic[BD[[rowIndex3,colIndex3]]],FieldSize\[Rule]fieldSize]],Range[cols3]],HoldAll]},Grid[Map[makeRow3,Range[rows3]]]],Initialization->(If[Not[initialized3===True],BD=Table[Null,{rows3},{cols3}]];
initialized3=True),LocalizeVariables\[Rule]False]]}
},Alignment\[Rule]Left]],

Xcell=Reap[Do[Sow[Subscript[Xcell, i]=DiagonalMatrix[Table[" X ",{mdimension-i}],i]],{i,mdimension-1}]]//Last;
Xcell=Total[Transpose[Flatten[Xcell,1]]]/.{0\[Rule]Null};BD[[2;;mdimension,1;;mdimension]]=Xcell[[1;;mdimension-1,1;;mdimension]];

BD[[All,1]]=B1;
BD[[1]]=B1;
Do[BD[[i,i]]=" X ",{i,1,numberofbodies+1,1}];,Row[{Button["Cancel",FrontEndExecute@FrontEndToken["EvaluatorQuit"]],Button["Next - Efficiency evaluation",DialogReturn[]]}]}];*)

(* Section 2 Normal Power *)
(*The evaluation of normal power.  The results include constraint matirx, Cq,dependent and independent selection matices, Si and Sd, constraint forces, lambda, external forces,  Qe, internal powers,  Pc and external powers,  Pe.*)

(*Normal power*)
CC1=Take[C1,{2,mdimension},{2,mdimension}];
CC2=Take[C2,{2,mdimension},{2,mdimension}];
(*BDM=Take[BD,{2,mdimension},{2,mdimension}];*)

(*\\\\\Form matrix qi/////*)
bodynumberofground=1;
(*Assigning input bodies*)
If[numberofinputs==1,
bodynumberofinput1=bodynumberofground+1;
,Null];
If[numberofinputs==2,
bodynumberofinput1=bodynumberofground+1;
bodynumberofinput2=bodynumberofinput1+1;
,Null];

(*overlap setting*)
If[overlap1==True&&numberofcarriers==1,
bodynumberofpc1=bodynumberofinput1;
,Null];
If[overlap1==False&&numberofcarriers==1,
bodynumberofpc1=bodynumberofground+numberofinputs+1;
,Null];

(* Assigning planetary carrier number *)
If[numberofcarriers==2,
bodynumberofpc1=bodynumberofground+numberofinputs+1;
bodynumberofpc2=bodynumberofpc1+1;
,Null];
If[numberofcarriers==3,
bodynumberofpc1=bodynumberofground+numberofinputs+1;
bodynumberofpc2=bodynumberofpc1+1;
bodynumberofpc3=bodynumberofpc2+1;
,Null];
numberofq=2 numberofbodies;

Pgear=Count[CC2,G,2];
(*Print["Number of gear meshes = ",Pgear];*)
Pbearing=Count[CC2,B,2];
(*Print["Number of bearing connections = ",Pbearing];*)
numberofconstraints=Pgear+Pbearing;


If[numberofinputs==1,
qi=ConstantArray[0,{numberofinputs+2,1}];
qi[[3,1]]=Subscript[w, bodynumberofinput1];(*input wi1*),
Null];
If[numberofinputs==2,
qi=ConstantArray[0,{numberofinputs+2,1}];
qi[[3,1]]=Subscript[w, bodynumberofinput1];(*input wi1*)
qi[[4,1]]=Subscript[w, bodynumberofinput2],(*input wi2*)
Null];
(*Print["qi = ",qi];*)

(*\\\\\Form matrix Qd/////*)
If[numberofoutputs==1,
Qd=ConstantArray[0,{numberofconstraints,1}];
Qd[[numberofconstraints,1]]=Subscript[T, numberofbodies]; 
,Null];
If[numberofoutputs==2,
Qd=ConstantArray[0,{numberofconstraints,1}];
Qd[[numberofconstraints,1]]=Subscript[T, numberofbodies]; 
Qd[[numberofconstraints-2,1]]=Subscript[T, numberofbodies-1]; 
,Null];

(*\\\\\Form Si & Sd matrices/////*)
IM=IdentityMatrix[numberofq];
Do[Subscript[S, m]=IM[[All,m]],{m,numberofq}];
If[numberofinputs==1,
Si=Transpose[Join[{Subscript[S, 1]},{Subscript[S, 2]},{Subscript[S, 2bodynumberofinput1]}]];
Sd=Delete[IM,{2bodynumberofinput1}];
Sd=Delete[Sd,{2}];
Sd=Delete[Sd,{1}];
Sd=Transpose[Sd],
Null];
If[numberofinputs==2,
Si=Transpose[Join[{Subscript[S, 1]},{Subscript[S, 2]},{Subscript[S, 2bodynumberofinput1]},{Subscript[S, 2bodynumberofinput2]}]];
Sd=Delete[IM,{2bodynumberofinput2}];
Sd=Delete[Sd,{2bodynumberofinput1}];
Sd=Delete[Sd,{2}];
Sd=Delete[Sd,{1}];
Sd=Transpose[Sd],
Null];
(*Print["Independent matrix Si = ",MatrixForm[Si]];
Print["Dependent matrix Sd = ",MatrixForm[Sd]];*)

(*Forming constraint equations*)
PT1=Reap[Do[Sow[Subscript[PG, i]=Position[CC2,G][[i]]],{i,Pgear}]]//Last(*i=1~3*);
PT2=Reap[Do[Sow[Subscript[PB, j]=Position[CC2,B][[j]]],{j,Pbearing}]]//Last(*j=1~5*);
PT=Flatten[Join[PT1,PT2],1];

Table[Subscript[CT, n][g_,h_,x_,y_]:=Subscript[q, 2g-1]+x*Subscript[q, 2g]-(Subscript[q, 2h-1]+y*Subscript[q, 2h]),{g,PT[[All,1]]},{h,{PT[[All,2]]}},{n,numberofconstraints}];
Do[Subscript[PT, n]=PT[[n]],{n,numberofconstraints}];
Con=Reap[Do[Sow[Subscript[Cons, n]=Subscript[CT, n][Subscript[PT, n][[1]],Subscript[PT, n][[2]],(CC1[[Subscript[PT, n][[1]],Subscript[PT, n][[2]]]]),(CC1[[Subscript[PT, n][[2]],Subscript[PT, n][[1]]]])]],{n,numberofconstraints}]]//Last;
Con=Flatten[Con,2];
MatrixForm[Con];

(*Form constraint matrix*)
qq=Reap[Do[Sow[Subscript[q, m]],{m,numberofq}]]//Last;
qq=Flatten[qq,2];
Cq=Coefficient[#,qq]&/@Con;
(*Print["Constraint matrix Cq = ",MatrixForm[Cq]];*)

Do[Subscript[Sb, n]=Transpose[Join[{Subscript[S, 2(Subscript[PT, n][[2]])-1]},{Subscript[S, 2(Subscript[PT, n][[2]])]}]],{n,numberofconstraints}];

(*Constraint Forces*)
AA=Transpose[Cq.Sd];
BB=Transpose[Cq.Si];
MatrixForm[AA];
MatrixForm[BB];
lambda=LinearSolve[AA,Qd];
(*Print["Constraint forces = ",MatrixForm[lambda]];*)

(*Total External Forces*)
H=Si.BB.Inverse[AA]+Sd;
MatrixForm[H];
Qe=H.Qd;
(*Print["Total external forces = ",MatrixForm[Qe]];*)

(*Internal Powers*)
Gn=-Sd.(Inverse[(Cq.Sd)]).Cq.Si+Si;
MatrixForm[Gn];
q=Gn.qi;
MatrixForm[q]; 
(*Print["The velocity and angular velocity, q = ",q//MatrixForm];*)
lamT=Transpose[lambda];
IMm=IdentityMatrix[numberofconstraints];
Do[Subscript[y, n]=Transpose[{IMm[[All,n]]}],{n,numberofconstraints}];
Do[Subscript[yT, n]=Transpose[Subscript[y, n]],{n,numberofconstraints}];
Do[Subscript[SbT, n]=Transpose[Subscript[Sb, n]],{n,numberofconstraints}];
Pc=Reap[Do[Sow[Subscript[Pc, n]=-lamT.Subscript[y, n].Subscript[yT, n].Cq.Subscript[Sb, n].Subscript[SbT, n].q],{n,numberofconstraints}]]//Last;
Pc=Simplify[Flatten[Pc,2]];
(*Print["Internal Powers, Pc = ",Pc=Simplify[Join[Pc]]//MatrixForm]*)

(*External Powers*)
Do[Subscript[x, m]=Transpose[{Subscript[S, m]}],{m,numberofq}];
Pe=Reap[Do[Sow[Subscript[Pe, m]=Transpose[Qe].Subscript[x, m].Transpose[Subscript[x, m]].q],{m,numberofq}]]//Last;
Pe=Simplify[Flatten[Pe,2]];
(*Print["The vector of all powers due to external forces, Pe = ",MatrixForm[Pe]];*)

(*B Matrix Bmatrix*)
C3=ConstantArray[0,{numberofbodies,2 numberofbodies}];
Do[C3[[q,2q-1]]=1,{q,numberofbodies}];
Do[C3[[q,2q]]=1,{q,numberofbodies}];
C3;
Peb=Simplify[C3.Pe];
(*Print["External Powers, Peb = ",MatrixForm[Simplify[Peb]]];*)


(* Section 3 Virtual Power *)
(*The evaluation of virtual power.  The results include Depend and independ selection matices,Virtual Internal Powers,  Pcv and Virtual External Powers Pebv. *)
(*--------------Virtual Power-----------------*)
(*\\\\\\\\\\\\\\\\\\\\\1 inputs & 1 carriers, same body/////////////////////*)
If[numberofinputs==1&&overlap1==True,
Do[Subscript[S, m]=IM[[All,m]],{m,numberofq}];
Siv=Transpose[Join[{Subscript[S, 2bodynumberofground]},{Subscript[S, 2bodynumberofpc1-1]},{Subscript[S, 2bodynumberofpc1]}]];
Sdv=Delete[IM,{2bodynumberofpc1}];
Sdv=Delete[Sdv,{2 bodynumberofpc1-1}];
Sdv=Delete[Sdv,{2bodynumberofground}];
Sdv=Transpose[Sdv];
qiv=ConstantArray[0,{numberofinputs+2,1}];
qiv[[1,1]]=(q[[2 bodynumberofground,1]])-(q[[2 bodynumberofpc1,1]]);(*input wi1*)
(*Vectors of the Coordinates q-dot*)
Gv=-Sdv.(Inverse[Cq.Sdv]).Cq.Siv+Siv;
MatrixForm[Gv];
qv=Gv.qiv;
MatrixForm[qv];
(*Virtual Internal Powers*)
Pcv=Reap[Do[Sow[Subscript[Pcv, n]=-lamT.Subscript[y, n].Subscript[yT, n].Cq.Subscript[Sb, n].Subscript[SbT, n].qv],{n,numberofconstraints}]]//Last;
Pcv=Simplify[Flatten[Pcv,2]];
(*Virtual External Powers*)
Pev=Reap[Do[Sow[Subscript[Pev, m]=Transpose[Qe].Subscript[x, m].Transpose[Subscript[x, m]].qv],{m,numberofq}]]//Last;
Pev=Flatten[Pev,2];
Pebv=Simplify[C3.Pev];
,Null];

(*\\\\\\\\\\\\\\\\1 inputs & 1 carriers, different body//////////////////////*)
If[numberofinputs==1&&overlap1==False,
Do[Subscript[S, m]=IM[[All,m]],{m,numberofq}];
Siv=Transpose[Join[{Subscript[S, 2bodynumberofinput1]},{Subscript[S, 2bodynumberofpc1-1]},{Subscript[S, 2bodynumberofpc1]}]];
Sdv=Delete[IM,{2bodynumberofpc1}];
Sdv=Delete[Sdv,{2 bodynumberofpc1-1}];
Sdv=Delete[Sdv,{2bodynumberofinput1}];
Sdv=Transpose[Sdv];
qiv=ConstantArray[0,{numberofinputs+2,1}];
qiv[[1,1]]=Subscript[w, bodynumberofinput1]-(q[[2 bodynumberofpc1,1]]);(*input wi1*)
(*Vectors of the Coordinates q-dot*)
Gv=-Sdv.(Inverse[Cq.Sdv]).Cq.Siv+Siv;
MatrixForm[Gv];
qv=Gv.qiv;
MatrixForm[qv];
(*Virtual Internal Powers*)
Pcv=Reap[Do[Sow[Subscript[Pcv, n]=-lamT.Subscript[y, n].Subscript[yT, n].Cq.Subscript[Sb, n].Subscript[SbT, n].qv],{n,numberofconstraints}]]//Last;
Pcv=Simplify[Flatten[Pcv,2]];
(*Virtual External Powers*)
Pev=Reap[Do[Sow[Subscript[Pev, m]=Transpose[Qe].Subscript[x, m].Transpose[Subscript[x, m]].qv],{m,numberofq}]]//Last;
Pev=Flatten[Pev,2];
Pebv=Simplify[C3.Pev];
,Null];

(*2 inputs & 1 carriers*)
If[numberofinputs==2&&numberofcarriers==1,
Do[Subscript[S, m]=IM[[All,m]],{m,numberofq}];
Siv=Transpose[Join[{Subscript[S, 2bodynumberofinput1]},{Subscript[S, 2bodynumberofinput2]},{Subscript[S, 2bodynumberofpc1-1]},{Subscript[S, 2bodynumberofpc1]}]];
Sdv=Delete[IM,{2bodynumberofpc1}];
Sdv=Delete[Sdv,{2 bodynumberofpc1-1}];
Sdv=Delete[Sdv,{2bodynumberofinput2}];
Sdv=Delete[Sdv,{2bodynumberofinput1}];
Sdv=Transpose[Sdv];
qiv=ConstantArray[0,{numberofinputs+2,1}];
qiv[[1,1]]=Subscript[w, bodynumberofinput1]-(q[[2 bodynumberofpc1,1]]);(*input wi1*)
qiv[[2,1]]=Subscript[w, bodynumberofinput2]-(q[[2 bodynumberofpc1,1]]);(*input wi2*)
(*Vectors of the Coordinates q-dot*)
Gv=-Sdv.(Inverse[Cq.Sdv]).Cq.Siv+Siv;
MatrixForm[Gv];
qv=Gv.qiv;
MatrixForm[qv];
(*Virtual Internal Powers*)
Pcv=Reap[Do[Sow[Subscript[Pcv, n]=-lamT.Subscript[y, n].Subscript[yT, n].Cq.Subscript[Sb, n].Subscript[SbT, n].qv],{n,numberofconstraints}]]//Last;
Pcv=Simplify[Flatten[Pcv,2]];
(*Virtual External Powers*)
Pev=Reap[Do[Sow[Subscript[Pev, m]=Transpose[Qe].Subscript[x, m].Transpose[Subscript[x, m]].qv],{m,numberofq}]]//Last;
Pev=Flatten[Pev,2];
Pebv=Simplify[C3.Pev];
,Null];

(*2 inputs & 2 carriers*)
If[numberofinputs==2&&numberofcarriers==2,
Do[Subscript[S, m]=IM[[All,m]],{m,numberofq}];
Siv1=Transpose[Join[{Subscript[S, 2bodynumberofinput1]},{Subscript[S, 2bodynumberofinput2]},{Subscript[S, 2bodynumberofpc1-1]},{Subscript[S, 2bodynumberofpc1]}]];
Sdv1=Delete[IM,{2bodynumberofpc1}];
Sdv1=Delete[Sdv1,{2 bodynumberofpc1-1}];
Sdv1=Delete[Sdv1,{2bodynumberofinput2}];
Sdv1=Delete[Sdv1,{2bodynumberofinput1}];
Sdv1=Transpose[Sdv1];
qiv1=ConstantArray[0,{numberofinputs+2,1}];
qiv1[[1,1]]=Subscript[w, bodynumberofinput1]-(q[[2 bodynumberofpc1,1]]);(*input wi1*)
qiv1[[2,1]]=Subscript[w, bodynumberofinput2]-(q[[2 bodynumberofpc1,1]]);(*input wi2*)
Gv1=-Sdv1.(Inverse[Cq.Sdv1]).Cq.Siv1+Siv1;
MatrixForm[Gv1];
qv1=Gv1.qiv1;
MatrixForm[qv1];
(*Virtual Internal Powers*)
Pcv1=Reap[Do[Sow[Subscript[Pcv1, n]=-lamT.Subscript[y, n].Subscript[yT, n].Cq.Subscript[Sb, n].Subscript[SbT, n].qv1],{n,numberofconstraints}]]//Last;
Pcv1=Simplify[Flatten[Pcv1,2]];
(*Virtual External Powers*)
Pev1=Reap[Do[Sow[Subscript[Pev1, m]=Transpose[Qe].Subscript[x, m].Transpose[Subscript[x, m]].qv1],{m,numberofq}]]//Last;
Pev1=Flatten[Pev1,2];
Pebv1=Simplify[C3.Pev1];

Do[Subscript[S, m]=IM[[All,m]],{m,numberofq}];
Siv2=Transpose[Join[{Subscript[S, 2bodynumberofinput1]},{Subscript[S, 2bodynumberofinput2]},{Subscript[S, 2bodynumberofpc2-1]},{Subscript[S, 2bodynumberofpc2]}]];
Sdv2=Delete[IM,{2bodynumberofpc2}];
Sdv2=Delete[Sdv2,{2 bodynumberofpc2-1}];
Sdv2=Delete[Sdv2,{2bodynumberofinput2}];
Sdv2=Delete[Sdv2,{2bodynumberofinput1}];
Sdv2=Transpose[Sdv2];
qiv2=ConstantArray[0,{numberofinputs+2,1}];
qiv2[[1,1]]=Subscript[w, bodynumberofinput1]-(q[[2 bodynumberofpc2,1]]);(*input wi1*)
qiv2[[2,1]]=Subscript[w, bodynumberofinput2]-(q[[2 bodynumberofpc2,1]]);(*input wi2*)
Gv2=-Sdv2.(Inverse[Cq.Sdv2]).Cq.Siv2+Siv2;
MatrixForm[Gv2];
qv2=Gv2.qiv2;
MatrixForm[qv2];
(*Virtual Internal Powers*)
Pcv2=Reap[Do[Sow[Subscript[Pcv2, n]=-lamT.Subscript[y, n].Subscript[yT, n].Cq.Subscript[Sb, n].Subscript[SbT, n].qv2],{n,numberofconstraints}]]//Last;
Pcv2=Simplify[Flatten[Pcv2,2]];
(*Print["Virtual Internal Powers, Pcv2 = ",MatrixForm[Pcv2]];*)
(*Virtual External Powers*)
Pev2=Reap[Do[Sow[Subscript[Pev2, m]=Transpose[Qe].Subscript[x, m].Transpose[Subscript[x, m]].qv2],{m,numberofq}]]//Last;
Pev2=Flatten[Pev2,2];
Pebv2=Simplify[C3.Pev2];
,Null];

(*2 inputs & 3 carriers*)
If[numberofinputs==2&&numberofcarriers==3,
Do[Subscript[S, m]=IM[[All,m]],{m,numberofq}];
Siv1=Transpose[Join[{Subscript[S, 2bodynumberofinput1]},{Subscript[S, 2bodynumberofinput2]},{Subscript[S, 2bodynumberofpc1-1]},{Subscript[S, 2bodynumberofpc1]}]];
Sdv1=Delete[IM,{2bodynumberofpc1}];
Sdv1=Delete[Sdv1,{2 bodynumberofpc1-1}];
Sdv1=Delete[Sdv1,{2bodynumberofinput2}];
Sdv1=Delete[Sdv1,{2bodynumberofinput1}];
Sdv1=Transpose[Sdv1];
qiv1=ConstantArray[0,{numberofinputs+2,1}];
qiv1[[1,1]]=Subscript[w, bodynumberofinput1]-(q[[2 bodynumberofpc1,1]]);(*input wi1*)
qiv1[[2,1]]=Subscript[w, bodynumberofinput2]-(q[[2 bodynumberofpc1,1]]);(*input wi2*)
Gv1=-Sdv1.(Inverse[Cq.Sdv1]).Cq.Siv1+Siv1;
MatrixForm[Gv1];
qv1=Gv1.qiv1;
MatrixForm[qv1];
(*Virtual Internal Powers*)
Pcv1=Reap[Do[Sow[Subscript[Pcv1, n]=-lamT.Subscript[y, n].Subscript[yT, n].Cq.Subscript[Sb, n].Subscript[SbT, n].qv1],{n,numberofconstraints}]]//Last;
Pcv1=Simplify[Flatten[Pcv1,3]];
(*Virtual External Powers*)
Pev1=Reap[Do[Sow[Subscript[Pev1, m]=Transpose[Qe].Subscript[x, m].Transpose[Subscript[x, m]].qv1],{m,numberofq}]]//Last;
Pev1=Flatten[Pev1,2];
Pebv1=Simplify[C3.Pev1];

Do[Subscript[S, m]=IM[[All,m]],{m,numberofq}];
Siv2=Transpose[Join[{Subscript[S, 2bodynumberofinput1]},{Subscript[S, 2bodynumberofinput2]},{Subscript[S, 2bodynumberofpc2-1]},{Subscript[S, 2bodynumberofpc2]}]];
Sdv2=Delete[IM,{2bodynumberofpc2}];
Sdv2=Delete[Sdv2,{2 bodynumberofpc2-1}];
Sdv2=Delete[Sdv2,{2bodynumberofinput2}];
Sdv2=Delete[Sdv2,{2bodynumberofinput1}];
Sdv2=Transpose[Sdv2];
qiv2=ConstantArray[0,{numberofinputs+2,1}];
qiv2[[1,1]]=Subscript[w, bodynumberofinput1]-(q[[2 bodynumberofpc2,1]]);(*input wi1*)
qiv2[[2,1]]=Subscript[w, bodynumberofinput2]-(q[[2 bodynumberofpc2,1]]);(*input wi2*)
Gv2=-Sdv2.(Inverse[Cq.Sdv2]).Cq.Siv2+Siv2;
MatrixForm[Gv2];
qv2=Gv2.qiv2;
MatrixForm[qv2];
(*Virtual Internal Powers*)
Pcv2=Reap[Do[Sow[Subscript[Pcv2, n]=-lamT.Subscript[y, n].Subscript[yT, n].Cq.Subscript[Sb, n].Subscript[SbT, n].qv2],{n,numberofconstraints}]]//Last;
Pcv2=Simplify[Flatten[Pcv2,3]];
(*Virtual External Powers*)
Pev2=Reap[Do[Sow[Subscript[Pev2, m]=Transpose[Qe].Subscript[x, m].Transpose[Subscript[x, m]].qv2],{m,numberofq}]]//Last;
Pev2=Flatten[Pev2,2];
Pebv2=Simplify[C3.Pev2];

Do[Subscript[S, m]=IM[[All,m]],{m,numberofq}];
Siv3=Transpose[Join[{Subscript[S, 2bodynumberofinput1]},{Subscript[S, 2bodynumberofinput2]},{Subscript[S, 2bodynumberofpc3-1]},{Subscript[S, 2bodynumberofpc3]}]];
Sdv3=Delete[IM,{2bodynumberofpc3}];
Sdv3=Delete[Sdv3,{2 bodynumberofpc3-1}];
Sdv3=Delete[Sdv3,{2bodynumberofinput2}];
Sdv3=Delete[Sdv3,{2bodynumberofinput1}];
Sdv3=Transpose[Sdv3];
qiv3=ConstantArray[0,{numberofinputs+2,1}];
qiv3[[1,1]]=Subscript[w, bodynumberofinput1]-(q[[2 bodynumberofpc3,1]]);(*input wi1*)
qiv3[[2,1]]=Subscript[w, bodynumberofinput2]-(q[[2 bodynumberofpc3,1]]);(*input wi2*)
Gv3=-Sdv3.(Inverse[Cq.Sdv3]).Cq.Siv3+Siv3;
MatrixForm[Gv3];
qv3=Gv3.qiv3;
MatrixForm[qv3];
(*Virtual Internal Powers*)
Pcv3=Reap[Do[Sow[Subscript[Pcv3, n]=-lamT.Subscript[y, n].Subscript[yT, n].Cq.Subscript[Sb, n].Subscript[SbT, n].qv3],{n,numberofconstraints}]]//Last;
Pcv3=Simplify[Flatten[Pcv3,3]];
(*Virtual External Powers*)
Pev3=Reap[Do[Sow[Subscript[Pev3, m]=Transpose[Qe].Subscript[x, m].Transpose[Subscript[x, m]].qv3],{m,numberofq}]]//Last;
Pev3=Flatten[Pev3,2];
Pebv3=Simplify[C3.Pev3];
,Null];

(* Section 4 Efficiency Evaluation *)
Lossfactors=Table[{Subscript[\[Mu], i]},{i,1,Pgear}];
Tout=1;
\[Mu]f=0;
(*new code for Gpowers*)
If[numberofcarriers==1,
PGPv=Position[PT,bodynumberofpc1+1];
Gpowers=Pc[[1;;Pgear]];
Gpowers=ReplacePart[Gpowers,{PGPv[[1,1]],1}->Pcv[[PGPv[[1,1]],1]]];
Gpowers=ReplacePart[Gpowers,{PGPv[[2,1]],1}->Pcv[[PGPv[[2,1]],1]]];
Gpowers=Flatten[Transpose[Abs[Gpowers]]],Null
];

If[numberofcarriers==2,
PGPv1=Position[PT,bodynumberofpc1+2];
PGPv2=Position[PT,bodynumberofpc2+2];
Gpowers=Pc[[1;;Pgear]];
Gpowers=ReplacePart[Gpowers,{PGPv1[[1,1]],1}->Pcv1[[PGPv1[[1,1]],1]]];
Gpowers=ReplacePart[Gpowers,{PGPv1[[2,1]],1}->Pcv1[[PGPv1[[2,1]],1]]];
Gpowers=ReplacePart[Gpowers,{PGPv2[[1,1]],1}->Pcv2[[PGPv2[[1,1]],1]]];
Gpowers=ReplacePart[Gpowers,{PGPv2[[2,1]],1}->Pcv2[[PGPv2[[2,1]],1]]];
Gpowers=Flatten[Transpose[Abs[Gpowers]]],Null
];

(*Bearing power loss model*)
(*Relative bearing angular velocities*)
(*Table[Subscript[VCT, n][g_,h_]:=q[[2g]]-q[[2h]],{g,PT[[All,1]]},{h,{PT[[All,2]]}},{n,numberofconstraints}];
Do[Subscript[PT, n]=PT[[n]],{n,numberofconstraints}];
VCT=Reap[Do[Sow[Subscript[VCTs, n]=Subscript[VCT, n][Subscript[PT, n][[1]],Subscript[PT, n][[2]]]],{n,numberofconstraints}]]//Last;
BearingAngluarVelocities=Abs[Flatten[VCT,1][[Pgear+1;;numberofconstraints]]];

(*Friction torques*)
BPT=Flatten[PT2,1];
dm=Flatten[Reap[Do[Sow[Subscript[DM, i]={BDM[[BPT[[i,1]],BPT[[i,2]]]]}],{i,Pbearing}]]//Last,1];
Cforces=Flatten[Abs[lambda[[Pgear+1;;numberofconstraints]]]];
Bearingpowerloss=\[Mu]f Total[Cforces.((dm/2)*BearingAngluarVelocities)];*)

(*Eff eqt*)
If[numberofinputs==1&&overlap1==True,

eff=1-(Total[Gpowers.Lossfactors])/((1/2)(Total[Abs[Flatten[Peb,2]]]))/.{Subscript[T, numberofbodies]->Tout,Subscript[w, bodynumberofinput1]->1};(*effb=1-(Total[Gpowers.Lossfactors]+Bearingpowerloss)/((1/2)(Total[Abs[Flatten[Peb,2]]]))/.{Subscript[T, numberofbodies]\[Rule]Tout,Subscript[w, bodynumberofinput1]\[Rule]1};*)
,Null];


If[numberofinputs==1&&overlap1==False,

eff=1-(Total[Gpowers.Lossfactors])/((1/2)(Total[Abs[Flatten[Peb,2]]]))/.{Subscript[T, numberofbodies]->Tout,Subscript[w, bodynumberofinput1]->1};
(*effb=1-(Total[Gpowers.Lossfactors]+Bearingpowerloss)/((1/2)(Total[Abs[Flatten[Peb,2]]]))/.{Subscript[T, numberofbodies]\[Rule]Tout,Subscript[w, bodynumberofinput1]\[Rule]1};*)
,Null];

If[numberofinputs==2&&numberofcarriers==1,

eff=1-(Total[Gpowers.Lossfactors])/((1/2)(Total[Abs[Flatten[Peb,2]]]))/.{Subscript[T, numberofbodies]->Tout,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]};(*effb=1-(Total[Gpowers.Lossfactors]+Bearingpowerloss)/((1/2)(Total[Abs[Flatten[Peb,2]]]))/.{Subscript[T, numberofbodies]\[Rule]Tout,Subscript[w, bodynumberofinput1]\[Rule]Cos[\[Theta]],Subscript[w, bodynumberofinput2]\[Rule]Sin[\[Theta]]};*)
,Null];

If[numberofinputs==2&&numberofcarriers==2,

eff=1-(Total[Gpowers.Lossfactors])/((1/2)(Total[Abs[Flatten[Peb,2]]]))/.{Subscript[T, numberofbodies-1]->k,Subscript[T, numberofbodies]->1,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]};(*effb=1-(Total[Gpowers.Lossfactors]+Bearingpowerloss)/((1/2)(Total[Abs[Flatten[Peb,2]]]))/.{Subscript[T, numberofbodies-1]\[Rule]k,Subscript[T, numberofbodies]\[Rule]1,Subscript[w, bodynumberofinput1]\[Rule]Cos[\[Theta]],Subscript[w, bodynumberofinput2]\[Rule]Sin[\[Theta]]};*)
,Null];

eff=Simplify[eff];
(*effb=Simplify[effb];*)

(* Section 5 Graphical representation for normal and virtual power *)
(*The determination of postitve power direction and the amount of power passing through each links. The size of each link will change according to the power passing through it. *)

Table[Subscript[CON, n][g_,h_]:=g\[UndirectedEdge]h,{g,PT[[All,1]]},{h,{PT[[All,2]]}},{n,numberofconstraints}];
Do[Subscript[PT, n]=PT[[n]],{n,numberofconstraints}];
Connection=Reap[Do[Sow[Subscript[CONS, n]=Subscript[CON, n][Subscript[PT, n][[1]],Subscript[PT, n][[2]]]],{n,numberofconstraints}]]//Last;
Table[Subscript[CON, n][g_,h_]:={g\[UndirectedEdge]h},{g,PT[[All,1]]},{h,{PT[[All,2]]}},{n,numberofconstraints}];
Do[Subscript[PT, n]=PT[[n]],{n,numberofconstraints}];
Connection=Reap[Do[Sow[Subscript[CONS, n]=Subscript[CON, n][Subscript[PT, n][[1]],Subscript[PT, n][[2]]]],{n,numberofconstraints}]]//Last;
ConnectionT=Flatten[Connection,1];
Connection=Flatten[Connection,2];

(*EdgeSize*)
Table[Subscript[SIZE, n][g_]:={g},{g,PT[[All,1]]},{n,numberofconstraints}];
Do[Subscript[PT, n]=PT[[n]],{n,numberofconstraints}];
Size=Reap[Do[Sow[Subscript[SIZES, n]=Subscript[SIZE, n][Subscript[PT, n][[1]]]],{n,numberofconstraints}]]//Last;
Size=Flatten[Size,1];
branchpower=Join[Size,Sqrt[Abs[Pc]],2];
branchpower=Sort@branchpower//.{a___,{x_,y1_},{x_,y2_},b___}:>{a,{x,y1+y2},b};
BP=branchpower[[1+numberofinputs;;-(numberofoutputs+1)]];

If[numberofcarriers==1,Table[Subscript[SIZE, n][g_]:={g},{g,PT[[All,1]]},{n,numberofconstraints}];
Do[Subscript[PT, n]=PT[[n]],{n,numberofconstraints}];
Size=Reap[Do[Sow[Subscript[SIZES, n]=Subscript[SIZE, n][Subscript[PT, n][[1]]]],{n,numberofconstraints}]]//Last;
Size=Flatten[Size,1];
branchpowerv=Join[Size,Sqrt[Pcv],2];
branchpowerv=Sort@branchpowerv//.{a___,{x_,y1_},{x_,y2_},b___}:>{a,{x,y1+y2},b};
BPv=branchpowerv[[1+numberofinputs;;-(numberofoutputs+1)]];
,Null];

If[numberofcarriers==2,Table[Subscript[SIZE, n][g_]:={g},{g,PT[[All,1]]},{n,numberofconstraints}];
Do[Subscript[PT, n]=PT[[n]],{n,numberofconstraints}];
Size=Reap[Do[Sow[Subscript[SIZES, n]=Subscript[SIZE, n][Subscript[PT, n][[1]]]],{n,numberofconstraints}]]//Last;
Size=Flatten[Size,1];
branchpowerv1=Join[Size,Sqrt[Pcv1],2];
branchpowerv1=Sort@branchpowerv1//.{a___,{x_,y1_},{x_,y2_},b___}:>{a,{x,y1+y2},b};
BPv1=branchpowerv1[[1+numberofinputs;;-(numberofoutputs+1)]];
branchpowerv2=Join[Size,Sqrt[Pcv2],2];
branchpowerv2=Sort@branchpowerv2//.{a___,{x_,y1_},{x_,y2_},b___}:>{a,{x,y1+y2},b};
BPv2=branchpowerv2[[1+numberofinputs;;-(numberofoutputs+1)]];
,Null];

(*ConnectionType*)
GPL=Table[Connection[[i]]->Subscript[G, i],{i,1,Pgear}];
BPL=Table[Connection[[i]]->Subscript[B, i-Pgear],{i,Pgear+1,numberofconstraints}];
ConnectionType=Flatten[Join[GPL,BPL],2];

(*EdgepStyle for genuine power*)
GM=Flatten[PT1,1];
Table[Subscript[GEN, n][g_,h_]:={g->h},{g,GM[[All,1]]},{h,{GM[[All,2]]}},{n,Pgear}];
Do[Subscript[GE, n]=GM[[n]],{n,Pgear}];
Geniunepowers=Reap[Do[Sow[Subscript[GENS, n]=Subscript[GEN, n][Subscript[GE, n][[1]],Subscript[GE, n][[2]]]],{n,Pgear}]]//Last;
Gencolour=Flatten[Geniunepowers,1];

Table[Subscript[COLR, n][g_]:=g->Red,{g,Gencolour[[1]]},{n,Pgear}];
Do[Subscript[Gencolour, n]=Gencolour[[n]],{n,Pgear}];
Colourstyle=Reap[Do[Sow[Subscript[COLRS, n]=Subscript[COLR, n][Subscript[Gencolour, n][[1]]]],{n,Pgear}]]//Last;
Colourstyle=Flatten[Colourstyle];

Table[Subscript[COLR, n][g_]:=g->Green,{g,Gencolour[[1]]},{n,Pgear}];
Do[Subscript[Gencolour, n]=Gencolour[[n]],{n,Pgear}];
ColourstyleV=Reap[Do[Sow[Subscript[COLRS, n]=Subscript[COLR, n][Subscript[Gencolour, n][[1]]]],{n,Pgear}]]//Last;
ColourstyleV=Flatten[ColourstyleV];
(**************************)
gearlinks=numberofbodies-numberofinputs-numberofoutputs-1;

BPB=Flatten[ConstantArray[0,{1,numberofbodies}],2];
Do[BPB=ReplacePart[BPB,BP[[i,1]]->BP[[i,2]]],{i,1,gearlinks}];

If[numberofcarriers==1,
BPBv=Flatten[ConstantArray[0,{1,numberofbodies}],2];
Do[BPBv=ReplacePart[BPBv,BPv[[i,1]]->BPv[[i,2]]],{i,1,gearlinks}];
bodies=Range[1,numberofbodies];
linkpowersn=(Abs[Flatten[Abs[Sqrt[Peb]],2]+Abs[BPB]])/(Total[Flatten[Abs[Sqrt[Peb]]]]/2);
linkpowersv=(Abs[Flatten[Abs[Pebv],2]+Abs[BPBv]])/(Total[Flatten[Abs[Pebv]]]/2);
(*linkpowersv=(Abs[Flatten[Abs[Sqrt[Pebv]],2]+Abs[BPBv]])/(Total[Flatten[Abs[Sqrt[Pebv]]]]/2);*)
matchn=Transpose@{bodies,linkpowersn}/.List[a_,b_]->Rule[a,b];
matchv=Transpose@{bodies,linkpowersv}/.List[a_,b_]->Rule[a,b];
,Null]

If[numberofcarriers==2,
BPBv1=Flatten[ConstantArray[0,{1,numberofbodies}],2];
Do[BPBv1=ReplacePart[BPBv1,BPv1[[i,1]]->BPv1[[i,2]]],{i,1,gearlinks}];
BPBv2=Flatten[ConstantArray[0,{1,numberofbodies}],2];
Do[BPBv2=ReplacePart[BPBv2,BPv2[[i,1]]->BPv2[[i,2]]],{i,1,gearlinks}];
bodies=Range[1,numberofbodies];
linkpowersn=(Abs[Flatten[Abs[Sqrt[Peb]],2]+Abs[BPB]])/(Total[Flatten[Abs[Sqrt[Peb]]]]/2);
matchn=Transpose@{bodies,linkpowersn}/.List[a_,b_]->Rule[a,b];
linkpowersv1=(Abs[Flatten[Abs[Sqrt[Pebv1]],2]+Abs[BPBv1]])/(Total[Flatten[Abs[Sqrt[Pebv1]]]]/2);
matchv1=Transpose@{bodies,linkpowersv1}/.List[a_,b_]->Rule[a,b];
linkpowersv2=(Abs[Flatten[Abs[Sqrt[Pebv2]],2]+Abs[BPBv2]])/(Total[Flatten[Abs[Sqrt[Pebv2]]]]/2);
matchv2=Transpose@{bodies,linkpowersv2}/.List[a_,b_]->Rule[a,b];
,Null]

(*arrow thicknesses, internal powers, and external powers*)
If[numberofinputs==1&&numberofoutputs==1,
normalsize=matchn/. {Subscript[T, numberofbodies]->1};
virtualsize=matchv/. {Subscript[T, numberofbodies]->1};
(*Thickness for the arrows*)
PcT=Pc/.{Subscript[T, numberofbodies]->1};
PcvT=Pcv/.{Subscript[T, numberofbodies]->1};
Athicknessn=(First@#->Directive[Thickness[Last@#]])&/@Transpose@{Connection,Flatten[PcT]*0.08};
Athicknessv=(First@#->Directive[Thickness[Last@#]])&/@Transpose@{Connection,Flatten[PcvT]*0.08};
Percentagen=(First@#->Last@#)&/@Transpose@{bodies,linkpowersn/. {Subscript[T, numberofbodies]->1}};
Percentagev=(First@#->Last@#)&/@Transpose@{bodies,linkpowersv/. {Subscript[T, numberofbodies]->1}};
(*display of internal powers*)
internalpowersn=MapThread[{#1\[Colon]#2}(*[[#2]]*)&,{Connection,Flatten[Pc/. {Subscript[T, numberofbodies]->Tout}]}];
internalpowersn=DeleteCases[Flatten[internalpowersn],_\[Colon]0,{1}];
internalpowersv=MapThread[{#1\[Colon]#2}(*[[#2]]*)&,{Connection,Flatten[Pcv/. {Subscript[T, numberofbodies]->Tout}]}];
internalpowersv=DeleteCases[Flatten[internalpowersv],_\[Colon]0,{1}];
(*display of external powers*)
externalpowersn=(First@#\[Colon]Last@#)&/@Transpose@{bodies,Flatten[Peb]/. {Subscript[T, numberofbodies]->Tout}};
externalpowersn=DeleteCases[externalpowersn,_\[Colon]0];
externalpowersv=(First@#\[Colon]Last@#)&/@Transpose@{bodies,Flatten[Pebv]/. {Subscript[T, numberofbodies]->Tout}};
externalpowersv=DeleteCases[externalpowersv,_\[Colon]0];
,Null];

If[numberofinputs==2&&numberofoutputs==1,
normalsize=matchn/. {Subscript[T, numberofbodies]->1,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]};
virtualsize=matchv/. {Subscript[T, numberofbodies]->1,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]};
(*Thickness for the arrows*)
PcT=Pc/.{Subscript[T, numberofbodies]->1,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]};
PcvT=Pcv/.{Subscript[T, numberofbodies]->1,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]};
Athicknessn=(First@#->Directive[Thickness[Last@#]])&/@Transpose@{Connection,Flatten[PcT]*0.08};
Athicknessv=(First@#->Directive[Thickness[Last@#]])&/@Transpose@{Connection,Flatten[PcvT]*0.08};
Percentagen=(First@#->Last@#)&/@Transpose@{bodies,linkpowersn/. {Subscript[T, numberofbodies]->1,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]}};
Percentagev=(First@#->Last@#)&/@Transpose@{bodies,linkpowersv/. {Subscript[T, numberofbodies]->1,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]}};
(*display of internal powers*)
internalpowersn=MapThread[{#1\[Colon]#2}(*[[#2]]*)&,{Connection,Flatten[Pc/. {Subscript[T, numberofbodies]->Tout,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]}]}];
internalpowersn=DeleteCases[Flatten[internalpowersn],_\[Colon]0,{1}];
internalpowersv=MapThread[{#1\[Colon]#2}(*[[#2]]*)&,{Connection,Flatten[Pcv/. {Subscript[T, numberofbodies]->Tout,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]}]}];
internalpowersv=DeleteCases[Flatten[internalpowersv],_\[Colon]0,{1}];
(*display of external powers*)
externalpowersn=(First@#\[Colon]Last@#)&/@Transpose@{bodies,Flatten[Peb]/. {Subscript[T, numberofbodies]->Tout,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]}};
externalpowersn=DeleteCases[externalpowersn,_\[Colon]0];
externalpowersv=(First@#\[Colon]Last@#)&/@Transpose@{bodies,Flatten[Pebv]/. {Subscript[T, numberofbodies]->Tout,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]}};
externalpowersv=DeleteCases[externalpowersv,_\[Colon]0];
,Null];

If[numberofinputs==2&&numberofoutputs==2,
normalsize=matchn/. {Subscript[T, numberofbodies]->1,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]};
virtualsize1=matchv1/. {Subscript[T, numberofbodies]->1,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]};
virtualsize2=matchv2/. {Subscript[T, numberofbodies]->1,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]};
(*Thickness for the arrows*)
PcT=Pc/.{Subscript[T, numberofbodies]->1,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]};
PcvT1=Pcv1/.{Subscript[T, numberofbodies]->1,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]};
PcvT2=Pcv2/.{Subscript[T, numberofbodies]->1,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]};
Athicknessn=(First@#->Directive[Thickness[Last@#]])&/@Transpose@{Connection,Flatten[PcT]*0.08};
Athicknessv1=(First@#->Directive[Thickness[Last@#]])&/@Transpose@{Connection,Flatten[PcvT1]*0.08};
Athicknessv2=(First@#->Directive[Thickness[Last@#]])&/@Transpose@{Connection,Flatten[PcvT2]*0.08};
Percentagen=(First@#->Last@#)&/@Transpose@{bodies,linkpowersn/. {Subscript[T, numberofbodies]->1,Subscript[T, numberofbodies-1]->k,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]}};
Percentagev1=(First@#->Last@#)&/@Transpose@{bodies,linkpowersv1/. {Subscript[T, numberofbodies]->1,Subscript[T, numberofbodies-1]->k,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]}};
Percentagev2=(First@#->Last@#)&/@Transpose@{bodies,linkpowersv2/. {Subscript[T, numberofbodies]->1,Subscript[T, numberofbodies-1]->k,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]}};
(*display of internal powers*)
internalpowersn=MapThread[{#1\[Colon]#2}(*[[#2]]*)&,{Connection,Flatten[Pc/. {Subscript[T, numberofbodies]->Tout,Subscript[T, numberofbodies-1]->k,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]}]}];
internalpowersn=DeleteCases[Flatten[internalpowersn],_\[Colon]0,{1}];
internalpowersv1=MapThread[{#1\[Colon]#2}(*[[#2]]*)&,{Connection,Flatten[Pcv1/. {Subscript[T, numberofbodies]->Tout,Subscript[T, numberofbodies-1]->k,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]}]}];
internalpowersv1=DeleteCases[Flatten[internalpowersv1],_\[Colon]0,{1}];
internalpowersv2=MapThread[{#1\[Colon]#2}(*[[#2]]*)&,{Connection,Flatten[Pcv2/. {Subscript[T, numberofbodies]->Tout,Subscript[T, numberofbodies-1]->k,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]}]}];
internalpowersv2=DeleteCases[Flatten[internalpowersv2],_\[Colon]0,{1}];
(*display of external powers*)
externalpowersn=(First@#\[Colon]Last@#)&/@Transpose@{bodies,Flatten[Peb]/. {Subscript[T, numberofbodies]->Tout,Subscript[T, numberofbodies-1]->k,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]}};
externalpowersn=DeleteCases[externalpowersn,_\[Colon]0];
externalpowersv1=(First@#\[Colon]Last@#)&/@Transpose@{bodies,Flatten[Pebv1]/. {Subscript[T, numberofbodies]->Tout,Subscript[T, numberofbodies-1]->k,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]}};
externalpowersv1=DeleteCases[externalpowersv1,_\[Colon]0];
externalpowersv2=(First@#\[Colon]Last@#)&/@Transpose@{bodies,Flatten[Pebv2]/. {Subscript[T, numberofbodies]->Tout,Subscript[T, numberofbodies-1]->k,Subscript[w, bodynumberofinput1]->Cos[\[Theta]],Subscript[w, bodynumberofinput2]->Sin[\[Theta]]}};
externalpowersv2=DeleteCases[externalpowersv2,_\[Colon]0];
,Null];

(* Section 6 Results display *)
CreateDialog[TabView[{
"Constraint matrix"->
Grid[{{"Constraint matrix Cq = ",MatrixForm[Cq]}}],
"Normal Power"->
Grid[{{"Internal Powers, Pc = ",MatrixForm[Pc],"; External Powers, Peb = ",MatrixForm[Simplify[Peb]]}}],
"Virtual Power"->Grid[{{If[numberofcarriers==1,Grid[{{"Virtual Internal Powers, Pcv = ",MatrixForm[Pcv],"; Virtual External Powers, Pebv = ",MatrixForm[Pebv]}}]]},{If[numberofcarriers==2,Grid[{{"Virtual Internal Powers, Pcv1 = ",MatrixForm[Pcv1],"; Virtual External Powers, Pebv1 = ",MatrixForm[Pebv1],"; Virtual Internal Powers, Pcv2 = ",MatrixForm[Pcv2],"Virtual External Powers, Pebv2 = ",MatrixForm[Pebv2]}}]]}}],
"Efficiency expression and plots"->Panel[Column[{Grid[{{"The efficiency expression: \[Eta] = "},{Simplify[eff]}}],
(*Grid[{{"The efficiency expression: \[Eta] = "},{Simplify[effb]}}],*)Grid[{{If[numberofinputs==1,With[{expr=eff},Manipulate[Grid[{{"Efficiency"},{expr}}],{Subscript[\[Mu], 1],0,1,0.001},{Subscript[\[Mu], 2],0,1,0.001}]]]},
{If[numberofinputs==2&&numberofoutputs==1,
With[{expr=eff},Manipulate[Plot[expr,{\[Theta],0, 2Pi}, PlotRange->All],{Subscript[\[Mu], 1],0,1,0.001},{Subscript[\[Mu], 2],0,1,0.001},{Subscript[\[Mu], 3],0,1,0.001}]]]},
{If[numberofinputs==2&&numberofoutputs==2,
With[{expr=eff},Manipulate[Plot[expr,{\[Theta],0, 2Pi}, PlotRange->All],{Subscript[\[Mu], 1],0,1,0.001},{Subscript[\[Mu], 2],0,1,0.001},{Subscript[\[Mu], 3],0,1,0.001},{Subscript[\[Mu], 4],0,1,0.001},{Subscript[\[Mu], 5],0,1,0.001},{Subscript[\[Mu], 6],0,1,0.001},{Subscript[\[Mu], 7],0,1,0.001},{Subscript[\[Mu], 8],0,1,0.001},{Subscript[\[Mu], 9],0,1,0.001},{Subscript[\[Mu], 10],0,1,0.001},{k,1,10}]]]}}]}]],
"Power flow diagrams"->
If[numberofinputs==1,With[{bodynumberofpc1=bodynumberofpc1,bodynumberofground=bodynumberofground,bodynumberofinput1=bodynumberofinput1,normalsize=normalsize,virtualsize=virtualsize,Connection=Connection,ConnectionType=ConnectionType,Colourstyle=Colourstyle,ColourstyleV=ColourstyleV,(*Athicknessn=Athicknessn,Athicknessv=Athicknessv,*)Percentagen=Percentagen,Percentagev=Percentagev,internalpowersn=internalpowersn,internalpowersv=internalpowersv,externalpowersn=externalpowersn,externalpowersv=externalpowersv},Manipulate[Grid[{{"Normal power flow diagram",Graph[Connection,VertexLabels->Placed["Name",Center],VertexLabelStyle->Directive[Black,20],VertexSize->Large(*,VertexSize\[Rule]normalsize*),VertexStyle->(Sequence@@{#->RGBColor@{#/.Percentagen,0.1,0.1,0.7}}&/@bodies),Properties->{bodynumberofground->{VertexShape->Graphics[Circle[{0, 0}], ImageSize -> {24., Automatic}]}},EdgeLabels->ConnectionType,EdgeStyle->(*Athicknessn*)Colourstyle,(*EdgeShapeFunction\[Rule]({Arrowheads[{{.03,.95}}],Arrow[#1,{.02 #2[[1]],.01#2[[2]]}]}&),*)(*EdgeShapeFunction\[Rule]({Arrowheads[{{.1,1}}],Arrow[#1,{.05 #2[[1]]/.normalsize,.05 #2[[2]]/.normalsize}]}&),*)ImageSize->{200,200}],Grid[{{"Normal internal powers"},{NumberForm[Column[internalpowersn],2]}},Frame->All],Grid[{{"Normal external powers"},{NumberForm[Column[externalpowersn],2]}},Frame->All]},{"Virtual power flow diagram",Graph[Connection,VertexLabels->Placed["Name",Center],VertexLabelStyle->Directive[Black,20],VertexSize->Large(*,VertexSize\[Rule]virtualsize*),VertexStyle->(Sequence@@{#->RGBColor@{#/.Percentagev,0.1,0.1,0.7}}&/@bodies),Properties->{bodynumberofpc1->{VertexShape->Graphics[Circle[{0, 0}], ImageSize -> {24., Automatic}]}},EdgeLabels->ConnectionType,EdgeStyle->(*Athicknessv*)ColourstyleV,(*EdgeShapeFunction\[Rule]({Arrowheads[{{.03,.95}}],Arrow[#1,{.02 #2[[1]],.01#2[[2]]}]}&),*)(*EdgeShapeFunction\[Rule]({Arrowheads[{{.1,1}}],Arrow[#1,{.05 #2[[1]]/.virtualsize,.05 #2[[2]]/.virtualsize}]}&),*)ImageSize->{200,200}],Grid[{{"Virtual internal powers"},{NumberForm[Column[internalpowersv],2]}},Frame->All],Grid[{{"Virtual external powers"},{NumberForm[Column[externalpowersv],2]}},Frame->All]}}],{Subscript[w, 2],-1,1}]]
,If[numberofinputs==2&&numberofcarriers==1,
With[{bodynumberofpc1=bodynumberofpc1,bodynumberofground=bodynumberofground,bodynumberofinput1=bodynumberofinput1,bodynumberofinput2=bodynumberofinput2,normalsize=normalsize,virtualsize=virtualsize,Connection=Connection,ConnectionType=ConnectionType,Colourstyle=Colourstyle,ColourstyleV=ColourstyleV,(*Athicknessn=Athicknessn,Athicknessv=Athicknessv,*)Percentagen=Percentagen,Percentagev=Percentagev,internalpowersn=internalpowersn,internalpowersv=internalpowersv,externalpowersn=externalpowersn,externalpowersv=externalpowersv},Manipulate[Grid[{{"Normal power flow diagram",Graph[Connection,VertexLabels->Placed["Name",Center],VertexLabelStyle->Directive[Black,20],VertexSize->Large(*,VertexSize\[Rule]normalsize*),VertexStyle->(Sequence@@{#->RGBColor@{#/.Percentagen,0.1,0.1,0.7}}&/@bodies),Properties->{bodynumberofground->{VertexShape->Graphics[Circle[{0, 0}], ImageSize -> {24., Automatic}]}},EdgeLabels->ConnectionType,EdgeStyle->(*Athicknessn*)Colourstyle,(*EdgeShapeFunction\[Rule]({Arrowheads[{{.03,.95}}],Arrow[#1,{.02 #2[[1]],.01#2[[2]]}]}&),*)(*EdgeShapeFunction\[Rule]({Arrowheads[{{.1,1}}],Arrow[#1,{.05 #2[[1]]/.normalsize,.05 #2[[2]]/.normalsize}]}&),*)ImageSize->{200,200}],Grid[{{"Normal internal powers"},{NumberForm[Column[internalpowersn],2]}}],Grid[{{"Normal external powers"},{NumberForm[Column[externalpowersn],2]}}]},{"Virtual power flow diagram",Graph[Connection,VertexLabels->Placed["Name",Center],VertexLabelStyle->Directive[Black,20],VertexSize->Large(*,VertexSize\[Rule]virtualsize*),VertexStyle->(Sequence@@{#->RGBColor@{#/.Percentagev,0.1,0.1,0.7}}&/@bodies),Properties->{bodynumberofpc1->{VertexShape->Graphics[Circle[{0, 0}], ImageSize -> {24., Automatic}]}},EdgeLabels->ConnectionType,EdgeStyle->(*Athicknessv*)ColourstyleV,(*EdgeShapeFunction\[Rule]({Arrowheads[{{.03,.95}}],Arrow[#1,{.02 #2[[1]],.01#2[[2]]}]}&),*)(*EdgeShapeFunction\[Rule]({Arrowheads[{{.1,1}}],Arrow[#1,{.05 #2[[1]]/.virtualsize,.05 #2[[2]]/.virtualsize}]}&),*)ImageSize->{200,200}],Grid[{{"Virtual internal powers"},{NumberForm[Column[internalpowersv],2]}}],Grid[{{"Virtual external powers"},{NumberForm[Column[externalpowersv],2]}}]}}],{\[Theta],0,2 Pi}]]
,If[numberofinputs==2&&numberofcarriers==2,
With[{bodies=bodies,bodynumberofpc1=bodynumberofpc1,bodynumberofpc2=bodynumberofpc2,bodynumberofground=bodynumberofground,bodynumberofinput1=bodynumberofinput1,bodynumberofinput2=bodynumberofinput2,normalsize=normalsize,virtualsize1=virtualsize1,virtualsize2=virtualsize2,Connection=Connection,ConnectionType=ConnectionType,Colourstyle=Colourstyle,ColourstyleV=ColourstyleV,(*Athicknessn=Athicknessn,Athicknessv1=Athicknessv1,Athicknessv2=Athicknessv2,*)Percentagen=Percentagen,Percentagev1=Percentagev1,Percentagev2=Percentagev2,internalpowersn=internalpowersn,internalpowersv1=internalpowersv1,internalpowersv2=internalpowersv2,externalpowersn=externalpowersn,externalpowersv1=externalpowersv1,externalpowersv2=externalpowersv2},Manipulate[Grid[{{"Normal power flow diagram",Graph[Connection,VertexLabels->Placed["Name",Center],Properties->{bodynumberofground->{VertexShape->Graphics[Circle[{0, 0}], ImageSize -> {24., Automatic}]}},EdgeStyle->Colourstyle,(*EdgeShapeFunction\[Rule]({Arrowheads[{{.03,.95}}],Arrow[#1,{.02 #2[[1]],.01#2[[2]]}]}&),*)EdgeLabels->ConnectionType,VertexLabelStyle->Directive[Black,15],VertexSize->Large,VertexStyle->(Sequence@@{#->CMYKColor@{#/.Percentagen,0.01,0.01,.1,0.5}(*0.1,0.8,0.7*)}&/@bodies),ImageSize->{200,200}],Grid[{{"Normal internal powers"},{NumberForm[Column[internalpowersn],2]}}],Grid[{{"Normal external powers"},{NumberForm[Column[externalpowersn],2]}}]},{"Virtual power flow diagram 1",Graph[Connection,VertexLabels->Placed["Name",Center],Properties->{bodynumberofpc1->{VertexShape->Graphics[Circle[{0, 0}], ImageSize -> {24., Automatic}]}},EdgeStyle->ColourstyleV,(*EdgeShapeFunction\[Rule]({Arrowheads[{{.03,.95}}],Arrow[#1,{.02 #2[[1]],.01#2[[2]]}]}&),*)EdgeLabels->ConnectionType,VertexLabelStyle->Directive[Black,15],VertexSize->Large,VertexStyle->(Sequence@@{#->CMYKColor@{#/.Percentagev1,0.01,0.01,.1,0.5}}&/@bodies),ImageSize->{200,200}],Grid[{{"Virtual internal powers 1"},{NumberForm[Column[internalpowersv1],2]}}],Grid[{{"Virtual external powers 1"},{NumberForm[Column[externalpowersv1],2]}}]},
{"Virtual power flow diagram 2",Graph[Connection,VertexLabels->Placed["Name",Center],Properties->{bodynumberofpc2->{VertexShape->Graphics[Circle[{0, 0}], ImageSize -> {24., Automatic}]}},EdgeStyle->ColourstyleV,(*EdgeShapeFunction\[Rule]({Arrowheads[{{.03,.95}}],Arrow[#1,{.02 #2[[1]],.01#2[[2]]}]}&),*)EdgeLabels->ConnectionType,VertexLabelStyle->Directive[Black,15],VertexSize->Large,VertexStyle->(Sequence@@{#->CMYKColor@{#/.Percentagev2,0.01,0.01,.1,0.5}}&/@bodies),ImageSize->{200,200}],Grid[{{"Virtual internal powers 2"},{NumberForm[Column[internalpowersv2],2]}}],Grid[{{"Virtual external powers 2"},{NumberForm[Column[externalpowersv2],2]}}]}}],{\[Theta],0,2 Pi},{k,0,10}]],Null]]]
}]]
]
