% Solution of a resistive network with constant sources.

% Element types are as follows:
%
%    A: Voltage Controlled Voltage Source
%    B: Voltage Controlled Current Source
%    D: Current Controlled Voltage Source
%    E: Current Controlled Current Source
%    G: Conductance
%    I: Independent Current Source
%    O: Open Circuit
%    R: Resistance
%    S: Short Circuit
%    V: Independent Voltage Source

% Element data depends on the type of element.
%
%    Independent Sources, Resistors and Conductances
%        Reference Designator, Node A, Node B, Value
%
%    Dependent Sources
%        Reference Designator, Node A, Node B, Value, Control Element
%
%    Short and Open Circuits
%        Reference Designator, Node A, Node B

% Circuit nodes are assumed to be numbered sequentially starting with Node 1.

% Reference Designators and Control Elements may be up to seven characters in length.

% Voltages are specified as measured at Node A with respect to Node B.

% Currents are specified as flowing from Node A to Node B.

% All node equations assume the direction of positive current flow is exiting the node.



% Initialization
clear all
close all

m = 1;
sources = 0;
sources1 = 0;
conductances = 0;
shorts = 0;
opens = 0;

Type = '';
Designator = '';
Control = '';

Value = [];
NodeA = [];
NodeB = [];



% Input element data
fprintf ('\nEnter Network Elements\n')

while (1 == 1)

  fprintf ('Element%3d:  ', m)
  String = input ('', 's');
  if (isempty (String))
    break
  end
  
  % Parse the element data string
  String(isspace (String)) = [];
  String = [String ','];
  if (String(1) == 'S' || String(1) == 'O')
    String = [String '0,'];
  end
  Commas = strfind (String, ',');
  
  % Determine element type
  if (String(1:4) == 'VCVS')
    Type(m) = 'A';
  elseif (String(1:4) == 'VCCS')
    Type(m) = 'B';
  elseif (String(1:4) == 'CCVS')
    Type(m) = 'D';
  elseif (String(1:4) == 'CCCS')
    Type(m) = 'E';
  else
    Type(m) = String(1);
  end

  % Determine common element parameters
  if (length (Commas) >= 4)
    Designator(m,1:7) = [String(1:Commas(1)-1) blanks(8-Commas(1))];
    NodeA(m) = str2num (String(Commas(1)+1:Commas(2)-1));
    NodeB(m) = str2num (String(Commas(2)+1:Commas(3)-1));
    Value(m) = str2num (String(Commas(3)+1:Commas(4)-1));
  else
    Type(m) = 'X';
  end
  
  % Determine element specific parameters and increment element counters
  Control(m,1:7) = blanks(7);
  
  switch (Type(m))
    
    % Dependent source
    case {'A', 'B', 'D', 'E'}
      if (length (Commas) == 5)
        Control(m,1:Commas(5)-Commas(4)-1) = String(Commas(4)+1:Commas(5)-1);
        sources = sources+1;
        if (Type(m) == 'D' || Type(m) == 'E')
          sources1 = sources1+1;
        end
      else
        Type(m) = 'X';
      end

    % Conductance or resistor
    case {'G', 'R'}
      conductances = conductances+1;
      
    % Independent source
    case {'I', 'V'}
      sources = sources+1;

    % Open circuit
    case {'O'}
      opens = opens+1;
      
    % Short circuit
    case {'S'}
      shorts = shorts+1;
      
    % Invalid element
    otherwise
      fprintf ('Invalid Element\n')
      m = m-1;
      
  end

  m = m+1;

end



% Set up the network coefficient matrix A and source vector B
elements = m-1;
nodes = max (max (NodeA), max (NodeB));
n = nodes + sources + sources1 + shorts;
A = zeros (n, n);
B = zeros (n, 1);
X = zeros (n, 1);
k = 0;

