eit_spice

function spice = eit_spice(img, [name])
 Converts an EIT FEM model with assigned conductivities (an EIDORS "img") to a
 model reduced, fully connected mesh of resistors in SPICE format.
 If the FEM model has complex valued conductivities, the mesh will be an RLC
 mesh network.

 An optional subcircuit 'name' can be provided.

 TODO complex value support
 TODO fix electrode ordering for mixed PEM/CEM electrodes

 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

0001 function spice = eit_spice(img, name)
0002 %function spice = eit_spice(img, [name])
0003 % Converts an EIT FEM model with assigned conductivities (an EIDORS "img") to a
0004 % model reduced, fully connected mesh of resistors in SPICE format.
0005 % If the FEM model has complex valued conductivities, the mesh will be an RLC
0006 % mesh network.
0007 %
0008 % An optional subcircuit 'name' can be provided.
0009 %
0010 % TODO complex value support
0011 % TODO fix electrode ordering for mixed PEM/CEM electrodes
0012 %
0013 % CITATION_REQUEST:
0014 % AUTHOR: A Boyle and A Adler
0015 % TITLE: Integrating Circuit Simulation with EIT FEM Models
0016 % JOURNAL: 19th International Conference on Biomedical Applications of Electrical Impedance Tomography, Edinburgh, UK
0017 % YEAR: 2018
0018 %
0019 
0020 %  (C) 2018 A. Boyle, License: GPL version 2 or version 3
0021 
0022    if ischar(img) & strcmp(img, 'UNIT_TEST') unit_test(); return; end
0023 
0024    if nargin == 1
0025       name = 'eit';
0026    end
0027 
0028    Dprime = model_reduce(img);
0029 %  disp(full(1./Dprime))
0030 %  disp(full(-1./(Dprime - tril(Dprime))));
0031    spice = netlist(Dprime,name);
0032 
0033    if nargout == 0
0034       filename = [ name '.s' ];
0035       FILE = fopen(filename, 'w');
0036       fprintf(FILE,'%s\n',spice{:});
0037       fclose(FILE);
0038       eidors_msg(['saved SPICE netlist to ' filename]);
0039       return
0040    end
0041 end
0042 
0043 function Dprime = model_reduce(img)
0044    Y = calc_system_mat(img); Y= Y.E;
0045    nn= num_nodes(img);
0046    % Decompose into blocks; assumes that the nn+1:end nodes are CEM electrodes
0047    rm = 1:nn; % nodes to fold
0048    kp = nn+1:size(Y,1); % nodes to keep
0049    % Now handle PEM electrodes, by transferring nodes between the rm and el sets
0050    for i = 1:length(img.fwd_model.electrode)
0051       el = img.fwd_model.electrode(i);
0052       if length(el.nodes) == 1
0053          rm(rm == el.nodes) = [];
0054          kp(end+1) = el.nodes;
0055       end
0056    end
0057    % Note: C = B' ... we don't need to calculate it for symmetric matrices
0058    A = Y(rm,rm); B= Y(rm,kp); D = Y(kp,kp);
0059    Dprime  = D - B'*inv(A)*B;
0060 end
0061 
0062 function out = netlist(Dprime, name)
0063    nn = size(Dprime,1);
0064    ndr = floor(log10(nn*(nn-1)/2))+1; % number of digits for resistors
0065    nde = floor(log10(nn))+1; % number of digits for electrodes
0066    str = ['.subckt ' name ];
0067    for ii = 1:nn;
0068       str = [ str sprintf([' e%0' num2str(nde) 'd'], ii) ];
0069    end
0070    out = { str };
0071    str = ['re%0' num2str(ndr) 'd e%0' num2str(nde) 'd e%0' num2str(nde) 'd %s'];
0072    rr = 1;
0073    for ii = 1:nn;
0074       for jj = (ii+1):nn;
0075          val = sprintf('%0.17g',-1/Dprime(ii,jj));
0076          out(end+1,1) = { strrep(sprintf(str,rr,ii,jj,val),'+','') }; % we strip '+'
0077          rr = rr +1;
0078       end
0079    end
0080    out(end+1,1) = { '.ends' };
0081 end
0082 
0083 function unit_test()
0084    imdl = mk_common_model('a2s',8);
0085    stim = mk_stim_patterns(8,1,'{ad}','{ad}',{'meas_current'},1);
0086    imdl.fwd_model.stimulation = stim(1);
0087    imdl.fwd_model = rmfield(imdl.fwd_model, 'meas_select');
0088    img = mk_image(imdl,1);
0089    v = fwd_solve(img);
0090    disp('stim');
0091    disp(full(stim(1).stim_pattern));
0092    disp(stim(1).stimulation)
0093    disp('meas');
0094    disp(full(stim(1).meas_pattern));
0095    disp('voltages');
0096    disp(v.meas)
0097    eit_spice(img)
0098 end