system_mat_instrument

 SYSTEM_MAT_INSTRUMENT: 
  Calculate systems matrix inclusing modelling of instrument impedance
 img.fwd_model.system_mat = @system_mat_instrument
 img.fwd_model.system_mat_instrument.clist = [CONNECTION LIST]
 
 where
  connection list =
     [node1a, node1b, admittance_1ab]
     [node2a, node2b, admittance_2ab]

 Instrument electrode can be added using
  fmdl.electrode(end+1) = struct('nodes','instrument','z_contact',NaN);
  Such electrodes must be last

 CITATION_REQUEST:
 AUTHOR: A Adler
 TITLE: Modelling instrument admittances with EIDORS
 CONFERENCE: EIT 2021
 YEAR: 2021
 PAGE: 74
 LINK: https://zenodo.org/record/4940249

0001 function s_mat = system_mat_instrument( img);
0002 % SYSTEM_MAT_INSTRUMENT:
0003 %  Calculate systems matrix inclusing modelling of instrument impedance
0004 % img.fwd_model.system_mat = @system_mat_instrument
0005 % img.fwd_model.system_mat_instrument.clist = [CONNECTION LIST]
0006 %
0007 % where
0008 %  connection list =
0009 %     [node1a, node1b, admittance_1ab]
0010 %     [node2a, node2b, admittance_2ab]
0011 %
0012 % Instrument electrode can be added using
0013 %  fmdl.electrode(end+1) = struct('nodes','instrument','z_contact',NaN);
0014 %  Such electrodes must be last
0015 %
0016 % CITATION_REQUEST:
0017 % AUTHOR: A Adler
0018 % TITLE: Modelling instrument admittances with EIDORS
0019 % CONFERENCE: EIT 2021
0020 % YEAR: 2021
0021 % PAGE: 74
0022 % LINK: https://zenodo.org/record/4940249
0023 
0024 % (C) 2021 Andy Adler. License: GPL version 2 or version 3
0025 % $Id: system_mat_instrument.m 6478 2022-12-26 00:26:34Z aadler $
0026 
0027    if ischar(img) && strcmp(img,'UNIT_TEST'); do_unit_test; return; end
0028 
0029    citeme(mfilename)
0030 
0031    new_c_list = img.fwd_model.system_mat_instrument.connect_list;
0032 
0033 % pre-delete all 'instrument' nodes
0034 % system matrix can't handle them
0035 % BUG ALERT: instrument nodes must be at end
0036    for i=num_elecs(img):-1:1
0037       enodes = img.fwd_model.electrode(i).nodes;
0038       if ischar(enodes) && strcmp(enodes,'instrument');
0039          img.fwd_model.electrode(i)= [];
0040       end
0041    end
0042    s_mat= system_mat_1st_order(img);
0043    E = s_mat.E;
0044 
0045    n_nodes = num_nodes(img);
0046    if ~isempty(new_c_list);
0047       n_elecs = max([max(new_c_list(:,1:2)),num_elecs(img)]);
0048    else
0049       n_elecs = num_elecs(img);
0050    end
0051    n_max  = n_nodes + n_elecs; 
0052    Eo = sparse(n_max, n_max);
0053    Eo(1:size(E,1),1:size(E,2)) = E;
0054    for i=1:size(new_c_list,1);
0055      x = new_c_list(i,1) + n_nodes;
0056      y = new_c_list(i,2) + n_nodes;
0057      c = new_c_list(i,3);
0058      Eo(x,y) = Eo(x,y) - c;
0059      Eo(y,x) = Eo(y,x) - c;
0060      Eo(x,x) = Eo(x,x) + c;
0061      Eo(y,y) = Eo(y,y) + c;
0062    end
0063    
0064    s_mat.E = Eo;
0065 
0066 function do_unit_test
0067    % Example from conference paper
0068    fmdl = eidors_obj('fwd_model','eg', ...
0069        'nodes',[0,0;0,1;2,0;2,1], ...
0070        'elems',[1,2,3;2,3,4], ...
0071        'gnd_node',1);
0072    fmdl = mdl_normalize(fmdl,1);
0073    fmdl.system_mat = @eidors_default;
0074    fmdl.solve      = @eidors_default;
0075    img = mk_image(fmdl,[1,2]);
0076    s_mat= calc_system_mat(img); 
0077    Eok= [ 5,-4,-1, 0;-4, 6, 0,-2;
0078          -1, 0, 9,-8; 0,-2,-8,10]/4;
0079    unit_test_cmp('basic',s_mat.E,Eok);
0080 
0081    img.fwd_model.electrode = [ ...
0082      struct('nodes',[1,2],'z_contact',5/30), ...
0083      struct('nodes',[3,4],'z_contact',5/60)];
0084    s_mat= calc_system_mat(img); 
0085 
0086    Eok(:,5:6) = 0; Eok(5:6,:) = 0;
0087    Eok = Eok + [ ...
0088      2 1 0 0 -3 0; 1 2 0 0 -3 0;
0089      0  0 4 2 0 -6; 0 0  2  4 0 -6
0090     -3 -3 0 0 6  0; 0 0 -6 -6 0 12];
0091    unit_test_cmp('with CEM',s_mat.E,Eok,1e-13);
0092 
0093    img.fwd_model.stimulation = ...
0094       stim_meas_list([1,2,1,2]); 
0095 
0096    imgs = img; imgs.elem_data = [1;1];
0097    vv= fwd_solve(imgs);
0098    Ztest =  2+5/30+5/60;
0099    unit_test_cmp('CEM solve',vv.meas,Ztest,1e-13);
0100 
0101 % NOTE: there is a bug in the paper,
0102 % actually, the connection is to newEl
0103    NewElectrode = 3; % add new electrode
0104    Y12 = 1/10; % Y=1/10Ohms between E#1 and new electrode
0105    Y13 = 1/ 5; % Y=1/ 5Ohms between E#1 and new elec
0106    c_list = [ 1,3,Y12;
0107               2,NewElectrode,Y13];
0108    % Add instrument electrode
0109    img.fwd_model.electrode(NewElectrode) = ... 
0110        struct('nodes','instrument','z_contact',NaN);
0111    img.fwd_model.stimulation = ...
0112       stim_meas_list([1,3,1,3]); 
0113 
0114    img.fwd_model.system_mat = @system_mat_instrument;
0115    img.fwd_model.system_mat_instrument.connect_list = c_list;
0116    s_mat= calc_system_mat(img); 
0117 
0118    Eok(:,7) = 0; Eok(7,:) = 0;
0119    Eok = Eok + [ 0 0 0 0 0 0 0;
0120 0 0 0 0   0   0   0; 0 0 0 0   0   0   0; 
0121 0 0 0 0   0   0   0; 0 0 0 0  .1   0 -.1;
0122 0 0 0 0   0  .2 -.2; 0 0 0 0 -.1 -.2  .3];
0123    unit_test_cmp('instrument',s_mat.E,Eok,1e-13);
0124 
0125    imgs = img; imgs.elem_data = [1;1];
0126    vv= fwd_solve(imgs);
0127    Ztest = 1/( Y12 + 1/(Ztest + 1/Y13));
0128    
0129    unit_test_cmp('instrument solve1',vv.meas, Ztest,1e-13);
0130 
0131    % apply current to ground
0132    imgs.fwd_model.stimulation = struct( ...
0133       'stim_pattern',[1;0;NaN], ...
0134       'volt_pattern',[0;0;0], ...
0135       'meas_pattern',[1,0,-1;0,1,-1]);
0136    vv= fwd_solve(imgs);
0137    Ztest2 = 10/(7.25+10)*5;
0138    unit_test_cmp('instrument solve2',vv.meas, [Ztest;Ztest2],1e-13);
0139 
0140    % apply voltage to ground
0141    imgs.fwd_model.stimulation = struct( ...
0142       'stim_pattern',[NaN;0;NaN], ...
0143       'volt_pattern',[10;0;0], ...
0144       'meas_pattern',[1,0,-1;0,1,-1]);
0145    vv= fwd_solve(imgs);
0146    Ztest = 10*[1;5/7.25];
0147    unit_test_cmp('instrument solve3',vv.meas, Ztest,1e-13);
0148 
0149 
0150 %  disp(Eok)
0151 %  disp(full(s_mat.E))