(*:Title: NodAl.m *)

(*:SYMAP Version 1.0 *)

(*:Copyright 1997 by Jian Sun *)

(*:All rights reserved *)

(*:History:
	Original Version by Jian Sun
*)

StateSpaceModel[NetList_List, InputVars_List, SwitchStates_List] := 
  Block[{CharNums,BranchCurrent, NodeVoltage, StateVars, Dim, ib, vn, x, 
      CircuitEqus, VarCor, i, j},
		CharNums = PreProcessing[NetList];
		Dim = Length[CharNums[[3]]] + Length[CharNums[[4]]];
		BranchCurrent = Table[ib[j][t], {j, 1, CharNums[[1]]}];
		NodeVoltage = Table[vn[j][t], {j, 1, CharNums[[2]]}];
		StateVars = Table[x[i][t], {i, 1, Dim}];
		VarCor = ReplaceRule[NetList, StateVars, BranchCurrent, NodeVoltage, 
                     CharNums[[3]], CharNums[[4]]];
		CircuitEqus = NodalAnalysis[NetList, BranchCurrent, NodeVoltage];
		FindAB[CircuitEqus, SwitchStates, VarCor,
		       StateVars,InputVars, BranchCurrent, NodeVoltage]];

PreProcessing[NetList_List] := Block[
	{BranchNum, NodeNum, IndBranch, CapNode, CompList, i, j},
	 
	 BranchNum = Length[NetList];
	 NodeNum = Max[Union[Table[NetList[[i, 2, 1]], 
	                           {i, 1, BranchNum}],
	                     Table[NetList[[i, 2, 2]], 
	                           {i, 1, BranchNum}]]];

	IndBranch = {};
	For[i=1, i<=BranchNum, i++,
		If[SameQ[NetList[[i, 1]], Inductor], 
	   	   IndBranch = Join[IndBranch, {i}]]];
	   	   
	CapNode = {};
	For[i=1, i<=BranchNum, i++,
		If[NetList[[i, 1]]==Capacitor, 
	   	   CapNode = Join[CapNode, {NetList[[i, 2]]}] ]];
	   	   
	CompList = {};
	For[i=1, i<=BranchNum, i++,
		If[!MemberQ[{Switch},
		           NetList[[i, 1]]],
		   CompList = Join[CompList, {NetList[[i, 3]]}]
		  ] ];
	   	   
	{BranchNum, NodeNum, IndBranch, CapNode, CompList}
	]

NodalAnalysis[NetList_List, BranchCurrent_List, 
              NodeVoltage_List] := Block[
	{KCLEqus, KVLEqus, OneEqu, BranchNum, NodeNum, 
	 ConnectBranch, t, i, u, j, j1, j2, Par, k, l},

	BranchNum = Length[NetList];
	NodeNum = Length[NodeVoltage];

	ConnectBranch = Table[0, {BranchNum}];

	KCLEqus = {};
	KVLEqus = {};
	
	For[j=1, j<=BranchNum, j++,
	    j1 = NetList[[j, 2, 1]];
	    j2 = NetList[[j, 2, 2]];
	    Par = NetList[[j, 3]];

		If[SameQ[NetList[[j, 1]], VoltageSource],
		   OneEqu = {If[j1==0, 
		                0, 
		                NodeVoltage[[j1]]] - 
		                If[j2==0, 0, 
		                          NodeVoltage[[j2]]
		                  ] == Par}];

		If[SameQ[NetList[[j, 1]], CurrentSource],
		   OneEqu = {BranchCurrent[[j]] == Par}];
		                  
		If[SameQ[NetList[[j, 1]], Capacitor],
		   OneEqu = {Par D[If[j1==0, 
		                      0, 
		                      NodeVoltage[[j1]]] - 
		                      If[j2==0, 
		                         0, 
		                         NodeVoltage[[j2]]],
		                   t] ==
		             BranchCurrent[[j]]}];

		If[SameQ[NetList[[j, 1]], Inductor],
		   OneEqu = {If[j1==0, 
		                0, 
		                NodeVoltage[[j1]]] - 
		                If[j2==0, 
		                   0, 
		                   NodeVoltage[[j2]]
		                  ] ==
		             Par D[BranchCurrent[[j]], t]}];

		If[SameQ[NetList[[j, 1]], Resistor],
		   OneEqu = {If[j1==0, 0, NodeVoltage[[j1]]] - 
		             If[j2==0, 0, NodeVoltage[[j2]]] ==
		             Par BranchCurrent[[j]]}];

		If[SameQ[NetList[[j, 1]], Switch],
		   OneEqu = {Par (If[j1==0, 0, NodeVoltage[[j1]]] - 
		                  If[j2==0, 0, NodeVoltage[[j2]]]) +
		             (1 - Par) BranchCurrent[[j]] == 0}];

		If[SameQ[NetList[[j, 1]], Diode],
		   OneEqu = {Par (If[j1==0, 0, NodeVoltage[[j1]]] - 
		                  If[j2==0, 0, NodeVoltage[[j2]]]) +
		             (1 - Par) BranchCurrent[[j]] == 0}];

		KVLEqus = Join[KVLEqus, OneEqu] ];

	For[j=1, j<=NodeNum, j++, 
	    For[k=1, k<=BranchNum, k++,
	    	ConnectBranch[[k]] = 0;
	    	If[NetList[[k, 2, 1]]==j, ConnectBranch[[k]]=-1];
	    	If[NetList[[k, 2, 2]]==j, ConnectBranch[[k]]=1]
	    	];

	    OneEqu = {Sum[BranchCurrent[[k]] ConnectBranch[[k]],
	                  {k, 1, BranchNum}] == 0};
	    KCLEqus = Join[KCLEqus, OneEqu] ];	
	
	Join[KCLEqus, KVLEqus]
	]


