Description of spice

0001 function [Vn, In, nn] = spice_eit(netlist, freq)
0002 % Vn = spice_eit(netlist, [freq])
0003 % Modified Nodal Analysis (MNA) circuit solver
0004 % for (almost/simple) spice netlists
0005 % input
0006 %   netlist - by cell array { [ part nodes...  args ], ... }
0007 %             where parts are spice types:
0008 %               "Rn" resistor; 2 nodes; args(1) = impedance
0009 %               "Ln" inductor; 2 nodes; args(1) = inductance
0010 %               "Cn" capacitor; 2 nodes; args(1) = capacitance
0011 %               "Vn" voltage source; 2 nodes (-,+); args(1) = voltage
0012 %               "In" current source; 2 nodes (-,+); args(1) = current
0013 %               "Hn" current controlled voltage source, 4 nodes (-,+); args(1) = Vsrc, args(2) = gain
0014 %               "En" voltage controlled voltage source, 4 nodes (-,+,d-,d+); args(1) = gain
0015 %               "Fn" current controlled current source, 2 nodes (-,+); args(1) = Vsrc, args(2) = gain
0016 %               "Gn" voltage controlled current source, 4 nodes (-,+,d-,d+); args(1) = gain
0017 %   freq    - simulation frequency (default: 0 Hz, DC)
0018 %   note: ground node is the lowest numbered node
0019 % output
0020 %   Vn      - nodal voltages by row [node voltage]
0021 %   nn      - node numbers
0022 %
0023 % TODO reverse nodes for spice netlist (-,+) --> (+,-)
0024 % TODO support DC voltage/current: [[DC] {value}]; currently supports {value}
0025 % TODO support AC voltage/current: [AC {mag} [{phase}]]
0026 %
0027 % CITATION_REQUEST:
0028 % AUTHOR: A Boyle and A Adler
0029 % TITLE: Integrating Circuit Simulation with EIT FEM Models
0030 % JOURNAL: 19th International Conference on Biomedical Applications of Electrical Impedance Tomography, Edinburgh, UK
0031 % YEAR: 2018
0032 %
0033 
0034 %  (C) 2017, 2018 A. Boyle, License: GPL version 2 or version 3
0035 
0036 % based on http://matlabbyexamples.blogspot.ca/2011/11/circuit-solver-using-matlab-programming.html
0037 % and https://www.swarthmore.edu/NatSci/echeeve1/Ref/mna/MNA5.html
0038 % and https://cseweb.ucsd.edu/classes/sp08/cse245/slides/let1.ppt
0039 % and http://www.ecircuitcenter.com/SPICEsummary.htm
0040 % and http://users.ecs.soton.ac.uk/mz/CctSim/chap1_4.htm
0041 % and https://github.com/nik1106/MNA-MAT
0042 %    (MNA-MAT is a different implementation of a spice-only linear simulator for DC/transient sims in matlab)
0043 % and https://global.oup.com/us/companion.websites/fdscontent/uscompanion/us/static/companion.websites/9780199339136/Appendices/Appendix_B.pdf
0044 % and http://www2.ensc.sfu.ca/~ljilja/papers/4C823d01.pdf
0045 %    (non-linear components to ideal component circuits)
0046 
0047 if ischar(netlist) && strcmp(netlist,'UNIT_TEST'); unittest(); return; end
0048 
0049 if nargin < 2
0050    freq = 0;
0051 end
0052 w = 2*pi*freq;
0053 s = j*w;
0054 
0055 % collect nodes
0056 assert(iscell(netlist), 'error: expected cell array for netlist');
0057 nodes = zeros(length(netlist),2);
0058 delM = zeros(length(netlist),1); % init: tracking L -> V converted inductors
0059 M = 0; % number of voltage sources
0060 P = 0;
0061 for ii = 1:length(netlist)
0062    tt = netlist{ii};
0063    % save node numbers for later
0064    assert(isnumeric([tt{2:3}]), 'error: expected second and third field netlist to be node numbers');
0065    nodes(ii,:) = uint32([tt{2:3}]);
0066    id = tt{1};
0067    assert(ischar(id),'error: expected first field of netlist to be a component identifier');
0068    if abs(s) == 0
0069       % convert inductors to shorts at f=0; avoids 'matrix singular' errors
0070       switch id(1)
0071       case 'L'
0072          id = ['Vshort_f0_' id];
0073          tt{1} = id; tt{4} = 0;
0074          netlist{ii,1} = tt;
0075          delM(M+1) = 1;
0076       end
0077    end
0078    % count voltage sources and components;
0079    % this is a count of how many currents we will calculate
0080    if any(strcmp(id(1), {'V','E','H'}))
0081       M = M +1;
0082       Vidx{M} = id;
0083    end
0084    % extra entries in the sparse indices for large component stamps (dep. sources)
0085    switch id(1)
0086    case 'H'
0087       P = P+1;
0088    case 'E'
0089       P = P+2;
0090    end
0091 end
0092 nodes = uint32(nodes);
0093 gnd = min(nodes(:));
0094 delM = find(delM);
0095 
0096 % remap node#s to be valid array indices
0097 [ii] = unique(sort(nodes(:)));
0098 jj = [1:length(ii)]';
0099 mn = min(nodes(:));
0100 map(ii-mn+1) = jj;
0101 rmap(jj) = ii;
0102 
0103 nodes = map(nodes - mn +1);
0104 nodes = nodes -1; % gnd -> node#0
0105 
0106 
0107 % save node numbers for output
0108 nn = rmap(2:end)';
0109 
0110 
0111 % build matrices
0112 L = size(nodes,1);
0113 N = length(nn); % number of nodes
0114 
0115 %print_netlist('exec',netlist);
0116 
0117 Aii = ones(L*4+P,1); Ajj=Aii; Avv=Aii*0; % init
0118 Bii = ones(L*2,1); Bjj=Bii; Bvv=Bii*0; % init
0119 Mn = N+1;
0120 Pn = L*4+1;
0121 for tt = 1:L
0122    row = netlist{tt}; val = row{4}; id = row{1};
0123    n1 = nodes(tt,1); n2 = nodes(tt,2);
0124    switch(id(1))
0125    case 'V' % independent voltage source
0126       idx = [ ((tt-1)*4+1):((tt-1)*4+4) ];
0127       % no off-diagonal entries for nodes connected to gnd
0128       vp = +1; vn = -1;
0129       if n1 == 0
0130          n1 = 1;
0131          vn = 0;
0132       elseif n2 == 0
0133          n2 = 1;
0134          vp = 0;
0135       end
0136       Aii(idx) = [ Mn Mn n1 n2 ];
0137       Ajj(idx) = [ n1 n2 Mn Mn ];
0138       Avv(idx) = [ vn vp vn vp ];
0139       idx = (tt-1)*2+1;
0140       Bii(idx) = [ Mn  ];
0141       Bvv(idx) = [ val ];
0142       Mn = Mn +1;
0143       continue;
0144    case 'F' % current controlled current source (CCCS)
0145       % from http://users.ecs.soton.ac.uk/mz/CctSim/chap1_4.htm
0146       Vsrc = row{4};
0147       Mi = find(strcmp(Vsrc, Vidx)) + N; % index of branch current
0148       assert(length(Mi) == 1,sprintf('failed to find %s for %s',Vsrc,id));
0149       gain = row{5}; assert(isnumeric(gain)); % gain
0150       idx = [ ((tt-1)*4+1):((tt-1)*4+2) ];
0151       % nodes connected to gnd?
0152       d1 = 1; d2 = 1;
0153       if n1 == 0
0154          n1 = 1;
0155          d1 = 0;
0156       end
0157       if n2 == 0
0158          n2 = 1;
0159          d2 = 0;
0160       end
0161       Aii(idx) = [  n1  n2 ];
0162       Ajj(idx) = [  Mi  Mi ];
0163       Avv(idx) = [ +d1 -d2 ]*gain;
0164       continue;
0165    case 'G' % voltage controlled current source (VCCS)
0166       % from https://cseweb.ucsd.edu/classes/sp08/cse245/slides/let1.ppt
0167       n1 = uint32(row{2}); % n-
0168       n2 = uint32(row{3}); % n+
0169       n3 = uint32(row{4}); % diff-
0170       n4 = uint32(row{5}); % diff+
0171       gain = row{6}; assert(isnumeric(gain)); % gain
0172       idx = [ ((tt-1)*4+1):((tt-1)*4+4) ];
0173       % no off-diagonal entries for nodes connected to gnd
0174       d1 = 1; d2 = 1; d3 = 1; d4 = 1;
0175       if n1 == 0
0176          n1 = 1;
0177          d1 = 0;
0178       end
0179       if n2 == 0
0180          n2 = 1;
0181          d2 = 0;
0182       end
0183       if n3 == 0
0184          n3 = 1;
0185          d3 = 0;
0186       end
0187       if n4 == 0
0188          n4 = 1;
0189          d4 = 0;
0190       end
0191       % A = [G B;C D] --> same as V; except C is changed
0192       %            C  C  B  B
0193       Aii(idx) = [  n1     n2     n1     n2 ];
0194       Ajj(idx) = [  n4     n3     n3     n4 ];
0195       Avv(idx) = [ -d1*d4 -d2*d3 +d1*d3 +d2*d4 ]*gain;
0196       continue;
0197    case 'E' % voltage controlled voltage source (VCVS)
0198       % from http://users.ecs.soton.ac.uk/mz/CctSim/chap1_4.htm
0199       n1 = uint32(row{2}); % n-
0200       n2 = uint32(row{3}); % n+
0201       n3 = uint32(row{4}); % diff-
0202       n4 = uint32(row{5}); % diff+
0203       gain = row{6}; assert(isnumeric(gain)); % gain
0204       idx = [ ((tt-1)*4+1):((tt-1)*4+4) Pn Pn+1 ];
0205       % no off-diagonal entries for nodes connected to gnd
0206       d1 = 1; d2 = 1; d3 = 1; d4 = 1;
0207       if n1 == 0
0208          n1 = 1;
0209          d1 = 0;
0210       end
0211       if n2 == 0
0212          n2 = 1;
0213          d2 = 0;
0214       end
0215       if n3 == 0
0216          n3 = 1;
0217          d3 = 0;
0218       end
0219       if n4 == 0
0220          n4 = 1;
0221          d4 = 0;
0222       end
0223       % A = [G B;C D] --> same as V; except C is changed
0224       %            C  C  B  B
0225       G = gain;
0226       Aii(idx) = [  Mn  Mn  n1  n2  Mn    Mn   ];
0227       Ajj(idx) = [  n1  n2  Mn  Mn  n3    n4   ];
0228       Avv(idx) = [ -d1 +d2 -d1 +d2 +d3*G -d4*G ];
0229       Pn = Pn +2;
0230       Mn = Mn +1;
0231       continue;
0232    case 'H' % current controlled voltage source (CCVS)
0233       % from http://users.ecs.soton.ac.uk/mz/CctSim/chap1_4.htm
0234       Vsrc = row{4};
0235       Mi = find(strcmp(Vsrc, Vidx)) + N; % index of branch current
0236       assert(length(Mi) == 1,sprintf('failed to find %s for %s',Vsrc,id));
0237       gain = row{5}; assert(isnumeric(gain)); % gain
0238       idx = [ ((tt-1)*4+1):((tt-1)*4+4) Pn ];
0239       % no off-diagonal entries for nodes connected to gnd
0240       d1 = 1; d2 = 1; d3 = 1; d4 = 1;
0241       if n1 == 0
0242          n1 = 1;
0243          d1 = 0;
0244       end
0245       if n2 == 0
0246          n2 = 1;
0247          d2 = 0;
0248       end
0249       Aii(idx) = [  Mn  Mn  n1  n2  Mn   ];
0250       Ajj(idx) = [  n1  n2  Mn  Mn  Mi   ];
0251       Avv(idx) = [ -d1 +d2 -d1 +d2  gain ];
0252       Pn = Pn +1;
0253       Mn = Mn +1;
0254       continue;
0255    case 'I' % independent current source
0256       idx = [ ((tt-1)*2+1):((tt-1)*2+2) ];
0257       ip = val; in = -val;
0258       if n1 == 0
0259          n1 = 1;
0260          in = 0;
0261       elseif n2 == 0
0262          n2 = 1;
0263          ip = 0;
0264       end
0265       Bii(idx) = [ n1 n2 ];
0266       Bvv(idx) = [ in ip ];
0267       continue; 
0268    case 'R' % resistor
0269       y = 1/val;
0270    case 'C' % capacitor
0271       y = s*val;
0272    case 'L' % inductor
0273       y = 1/(s*val);
0274    otherwise
0275       error(['unhandled netlist component type: ' id]);
0276    end
0277    if n1 == n2; error(sprintf('bad netlist @ line#%d: n1(%d)==n2(%d)',tt,n1,n2)); end
0278 
0279    % passive elements
0280    yn = -y; y1 = y; y2 = y;
0281    % no off-diagonal entries for nodes connected to gnd
0282    if n1 == 0
0283       n1 = 1;
0284       yn = 0;
0285       y1 = 0;
0286    elseif n2 == 0
0287       n2 = 1;
0288       yn = 0;
0289       y2 = 0;
0290    end
0291    idx = [ ((tt-1)*4+1):((tt-1)*4+4) ];
0292    Aii(idx) = [n1 n2 n1 n2];
0293    Ajj(idx) = [n2 n1 n1 n2];
0294    Avv(idx) = [yn yn y1 y2];
0295    % ... sparse-ified
0296    %A(n1,n2) = A(n1,n2)-y;
0297    %A(n2,n1) = A(n2,n1)-y;
0298    %A(n1,n1) = A(n1,n1)+y;
0299    %A(n2,n2) = A(n2,n2)+y;
0300 end
0301 % delete ground node
0302 %idx = find(ii==gnd); ii(idx) = []; jj(idx) = []; vv(idx) = [];
0303 %idx = find(jj==gnd); ii(idx) = []; jj(idx) = []; vv(idx) = [];
0304 % build matrix
0305 A = sparse(Aii,Ajj,Avv,N+M,N+M);
0306 B = sparse(Bii,Bjj,Bvv,N+M,1);
0307 
0308 %AA=full(A)
0309 %BB=full(B)
0310 
0311 % should be sparse and symmetric! check we use the right solver
0312 X = A\B;
0313 Vn = full(X(1:N));
0314 In = full(X(N+1:N+M));
0315 In(delM) = [];
0316 
0317 function unittest()
0318    pass = 1;
0319    netlist_rlc = {{'R1', 0, 1, 1e3},
0320                   {'L1', 0, 1, 1},
0321                   {'C1', 0, 1, 5e-6},
0322                   {'RC2', 2, 1, 1},
0323                   {'R2', 3, 2, 1},
0324                   {'V1', 0, 3, 1}};
0325    tol = eps*10;
0326    for ii = 1:9
0327       switch(ii)
0328       case 1
0329          desc = 'voltage divider';
0330          netlist = {{'R1', 1, 2, 1},
0331                     {'R2', 0, 2, 1},
0332                     {'V0', 0, 1, 2}};
0333          Ev = [2,1]';
0334          Ei = [-1]';
0335          Enn = [1 2]';
0336          freq = 0;
0337       case 2
0338          desc = 'resistor network; gnd = n#4';
0339          % from http://matlabbyexamples.blogspot.ca/2011/11/circuit-solver-using-matlab-programming.html
0340          netlist = {{'R1', 1, 2, 2},
0341                     {'R2', 1, 3, 4},
0342                     {'R3', 2, 3, 5.2},
0343                     {'R4', 3, 4, 6},
0344                     {'R5', 2, 4, 3},
0345                     {'V1', 4, 1, 6},
0346                     {'V0', 0, 4, 0}};
0347          Ev = [6,3.6,3.6,0]';
0348          Ei = [-1.8,0]';
0349          Enn = [1 2 3 4]';
0350          freq = 0;
0351       case 3
0352          desc = 'reactive elements, f=0';
0353          netlist = netlist_rlc;
0354          Ev = [0.0,0.5,1.0]';
0355          Ei = [-0.5]';
0356          Enn = [1 2 3]';
0357          freq = 0;
0358       case 4
0359          desc = 'reactive elements, f=1';
0360          netlist = netlist_rlc;
0361          Ev = [0.90655 + 0.28793i;
0362                0.95328 + 0.14397i;
0363                1];
0364          Ei = -0.04672 + 0.14397i;
0365          Enn = [1 2 3]';
0366          freq = 1;
0367          tol = 6e-6;
0368       case 5
0369          desc = 'reactive elements, f=10';
0370          netlist = netlist_rlc;
0371          Ev = [0.99704 + 0.03105i;
0372                0.99852 + 0.01552i;
0373                1.0];
0374          Ei = -0.00148 + 0.01552i;
0375          Enn = [1 2 3]';
0376          freq = 10;
0377       case 6
0378          desc = 'voltage controlled current source VCCS G';
0379          netlist = {{'V1', 0, 1, 1},
0380                     {'R1', 0, 1, 1},
0381                     {'RL', 0, 2, 1},
0382                     {'G1', 0, 2, 0, 1, 5}};
0383          Ev = [1,-5]';
0384          Ei = [-1]';
0385          Enn = [1 2]';
0386          freq = 0;
0387       case 7
0388          desc = 'current controlled current source CCCS F';
0389          netlist = {{'V1', 0, 1, 1},
0390                     {'R1', 0, 1, 1},
0391                     {'RL', 0, 2, 1},
0392                     {'F1', 0, 2, 'V1', 5}};
0393          Ev = [1,-5]';
0394          Ei = [-1]';
0395          Enn = [1 2]';
0396          freq = 0;
0397       case 8
0398          desc = 'voltage controlled voltage source VCVS E';
0399          netlist = {{'V1', 0, 1, 1},
0400                     {'R1', 0, 1, 1},
0401                     {'RL', 0, 2, 1},
0402                     {'E1', 0, 2, 0, 1, 5}};
0403          Ev = [1,+5]';
0404          Ei = [-1,-5]';
0405          Enn = [1 2]';
0406          freq = 0;
0407       case 9
0408          desc = 'current controlled voltage source CCVS H';
0409          netlist = {{'V1', 0, 1, 1},
0410                     {'R1', 0, 1, 1},
0411                     {'RL', 0, 2, 1},
0412                     {'H1', 0, 2, 'V1', 5}};
0413          Ev = [1,+5]';
0414          Ei = [-1,-5]';
0415          Enn = [1 2]';
0416          freq = 0;
0417       otherwise
0418          error('oops');
0419       end
0420       desc = [sprintf('#%d ',ii) desc];
0421 
0422       disp('-----------------------------------------');
0423       print_netlist(desc,netlist);
0424       [v,i,nn]=spice_eit(netlist,freq);
0425       pass=cmp(pass,'voltages',v,Ev,tol);
0426       pass=cmp(pass,'currents',i,Ei,tol);
0427       pass=cmp(pass,'node#s',nn,Enn);
0428    end
0429    disp('-----------------------------------------');
0430    if pass
0431       disp('overall: PASS');
0432    else
0433       disp('overall: FAIL');
0434    end
0435 
0436 function print_netlist(desc,n)
0437    fprintf('NETLIST: %s\n',desc);
0438    for ii =1:length(n)
0439       nn = n{ii}; id = nn{1};
0440       if any(strcmp(id(1),{'G','E'}))
0441          assert(length(nn) == 6,sprintf('bad netlist, row#%d has %d elements',ii,length(nn)));
0442          fprintf('  %-5s (%2d,%2d)<-(%2d,%2d)  %0.3f\n',nn{1},nn{2},nn{3},nn{4},nn{5},nn{6});
0443       elseif any(strcmp(id(1),{'F','H'}))
0444          assert(length(nn) == 5,sprintf('bad netlist, row#%d has %d elements',ii,length(nn)));
0445          fprintf('  %-5s (%2d,%2d)<-(%-5s)  %0.3f\n',nn{1},nn{2},nn{3},nn{4},nn{5});
0446       else
0447          assert(length(nn) == 4,sprintf('bad netlist, row#%d has %d elements',ii,length(nn)));
0448          fprintf('  %-5s (%2d,%2d)           %0.3f\n',nn{1},nn{2},nn{3},nn{4});
0449       end
0450    end
0451 
0452 function pass=cmp(pass,str,X,Y,tol)
0453    if nargin < 5
0454       tol = 0;
0455    end
0456    if any(size(X) ~= size(Y))
0457       fprintf('TEST: %-30s fail; size mismatch (%d,%d)!=(%d,%d)\n',str,size(X),size(Y));
0458       pass=0;
0459       return
0460    end
0461    if tol == 0
0462       if any(X ~= Y)
0463          [ X Y ]
0464          fprintf('TEST: %-30s fail (%g)\n',str,abs(max(X(:)-Y(:))));
0465          pass=0;
0466          return;
0467       end
0468    else
0469       err=abs(X-Y);
0470       err(err<tol)=[];
0471       if length(err) > 0
0472          [ X Y ]
0473          fprintf('TEST: %-30s fail (%g, tol=%g)\n',str,norm(err),tol);
0474          pass=0;
0475          return;
0476       end
0477    end
0478    fprintf('TEST: %-30s pass\n',str);
spice_eit

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
