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 6926 2024-05-31 15:34:13Z 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 % standard field order
0414 imdl = eidors_obj('set', imdl);
0415 
0416 if save_model_to_disk
0417    save(filename,'imdl','SALT');
0418 end
0419 
0420 % convert 2x3 array to "x1,y1,z1;x2,y2,z2" string
0421 % convert 1x3 array to "x1,y1,z1" string
0422 % if 1x1 array, then use xy_ctr=x1,y1 and [0 0 1]=x2,y2,z2 (z+)
0423 function s = a2s(a)
0424 if length(a(:)) == 3
0425    s = sprintf('%f,%f,%f', ...
0426                a(1), a(2), a(3));
0427 else
0428    s = sprintf('%f,%f,%f;%f,%f,%f', ...
0429                 a(1,1), a(1,2), a(1,3), ...
0430                a(2,1), a(2,2), a(2,3));
0431 end
0432 
0433 % returns R rotation/scaling and X0 offset
0434 % xyz1 = R * xyz + X; % rotate and scale to +/- 1
0435 function [R,X,var] = rot_line_to_xaxis(xyz)
0436 x = xyz(:,1); y=xyz(:,2); z=xyz(:,3);
0437 
0438 % fit line to points
0439 % http://www.mathworks.com/matlabcentral/newsreader/view_thread/32502
0440 p = mean(xyz);
0441 [U,S,V] = svd([x-p(1), y-p(2), z-p(3)]);
0442 if V(end,1) ~= 0
0443    N=1/V(end,1)*V(:,1);
0444 else
0445    N=V(:,1);
0446 end
0447 A=p' + dot( xyz(1,  :) - p, N ) * N/norm(N)^2;
0448 B=p' + dot( xyz(end,:) - p, N ) * N/norm(N)^2;
0449 
0450 % rotate line to +y-axis
0451 a = N/norm(N);
0452 b = [ 1 0 0 ]'; % +x
0453 % http://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d
0454 v = cross(a, b);
0455 s = norm(v); c = dot(a, b); % sin, cos
0456 xt = @(v) [   0  -v(3)  v(2); ... % skew symmetric cross-product of v
0457             v(3)    0  -v(1); ... % diagnoal is the scaling identity matrix
0458            -v(2)  v(1)    0];
0459 if abs(s) < eps*1e3
0460    R = eye(3);
0461 else
0462    R = (eye(3) + xt(v) + xt(v)^2 * (1-c)/s^2);
0463    R=R'; % for right multiply: xyz*R
0464 end
0465 
0466 X = repmat(p, size(xyz,1), 1);
0467 R = R / max(range((xyz-X)*R)) * 2;
0468 
0469 DEBUG=0;
0470 if DEBUG
0471    clf;
0472    subplot(211);
0473    plot3(x,y,z,'bo');
0474    hold on;
0475    plot3(x(1),y(1),z(1),'go');
0476    plot3(p(1),p(2),p(3),'ro');
0477    plot3([A(1),B(1)],[A(2),B(2)],[A(3),B(3)]);
0478    grid; axis equal; hold off;
0479    xlabel('x'); ylabel('y'); zlabel('z');
0480    title('pre-rotation');
0481    subplot(212);
0482    xx = (xyz-X)*R;
0483    plot3(xx(:,1), xx(:,2), xx(:,3),'bo');
0484    hold on;
0485    plot3(xx(1,1), xx(1,2), xx(1,3),'go');
0486    plot3(0,0,0,'ro');
0487    grid; axis equal; hold off;
0488    xlabel('x'); ylabel('y'); zlabel('z');
0489    title('post-rotation');
0490 end
0491 
0492 % the (max-min) range of a variable's values
0493 function r = range(a)
0494 r = max(a(:))-min(a(:));
0495 
0496 function [mdl2,idx2] = copy_mdl2d_from3d(mdl3,idx3,sel);
0497 % AB: taken from EIDORS function ng_mk_gen_models() subfunction of the same name
0498 % AB: NEW: sel = 'x', 'y' or 'z' -- default was Z, we want X
0499    if sel == 'x'
0500       % swap Z and X
0501       T = [ 0 0 1; 0 1 0; 1 0 0 ];
0502    elseif sel == 'y'
0503       % swap Z and Y
0504       T = [ 1 0 0; 0 0 1; 0 1 0 ];
0505    elseif sel == 'z'
0506       T = eye(3);
0507    else
0508       error('sel must be "x", "y" or "z"');
0509    end
0510    mdl3.nodes = mdl3.nodes * T; % AB: SWAP axes
0511 
0512    % set name
0513    mdl2 = eidors_obj('fwd_model',sprintf('%s 2D',mdl3.name));
0514 
0515    % set nodes
0516    [bdy,idx] = find_boundary(mdl3.elems);
0517    vtx = mdl3.nodes;
0518    z_vtx = reshape(vtx(bdy,3), size(bdy) );
0519    z_vtx_thres = max(z_vtx(:))-10*eps*range(z_vtx(:));
0520    lay0  = find( all(z_vtx >= z_vtx_thres, 2) );
0521    bdy0  = bdy( lay0, :);
0522 
0523    vtx0  = unique(bdy0(:));
0524    mdl2.nodes = vtx(vtx0,1:2);
0525 
0526    % set elems
0527    nmap  = zeros(size(vtx,1),1); nmap(vtx0) = 1:length(vtx0);
0528    bdy0  = reshape(nmap(bdy0), size(bdy0) ); % renumber to new scheme
0529    mdl2.elems = bdy0;
0530 
0531    % set boundary
0532    mdl2.boundary = find_boundary( mdl2.elems);
0533 
0534    % set gnd_node
0535    mdl2.gnd_node = nmap(mdl3.gnd_node);
0536 
0537    % set material indices
0538    % TODO: vectorize code
0539    idx2 = {};
0540    idx0  = idx( lay0, :);
0541    for i=1:size(idx3,2)
0542       idx2{i} = [];
0543       ii = 1;
0544       for j=1:size(idx3{i},1)
0545          idx_tmp = find( idx0==idx3{i}(j) );
0546          if not(isempty(idx_tmp))
0547             idx2{i}(ii,1) = idx_tmp(1,1);
0548             ii = ii + 1;
0549          end
0550       end
0551    end
0552 
0553    % set electrode
0554    if isfield(mdl3,'electrode')
0555       mdl2.electrode = mdl3.electrode;
0556       for i=1:length(mdl2.electrode);
0557          nn = mdl3.nodes(mdl3.electrode(i).nodes,:);
0558          enodes = nmap( mdl2.electrode(i).nodes );
0559          enodes(enodes==0) = []; % Remove 3D layers
0560          mdl2.electrode(i).nodes = enodes(:)';
0561       end
0562    end
0563 
0564    ignore = {'electrode', 'nodes', 'boundary', 'elems', 'gnd_node', 'boundary_numbers', 'mat_idx'};
0565    for n=fieldnames(mdl3)'
0566       if ~any(strcmp(n,ignore))
0567          mdl2.(n{:}) = mdl3.(n{:});
0568       end
0569    end
0570 
0571 function mdl = convert_cem2pem(mdl, xyzc)
0572    if ~isfield(mdl, 'electrode')
0573       return;
0574    end
0575    nd = size(mdl.nodes,2); % number of dimensions
0576    for i=1:length(mdl.electrode)
0577       en = mdl.electrode(i).nodes;
0578       nn = mdl.nodes(en,:); % all nodes for this electrode
0579       if nd == 2
0580          np = xyzc(i,[1 3]); % true elec location (2d)
0581       else
0582          np = xyzc(i,:); % true elec location (3d)
0583       end
0584       D = pdist([np; nn]);
0585       [~,idx] = min(D(1,2:end)); % closest CEM node to ideal PEM location
0586       mdl.electrode(i).nodes = en(idx);
0587       % NOTE: used to bump node's location but now we place a corner of the
0588       % square electrode at precisely the correct location
0589       %mdl.nodes(en,:) = np; % shift PEM node to true location
0590    end
0591 
0592 function D2 = pdist(X) % row vectors
0593    if nargin == 2; error('only supports Euclidean distances'); end
0594    %D2 = bsxfun(@plus, dot(X, X, 1)', dot(Y, Y, 1)) - 2*(X'*Y) % 1d
0595    D2 = bsxfun(@minus, X, permute(X,[3 2 1]));
0596    D2 = squeeze(sqrt(sum(D2.^2,2)));
0597 
0598 function [mdl, c] = shift_electrode_positions(mdl, dx)
0599    for i = 1:length(mdl.electrode)
0600       en = mdl.electrode(i).nodes;
0601       ex = mdl.nodes(en,:);
0602       ep(i,:) = (max(ex,[],1) + min(ex,[],1))/2; % mid-point
0603    end
0604    [mdl, c] = correct_electrode_positions(mdl, ep + dx);
0605 
0606 function [mdl, c] = correct_electrode_positions(mdl, xyzc)
0607    %eidors_msg('correct_electrode_positions');
0608    nd = size(mdl.nodes,2);
0609    c = 0; err = 1;
0610    while max(err) > eps
0611       for n = 1:nd
0612          switch(n)
0613             case 1
0614                mdl = mash_nodes(mdl, 'shift_all',     1, 1, xyzc); % X (downslope)
0615             case 2 % 2d: Y, 3d: Z
0616                mdl = mash_nodes(mdl, 'shift_surface', 1, nd, xyzc); % Y or Z (vertical)
0617             case 3 % note: we only do this for 3d
0618                mdl = mash_nodes(mdl, 'shift_middle',  1, 2, xyzc); % Y (cross-slope)
0619             otherwise
0620                error('duh!');
0621          end
0622       end
0623       err = sqrt(sum(mdl_elec_err(mdl, xyzc).^2,2));
0624       c=c+1;
0625       if c >= 100
0626          break;
0627       end
0628    end
0629 
0630 function   mdl = mash_nodes(mdl, method, idm, dim, elec_true)
0631    elec_err = mdl_elec_err(mdl, elec_true);
0632    err = elec_err(:,dim);
0633 
0634    % add borders for electrode positions at
0635    % the volume boundary and 50% of the edge to electrode distance
0636    xq = mdl.nodes(:,idm);
0637    x = [min(xq); ...
0638         mean([min(xq) min(elec_true(:,idm))]); ...
0639         elec_true(:,idm);
0640         mean([max(xq) max(elec_true(:,idm))]); ...
0641         max(xq)];
0642    v = [0; 0; err; 0; 0];
0643 
0644    % scale error to match the electrode locations
0645    interp_method = 'linear';
0646    if strcmp(method, 'shift_surface')
0647          interp_method = 'pchip';
0648    end
0649    vq = interp1(x, v, xq, interp_method, 'extrap');
0650 
0651    switch method
0652       case 'shift_all'
0653          vqs = 1;
0654       case 'shift_middle'
0655          yq = mdl.nodes(:,dim);
0656          yqr = max(yq)-min(yq); % range
0657          yqm = (max(yq) + min(yq))/2; % middle = (max + min)/2
0658          vqs = 1-abs(yq-yqm)./(yqr/2); % scale the shift depending on x's distance from midline
0659       case 'shift_surface' % assumes positive surface
0660          yq = mdl.nodes(:,dim);
0661          yqr = max(yq) - min(yq);
0662          yqm = min(yq); % min
0663          vqs = abs(yq - yqm)./yqr; % scale by distance from surface
0664       otherwise
0665          error(['unrecognized method: ',method]);
0666    end
0667    mdl.nodes(:,dim) = mdl.nodes(:,dim) + (vq .* vqs);
0668 
0669 % calculate the error in electrode position for the fwd_model
0670 function err = mdl_elec_err(mdl, xyzc)
0671    if ~isfield(mdl, 'electrode')
0672       error('electrodes not available on this model, must supply positional error');
0673    end
0674 
0675    nel=length(mdl.electrode); % number of electrodes
0676    nd=size(mdl.nodes,2); % number of dimensions
0677 
0678    eu = ones(nel,nd)*NaN; % init
0679    for i=1:length(mdl.electrode)
0680       nn = mdl.nodes(mdl.electrode(i).nodes,:); % nodes per electrode
0681       eu(i,:) = (max(nn,[],1) + min(nn,[],1))./2; % approx centre of each electrode
0682    end
0683    err = xyzc(:,1:nd) - eu;
0684 
0685 function test_fwd_rec_match(imdl,idx,str)
0686    fwd = imdl.fwd_model.nodes(:,idx);
0687    mxf = max(fwd); mnf = min(fwd);
0688    rec = imdl.rec_model.nodes;
0689    mxr = max(rec); mnr = min(fwd);
0690    unit_test_cmp(['fwd-rec match: ',str],[mxf,mnf],[mxr,mnr],10*eps);
0691 
0692 function do_unit_test
0693    ne = 16;
0694    imdl = mk_geophysics_model('h2p5a', ne);
0695    imdl.fwd_model.stimulation = stim_pattern_geophys(ne, 'Wenner');
0696    img = mk_image(imdl.fwd_model, 1);
0697    img.fwd_solve.get_all_meas = 1;
0698    vh = fwd_solve_halfspace(img);
0699    vd = fwd_solve(img);
0700 clf; h=plot([vh.meas vd.meas],'o--'); legend('analytic','FEM'); set(gca,'box','off'); set(h,'LineWidth',2);
0701 
0702    imdl1 = mk_geophysics_model('h2a',[1:6]);
0703    imdl2 = mk_geophysics_model('h2a',[1:6]');
0704    imdl3 = mk_geophysics_model('h2a',[1:6]'*[1 0]);
0705    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
0706    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
0707    imdl4 = mk_geophysics_model('h2a',[1:6]'*[1 0] + ([1:6]*0+2)'*[0 1]);
0708    R = @(x) [cosd(x) -sind(x); sind(x) cosd(x)]; % rotation matrix
0709    X = [0 2];
0710    imdl5 = mk_geophysics_model('h2a',([1:6]'*[1 0] + ([1:6]*0+1)'*X)*R(-135));
0711    elec_pos_2d = [1 1; 2 2; 3 1; 4 1.5];
0712    elec_pos_3d = [1 0 0; 2 0.5 1; 3 -0.5 2.5; 10 0 3];
0713    imdl2dc = mk_geophysics_model('h2a', elec_pos_2d); % 2D complete electrode model
0714    imdl3dc = mk_geophysics_model('h3a', elec_pos_3d); % 3D complete electrode model
0715    imdl2dp = mk_geophysics_model('H2a', elec_pos_2d); % 2D point electrode model
0716    imdl3dp = mk_geophysics_model('H3a', elec_pos_3d); % 3D point electrode model
0717    imdl2df = mk_geophysics_model('H2a', elec_pos_2d, {'threshold', Inf}); % 2D point electrode model... without mashing
0718    imdl3df = mk_geophysics_model('H3a', elec_pos_3d, {'threshold', Inf}); % 3D point electrode model... without mashing
0719 
0720    % dual meshes
0721    imdlh32 = mk_geophysics_model('h32a', elec_pos_3d);
0722    imdlh22 = mk_geophysics_model('h22a', elec_pos_2d);
0723    imdlH32 = mk_geophysics_model('H32a', elec_pos_3d);
0724    imdlH22 = mk_geophysics_model('H22a', elec_pos_2d);
0725 
0726    % std dual meshes w/ 16 elec
0727    imdlh32_16 = mk_geophysics_model('h32a', 16);
0728    test_fwd_rec_match(imdlh32_16,[1,3],'h32a')
0729    imdlh22_16 = mk_geophysics_model('h22a', 16);
0730    test_fwd_rec_match(imdlh22_16,[1,2],'h22a')
0731    imdlH32_16 = mk_geophysics_model('H32a', 16);
0732    test_fwd_rec_match(imdlH32_16,[1,3],'h22a')
0733    imdlH22_16 = mk_geophysics_model('H22a', 16);
0734    test_fwd_rec_match(imdlH22_16,[1,2],'h22a')
0735 
0736 if 0
0737    % Nolwenn's grid
0738    [yy,xx] = meshgrid(0:3:24, [0 4:2:20*2 20*2+4]); % 9 x 21 grid
0739    xyz = [xx(:) yy(:) zeros(length(yy(:)),1)];
0740 
0741    % test grid
0742    % 1. bad electrodes ... ? fixme?
0743    % 2. mash_nodes breaks for co-located nodes in any particular dimension... need 2d interp?
0744    [yy,xx]=meshgrid(1:3,1:3); xyz = [xx(:) yy(:) zeros(length(yy(:)),1)]; % 3 x 3 grid
0745    imdl = mk_geophysics_model('H3a',xyz);
0746    % TODO add 'g3a' and 'G3a' types?
0747 end
0748 
0749    unit_test_cmp('h2a halfspace vs default TEST', norm(vh.meas - vd.meas), 0, 4e-3);
0750    unit_test_cmp('1d elec list equivalence (row/col)',unit_test_elec_pos(imdl1), unit_test_elec_pos(imdl2));
0751    unit_test_cmp('1d vs. 2d elec list equivalence',unit_test_elec_pos(imdl1), unit_test_elec_pos(imdl3));
0752    unit_test_cmp('2D PEM w/o node mashing, no vertical relief',[1:6]'*[1 0], unit_test_elec_pos(imdl3Hnm2d), eps);
0753    unit_test_cmp('2D PEM *is* PEM',length(imdl3Hnm2d.fwd_model.electrode(1).nodes),1)
0754    unit_test_cmp('3D PEM w/o node mashing, no vertical relief',[1:6]'*[1 0 0], unit_test_elec_pos(imdl3Hnm3d), eps);
0755    unit_test_cmp('3D PEM *is* PEM',length(imdl3Hnm3d.fwd_model.electrode(1).nodes),1)
0756    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]);
0757    unit_test_cmp('1d vs. 2d + y=2 - 135 deg elec eq', ...
0758                  unit_test_elec_pos(imdl1), ...
0759                  unit_test_elec_pos(imdl5, R(135), -X), eps*10);
0760    unit_test_cmp('2d with vertical geometry (mash nodes) CEM', ...
0761                  elec_pos_2d, ...
0762                  unit_test_elec_pos(imdl2dc), 0.01);
0763    unit_test_cmp('3d with vertical geometry (mash nodes) CEM', ...
0764                  elec_pos_3d, ...
0765                  unit_test_elec_pos(imdl3dc), 0.001);
0766    unit_test_cmp('2d with vertical geometry (mash nodes) PEM', ...
0767                  elec_pos_2d, ...
0768                  unit_test_elec_pos(imdl2dp), 10*eps);
0769    unit_test_cmp('3d with vertical geometry (mash nodes) PEM', ...
0770                  elec_pos_3d, ...
0771                  unit_test_elec_pos(imdl3dp), 10*eps);
0772    unit_test_cmp('2d with vertical geometry (w/o mash nodes) PEM', ...
0773                  elec_pos_2d, ...
0774                  unit_test_elec_pos(imdl2df), -10*eps);
0775    unit_test_cmp('3d with vertical geometry (w/o mash nodes) PEM', ...
0776                  elec_pos_3d, ...
0777                  unit_test_elec_pos(imdl3df), -10*eps);
0778 
0779    unit_test_cmp('h32a dual 3d/2d', ...
0780                  elec_pos_3d, ...
0781                  unit_test_elec_pos(imdlh32), 0.001);
0782    unit_test_cmp('H32a dual 3d/2d', ...
0783                  elec_pos_3d, ...
0784                  unit_test_elec_pos(imdlH32), 10*eps);
0785    unit_test_cmp('h22a dual 2d/2d', ...
0786                  elec_pos_2d, ...
0787                  unit_test_elec_pos(imdlh22), 0.002);
0788    unit_test_cmp('H22a dual 2d/2d', ...
0789                  elec_pos_2d, ...
0790                  unit_test_elec_pos(imdlH22), 10*eps);
0791 
0792 clf; subplot(221); show_fem(imdl1.fwd_model); title('models match? A');
0793      subplot(222); show_fem(imdl5.fwd_model); title('models match? C');
0794      subplot(223); show_fem(imdl2dc.fwd_model); title('2d deformations');
0795      subplot(224); show_fem(imdl3dc.fwd_model); title('3d deformations'); view([0 -1 0.01]);
0796 
0797 if 1
0798    clf; subplot(221); show_fem(imdlh22.fwd_model);
0799         subplot(222); show_fem(imdlh32.fwd_model);
0800         subplot(223); show_fem(imdlH22.fwd_model);
0801         subplot(224); show_fem(imdlH32.fwd_model);
0802 end
0803 
0804 if 0
0805    clf; subplot(221); show_fem(imdlh22.rec_model);
0806         subplot(222); show_fem(imdlh32.rec_model);
0807         subplot(223); show_fem(imdlH22.rec_model);
0808         subplot(224); show_fem(imdlH32.rec_model);
0809 end
0810 
0811    imdl = mk_geophysics_model('h2p5a', ne, {'extra_ng_code', 'solid tt = orthobrick(-1,-1,-1;-0.5,0,-0.5);\ntlo tt;\n'});
0812    clf; show_fem(imdl.fwd_model);
0813    % Note: this just has to be swallowed safely...
0814 
0815 function xyz = unit_test_elec_pos(imdl, R, X)
0816    if nargin < 2; R = 1; end
0817    if nargin < 3; X = 0; end
0818    fmdl = imdl.fwd_model;
0819    xyz = zeros(length(fmdl.electrode),size(fmdl.nodes,2))*NaN;
0820    for i = 1:length(fmdl.electrode)
0821       nn = fmdl.nodes(fmdl.electrode(i).nodes,:);
0822       xyz(i,:) = (max(nn,[],1) + min(nn,[],1))/2;
0823       xyz(i,:) = xyz(i,:)*R + X;
0824    end