ReplaceRule[NetList_List, StateVars_List, BranchCurrent_List, 
            NodeVoltage_List, IndBranch_List,
            CapNode_List] := Block[{i, VarCor},

	VarCor =
	Join[Table[BranchCurrent[[IndBranch[[i]]]] ->
                   StateVars[[i]], 
                   {i, 1, Length[IndBranch]}],
             Table[If[CapNode[[i, 1]]==0, 0,
                      NodeVoltage[[CapNode[[i, 1]]]]] ->
                   If[CapNode[[i, 2]]==0, 0,
                      NodeVoltage[[CapNode[[i, 2]]]]] +
                      StateVars[[Length[IndBranch]+i]],
                   {i, 1, Length[CapNode]}] ];

	For[i=1, i<=Length[VarCor], i++,
	    If[VarCor[[i, 1]]==0, VarCor =
	       ReplacePart[VarCor, VarCor[[i, 2, 1]] -> -VarCor[[i, 2, 2]], i]]
	   ];
		
	StylePrint[VarCor, "Output"];
	VarCor
    ]

StateSpace[GenCircuitEqus_List, Topology_List,
           VarCor_List, StateVars_List,
           BranchCurrent_List, NodeVoltage_List] := Block[
	{CircuitEqus, IntermEqus, StateEqus, AnciEqus, i, j},
	
	CircuitEqus = Evaluate[GenCircuitEqus /. Topology];
	IntermEqus = CircuitEqus /. Join[VarCor, D[VarCor, t]];

	AnciEqus = AncillaryEqus[IntermEqus, StateVars,
	                         BranchCurrent, NodeVoltage];

	StateEqus = Eliminate[IntermEqus,
	                      Join[BranchCurrent, NodeVoltage]];

	StateEqus = DeTrivEqus[StateEqus, StateVars];
	AnciEqus = Join[AnciEqus,
			        AncillaryEqus[StateEqus, StateVars,
	                              BranchCurrent, NodeVoltage]];

	Join[StateEqus, AnciEqus, D[AnciEqus, t]]
	]		

