Home > eidors > solvers > forward > spice_eit.m

spice_eit

PURPOSE ^

Vn = spice_eit(netlist, [freq])

SYNOPSIS ^

function [Vn, In, nn] = spice_eit(netlist, freq)

DESCRIPTION ^

 Vn = spice_eit(netlist, [freq])
 Modified Nodal Analysis (MNA) circuit solver
 for (almost/simple) spice netlists
 input
   netlist - by cell array { [ part nodes...  args ], ... }
             where parts are spice types:
               "Rn" resistor; 2 nodes; args(1) = impedance
               "Ln" inductor; 2 nodes; args(1) = inductance
               "Cn" capacitor; 2 nodes; args(1) = capacitance
               "Vn" voltage source; 2 nodes (-,+); args(1) = voltage
               "In" current source; 2 nodes (-,+); args(1) = current
               "Hn" current controlled voltage source, 4 nodes (-,+); args(1) = Vsrc, args(2) = gain
               "En" voltage controlled voltage source, 4 nodes (-,+,d-,d+); args(1) = gain
               "Fn" current controlled current source, 2 nodes (-,+); args(1) = Vsrc, args(2) = gain
               "Gn" voltage controlled current source, 4 nodes (-,+,d-,d+); args(1) = gain
   freq    - simulation frequency (default: 0 Hz, DC)
   note: ground node is the lowest numbered node
 output
   Vn      - nodal voltages by row [node voltage]
   nn      - node numbers

 TODO reverse nodes for spice netlist (-,+) --> (+,-)
 TODO support DC voltage/current: [[DC] {value}]; currently supports {value}
 TODO support AC voltage/current: [AC {mag} [{phase}]]

 CITATION_REQUEST:
 AUTHOR: A Boyle and A Adler
 TITLE: Integrating Circuit Simulation with EIT FEM Models 
 JOURNAL: 19th International Conference on Biomedical Applications of Electrical Impedance Tomography, Edinburgh, UK
 YEAR: 2018

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Fri 01-Jun-2018 15:59:55 by m2html © 2005