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 6797 2024-04-20 19:37:13Z 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 
0034    img.fwd_model = remove_instrument_electrodes(img.fwd_model);
0035    s_mat= system_mat_1st_order(img);
0036    E = s_mat.E;
0037 
0038    n_nodes = num_nodes(img);
0039    if ~isempty(new_c_list);
0040       n_elecs = max([max(new_c_list(:,1:2)),num_elecs(img)]);
0041    else
0042       n_elecs = num_elecs(img);
0043    end
0044    n_max  = n_nodes + n_elecs; 
0045    Eo = sparse(n_max, n_max);
0046    Eo(1:size(E,1),1:size(E,2)) = E;
0047    for i=1:size(new_c_list,1);
0048      x = new_c_list(i,1) + n_nodes;
0049      y = new_c_list(i,2) + n_nodes;
0050      c = new_c_list(i,3);
0051      Eo(x,y) = Eo(x,y) - c;
0052      Eo(y,x) = Eo(y,x) - c;
0053      Eo(x,x) = Eo(x,x) + c;
0054      Eo(y,y) = Eo(y,y) + c;
0055    end
0056    
0057    s_mat.E = Eo;
0058 
0059 function do_unit_test
0060    do_unit_test_x;
0061    do_unit_test_paper
0062    do_unit_test_2d
0063    do_unit_test_3d
0064    do_cache_test
0065 
0066 % This is to test a bug in dev/todo.txt for 3.12. Bug was not seen: 2024-04-20
0067 function do_cache_test
0068    fmdl = getfield( ...
0069     mk_common_model('a2C2',2),'fwd_model');
0070    fmdl.electrode(end+1) = struct( ...
0071           'nodes','instrument','z_contact',NaN);
0072    fmdl.stimulation = stim_meas_list([1,3,1,3],3);
0073    fmdl.system_mat = @system_mat_instrument;
0074    fmdl.system_mat_instrument.connect_list = ...
0075         [2,3,10];
0076    vA = fwd_solve(mk_image(fmdl,1));
0077    fmdl.system_mat_instrument.connect_list = ...
0078         [2,3,20];
0079    vB = fwd_solve(mk_image(fmdl,1)) ;
0080    unit_test_cmp('check different (cache)', vA.meas > vB.meas, true);
0081 
0082 function do_unit_test_x
0083    % Example from conference paper
0084    fmdl = mk_circ_tank(1,[],2);
0085    zC = 0.1;
0086    for i=1:num_elecs(fmdl);
0087       fmdl.electrode(i).nodes = num_nodes(fmdl) - 2*i + [1,2];
0088       fmdl.electrode(i).z_contact = zC;
0089    end
0090 
0091    fmdl.system_mat = @system_mat_instrument;
0092    fmdl.system_mat_instrument.connect_list = [];
0093    fmdl.stimulation = stim_meas_list([1,2,1,2], ...
0094                                 num_elecs(fmdl) ); 
0095    img= mk_image(fmdl,1);
0096    s_mat= calc_system_mat(img);
0097    vA = fwd_solve(img);
0098    % Solve with no resistors vA = 1+len*zC;
0099    unit_test_cmp('square:', vA.meas, 1+sqrt(2)*zC,1e-14);
0100 
0101    NewElectrode = num_elecs(fmdl)+1;
0102    fmdl.electrode(NewElectrode) = struct( ...
0103           'nodes','instrument','z_contact',NaN);
0104    zC = 1000;
0105    fmdl.system_mat_instrument.connect_list = ...
0106         [1,3,1/zC];
0107    fmdl.stimulation = stim_meas_list([1,3,1,3], ...
0108                                 num_elecs(fmdl) ); 
0109    img = mk_image(fmdl,1);
0110 
0111 %  show_fem(img)
0112    vB = fwd_solve(img);
0113    unit_test_cmp('square2:', vB.meas, zC,1e-12);
0114 
0115    img.fwd_model.system_mat_instrument.connect_list = ...
0116         [2,3,2/zC;
0117          1,3,1/zC];
0118    vC = fwd_solve(img);
0119    test = 1/(1/zC + 1/(zC/2 + 1 + 0.1*sqrt(2)));
0120    unit_test_cmp('square3:', vC.meas, test,1e-12);
0121 
0122 function do_unit_test_2d
0123    img = mk_image(mk_common_model('a2C2',2));
0124    s_mat = calc_system_mat(img);
0125    test= zeros(1,num_nodes(img)+2);
0126    test(1,1:5)= [4,-1,-1,-1,-1];
0127    unit_test_cmp('basic',s_mat.E(1,:),test,10*eps);
0128 
0129    img.fwd_model.stimulation = stim_meas_list([1,2,1,2],2);
0130    v0 = fwd_solve(img);
0131    Y12 = 1/10; % Y=1/10Ohms between E#1 and new electrode
0132    c_list = [ 1,2,Y12];
0133    img.fwd_model.system_mat = @system_mat_instrument;
0134    img.fwd_model.system_mat_instrument.connect_list = c_list;
0135    vv = fwd_solve(img);
0136    test = 1/(1/v0.meas + Y12);
0137    unit_test_cmp('Resistor',vv.meas,test,1e3*eps);
0138   
0139    NewElectrode = num_elecs(img)+1; % add new electrode
0140    Y13 = 1/ 5; % Y=1/ 5Ohms between E#1 and new elec
0141    c_list = [ 1,2,Y12;
0142               2,NewElectrode,Y13];
0143    img.fwd_model.system_mat_instrument.connect_list = c_list;
0144    % Add instrument electrode
0145    img.fwd_model.electrode(NewElectrode) = ... 
0146        struct('nodes','instrument','z_contact',NaN);
0147    img.fwd_model.stimulation = stim_meas_list([1,2,1,3],NewElectrode);
0148    s_mat= calc_system_mat(img);
0149    vv = fwd_solve(img);
0150 
0151 function do_unit_test_paper
0152    % Example from conference paper
0153    fmdl = eidors_obj('fwd_model','eg', ...
0154        'nodes',[0,0;0,1;2,0;2,1], ...
0155        'elems',[1,2,3;2,3,4], ...
0156        'gnd_node',1);
0157    fmdl = mdl_normalize(fmdl,1);
0158    fmdl.system_mat = @eidors_default;
0159    fmdl.solve      = @eidors_default;
0160    img = mk_image(fmdl,[1,2]);
0161    s_mat= calc_system_mat(img); 
0162    Eok= [ 5,-4,-1, 0;-4, 6, 0,-2;
0163          -1, 0, 9,-8; 0,-2,-8,10]/4;
0164    unit_test_cmp('basic',s_mat.E,Eok);
0165 
0166    img.fwd_model.electrode = [ ...
0167      struct('nodes',[1,2],'z_contact',5/30), ...
0168      struct('nodes',[3,4],'z_contact',5/60)];
0169    s_mat= calc_system_mat(img); 
0170 
0171    Eok(:,5:6) = 0; Eok(5:6,:) = 0;
0172    Eok = Eok + [ ...
0173      2 1 0 0 -3 0; 1 2 0 0 -3 0;
0174      0  0 4 2 0 -6; 0 0  2  4 0 -6
0175     -3 -3 0 0 6  0; 0 0 -6 -6 0 12];
0176    unit_test_cmp('with CEM',s_mat.E,Eok,1e-13);
0177 
0178    img.fwd_model.stimulation = ...
0179       stim_meas_list([1,2,1,2], ...
0180       num_elecs(img) ); 
0181 
0182    imgs = img; imgs.elem_data = [1;1];
0183    vv= fwd_solve(imgs);
0184    Ztest =  2+5/30+5/60;
0185    unit_test_cmp('CEM solve',vv.meas,Ztest,1e-13);
0186 
0187 % NOTE: there is a bug in the paper,
0188 % actually, the connection is to newEl
0189    NewElectrode = num_elecs(img)+1;
0190    Y12 = 1/10; % Y=1/10Ohms between E#1 and new electrode
0191    Y13 = 1/ 5; % Y=1/ 5Ohms between E#1 and new elec
0192    c_list = [ 1,3,Y12;
0193               2,NewElectrode,Y13];
0194    % Add instrument electrode
0195    img.fwd_model.electrode(NewElectrode) = ... 
0196        struct('nodes','instrument','z_contact',NaN);
0197    img.fwd_model.stimulation = ...
0198       stim_meas_list([1,3,1,3], ...
0199       num_elecs(img) ); 
0200 
0201    img.fwd_model.system_mat = @system_mat_instrument;
0202    img.fwd_model.system_mat_instrument.connect_list = c_list;
0203    s_mat= calc_system_mat(img); 
0204    save thisimg img
0205 
0206    Eok(:,7) = 0; Eok(7,:) = 0;
0207    Eok = Eok + [ 0 0 0 0 0 0 0;
0208 0 0 0 0   0   0   0; 0 0 0 0   0   0   0; 
0209 0 0 0 0   0   0   0; 0 0 0 0  .1   0 -.1;
0210 0 0 0 0   0  .2 -.2; 0 0 0 0 -.1 -.2  .3];
0211    unit_test_cmp('instrument',s_mat.E,Eok,1e-13);
0212 
0213    imgs = img; imgs.elem_data = [1;1];
0214    vv= fwd_solve(imgs);
0215    Ztest = 1/( Y12 + 1/(Ztest + 1/Y13));
0216    
0217    unit_test_cmp('instrument solve1',vv.meas, Ztest,1e-13);
0218 
0219    % apply current to ground
0220    imgs.fwd_model.stimulation = struct( ...
0221       'stim_pattern',[1;0;NaN], ...
0222       'volt_pattern',[0;0;0], ...
0223       'meas_pattern',[1,0,-1;0,1,-1]);
0224    pp = fwd_model_parameters(img.fwd_model);
0225    vv= fwd_solve(imgs);
0226    Ztest2 = 10/(7.25+10)*5;
0227    unit_test_cmp('instrument solve2',vv.meas, [Ztest;Ztest2],1e-13);
0228 
0229    % apply voltage to ground
0230    imgs.fwd_model.stimulation = struct( ...
0231       'stim_pattern',[NaN;0;NaN], ...
0232       'volt_pattern',[10;0;0], ...
0233       'meas_pattern',[1,0,-1;0,1,-1]);
0234    vv= fwd_solve(imgs);
0235    Ztest = 10*[1;5/7.25];
0236    unit_test_cmp('instrument solve3',vv.meas, Ztest,1e-13);
0237 
0238 
0239 %  disp(Eok)
0240 %  disp(full(s_mat.E))
0241 
0242 function do_unit_test_3d
0243    fmdl = mk_common_model('n3r2',[16,2]);
0244    % get the fwd_model. Remove stims to replace
0245    fmdl = rmfield(fmdl.fwd_model,'stimulation');
0246    gndE= num_elecs(fmdl) + 1;
0247    fmdl.electrode(gndE) = struct('nodes','instrument','z_contact',NaN);
0248    Zgnd = 1e-3;% Ohms -- impedance to ground
0249    VoltIn = 10;
0250    stim.meas_pattern = [eye(gndE-1), -Zgnd*ones(gndE-1,1)];
0251    stim.meas_pattern(1:16,:) = [];
0252 
0253    c_list =  [17:32]'*[1,0,0] + [0,gndE, 1/Zgnd];
0254    fmdl.system_mat = @system_mat_instrument;
0255    fmdl.system_mat_instrument.connect_list = c_list;
0256    for i=1:16
0257       stim.volt_pattern = [zeros(gndE,1)];
0258       stim.stim_pattern = stim.volt_pattern;
0259       stim.volt_pattern(i) = VoltIn;
0260       stim.stim_pattern([i,gndE]) = NaN;
0261       fmdl.stimulation(i) = stim;
0262    end
0263    img = mk_image(fmdl,1);
0264    vv = fwd_solve(img);
0265    unit_test_cmp('3D shape',[max(vv.meas) min(vv.meas)], ...
0266     [ 1.152029995247298e-06, 1.150256267849528e-06], 1e-14);