0001 function [Vn, In, nn] = spice_eit(netlist, freq)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
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
0051 assert(iscell(netlist), 'error: expected cell array for netlist');
0052 nodes = zeros(length(netlist),2);
0053 delM = zeros(length(netlist),1);
0054 M = 0;
0055 P = 0;
0056 for ii = 1:length(netlist)
0057 tt = netlist{ii};
0058
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
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
0074
0075 if any(strcmp(id(1), {'V','E','H'}))
0076 M = M +1;
0077 Vidx{M} = id;
0078 end
0079
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
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;
0100
0101
0102
0103 nn = rmap(2:end)';
0104
0105
0106
0107 L = size(nodes,1);
0108 N = length(nn);
0109
0110
0111
0112 Aii = ones(L*4+P,1); Ajj=Aii; Avv=Aii*0;
0113 Bii = ones(L*2,1); Bjj=Bii; Bvv=Bii*0;
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'
0121 idx = [ ((tt-1)*4+1):((tt-1)*4+4) ];
0122
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'
0140
0141 Vsrc = row{4};
0142 Mi = find(strcmp(Vsrc, Vidx)) + N;
0143 assert(length(Mi) == 1,sprintf('failed to find %s for %s',Vsrc,id));
0144 gain = row{5}; assert(isnumeric(gain));
0145 idx = [ ((tt-1)*4+1):((tt-1)*4+2) ];
0146
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'
0161
0162 n1 = uint32(row{2});
0163 n2 = uint32(row{3});
0164 n3 = uint32(row{4});
0165 n4 = uint32(row{5});
0166 gain = row{6}; assert(isnumeric(gain));
0167 idx = [ ((tt-1)*4+1):((tt-1)*4+4) ];
0168
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
0187
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'
0193
0194 n1 = uint32(row{2});
0195 n2 = uint32(row{3});
0196 n3 = uint32(row{4});
0197 n4 = uint32(row{5});
0198 gain = row{6}; assert(isnumeric(gain));
0199 idx = [ ((tt-1)*4+1):((tt-1)*4+4) Pn Pn+1 ];
0200
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
0219
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'
0228
0229 Vsrc = row{4};
0230 Mi = find(strcmp(Vsrc, Vidx)) + N;
0231 assert(length(Mi) == 1,sprintf('failed to find %s for %s',Vsrc,id));
0232 gain = row{5}; assert(isnumeric(gain));
0233 idx = [ ((tt-1)*4+1):((tt-1)*4+4) Pn ];
0234
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'
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'
0264 y = 1/val;
0265 case 'C'
0266 y = s*val;
0267 case 'L'
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
0275 yn = -y; y1 = y; y2 = y;
0276
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
0291
0292
0293
0294
0295 end
0296
0297
0298
0299
0300 A = sparse(Aii,Ajj,Avv,N+M,N+M);
0301 B = sparse(Bii,Bjj,Bvv,N+M,1);
0302
0303
0304
0305
0306
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
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);