mk_geophysics_model

 imdl = mk_geophysics_model(str, ne, [option])

 ne  - number of electrodes, 5 metre spacing (+5,+10,...)
       and 0.1 metre diameter
         OR
       a list of electrode locations in the x-dimension or a 2- or
       3-dimensional array, one electrode per row, missing columns
       will be set to zero
       ne = 16
       ne = [4 6 10 20] % 1d: (x)
       ne = [0 1; 2 1.5; 3 1.2; 7 2.5] % 2d: (x,y)
       ne = [0 0.1 1; 2 -0.1 1.5; 3 -0.15 1.2; 7 0 2.5] % 3d: (x,y,z)
 str - model, x = see hmax_rec
       h2x -   2D half-space, linear CEM array (2d fwd)
       h2p5x - 2.5D half-space, linear CEM array (2d fwd + Fourier y-dimension)
       h3x   - 3D half-space, linear CEM array (3d fwd)
       h22x  - 2D half-space, linear CEM array (2d fwd, 2d rec)
       h32x  - 2D half-space, linear CEM array (3d fwd, 2d rec)
       h33x  - 2D half-space, linear CEM array (3d fwd, 3d rec)
       H2x -   2D half-space, linear PEM array (2d fwd)
       H2p5x - 2.5D half-space, linear PEM array (2d fwd + Fourier y-dimension)
       H3x   - 3D half-space, linear PEM array (3d fwd)
       H22x  - 2D half-space, linear PEM array (2d fwd, 2d rec)
       H32x  - 2D half-space, linear PEM array (3d fwd, 2d rec)
       H33x  - 2D half-space, linear PEM array (3d fwd, 3d rec)
 opt - override default configuration options (optional cell array)
       'hmax_fwd' - fine reconstruction mesh density, given an
                    array width xw, and electrode spacing es
                    Note that for ne=16, 'A' and 'a' are equivalent.
                ['a' : hmax_fwd=xw/1;    ['A' : hmax_fwd=es*16;
                 'b' : hmax_fwd=xw/2;     'B' : hmax_fwd=es*8;
                 'c' : hmax_fwd=xw/4;     'C' : hmax_fwd=es*4;
                 ...                       ...
                 'z' : hmax_fwd=xw/2^25]  'Z' : hmax_fwd=es*2^-21]
       'hmax_fwd_inner'
                  - reconstruction model mesh density for the inner region
                    [default: 1/2 of the outer region density hmax_fwd]
       'hmax_rec' - reconstruction model mesh density [hmax_fwd*2]
       'elec_width' - electrode width [0.1 m]
                    width = 0 requests a PEM, rather than CEM
       'z_contact' - electrode contact impedance [0.01 \Ohm.m]
       'elec_spacing' - distance between electrode centers [5 m]
       'extend_x' - extra mesh in the principle axis of the
                    electrode array, multiple of array width [1]
       'extend_y' - extra mesh in the minor axis of the
                    electrode array, multiple of array width
                    (3D models only) [1]
       'extend_z' - extra depth of model, multiple of array width [1]
       'extend_inner_x'
                  - inner (denser) mesh, multiple of array width [3/5]
       'extend_inner_y'
                  - inner (denser) mesh, multiple of array width [3/5]
                    (3D models only)
       'extend_inner_z'
                  - inner (denser) mesh, multiple of array width [2/5]
       'skip_c2f' - skip building the rec_model to fwd_model mapping [0]
       'threshold' - threshold for electrode placement errors [1e-12]
                    if error exceeds threshold, mash nodes by cubic interp
                    maintaining the side and lower boundaries of the
                    model... the nodes will be mashed until the electrodes
                    conform to the requested electrode positions,
                    correcting for Netgen inaccuracies and lack of
                    vertical profile (Inf = disable node mashing)
       'build_stim' - use stim_pattern_geophys to build a standard geophysics
                      stim/meas sequence and add it to the model, based on the
                      number of electrodes, and assuming a co-linear array;
                      set to 'none' to skip this step
                      [default: 'Wenner']
       'extra_ng_code' - pass extra netgen code, any 'tlo' are subtracted from
                      the inner region, units are scaled, so left to right-most
                      electrodes are -1 to +1
                      [default: {}]

 The linear electrode array runs in the +X direction at Z=0. For
 the 3D model, the Y-axis is perpendicular to the electrode array.

 (C) 2015--2017 A. Boyle
 License: GPL version 2 or version 3

