function [f,g,h,e] = mknw(desc,mode)

% MKNW	  Build an electrical network representation
%   [F,G,H,E] = MKNW(DESC) The electrical network is described in its natural
%   variables: the voltages at its nodes and the currents through the
%   components. If x is the column vector formed with these voltages and these
%   currents, u the inputs (voltage or current sources) and y the outputs, then
%   the equations describing the circuit are:
%       Fx + G sx = Hu
%               y = Ex
%   [A,B,C,D] = MKNW(DESC,'SS') returns the state-space representation:
%      sx = Ax + Bu
%       y = Cx + Du
%
% DESC is a matrix whose raws have the following form :
% {'$' node1 node2 node3 value}
%  'r' m n 0 value  resistance between node m and n
%  'c' m n 0 value  capacitance between node m and n
%  'l' m n 0 value  inductance between node m and n
%  'm' 0 0 0 value  mutual between the two previous 'l'
%  'k' 0 0 0 value  coupling factor between the two previous 'l'
%  '-' m n 0 0      connection between node m and n
%  'u' m n 0 0      voltage source between node m and n
%  'i' m n 0 0      current source between node m and n
%  'g' m 0 0 0      node m is grounded
%  'a' m n l gain   opamp, + at m, - at n, output at l, ideal if gain = Inf
%  'o' m n 0 0      voltage output between node m and n
%  'p' 0 0 0 0      current output on the previous component
%
%   Mutual is limited to the coupling between two inductances. 'p' after 'a' it
%   gives the output current. Except for 'm','k' and 'p', the declaration order
%   is free. The input/output numbering is given by the declaration order.
%
% [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

desc = [abs(char(desc(:,1))) cell2mat(desc(:,2:end))];
nc =  sum(desc(:,1) == 'r') + ...
      sum(desc(:,1) == 'c') + ...
      sum(desc(:,1) == 'l') + ...
      sum(desc(:,1) == '-') + ...
      sum(desc(:,1) == 'a');
nu =  sum(desc(:,1) == 'u');
ni =  sum(desc(:,1) == 'i');
ig =  desc(desc(:,1) == 'g',2);
ng =  length(ig);
nO =  sum(desc(:,1) == 'o') + ...
      sum(desc(:,1) == 'p');

nV = max(max(desc(:,2:4)));
nC = nc + nu + ni;
nI = nu + ni;

fV = zeros(nC,nV);
fC = zeros(nC,nC);
gV = zeros(nC,nV);
gC = zeros(nC,nC);
k  = zeros(nV,nC);
h  = zeros(nC+nV,nI);
eV = zeros(nO,nV);
eC = zeros(nO,nC);

i  = 0;
iI = 0;
iO = 0;
for kd = 1:size(desc,1)
 d = desc(kd,:);
 switch char(d(1))
  case 'r' % resistance
   i = i + 1;
   fV(i,d(2)) = +1;
   fV(i,d(3)) = -1;
   fC(i,i)    = d(5);
   k (d(2),i) = +1;
   k (d(3),i) = -1;
  case 'c' % capacitance
   i = i + 1;
   gV(i,d(2)) = +d(5);
   gV(i,d(3)) = -d(5);
   fC(i,i)    = +1;
   k (d(2),i) = +1;
   k (d(3),i) = -1;
  case 'l' % inductance
   i = i + 1;
   fV(i,d(2)) = +1;
   fV(i,d(3)) = -1;
   gC(i,i)    = d(5);
   k (d(2),i) = +1;
   k (d(3),i) = -1;
  case 'm' % mutual inductance
   gC(i,i-1)  = d(5);
   gC(i-1,i)  = d(5);
  case 'k' % coupling
   gC(i,i-1)  = d(5)*sqrt(gC(i,i)*gC(i-1,i-1));
   gC(i-1,i)  = gC(i,i-1);
  case '-' % simple connection
   i = i + 1;
   fV(i,d(2)) = +1;
   fV(i,d(3)) = -1;
   k (d(2),i) = +1;
   k (d(3),i) = -1;
  case 'a' % amplifier
   i = i + 1;
   if finite(d(5))
     fV(i,d(2)) = +d(5);
     fV(i,d(3)) = -d(5);
     fV(i,d(4)) = -1;
   else
     fV(i,d(2)) = +1;
     fV(i,d(3)) = -1;
   end
   k (d(4),i) =  1;
  case 'u' % voltage source
   i = i + 1;
   iI = iI + 1;
   fV(i,d(2)) = +1;
   fV(i,d(3)) = -1;
   h (i,iI)   = +1;
   k (d(2),i) = -1;
   k (d(3),i) = +1;
  case 'i' % current source
   i  = i  + 1;
   iI = iI + 1;
   fC(i,i)    = +1;
   h (i,iI)   = +1;
   k (d(2),i) = +1;
   k (d(3),i) = -1;
  case 'o' % output node
   iO = iO + 1;
   eV(iO,d(2)) = +1;
   eV(iO,d(3)) = -1;
  case 'p'
   iO = iO + 1;
   eC(iO,i   ) = +1;
 end
end

% remove Kirchhof law and add constraint for ground nodes
k(ig,:) = [];
gnd = zeros(ng,nV);
gnd(1:ng,ig) = eye(ng,ng);

% build the large matrices
f = [fV -fC;zeros(nV-ng,nV) k;gnd zeros(ng,nC)];
g = [gV -gC;zeros(nV,nV+nC)];
e = [eV eC];

% state-space representation
if nargin < 2 | ~strcmp(upper(mode),'SS'), return, end
[m,d,n] = svd(g);
i  = diag(d) > d(1,1)*length(d)*eps;
in = find(~i);
i  = find(i);
f  = m'*f*n;
h  = m'*h;
e  = e*n;
f1 = f( i, i);
f2 = f( i,in);
f3 = f(in, i);
f4 = f(in,in);
g  = diag(1./diag(d(i,i)));
h1 = h( i, :);
h2 = h(in, :);
e1 = e( :, i);
e2 = e( :,in);
[m,d,n] = svd(f4);
i  = find(diag(d) > d(1,1)*length(d)*eps);
f2 = f2*n;
f3 = m'*f3;
h2 = m'*h2;
e2 = e2*n;
f2 = f2(:,i);
f3 = f3(i,:);
f4 = diag(1./diag(d(i,i)));
h2 = h2(i,:);
e2 = e2(:,i);
[f,g,h,e] = minreal(g*(f2*f4*f3-f1),g*(h1-f2*f4*h2),e1-e2*f4*f3,e2*f4*h2);
 
