mk_geophysics_model

PURPOSE ^

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

SYNOPSIS ^

function imdl = mk_geophysics_model(str, ne, opt);

DESCRIPTION ^

 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
       ne = [0 1; 2 1.5; 3 1.2; 7 2.5] % 2d
       ne = [0 0.1 1; 2 -0.1 1.5; 3 -0.15 1.2; 7 0 2.5] % 3d
 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_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]
       '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)

 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, 2016 A. Boyle
 License: GPL version 2 or version 3

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005