0001 function imdl = mk_geophysics_model(str, ne, opt);
0002 % imdl = mk_geophysics_model(str, ne, [option])
0003 %
0004 % ne  - number of electrodes, 5 metre spacing (+5,+10,...)
0005 %       and 0.1 metre diameter
0006 %         OR
0007 %       a list of electrode locations in the x-dimension or a 2- or
0008 %       3-dimensional array, one electrode per row, missing columns
0009 %       will be set to zero
0010 %       ne = 16
0011 %       ne = [4 6 10 20] % 1d: (x)
0012 %       ne = [0 1; 2 1.5; 3 1.2; 7 2.5] % 2d: (x,y)
0013 %       ne = [0 0.1 1; 2 -0.1 1.5; 3 -0.15 1.2; 7 0 2.5] % 3d: (x,y,z)
0014 % str - model, x = see hmax_rec
0015 %       h2x -   2D half-space, linear CEM array (2d fwd)
0016 %       h2p5x - 2.5D half-space, linear CEM array (2d fwd + Fourier y-dimension)
0017 %       h3x   - 3D half-space, linear CEM array (3d fwd)
0018 %       h22x  - 2D half-space, linear CEM array (2d fwd, 2d rec)
0019 %       h32x  - 2D half-space, linear CEM array (3d fwd, 2d rec)
0020 %       h33x  - 2D half-space, linear CEM array (3d fwd, 3d rec)
0021 %       H2x -   2D half-space, linear PEM array (2d fwd)
0022 %       H2p5x - 2.5D half-space, linear PEM array (2d fwd + Fourier y-dimension)
0023 %       H3x   - 3D half-space, linear PEM array (3d fwd)
0024 %       H22x  - 2D half-space, linear PEM array (2d fwd, 2d rec)
0025 %       H32x  - 2D half-space, linear PEM array (3d fwd, 2d rec)
0026 %       H33x  - 2D half-space, linear PEM array (3d fwd, 3d rec)
0027 % opt - override default configuration options (optional cell array)
0028 %       'hmax_fwd' - fine reconstruction mesh density, given an
0029 %                    array width xw, and electrode spacing es
0030 %                    Note that for ne=16, 'A' and 'a' are equivalent.
0031 %                ['a' : hmax_fwd=xw/1;    ['A' : hmax_fwd=es*16;
0032 %                 'b' : hmax_fwd=xw/2;     'B' : hmax_fwd=es*8;
0033 %                 'c' : hmax_fwd=xw/4;     'C' : hmax_fwd=es*4;
0034 %                 ...                       ...
0035 %                 'z' : hmax_fwd=xw/2^25]  'Z' : hmax_fwd=es*2^-21]
0036 %       'hmax_fwd_inner'
0037 %                  - reconstruction model mesh density for the inner region
0038 %                    [default: 1/2 of the outer region density hmax_fwd]
0039 %       'hmax_rec' - reconstruction model mesh density [hmax_fwd*2]
0040 %       'elec_width' - electrode width [0.1 m]
0041 %                    width = 0 requests a PEM, rather than CEM
0042 %       'z_contact' - electrode contact impedance [0.01 \Ohm.m]
0043 %       'elec_spacing' - distance between electrode centers [5 m]
0044 %       'extend_x' - extra mesh in the principle axis of the
0045 %                    electrode array, multiple of array width [1]
0046 %       'extend_y' - extra mesh in the minor axis of the
0047 %                    electrode array, multiple of array width
0048 %                    (3D models only) [1]
0049 %       'extend_z' - extra depth of model, multiple of array width [1]
0050 %       'extend_inner_x'
0051 %                  - inner (denser) mesh, multiple of array width [3/5]
0052 %       'extend_inner_y'
0053 %                  - inner (denser) mesh, multiple of array width [3/5]
0054 %                    (3D models only)
0055 %       'extend_inner_z'
0056 %                  - inner (denser) mesh, multiple of array width [2/5]
0057 %       'skip_c2f' - skip building the rec_model to fwd_model mapping [0]
0058 %       'threshold' - threshold for electrode placement errors [1e-12]
0059 %                    if error exceeds threshold, mash nodes by cubic interp
0060 %                    maintaining the side and lower boundaries of the
0061 %                    model... the nodes will be mashed until the electrodes
0062 %                    conform to the requested electrode positions,
0063 %                    correcting for Netgen inaccuracies and lack of
0064 %                    vertical profile (Inf = disable node mashing)
0065 %       'build_stim' - use stim_pattern_geophys to build a standard geophysics
0066 %                      stim/meas sequence and add it to the model, based on the
0067 %                      number of electrodes, and assuming a co-linear array;
0068 %                      set to 'none' to skip this step
0069 %                      [default: 'Wenner']
0070 %       'extra_ng_code' - pass extra netgen code, any 'tlo' are subtracted from
0071 %                      the inner region, units are scaled, so left to right-most
0072 %                      electrodes are -1 to +1
0073 %                      [default: {}]
0074 %
0075 % The linear electrode array runs in the +X direction at Z=0. For
0076 % the 3D model, the Y-axis is perpendicular to the electrode array.
0077 %
0078 % (C) 2015--2017 A. Boyle
0079 % License: GPL version 2 or version 3
0080 
0081 % model: 64 electrode, 2d half-space
0082 % Once upon a time, this code started out from the following tutorial.
0083 % model from http://eidors3d.sourceforge.net/tutorial/other_models/square_mesh.shtml
0084 
0085 if ischar(str) && strcmp(str,'UNIT_TEST'); do_unit_test; return; end
0086 copt.fstr = 'mk_geophysics_model';
0087 if nargin < 3
0088    opt = {};
0089 end
0090 SALT='z$Id: mk_geophysics_model.m 6521 2022-12-30 21:33:17Z bgrychtol $'; % stick a key in the model 'save' file, so we can expire them when the model definitions age
0091 imdl = eidors_cache(@mk_model,{str, ne, opt, SALT}, copt);
0092 
0093 function imdl=mk_model(str,ne,opt,SALT);
0094 if str(1) ~= 'h' && str(1) ~= 'H'
0095    error([str ': I only know how to build linear half-space models: h***']);
0096 end
0097 
0098 MDL_2p5D_CONFIG = 0;
0099 switch str(2:end-1)
0100    case {'2', '3'} % simple meshes
0101       FMDL_DIM = str(2) - '0';
0102       CMDL_DIM = 0; % no cmdl
0103    case '2p5' % 2.5D Fourier transformed
0104       FMDL_DIM = 2;
0105       CMDL_DIM = 0; % no cmdl
0106       MDL_2p5D_CONFIG = 1;
0107    case {'22', '33', '32'} % dual meshes
0108       FMDL_DIM = str(2) - '0';
0109       CMDL_DIM = str(3) - '0';
0110    otherwise
0111       error([str ': unrecognized model type']);
0112 end
0113 assert(CMDL_DIM ~= 3, '3d rec_model not yet tested');
0114 
0115 skip_c2f = 0;
0116 if str(1) == 'h'
0117    elec_width = 0.1;
0118 else % str(1) == 'H'
0119    elec_width = 0;
0120 end
0121 z_contact= 0.01;
0122 nodes_per_elec= 3; %floor(elec_width/hmax_rec*10);
0123 elec_spacing= 5.0;
0124 threshold = 1e-12;
0125 save_model_to_disk = (length(ne) == 1) && (length(opt) == 0);
0126 
0127 extend_x = 1;
0128 extend_y = 1;
0129 extend_z = 1;
0130 extend_inner_x = 3/5;
0131 extend_inner_y = 3/5;
0132 extend_inner_z = 2/5;
0133 build_stim = 'Wenner';
0134 extra_ng_code = '';
0135 if length(opt) > 0 % allow overriding the default values
0136    assert(round(length(opt)/2)*2 == length(opt),'option missing value?');
0137    expect = {'hmax_rec','hmax_fwd', 'hmax_fwd_inner', ...
0138              'elec_width','z_contact','elec_spacing',...
0139              'extend_x', 'extend_y', 'extend_z', ...
0140              'extend_inner_x', 'extend_inner_y', 'extend_inner_z', ...
0141              'skip_c2f', 'threshold', 'build_stim', 'extra_ng_code'};
0142    opts = struct(opt{:})
0143    for i = fieldnames(opts)'
0144       assert(any(strcmp(i,expect)), ['unexpected option: ',i{:}]);
0145       eval([i{:} ' = opts.(i{:});']);
0146    end
0147    if (str(1) == 'H') && isfield(opts, 'elec_width')
0148       error('requested "H" PEM model but configured "elec_width" option');
0149    end
0150 end
0151 if length(ne) == 1 % ne: number of electrodes
0152    xw=(ne-1)*elec_spacing; % array width
0153    %xs=-(ne-1)*elec_spacing/2; % array centered
0154    xs=+5; % array at left-most at +5
0155    xyz = xs+([1:ne]'-1)*elec_spacing;
0156 else
0157    xyz = ne; % must be a set of coordinates for the electrodes...
0158 end
0159 if size(xyz,1) == 1
0160    xyz = xyz'; % flip to column
0161 end
0162 xyz = [xyz zeros(size(xyz,1),3-size(xyz,2))]; % [x 0 0] or [ x y 0 ] or [ x y z ]
0163 ne=size(xyz,1);
0164 [R, X] = rot_line_to_xaxis(xyz);
0165 % % rescale, centre electrodes so NetGen can be happy
0166 xyzc = (xyz - X)*R; % centre and scale electrodes: -1 to +1 y-axis
0167 % xyzc = xyz; % centre and scale electrodes: -1 to +1 y-axis
0168 xw=max(xyzc(:,1))-min(xyzc(:,1));
0169 xs=min(xyzc(:,1));
0170 elec_spacing = min(min(pdist(xyzc) + diag(inf*(1:size(xyzc,1))))); % min spacing btw elec
0171 
0172 if ~exist('hmax_fwd','var')
0173    if str(end)-'a' >= 0
0174       hmax_fwd = xw*2^-(str(end)-'a');
0175    else
0176       hmax_fwd = elec_spacing*2^-(str(end)-'A'-4);
0177    end
0178 end
0179 if ~exist('hmax_rec','var') % allow hmax_rec to depend on configured hmax_fwd
0180    hmax_rec=hmax_fwd*2.01; % avoid parametrization aliasing
0181 end
0182 if ~exist('hmax_fwd_inner','var') % allow hmax_rec to depend on configured hmax_fwd
0183    hmax_fwd_inner=hmax_fwd/2.0;
0184 end
0185 
0186 if save_model_to_disk
0187    isOctave = exist('OCTAVE_VERSION');
0188    if isOctave
0189       octavestr='-octave-';
0190    else
0191       octavestr='-';
0192    end
0193    % NOTE models are stored in the directory specified by eidors_cache('cache_path')
0194    filename=fullfile(eidors_cache('cache_path'),sprintf('mk_geophysics_model%simdl-%s-%03del.mat',octavestr,str,ne));
0195    if exist(filename, 'file') == 2
0196       tmp = whos('-file',filename,'SALT');
0197       if length(tmp) > 0
0198          tmp = load(filename,'SALT');
0199          fSALT = tmp.SALT;
0200       else
0201          fSALT = 'deadbeef';
0202       end
0203       if strcmp(fSALT, SALT)
0204          tmp=load(filename,'imdl');
0205          imdl = tmp.imdl;
0206          eidors_msg(sprintf('%s: %s, %d electrode model loaded from file',filename,str,ne));
0207          return
0208       end % hmm, the SALT doesn't match so we go back to generating a new model from scratch
0209    end
0210 end
0211 
0212 assert(extend_x>0,'extend_x must be > 0');
0213 assert(extend_y>0,'extend_y must be > 0');
0214 assert(extend_z>-1,'extend_z must be > -1');
0215 
0216 % Calculate cmdl, fmdl and inner fmdl (fmdlin) min/max coordinates.
0217 % After rescaling (normalizing), electrode major axis is along the x-axis, and
0218 % the electrode array will be 2 units long (-1,+1).
0219 fmdlbox =   [-(1+2*extend_x)        +(1+2*extend_x);
0220              -(1+2*extend_y)        +(1+2*extend_y);
0221              -(2+2*extend_z)        +3];
0222 fmdlinbox = [-(1+2*extend_inner_x)  +(1+2*extend_inner_x);
0223              -(1+2*extend_inner_y)  +(1+2*extend_inner_y);
0224              -(2+2*extend_inner_z)  +2];
0225 cmdlbox = fmdlbox;
0226 
0227 assert(all(fmdlbox(:,1) <= fmdlinbox(:,1)), 'oops, inner mesh must be smaller than outer mesh');
0228 assert(all(fmdlbox(:,2) >= fmdlinbox(:,2)), 'oops, inner mesh must be smaller than outer mesh');
0229 assert(all(fmdlbox(:,1) <= cmdlbox(:,1)), 'oops, reconstruction mesh must be smaller than forward mesh');
0230 assert(all(fmdlbox(:,2) >= cmdlbox(:,2)), 'oops, reconstruction mesh must be smaller than forward mesh');
0231 
0232 % build shape string for NetGen
0233 if elec_width ~= 0
0234    elec_shape = [elec_width*norm(R)/2, 0, elec_width*norm(R)/2/(nodes_per_elec-1)]; % CEM, circular, maxh
0235    elec_pos   = [ xyzc(:,1:FMDL_DIM), repmat([zeros(1,3-FMDL_DIM+2) 1],ne,1) ]; % p(x,y,z=0), n(0,0,1)
0236    cem2pem = 0;
0237 else
0238    if FMDL_DIM == 3
0239       elec_width = elec_spacing/2;
0240       elec_shape = [elec_width, elec_width, hmax_fwd_inner]; % PEM, sz, maxh
0241       elec_pos   = [ xyzc, repmat([0 0 1],ne,1) ]; % p(x,y,z=0), n(0,0,1)
0242       elec_pos(:,1) = elec_pos(:,1) + elec_width/2;
0243       elec_pos(:,2) = elec_pos(:,2) + elec_width/2;
0244       cem2pem = 1;
0245    else % FMDL_DIM == 2
0246       elec_width = elec_spacing/2;
0247       elec_shape = [elec_width, elec_width, hmax_fwd_inner]; % rectangular CEM->PEM, maxh
0248       elec_pos   = [ xyzc(:,1:2), repmat([0 0 0 1],ne,1) ]; % p(x,y,z=0), n(0,0,1)
0249       elec_pos(:,1) = elec_pos(:,1) + elec_width/2;
0250       cem2pem = 1;
0251    end
0252 end
0253 skip_tlo = '';
0254 for idx = strfind(extra_ng_code, 'tlo ')
0255    sc = strfind(extra_ng_code, ';'); % semicolons
0256    sc = min(sc + (sc<idx)*1e10);
0257    tlo = extra_ng_code(idx+4:sc-1);
0258    skip_tlo = [skip_tlo ' and (not ' tlo ')'];
0259 end
0260 elec_pos(:,FMDL_DIM) = 0; % flatten electrode positions onto the "ps" plane
0261 if FMDL_DIM == 3
0262    shape_str = [...
0263                 sprintf('solid ps = plane(%s);\n', a2s([0 0 0; 0 0 1])), ...
0264                 sprintf('solid bi = orthobrick(%s);\n', a2s(fmdlinbox')), ...
0265                 sprintf('solid bo = orthobrick(%s);\n', a2s(fmdlbox')), ...
0266                 extra_ng_code, ...
0267                 sprintf('solid ri = bi and ps%s -maxh=%f;\n', skip_tlo, hmax_fwd_inner), ...
0268                 sprintf('solid ro = bo and ps and (not bi) -maxh=%f;\n', hmax_fwd), ...
0269                 sprintf('solid mainobj = ri;\n')];
0270    % Note that ri must be the 'mainobj' so that it can intersect with the electrodes
0271    % additional top level objects for netgen
0272 else % netgen 2d model
0273    % to create the 2D slice we need to give NetGen something to work with
0274    %Need some width to let netgen work, but not too much so
0275    % that it meshes the entire region
0276    sw = range(fmdlinbox(1,:)) / 5; % initial estimate
0277    sw = min(sw,2*hmax_fwd); % coarse model slice width
0278    ri2d = fmdlinbox'; ri2d(:,2) = [-sw 0 ];
0279    ro2d = fmdlbox'; ro2d(:,2) = [-sw 0 ];
0280    shape_str = [...
0281                 sprintf('solid ps = plane(%s);\n', a2s([0 0 0; 0 0 1])), ...
0282                 sprintf('solid bi = orthobrick(%s);\n', a2s(ri2d)), ...
0283                 sprintf('solid bo = orthobrick(%s);\n', a2s(ro2d)), ...
0284                 extra_ng_code, ...
0285                 sprintf('solid ri = bi and ps%s -maxh=%f;\n', skip_tlo, hmax_fwd_inner), ...
0286                 sprintf('solid ro = bo and ps and (not bi) -maxh=%f;\n', hmax_fwd), ...
0287                 sprintf('solid mainobj = ri;\n')];
0288 end
0289 % fprintf('SHAPE_STR: %s', shape_str); elec_pos
0290 elec_obj = 'ps';
0291 [fmdl, mat_idx] = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj, 'tlo ro');
0292 if FMDL_DIM == 2 % 2D
0293    % now convert the roughly 2D slice into a true 2D plane
0294    [fmdl, mat_idx] = copy_mdl2d_from3d(fmdl, mat_idx, 'y');
0295 else % 3D
0296    if CMDL_DIM ~= 0
0297       % c2f
0298       cmdl.mk_coarse_fine_mapping.f2c_offset  = [0 0 0];
0299       cmdl.mk_coarse_fine_mapping.f2c_project = [1 0 0; 0 0 1; 0 1 0];
0300    end
0301 end
0302 
0303 if cem2pem
0304    fmdl = convert_cem2pem(fmdl, xyzc);
0305 end
0306 
0307 % 2d cmdl
0308 xllim = fmdlbox(1,1);
0309 xrlim = fmdlbox(1,2);
0310 zdepth = fmdlbox(3,1);
0311 % AA - 2022-05 Updated from hmax_rec/2
0312 xr=max(floor((xrlim-xllim)/hmax_rec),1)*2+1; % odd number
0313 yr=max(floor(-zdepth/hmax_rec),1)*2+1; % odd number
0314 [x,y] = meshgrid( linspace(xllim,xrlim,xr), linspace(zdepth,0,yr));
0315 vtx= [x(:),y(:)];
0316 if CMDL_DIM ~= 0
0317    cmdl= mk_fmdl_from_nodes( vtx,{vtx(1,:)}, z_contact, 'sq_m2');
0318 end
0319 
0320 % stick electrode nodes into cmdl so that show_fem will plot them
0321 for i=1:ne
0322    n=fmdl.electrode(i).nodes;
0323    nn=length(n);
0324    nx=fmdl.nodes(n,:);
0325 
0326    fmdl.electrode(i).z_contact = z_contact;
0327    if CMDL_DIM ~= 0
0328       nnc = length(cmdl.nodes);
0329       cmdl.nodes = [cmdl.nodes; nx(:,[1 FMDL_DIM])];
0330       cmdl.electrode(i).nodes = (nnc+1):(nnc+nn);
0331       cmdl.electrode(i).z_contact = z_contact;
0332    end
0333 end
0334 
0335 % fix electrode locations if necessary
0336 elec_err = sqrt(sum(mdl_elec_err(fmdl, xyzc).^2,2));
0337 if max(elec_err) > threshold % put electrodes in the right place
0338    [fmdl, cf] = correct_electrode_positions(fmdl, xyzc);
0339 
0340    if CMDL_DIM ~= 0
0341       [cmdl, cc] = correct_electrode_positions(cmdl, xyzc);
0342    end
0343 end
0344 % save functions for later use
0345 fmdl.mv_elec = @shift_electrode_positions;
0346 
0347 % reverse the centre and scaling
0348 nn = size(fmdl.nodes,1);
0349 Xn = repmat(X(1,:), nn, 1);
0350 fmdl.nodes = ([fmdl.nodes zeros(nn,3-FMDL_DIM)]/ R) + Xn;
0351 fmdl.nodes = fmdl.nodes(:,1:FMDL_DIM);
0352 
0353 if CMDL_DIM ~= 0
0354    nn = size(cmdl.nodes,1);
0355    Xn = repmat(X(1,:), nn, 1);
0356    if CMDL_DIM ~= FMDL_DIM % 2D
0357       cmdl.nodes = ([cmdl.nodes(:,1) zeros(nn,1) cmdl.nodes(:,2:end)]/ R) + Xn;
0358       cmdl.nodes = cmdl.nodes(:,[1 3]);
0359 %     cmdl.nodes = cmdl.nodes(:,[1 3])/R([1 3],[1 3]); % ONLY DIVIDE ONCE
0360    else % CMDL_DIM == FMDL_DIM
0361       cmdl.nodes = ([cmdl.nodes zeros(nn,3-CMDL_DIM)]/ R) + Xn;
0362       cmdl.nodes = cmdl.nodes(:,1:CMDL_DIM);
0363    end
0364 end
0365 
0366 % check that all electrodes were found
0367 for i = 1:length(fmdl.electrode)
0368    nn = fmdl.electrode(i).nodes;
0369    assert(length(nn(:)) > 0, sprintf('electrode#%d: failed to find nodes',i));
0370 end
0371 
0372 % show_fem(fmdl); xlabel('x'); ylabel('y'); zlabel('z');
0373 
0374 % TODO ... Actually we want a node in the middle: the boundary is bad too...
0375 [~, gn] = min(fmdl.nodes(:,end));
0376 fmdl.gnd_node = gn; % make sure the ground node is away from surface electrodes
0377 if CMDL_DIM ~= 0
0378    [~, gn] = min(cmdl.nodes(:,end));
0379    cmdl.gnd_node = gn; % make sure the ground node is away from surface electrodes
0380 end
0381 
0382 
0383 imdl= mk_common_model('a2d0c',4); % 2d model
0384 imdl.fwd_model = fmdl;
0385 imdl.name = ['EIDORS mk_geophysics_model ' str];
0386 imdl.fwd_model.name = ['EIDORS mk_geophysics_model fwd model ' str];
0387 if CMDL_DIM ~= 0
0388    imdl.rec_model.name = ['EIDORS mk_geophysics_model rec model ' str];
0389    imdl.rec_model = cmdl;
0390    % EIDORS "analytic_c2f" gets stuck, do an approximate one
0391    %eidors_default('set','mk_coarse_fine_mapping','mk_analytic_c2f');
0392    %eidors_default('set','mk_coarse_fine_mapping','mk_approx_c2f');
0393    %[c2f,out] = mk_coarse_fine_mapping(fmdl,cmdl);
0394    if ~skip_c2f
0395       [c2f,out] = mk_approx_c2f(fmdl,cmdl);
0396       imdl.fwd_model.coarse2fine = c2f;
0397       imdl.fwd_model.background = out;
0398    end
0399 end
0400 
0401 if MDL_2p5D_CONFIG
0402    imdl.fwd_model.jacobian   = @jacobian_adjoint_2p5d_1st_order;
0403    imdl.fwd_model.solve      = @fwd_solve_2p5d_1st_order;
0404    imdl.fwd_model.system_mat = @system_mat_2p5d_1st_order;
0405    imdl.fwd_model.jacobian_adjoint_2p5d_1st_order.k = [0 3];
0406    imdl.fwd_model.fwd_solve_2p5d_1st_order.k = [0 3];
0407 end
0408 
0409 if ~strcmp(build_stim,'none')
0410    imdl.fwd_model.stimulation = stim_pattern_geophys(length(imdl.fwd_model.electrode), build_stim);
0411 end
0412 
0413 if save_model_to_disk
0414    save(filename,'imdl','SALT');
0415 end
0416 
0417 % convert 2x3 array to "x1,y1,z1;x2,y2,z2" string
0418 % convert 1x3 array to "x1,y1,z1" string
0419 % if 1x1 array, then use xy_ctr=x1,y1 and [0 0 1]=x2,y2,z2 (z+)
0420 function s = a2s(a)
0421 if length(a(:)) == 3
0422    s = sprintf('%f,%f,%f', ...
0423                a(1), a(2), a(3));
0424 else
0425    s = sprintf('%f,%f,%f;%f,%f,%f', ...
0426                 a(1,1), a(1,2), a(1,3), ...
0427                a(2,1), a(2,2), a(2,3));
0428 end
0429 
0430 % returns R rotation/scaling and X0 offset
0431 % xyz1 = R * xyz + X; % rotate and scale to +/- 1
0432 function [R,X,var] = rot_line_to_xaxis(xyz)
0433 x = xyz(:,1); y=xyz(:,2); z=xyz(:,3);
0434 
0435 % fit line to points
0436 % http://www.mathworks.com/matlabcentral/newsreader/view_thread/32502
0437 p = mean(xyz);
0438 [U,S,V] = svd([x-p(1), y-p(2), z-p(3)]);
0439 if V(end,1) ~= 0
0440    N=1/V(end,1)*V(:,1);
0441 else
0442    N=V(:,1);
0443 end
0444 A=p' + dot( xyz(1,  :) - p, N ) * N/norm(N)^2;
0445 B=p' + dot( xyz(end,:) - p, N ) * N/norm(N)^2;
0446 
0447 % rotate line to +y-axis
0448 a = N/norm(N);
0449 b = [ 1 0 0 ]'; % +x
0450 % http://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d
0451 v = cross(a, b);
0452 s = norm(v); c = dot(a, b); % sin, cos
0453 xt = @(v) [   0  -v(3)  v(2); ... % skew symmetric cross-product of v
0454             v(3)    0  -v(1); ... % diagnoal is the scaling identity matrix
0455            -v(2)  v(1)    0];
0456 if abs(s) < eps*1e3
0457    R = eye(3);
0458 else
0459    R = (eye(3) + xt(v) + xt(v)^2 * (1-c)/s^2);
0460    R=R'; % for right multiply: xyz*R
0461 end
0462 
0463 X = repmat(p, size(xyz,1), 1);
0464 R = R / max(range((xyz-X)*R)) * 2;
0465 
0466 DEBUG=0;
0467 if DEBUG
0468    clf;
0469    subplot(211);
0470    plot3(x,y,z,'bo');
0471    hold on;
0472    plot3(x(1),y(1),z(1),'go');
0473    plot3(p(1),p(2),p(3),'ro');
0474    plot3([A(1),B(1)],[A(2),B(2)],[A(3),B(3)]);
0475    grid; axis equal; hold off;
0476    xlabel('x'); ylabel('y'); zlabel('z');
0477    title('pre-rotation');
0478    subplot(212);
0479    xx = (xyz-X)*R;
0480    plot3(xx(:,1), xx(:,2), xx(:,3),'bo');
0481    hold on;
0482    plot3(xx(1,1), xx(1,2), xx(1,3),'go');
0483    plot3(0,0,0,'ro');
0484    grid; axis equal; hold off;
0485    xlabel('x'); ylabel('y'); zlabel('z');
0486    title('post-rotation');
0487 end
0488 
0489 % the (max-min) range of a variable's values
0490 function r = range(a)
0491 r = max(a(:))-min(a(:));
0492 
0493 function [mdl2,idx2] = copy_mdl2d_from3d(mdl3,idx3,sel);
0494 % AB: taken from EIDORS function ng_mk_gen_models() subfunction of the same name
0495 % AB: NEW: sel = 'x', 'y' or 'z' -- default was Z, we want X
0496    if sel == 'x'
0497       % swap Z and X
0498       T = [ 0 0 1; 0 1 0; 1 0 0 ];
0499    elseif sel == 'y'
0500       % swap Z and Y
0501       T = [ 1 0 0; 0 0 1; 0 1 0 ];
0502    elseif sel == 'z'
0503       T = eye(3);
0504    else
0505       error('sel must be "x", "y" or "z"');
0506    end
0507    mdl3.nodes = mdl3.nodes * T; % AB: SWAP axes
0508 
0509    % set name
0510    mdl2 = eidors_obj('fwd_model',sprintf('%s 2D',mdl3.name));
0511 
0512    % set nodes
0513    [bdy,idx] = find_boundary(mdl3.elems);
0514    vtx = mdl3.nodes;
0515    z_vtx = reshape(vtx(bdy,3), size(bdy) );
0516    z_vtx_thres = max(z_vtx(:))-10*eps*range(z_vtx(:));
0517    lay0  = find( all(z_vtx >= z_vtx_thres, 2) );
0518    bdy0  = bdy( lay0, :);
0519 
0520    vtx0  = unique(bdy0(:));
0521    mdl2.nodes = vtx(vtx0,1:2);
0522 
0523    % set elems
0524    nmap  = zeros(size(vtx,1),1); nmap(vtx0) = 1:length(vtx0);
0525    bdy0  = reshape(nmap(bdy0), size(bdy0) ); % renumber to new scheme
0526    mdl2.elems = bdy0;
0527 
0528    % set boundary
0529    mdl2.boundary = find_boundary( mdl2.elems);
0530 
0531    % set gnd_node
0532    mdl2.gnd_node = nmap(mdl3.gnd_node);
0533 
0534    % set material indices
0535    % TODO: vectorize code
0536    idx2 = {};
0537    idx0  = idx( lay0, :);
0538    for i=1:size(idx3,2)
0539       idx2{i} = [];
0540       ii = 1;
0541       for j=1:size(idx3{i},1)
0542          idx_tmp = find( idx0==idx3{i}(j) );
0543          if not(isempty(idx_tmp))
0544             idx2{i}(ii,1) = idx_tmp(1,1);
0545             ii = ii + 1;
0546          end
0547       end
0548    end
0549 
0550    % set electrode
0551    if isfield(mdl3,'electrode')
0552       mdl2.electrode = mdl3.electrode;
0553       for i=1:length(mdl2.electrode);
0554          nn = mdl3.nodes(mdl3.electrode(i).nodes,:);
0555          enodes = nmap( mdl2.electrode(i).nodes );
0556          enodes(enodes==0) = []; % Remove 3D layers
0557          mdl2.electrode(i).nodes = enodes(:)';
0558       end
0559    end
0560 
0561    ignore = {'electrode', 'nodes', 'boundary', 'elems', 'gnd_node', 'boundary_numbers', 'mat_idx'};
0562    for n=fieldnames(mdl3)'
0563       if ~any(strcmp(n,ignore))
0564          mdl2.(n{:}) = mdl3.(n{:});
0565       end
0566    end
0567 
0568 function mdl = convert_cem2pem(mdl, xyzc)
0569    if ~isfield(mdl, 'electrode')
0570       return;
0571    end
0572    nd = size(mdl.nodes,2); % number of dimensions
0573    for i=1:length(mdl.electrode)
0574       en = mdl.electrode(i).nodes;
0575       nn = mdl.nodes(en,:); % all nodes for this electrode
0576       if nd == 2
0577          np = xyzc(i,[1 3]); % true elec location (2d)
0578       else
0579          np = xyzc(i,:); % true elec location (3d)
0580       end
0581       D = pdist([np; nn]);
0582       [~,idx] = min(D(1,2:end)); % closest CEM node to ideal PEM location
0583       mdl.electrode(i).nodes = en(idx);
0584       % NOTE: used to bump node's location but now we place a corner of the
0585       % square electrode at precisely the correct location
0586       %mdl.nodes(en,:) = np; % shift PEM node to true location
0587    end
0588 
0589 function D2 = pdist(X) % row vectors
0590    if nargin == 2; error('only supports Euclidean distances'); end
0591    %D2 = bsxfun(@plus, dot(X, X, 1)', dot(Y, Y, 1)) - 2*(X'*Y) % 1d
0592    D2 = bsxfun(@minus, X, permute(X,[3 2 1]));
0593    D2 = squeeze(sqrt(sum(D2.^2,2)));
0594 
0595 function [mdl, c] = shift_electrode_positions(mdl, dx)
0596    for i = 1:length(mdl.electrode)
0597       en = mdl.electrode(i).nodes;
0598       ex = mdl.nodes(en,:);
0599       ep(i,:) = (max(ex,[],1) + min(ex,[],1))/2; % mid-point
0600    end
0601    [mdl, c] = correct_electrode_positions(mdl, ep + dx);
0602 
0603 function [mdl, c] = correct_electrode_positions(mdl, xyzc)
0604    %eidors_msg('correct_electrode_positions');
0605    nd = size(mdl.nodes,2);
0606    c = 0; err = 1;
0607    while max(err) > eps
0608       for n = 1:nd
0609          switch(n)
0610             case 1
0611                mdl = mash_nodes(mdl, 'shift_all',     1, 1, xyzc); % X (downslope)
0612             case 2 % 2d: Y, 3d: Z
0613                mdl = mash_nodes(mdl, 'shift_surface', 1, nd, xyzc); % Y or Z (vertical)
0614             case 3 % note: we only do this for 3d
0615                mdl = mash_nodes(mdl, 'shift_middle',  1, 2, xyzc); % Y (cross-slope)
0616             otherwise
0617                error('duh!');
0618          end
0619       end
0620       err = sqrt(sum(mdl_elec_err(mdl, xyzc).^2,2));
0621       c=c+1;
0622       if c >= 100
0623          break;
0624       end
0625    end
0626 
0627 function   mdl = mash_nodes(mdl, method, idm, dim, elec_true)
0628    elec_err = mdl_elec_err(mdl, elec_true);
0629    err = elec_err(:,dim);
0630 
0631    % add borders for electrode positions at
0632    % the volume boundary and 50% of the edge to electrode distance
0633    xq = mdl.nodes(:,idm);
0634    x = [min(xq); ...
0635         mean([min(xq) min(elec_true(:,idm))]); ...
0636         elec_true(:,idm);
0637         mean([max(xq) max(elec_true(:,idm))]); ...
0638         max(xq)];
0639    v = [0; 0; err; 0; 0];
0640 
0641    % scale error to match the electrode locations
0642    interp_method = 'linear';
0643    if strcmp(method, 'shift_surface')
0644          interp_method = 'pchip';
0645    end
0646    vq = interp1(x, v, xq, interp_method, 'extrap');
0647 
0648    switch method
0649       case 'shift_all'
0650          vqs = 1;
0651       case 'shift_middle'
0652          yq = mdl.nodes(:,dim);
0653          yqr = max(yq)-min(yq); % range
0654          yqm = (max(yq) + min(yq))/2; % middle = (max + min)/2
0655          vqs = 1-abs(yq-yqm)./(yqr/2); % scale the shift depending on x's distance from midline
0656       case 'shift_surface' % assumes positive surface
0657          yq = mdl.nodes(:,dim);
0658          yqr = max(yq) - min(yq);
0659          yqm = min(yq); % min
0660          vqs = abs(yq - yqm)./yqr; % scale by distance from surface
0661       otherwise
0662          error(['unrecognized method: ',method]);
0663    end
0664    mdl.nodes(:,dim) = mdl.nodes(:,dim) + (vq .* vqs);
0665 
0666 % calculate the error in electrode position for the fwd_model
0667 function err = mdl_elec_err(mdl, xyzc)
0668    if ~isfield(mdl, 'electrode')
0669       error('electrodes not available on this model, must supply positional error');
0670    end
0671 
0672    nel=length(mdl.electrode); % number of electrodes
0673    nd=size(mdl.nodes,2); % number of dimensions
0674 
0675    eu = ones(nel,nd)*NaN; % init
0676    for i=1:length(mdl.electrode)
0677       nn = mdl.nodes(mdl.electrode(i).nodes,:); % nodes per electrode
0678       eu(i,:) = (max(nn,[],1) + min(nn,[],1))./2; % approx centre of each electrode
0679    end
0680    err = xyzc(:,1:nd) - eu;
0681 
0682 function test_fwd_rec_match(imdl,idx,str)
0683    fwd = imdl.fwd_model.nodes(:,idx);
0684    mxf = max(fwd); mnf = min(fwd);
0685    rec = imdl.rec_model.nodes;
0686    mxr = max(rec); mnr = min(fwd);
0687    unit_test_cmp(['fwd-rec match: ',str],[mxf,mnf],[mxr,mnr],10*eps);
0688 
0689 function do_unit_test
0690    ne = 16;
0691    imdl = mk_geophysics_model('h2p5a', ne);
0692    imdl.fwd_model.stimulation = stim_pattern_geophys(ne, 'Wenner');
0693    img = mk_image(imdl.fwd_model, 1);
0694    img.fwd_solve.get_all_meas = 1;
0695    vh = fwd_solve_halfspace(img);
0696    vd = fwd_solve(img);
0697 clf; h=plot([vh.meas vd.meas],'o--'); legend('analytic','FEM'); set(gca,'box','off'); set(h,'LineWidth',2);
0698 
0699    imdl1 = mk_geophysics_model('h2a',[1:6]);
0700    imdl2 = mk_geophysics_model('h2a',[1:6]');
0701    imdl3 = mk_geophysics_model('h2a',[1:6]'*[1 0]);
0702    imdl3Hnm2d = mk_geophysics_model('H2a',[1:6],{'threshold',Inf}); % try without mashing nodes, no veritcal geometry... electrodes should be precisely located if the electrodes were correctly placed
0703    imdl3Hnm3d = mk_geophysics_model('H3a',[1:6],{'threshold',Inf}); % try without mashing nodes, no veritcal geometry... electrodes should be precisely located if the electrodes were correctly placed
0704    imdl4 = mk_geophysics_model('h2a',[1:6]'*[1 0] + ([1:6]*0+2)'*[0 1]);
0705    R = @(x) [cosd(x) -sind(x); sind(x) cosd(x)]; % rotation matrix
0706    X = [0 2];
0707    imdl5 = mk_geophysics_model('h2a',([1:6]'*[1 0] + ([1:6]*0+1)'*X)*R(-135));
0708    elec_pos_2d = [1 1; 2 2; 3 1; 4 1.5];
0709    elec_pos_3d = [1 0 0; 2 0.5 1; 3 -0.5 2.5; 10 0 3];
0710    imdl2dc = mk_geophysics_model('h2a', elec_pos_2d); % 2D complete electrode model
0711    imdl3dc = mk_geophysics_model('h3a', elec_pos_3d); % 3D complete electrode model
0712    imdl2dp = mk_geophysics_model('H2a', elec_pos_2d); % 2D point electrode model
0713    imdl3dp = mk_geophysics_model('H3a', elec_pos_3d); % 3D point electrode model
0714    imdl2df = mk_geophysics_model('H2a', elec_pos_2d, {'threshold', Inf}); % 2D point electrode model... without mashing
0715    imdl3df = mk_geophysics_model('H3a', elec_pos_3d, {'threshold', Inf}); % 3D point electrode model... without mashing
0716 
0717    % dual meshes
0718    imdlh32 = mk_geophysics_model('h32a', elec_pos_3d);
0719    imdlh22 = mk_geophysics_model('h22a', elec_pos_2d);
0720    imdlH32 = mk_geophysics_model('H32a', elec_pos_3d);
0721    imdlH22 = mk_geophysics_model('H22a', elec_pos_2d);
0722 
0723    % std dual meshes w/ 16 elec
0724    imdlh32_16 = mk_geophysics_model('h32a', 16);
0725    test_fwd_rec_match(imdlh32_16,[1,3],'h32a')
0726    imdlh22_16 = mk_geophysics_model('h22a', 16);
0727    test_fwd_rec_match(imdlh22_16,[1,2],'h22a')
0728    imdlH32_16 = mk_geophysics_model('H32a', 16);
0729    test_fwd_rec_match(imdlH32_16,[1,3],'h22a')
0730    imdlH22_16 = mk_geophysics_model('H22a', 16);
0731    test_fwd_rec_match(imdlH22_16,[1,2],'h22a')
0732 
0733 if 0
0734    % Nolwenn's grid
0735    [yy,xx] = meshgrid(0:3:24, [0 4:2:20*2 20*2+4]); % 9 x 21 grid
0736    xyz = [xx(:) yy(:) zeros(length(yy(:)),1)];
0737 
0738    % test grid
0739    % 1. bad electrodes ... ? fixme?
0740    % 2. mash_nodes breaks for co-located nodes in any particular dimension... need 2d interp?
0741    [yy,xx]=meshgrid(1:3,1:3); xyz = [xx(:) yy(:) zeros(length(yy(:)),1)]; % 3 x 3 grid
0742    imdl = mk_geophysics_model('H3a',xyz);
0743    % TODO add 'g3a' and 'G3a' types?
0744 end
0745 
0746    unit_test_cmp('h2a halfspace vs default TEST', norm(vh.meas - vd.meas), 0, 4e-3);
0747    unit_test_cmp('1d elec list equivalence (row/col)',unit_test_elec_pos(imdl1), unit_test_elec_pos(imdl2));
0748    unit_test_cmp('1d vs. 2d elec list equivalence',unit_test_elec_pos(imdl1), unit_test_elec_pos(imdl3));
0749    unit_test_cmp('2D PEM w/o node mashing, no vertical relief',[1:6]'*[1 0], unit_test_elec_pos(imdl3Hnm2d), eps);
0750    unit_test_cmp('2D PEM *is* PEM',length(imdl3Hnm2d.fwd_model.electrode(1).nodes),1)
0751    unit_test_cmp('3D PEM w/o node mashing, no vertical relief',[1:6]'*[1 0 0], unit_test_elec_pos(imdl3Hnm3d), eps);
0752    unit_test_cmp('3D PEM *is* PEM',length(imdl3Hnm3d.fwd_model.electrode(1).nodes),1)
0753    unit_test_cmp('1d vs. 2d + y=2 elec list equivalence',unit_test_elec_pos(imdl1), unit_test_elec_pos(imdl4)-([1:6]*0+2)'*[0 1]);
0754    unit_test_cmp('1d vs. 2d + y=2 - 135 deg elec eq', ...
0755                  unit_test_elec_pos(imdl1), ...
0756                  unit_test_elec_pos(imdl5, R(135), -X), eps*10);
0757    unit_test_cmp('2d with vertical geometry (mash nodes) CEM', ...
0758                  elec_pos_2d, ...
0759                  unit_test_elec_pos(imdl2dc), 0.01);
0760    unit_test_cmp('3d with vertical geometry (mash nodes) CEM', ...
0761                  elec_pos_3d, ...
0762                  unit_test_elec_pos(imdl3dc), 0.001);
0763    unit_test_cmp('2d with vertical geometry (mash nodes) PEM', ...
0764                  elec_pos_2d, ...
0765                  unit_test_elec_pos(imdl2dp), 10*eps);
0766    unit_test_cmp('3d with vertical geometry (mash nodes) PEM', ...
0767                  elec_pos_3d, ...
0768                  unit_test_elec_pos(imdl3dp), 10*eps);
0769    unit_test_cmp('2d with vertical geometry (w/o mash nodes) PEM', ...
0770                  elec_pos_2d, ...
0771                  unit_test_elec_pos(imdl2df), -10*eps);
0772    unit_test_cmp('3d with vertical geometry (w/o mash nodes) PEM', ...
0773                  elec_pos_3d, ...
0774                  unit_test_elec_pos(imdl3df), -10*eps);
0775 
0776    unit_test_cmp('h32a dual 3d/2d', ...
0777                  elec_pos_3d, ...
0778                  unit_test_elec_pos(imdlh32), 0.001);
0779    unit_test_cmp('H32a dual 3d/2d', ...
0780                  elec_pos_3d, ...
0781                  unit_test_elec_pos(imdlH32), 10*eps);
0782    unit_test_cmp('h22a dual 2d/2d', ...
0783                  elec_pos_2d, ...
0784                  unit_test_elec_pos(imdlh22), 0.002);
0785    unit_test_cmp('H22a dual 2d/2d', ...
0786                  elec_pos_2d, ...
0787                  unit_test_elec_pos(imdlH22), 10*eps);
0788 
0789 clf; subplot(221); show_fem(imdl1.fwd_model); title('models match? A');
0790      subplot(222); show_fem(imdl5.fwd_model); title('models match? C');
0791      subplot(223); show_fem(imdl2dc.fwd_model); title('2d deformations');
0792      subplot(224); show_fem(imdl3dc.fwd_model); title('3d deformations'); view([0 -1 0.01]);
0793 
0794 if 1
0795    clf; subplot(221); show_fem(imdlh22.fwd_model);
0796         subplot(222); show_fem(imdlh32.fwd_model);
0797         subplot(223); show_fem(imdlH22.fwd_model);
0798         subplot(224); show_fem(imdlH32.fwd_model);
0799 end
0800 
0801 if 0
0802    clf; subplot(221); show_fem(imdlh22.rec_model);
0803         subplot(222); show_fem(imdlh32.rec_model);
0804         subplot(223); show_fem(imdlH22.rec_model);
0805         subplot(224); show_fem(imdlH32.rec_model);
0806 end
0807 
0808    imdl = mk_geophysics_model('h2p5a', ne, {'extra_ng_code', 'solid tt = orthobrick(-1,-1,-1;-0.5,0,-0.5);\ntlo tt;\n'});
0809    clf; show_fem(imdl.fwd_model);
0810    % Note: this just has to be swallowed safely...
0811 
0812 function xyz = unit_test_elec_pos(imdl, R, X)
0813    if nargin < 2; R = 1; end
0814    if nargin < 3; X = 0; end
0815    fmdl = imdl.fwd_model;
0816    xyz = zeros(length(fmdl.electrode),size(fmdl.nodes,2))*NaN;
0817    for i = 1:length(fmdl.electrode)
0818       nn = fmdl.nodes(fmdl.electrode(i).nodes,:);
0819       xyz(i,:) = (max(nn,[],1) + min(nn,[],1))/2;
0820       xyz(i,:) = xyz(i,:)*R + X;
0821    end