for m = 1:elements

  % Locate the control element for a dependent source
  if (ismember (Type(m), 'ABDE'))
    k1 = 0;
    for m1 = 1:elements
      if (Type(m1) == 'D' || Type(m1) == 'E')
        k1 = k1+2;
      elseif (Type(m1) == 'G' || Type(m1) == 'O')
        G = Value(m1);
      elseif (Type(m1) == 'R')
        G = 1/Value(m1);
      else
        k1 = k1+1;
      end
      if (Designator(m1,:) == Control(m,:))
        break
      end
    end
  end    

  % All sources and short circuits
  if (ismember (Type(m), 'ABDEISV'))
    k = k+1;
    if (NodeA(m) ~= 0)
      A(NodeA(m), nodes+k) = 1;
    end
    if (NodeB(m) ~= 0)
      A(NodeB(m), nodes+k) = -1;
    end
  end

  switch (Type(m))

    % Voltage controlled voltage source
    case {'A'}
      if (NodeA(m) ~= 0)
        A(nodes+k, NodeA(m)) = 1;
      end
      if (NodeB(m) ~= 0)
        A(nodes+k, NodeB(m)) = -1;
      end
      if (NodeA(m1) ~= 0)
        A(nodes+k, NodeA(m1)) = A(nodes+k, NodeA(m1)) - Value(m);
      end
      if (NodeB(m1) ~= 0)
        A(nodes+k, NodeB(m1)) = A(nodes+k, NodeB(m1)) + Value(m);
      end

    % Voltage controlled current source
    case {'B'}
      if (NodeA(m1) ~= 0)
        A(nodes+k, NodeA(m1)) = -Value(m);
      end
      if (NodeB(m1) ~= 0)
        A(nodes+k, NodeB(m1)) = Value(m);
      end
      A(nodes+k, nodes+k) = 1;

    % Current controlled voltage source
    case {'D'}
      if (NodeA(m) ~= 0)
        A(nodes+k, NodeA(m)) = 1;
      end
      if (NodeB(m) ~= 0)
        A(nodes+k, NodeB(m)) = -1;
      end
      A(nodes+k, nodes+k+1) = -Value(m);

    % Current controlled current source
    case {'E'}
      A(nodes+k, nodes+k) = 1;
      A(nodes+k, nodes+k+1) = -Value(m);
    
    % Conductance or resistor
    case {'G', 'R'}
      if (Type(m) == 'G')
        G = Value(m);
      else
        G = 1/Value(m);
      end
      if (NodeA(m) ~=0)
        A(NodeA(m), NodeA(m)) = A(NodeA(m), NodeA(m)) + G;
        if (NodeB(m) ~= 0)
          A(NodeA(m), NodeB(m)) = A(NodeA(m), NodeB(m)) - G;
        end
      end
      if (NodeB(m) ~= 0)
        A(NodeB(m), NodeB(m)) = A(NodeB(m), NodeB(m)) + G;
        if (NodeA(m) ~= 0)
          A(NodeB(m), NodeA(m)) = A(NodeB(m), NodeA(m)) - G;
        end
      end

    % Current source
    case {'I'}
      A(nodes+k, nodes+k) = 1;
      B(nodes+k,1) = Value(m);

    % Voltage source or short circuit
    case {'V', 'S'}
      if (NodeA(m) ~= 0)
        A(nodes+k, NodeA(m)) = 1;
      end
      if (NodeB(m) ~= 0)
        A(nodes+k, NodeB(m)) = -1;
      end
      if (Type(m) == 'V')
        B(nodes+k,1) = Value(m);
      else
        B(nodes+k,1) = 0;
      end

  end

  % Dependent current sources
  if (Type(m) == 'D' || Type(m) == 'E')
    k = k+1;
    if (Type(m1) == 'G' || Type(m1) == 'R' || Type(m1) == 'O')
      if (NodeA(m1) ~= 0)
        A(nodes+k, NodeA(m1)) = -G;
      end
      if (NodeB(m1) ~= 0)
        A(nodes+k, NodeB(m1)) = G;
      end
    else
      A(nodes+k, nodes+k1) = -1;
    end
    A(nodes+k, nodes+k) = 1;
  end

end



% Solve the network equations for the steady-state
X(:,1) = A \ B(:,1);



% Display steady-state conditions
fprintf ('\n\nNODE VOLTAGES\n')

% Display node voltages
fprintf ('Node           Voltage\n')
fprintf ('   0         Reference\n')
for m = 1:nodes
  fprintf ('%4d%18s\n', m, [num2metric(X(m)) 'V'])
end  