AncillaryEqus[ISEqus_List, StateVars_List,
              BranchCurrent_List, NodeVoltage_List] := Block[
	{i, j, CircuitVars, AnciEqus, NoStateDeriv, NoCirVars},

	CircuitVars = Join[BranchCurrent, NodeVoltage];
	AnciEqus = {};
	For[i=1, i<=Length[ISEqus], i++,
		NoStateDeriv = !MemberQ[Join[
			  Table[MemberQ[ISEqus[[i, 1]],
		                    D[StateVars, t][[j]], 
		                    Infinity] ||
		            SameQ[ISEqus[[i, 1]],
		                  D[StateVars, t][[j]]],
		           {j, 1, Length[StateVars]}],

		      Table[MemberQ[ISEqus[[i, 2]],
		                    D[StateVars, t][[j]],
		                    Infinity] ||
		            SameQ[ISEqus[[i, 2]],
		                  D[StateVars, t][[j]]],
		           {j, 1, Length[StateVars]}]], True];

		NoCirVars = !MemberQ[Join[
			  Table[MemberQ[ISEqus[[i, 1]],
		                    CircuitVars[[j]], Infinity] ||
		            SameQ[ISEqus[[i, 1]], CircuitVars[[j]]],
		           {j, 1, Length[CircuitVars]}],

		      Table[MemberQ[ISEqus[[i, 2]],
		                    CircuitVars[[j]], Infinity] ||
		            SameQ[ISEqus[[i, 2]], CircuitVars[[j]]],
		           {j, 1, Length[CircuitVars]}]], True];
		
		If[And[NoStateDeriv, NoCirVars],
		   AnciEqus = Join[AnciEqus, {ISEqus[[i]]}]];
	  ];

	AnciEqus]		

DeTrivEqus[StateEqus_, StateVars_List] := Block[
	{i, j, NewStateEqus},

	NewStateEqus = {};
	For[i=1, i<=Length[StateEqus], i++,
		If[!SameQ[StateEqus[[i]], True],
	       NewStateEqus = Join[NewStateEqus, 
	                           {StateEqus[[i]]}]
	   ]];

	NewStateEqus
]

FindAB[CircuitEqus_List, SwitchStates_List,VarCor_List,
	StateVars_List,InputVars_List, BranchCurrent_List,
	NodeVoltage_List] := Block[{i, A, B, StateEqus, StandEqus},
	A = {};
	B = {};
	For[i=1, i<=Length[SwitchStates], i++,
          StateEqus = StateSpace[CircuitEqus, SwitchStates[[i]], 
                               VarCor, StateVars, BranchCurrent, 
                               NodeVoltage];
        StandEqus = StandardMatrixForm[StateEqus, StateVars];
        A = Join[A, {SystemMatrix[StandEqus, StateVars]}];
        B = Join[B, {InputMatrix[StandEqus,InputVars]}];
	];
	{A, B, StateVars}]

StandardMatrixForm[StateEqus_List, StateVars_List] := Block[
	{StandEqus, SoluOrder, ExtraCons, OneEqu,
	 NoDeriv, SingleIterm, i, j},
	StandEqus = Solve[StateEqus, D[StateVars, t]][[1]];
	
	ExtraCons = {};
	For[i=1, i<=Length[StateEqus], i++,
		OneEqu = Simplify[StateEqus[[i, 1]] -
		                  StateEqus[[i, 2]]];
		NoDeriv = !MemberQ[Table[MemberQ[OneEqu, 
		                                 D[StateVars, t][[j]],
		                                 Infinity],
		                         {j, 1, Length[StateVars]}],
		                   True];
		SingleTerm = Or[SameQ[Head[OneEqu], Times],
		                SameQ[Head[OneEqu], Symbol]];	

		If[NoDeriv && SingleTerm,
		   For[j=1, j<=Length[StateVars], j++,
		       If[MemberQ[OneEqu, StateVars[[j]]],
		          ExtraCons = Join[ExtraCons,
		                           {StateVars[[j]] -> 0}]
		         ];
		      ];
		   ];
	    ];

	StandEqus = Evaluate[StandEqus /. ExtraCons];

	SoluOrder = Table[Position[
	            		Table[StandEqus[[i, 1]],
	                  	      {i, 1, Length[StateVars]}],
	                    D[StateVars, t][[j]]][[1, 1]],
                      {j, 1, Length[StandEqus]}];
                      
   StandEqus = Table[StandEqus[[SoluOrder[[i]]]],
                     {i, 1, Length[StateVars]}];
   StandEqus
   ]

SystemMatrix[Equs_List, StateVars_List] := Block[
	{A, i, j},

	A = Table[D[Equs[[i, 2]], StateVars[[j]]],
                {i, 1, Length[Equs]}, 
                {j, 1, Length[StateVars]}];

   Simplify[A]
  ]

InputMatrix[Equs_List, InputVars_List] := Block[
	{B, i, j},

	B = Table[D[Equs[[i, 2]], InputVars[[j]]],
                {i, 1, Length[Equs]}, 
                {j, 1, Length[InputVars]}];

   Simplify[B]
  ]