fwd_solve_halfspace

 FWD_SOLVE_HALFSPACE: data = fwd_solve_halfspace(img)
 Fwd solver for an analytic halfspace
 Input:
    img       = image struct
 Output:
    data = measurements struct
 Options: (to return internal node information)
    img.fwd_solve.get_all_meas = 1 (data.volt = all FEM nodes, but not CEM)
    img.fwd_solve.get_all_nodes= 1 (data.volt = all nodes, including CEM)

 *** Only supports point electrode models (PEM), a CEM will be approximated as
     a PEM at the center of the electrode geometry.

 *** Surface topology is NOT modelled. It is assumed that the electrodes are
     approximately coplanar and on the surface.

 *** The domain is assumed to be homogeneous. The mean conductivity will be used.

 *** Contact impedances z_c are ignored. (TODO)

 (C) 2016 A. Boyle. License: GPL version 2 or version 3
 dipole model: [Loke , "Tutorial: 2-D and 3-D electrical imaging surveys", 2004]
 Ref: TODO

0001 function data = fwd_solve_halfspace(fwd_model, img)
0002 % FWD_SOLVE_HALFSPACE: data = fwd_solve_halfspace(img)
0003 % Fwd solver for an analytic halfspace
0004 % Input:
0005 %    img       = image struct
0006 % Output:
0007 %    data = measurements struct
0008 % Options: (to return internal node information)
0009 %    img.fwd_solve.get_all_meas = 1 (data.volt = all FEM nodes, but not CEM)
0010 %    img.fwd_solve.get_all_nodes= 1 (data.volt = all nodes, including CEM)
0011 %
0012 % *** Only supports point electrode models (PEM), a CEM will be approximated as
0013 %     a PEM at the center of the electrode geometry.
0014 %
0015 % *** Surface topology is NOT modelled. It is assumed that the electrodes are
0016 %     approximately coplanar and on the surface.
0017 %
0018 % *** The domain is assumed to be homogeneous. The mean conductivity will be used.
0019 %
0020 % *** Contact impedances z_c are ignored. (TODO)
0021 %
0022 % (C) 2016 A. Boyle. License: GPL version 2 or version 3
0023 % dipole model: [Loke , "Tutorial: 2-D and 3-D electrical imaging surveys", 2004]
0024 % Ref: TODO
0025 
0026 % correct input paralemeters if function was called with only img
0027 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0028 
0029 if nargin == 1
0030    img= fwd_model;
0031 elseif  strcmp(getfield(warning('query','EIDORS:DeprecatedInterface'),'state'),'on')
0032    error('EIDORS:DeprecatedInterface', ...
0033       'Calling FWD_SOLVE_HALFSPACE with two arguments is deprecated');
0034 end
0035 fwd_model= img.fwd_model;
0036 
0037 img = data_mapper(img);
0038 if ~ismember(img.current_params, supported_params)
0039     error('EIDORS:PhysicsNotSupported', '%s does not support %s', ...
0040     'FWD_SOLVE_HALFSPACE',img.current_params);
0041 end
0042 pp= fwd_model_parameters( fwd_model, 'skip_VOLUME' );
0043 
0044 % all calcs use resistivity
0045 img = convert_img_units(img, 'conductivity');
0046 rho = 1/mean(img.elem_data);
0047 
0048 % create a data structure to return
0049 data.meas= analytic_meas(fwd_model, rho);
0050 data.time= NaN; % unknown
0051 data.name= 'solved by fwd_solve_halfspace';
0052 if isfield(img, 'fwd_solve') && isfield(img.fwd_solve, 'get_all_meas') && img.fwd_solve.get_all_meas
0053    v = analytic_nodal(fwd_model, rho);
0054    data.volt = v(1:pp.n_node,:); % but not on CEM nodes
0055 end
0056 try; if img.fwd_solve.get_all_nodes
0057    error('does not support CEM: can''t get_all_nodes; try get_all_meas');
0058 end; end
0059 
0060 
0061 function do_unit_test
0062    ne = 16;
0063    imdl = mk_geophysics_model('h2p5a', ne);
0064    imdl.fwd_model.stimulation = stim_pattern_geophys(ne, 'Wenner');
0065    img = mk_image(imdl.fwd_model, 1);
0066    img.fwd_solve.get_all_meas = 1;
0067    vh = fwd_solve_halfspace(img);
0068    vd = fwd_solve(img);
0069    scale=1e3; uvh=unique(round(vh.meas*scale)/scale);
0070    unit_test_cmp('h2a halfspace values TEST', uvh, ...
0071   [ 0.0060; 0.0080; 0.0110; 0.0160; 0.0320 ], 1/scale);
0072    unit_test_cmp('h2a halfspace vs default TEST', norm(vh.meas - vd.meas), 0, 4e-3);
0073 clf; h=plot([vh.meas vd.meas],'o--'); legend('analytic','FEM'); set(gca,'box','off'); set(h,'LineWidth',2);
0074 %   img.fwd_solve.get_all_meas = 1;
0075 %   vh = fwd_solve_halfspace(img);
0076 %   gnd_node = img.fwd_model.nodes(img.fwd_model.gnd_node,:);
0077 %clf; show_fem(img); hold on; plot3(gnd_node(1),gnd_node(2),gnd_node(3),'o');
0078 %clf; plot(vh.volt);
0079 
0080 % returns the potential difference
0081 % half-space with PEM electrodes
0082 % [Loke , "Tutorial: 2-D and 3-D electrical imaging surveys", 2004]
0083 function dv = analytic_meas(mdl, rho)
0084    ne = length(mdl.electrode); % number of electrodes
0085    nd = size(mdl.nodes, 2); % dimensions in model
0086    ep = zeros(ne,nd); % init electrode position
0087    for i = 1:ne
0088       en = mdl.electrode(i).nodes;
0089       ep(i,:) = mean(mdl.nodes(en), 2);
0090    end
0091    list = stim_meas_list(mdl.stimulation);
0092    a = list(:,1); % s+
0093    b = list(:,2); % s-
0094    m = list(:,3); % m+
0095    n = list(:,4); % m-
0096    AM = dist(ep(a), ep(m));
0097    BM = dist(ep(b), ep(m));
0098    AN = dist(ep(a), ep(n));
0099    BN = dist(ep(b), ep(n));
0100 % TODO  zc = [mdl.electrode(:).z_contact]';
0101 % TODO  MNzc = zc(m) + zc(n);
0102 
0103    k = geometric_factor(AM, BM, AN, BN);
0104    I = stim_current(mdl);
0105    dv = (rho*I)./k; % TODO AB missing zc
0106 
0107 % returns the potential everywhere relative to the ground node
0108 % half-space with PEM electrodes
0109 % [Loke , "Tutorial: 2-D and 3-D electrical imaging surveys", 2004]
0110 function v = analytic_nodal(mdl, rho)
0111    ne = length(mdl.electrode); % number of electrodes
0112    ns = length(mdl.stimulation); % number of stim patterns
0113    nn = size(mdl.nodes, 1); % nodes in the model
0114    nd = size(mdl.nodes, 2); % dimensions in model
0115    ep = zeros(ne,nd); % init electrode position
0116    for i = 1:ne
0117       en = mdl.electrode(i).nodes;
0118       ep(i,:) = mean(mdl.nodes(en), 2);
0119    end
0120    a = ones(ns,1)*nan; b = a; I = a;
0121    for i=1:ns % same order as pp.QQ
0122       s = mdl.stimulation(i).stim_pattern;
0123       assert_sane_stim(mdl.stimulation(i));
0124       a(i) = find(s > 0);
0125       b(i) = find(s < 0);
0126       I(i) = sum(s(s > 0));
0127    end
0128    epm = mdl.nodes(mdl.gnd_node); % ground node
0129    epn = mdl.nodes; % all nodes
0130    AM = dist(ep(a), epm);
0131    BM = dist(ep(b), epm);
0132    abN = dist(kron(ep,ones(nn,1)), repmat(epn,ne,1)); % every electrode to every node
0133    abN = reshape(abN, nn, ne); % cols = electrode#, rows = node#
0134 % TODO  zc = [mdl.electrode(:).z_contact]';
0135 % TODO  MNzc = zc(m) + zc(n);
0136 
0137    v = ones(nn,length(I));
0138    for i=1:ne
0139       ai = a(i); bi = b(i);
0140       k = geometric_factor(AM(ai), BM(bi), abN(:,ai), abN(:,bi));
0141       v(:,i) = rho*I(i)./k; % TODO AB missing zc
0142    end
0143 
0144 % transfer resistance:" stimulus current over measured voltage
0145 function [I, a, b] = stim_current(mdl)
0146    stim = mdl.stimulation;
0147    I = [];
0148    N = 0;
0149    for i = 1:length(stim)
0150       s = stim(i).stim_pattern;
0151       assert_sane_stim(stim(i));
0152       n = size(stim(i).meas_pattern, 1);
0153       I(N+[1:n]) = sum(s(s > 0));
0154       N = N + n;
0155    end
0156    I=I';
0157 
0158 function assert_sane_stim(stim)
0159    assert(length(stim) == 1); % only one stim pattern per "assert_sane" check
0160    s = stim.stim_pattern;
0161    assert(strcmp(stim.stimulation, 'Amp')); % scaling?
0162    assert(length(find(abs(s) > 0)) == 2); % pairwise stimulus only
0163    assert(sum(s) == 0); % dipole
0164    assert(size(s,2) == 1); % single dipole per stimulation(i)
0165 
0166 % distance from A to B
0167 function AB = dist(A, B)
0168    AB = sqrt(sum((B - A).^2,2));
0169 
0170 % analytic model, for half-space with point electrodes (PEM)
0171 function k = geometric_factor(AM, BM, AN, BN)
0172    k = 2*pi./(1./AM-1./BM-1./AN+1./BN);