% Display sources data
if (sources > 0)
  k = 0;
  Psupplied = 0;
  fprintf ('\nSOURCES')
  fprintf ('\nElement          Value   NodeA   NodeB   Control        Voltage        Current     P supplied\n')
  
  for m = 1:elements
    if (ismember (Type(m), 'ABDEISV'))
      k = k+1;
    end
    
    if (ismember (Type(m), 'ABDEIV'))
      if (NodeA(m) == 0)
        VAB = -X(NodeB(m),1);
      elseif (NodeB(m) == 0)
        VAB = X(NodeA(m),1);
      else
        VAB = X(NodeA(m),1) - X(NodeB(m),1);
      end
      
      IAB = X(nodes+k,1);
      
      if (Type(m) == 'B')
	    String = 'S';
	  elseif (Type(m) == 'D')
	    String = 'O';
      elseif (Type(m) == 'I')
        String = 'A';
      elseif (Type(m) == 'V')
        String = 'V';
      else
        String = ' ';
      end
      
      fprintf ('%7s%15s%8d%8d   %7s%15s%15s%15s\n', Designator(m,:), [num2metric(Value(m)) String], NodeA(m), NodeB(m), ...
              Control(m,:), [num2metric(VAB) 'V'], [num2metric(IAB) 'A'], [num2metric(-VAB*IAB) 'W'])
      
      Psupplied = Psupplied - VAB*IAB;
    end

    if (Type(m) == 'D' || Type(m) == 'E')
      k = k+1;
    end
  end

  fprintf ('                                                                                 ------------')
  fprintf ('\n                                                                              %15s\n', ...
          [num2metric(Psupplied) 'W'])
end

% Display passive elements data
if (conductances > 0)
  k = 0;
  Pabsorbed = 0;
  fprintf ('\nPASSIVE ELEMENTS')
  fprintf ('\nElement          Value   NodeA   NodeB                  Voltage        Current     P absorbed\n')

  for m = 1:elements
    if (ismember (Type(m), 'ABDEISV'))
      k = k+1;
    end
    
    if (Type(m) == 'G' || Type(m) == 'R')
      if (NodeA(m) == 0)
        VAB = -X(NodeB(m),1);
      elseif (NodeB(m) == 0)
        VAB = X(NodeA(m),1);
      else
        VAB = X(NodeA(m),1) - X(NodeB(m),1);
      end
      
      if (Type(m) == 'G')
        IAB = VAB * Value(m);
        fprintf ('%7s%15s%8d%8d   %7s%15s%15s%15s\n', Designator(m,:), [num2metric(Value(m)) 'S'], NodeA(m), NodeB(m), ...
                Control(m,:), [num2metric(VAB) 'V'], [num2metric(IAB) 'A'], [num2metric(VAB*IAB) 'W'])
      elseif (Type(m) == 'R')
        IAB = VAB / Value(m);
        fprintf ('%7s%15s%8d%8d   %7s%15s%15s%15s\n', Designator(m,:), [num2metric(Value(m)) 'O'], NodeA(m), NodeB(m), ...
                Control(m,:), [num2metric(VAB) 'V'], [num2metric(IAB) 'A'], [num2metric(VAB*IAB) 'W'])
      end
      
      Pabsorbed = Pabsorbed + VAB*IAB;
    end
    
    if (Type(m) == 'D' || Type(m) == 'E')
      k = k+1;
    end
  end
  
  fprintf ('                                                                                 ------------')
  fprintf ('\n                                                                              %15s\n', ...
          [num2metric(Pabsorbed), 'W'])
end

% Display shorts and opens data
if (shorts+opens > 0)
  k = 0;
  fprintf ('\nSHORTS and OPENS')
  fprintf ('\nElement                  NodeA   NodeB                  Voltage        Current\n')

  for m = 1:elements
    if (ismember (Type(m), 'ABDEISV'))
      k = k+1;
    end
    
    if (Type(m) == 'S')
      IAB = X(nodes+k,1);
      fprintf ('%7s               %8d%8d                         %15s\n', Designator(m,:), NodeA(m), NodeB(m), [num2metric(IAB) 'A'])
    end
    
    if (Type(m) == 'O')
      if (NodeA(m) == 0)
        VAB = -X(NodeB(m),1);
      elseif (NodeB(m) == 0)
        VAB = X(NodeA(m),1);
      else
        VAB = X(NodeA(m),1) - X(NodeB(m),1);
      end    
      fprintf ('%7s               %8d%8d          %15s\n', Designator(m,:), NodeA(m), NodeB(m), [num2metric(VAB) 'V'])
    end
    
    if (Type(m) == 'D' || Type(m) == 'E')
      k = k+1;
    end
  end

end

fprintf ('\n\n')