ng_mk_ellip_models

 NG_MAKE_ELLIP_MODELS: create elliptical models using netgen
[fmdl,mat_idx] = ng_mk_ellip_models(ellip_shape, elec_pos, ...
                 elec_shape, extra_ng_code);
 INPUT:
 ellip_shape = {height, [x_radius, y_radius, [maxsz]]}
    if height = 0 -> calculate a 2D shape
    x_radius, y_radius (OPT)  -> elliptical eccentricity in x,y directions(default = 1)
    maxsz  (OPT)  -> max size of mesh elems (default = course mesh)

 ELECTRODE POSITIONS:
  elec_pos = [n_elecs_per_plane,z_planes] 
     OR
  elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)

 ELECTRODE SHAPES::
  elec_shape = [width,height, maxsz]  % Rectangular elecs
     OR
  elec_shape = [radius, 0, maxsz ]    % Circular elecs
     OR 
  elec_shape = [0, 0, maxsz ]         % Point electrodes
    (point elecs does some tricks with netgen, so the elecs aren't exactly where you ask)

 Specify either a common electrode shape or for each electrode

 EXTRA_NG_CODE
   string of extra code to put into netgen geo file. Normally this
   would be to insert extra materials into the space

 OUTPUT:
  fmdl    - fwd_model object
  mat_idx - indices of materials (if extra_ng_code is used)
    Note mat_idx does not work in 2D. Netgen does not provide it.


 USAGE EXAMPLES:
 Simple 3D ellipse. Major, minor axes = [1.5, 0.8]. No electrodes
     fmdl= ng_mk_ellip_models([1,1.5,0.8],[0],[]);  show_fem(fmdl);
 
 Simple 2D cylinder. Axes = [1.5,0.8]. Refined to 0.1
     fmdl= ng_mk_ellip_models([0,1.5,0.8,0.1],[],[]); show_fem(fmdl);
 
 3D cylinder. Axes = [1.5,0.8]. 2 planes of 8 elecs with radius 0.1
     fmdl= ng_mk_ellip_models([1,1.2,0.8],[8,0.3,0.7],[0.1]); show_fem(fmdl);
 
 3D cylinder. Axes= [1.3,1] = 1. 7 rect elecs with no refinement
     fmdl= ng_mk_ellip_models([3,1.3],[7,1],[0.2,0.3]); show_fem(fmdl);
 
 2D cylinder. Axes = [1.2,0.8]. 11 rect elecs with refinement. Rotated 45 degrees
     fmdl= ng_mk_ellip_models([0,1.2,0.8],[11],[0.2,0,0.05]); 
     th = 45* [2*pi/360];
     fmdl.nodes = fmdl.nodes*[cos(th),sin(th);-sin(th),cos(th)]; show_fem(fmdl);
 
 2D cylinder. elecs at 0, 90 and 120 degrees
     fmdl= ng_mk_ellip_models([0,1.2,0.8],[0;90;120],[0.2,0,0.03]); show_fem(fmdl);
 
 3D cylinder. Various elecs at 0, 30, 60, 90 in planes
     el_pos = [0,0.5;30,1;60,1.5;90,2.0];
     el_sz  = [0.2,0,0.1;0.1,0,0.05;0.2,0.2,0.02;0.2,0.4,0.5];
     fmdl= ng_mk_ellip_models([3,0.8,1.2],el_pos,el_sz); show_fem(fmdl);
 
 Simple 3D cylinder with a ball
     extra={'ball','solid ball = sphere(0.5,0.5,1;0.4);'};
     [fmdl,mat_idx]= ng_mk_ellip_models([2,1.2,0.8],[8,1],[0.1],extra); 
     img= mk_image(fmdl, 1);
     img.elem_data(mat_idx{2}) = 2;   show_fem(img);

0001 function [fmdl,mat_idx] = ng_mk_ellip_models(ellip_shape, elec_pos, ...
0002                   elec_shape, extra_ng_code);
0003 % NG_MAKE_ELLIP_MODELS: create elliptical models using netgen
0004 %[fmdl,mat_idx] = ng_mk_ellip_models(ellip_shape, elec_pos, ...
0005 %                 elec_shape, extra_ng_code);
0006 % INPUT:
0007 % ellip_shape = {height, [x_radius, y_radius, [maxsz]]}
0008 %    if height = 0 -> calculate a 2D shape
0009 %    x_radius, y_radius (OPT)  -> elliptical eccentricity in x,y directions(default = 1)
0010 %    maxsz  (OPT)  -> max size of mesh elems (default = course mesh)
0011 %
0012 % ELECTRODE POSITIONS:
0013 %  elec_pos = [n_elecs_per_plane,z_planes]
0014 %     OR
0015 %  elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)
0016 %
0017 % ELECTRODE SHAPES::
0018 %  elec_shape = [width,height, maxsz]  % Rectangular elecs
0019 %     OR
0020 %  elec_shape = [radius, 0, maxsz ]    % Circular elecs
0021 %     OR
0022 %  elec_shape = [0, 0, maxsz ]         % Point electrodes
0023 %    (point elecs does some tricks with netgen, so the elecs aren't exactly where you ask)
0024 %
0025 % Specify either a common electrode shape or for each electrode
0026 %
0027 % EXTRA_NG_CODE
0028 %   string of extra code to put into netgen geo file. Normally this
0029 %   would be to insert extra materials into the space
0030 %
0031 % OUTPUT:
0032 %  fmdl    - fwd_model object
0033 %  mat_idx - indices of materials (if extra_ng_code is used)
0034 %    Note mat_idx does not work in 2D. Netgen does not provide it.
0035 %
0036 %
0037 % USAGE EXAMPLES:
0038 % Simple 3D ellipse. Major, minor axes = [1.5, 0.8]. No electrodes
0039 %     fmdl= ng_mk_ellip_models([1,1.5,0.8],[0],[]);  show_fem(fmdl);
0040 %
0041 % Simple 2D cylinder. Axes = [1.5,0.8]. Refined to 0.1
0042 %     fmdl= ng_mk_ellip_models([0,1.5,0.8,0.1],[],[]); show_fem(fmdl);
0043 %
0044 % 3D cylinder. Axes = [1.5,0.8]. 2 planes of 8 elecs with radius 0.1
0045 %     fmdl= ng_mk_ellip_models([1,1.2,0.8],[8,0.3,0.7],[0.1]); show_fem(fmdl);
0046 %
0047 % 3D cylinder. Axes= [1.3,1] = 1. 7 rect elecs with no refinement
0048 %     fmdl= ng_mk_ellip_models([3,1.3],[7,1],[0.2,0.3]); show_fem(fmdl);
0049 %
0050 % 2D cylinder. Axes = [1.2,0.8]. 11 rect elecs with refinement. Rotated 45 degrees
0051 %     fmdl= ng_mk_ellip_models([0,1.2,0.8],[11],[0.2,0,0.05]);
0052 %     th = 45* [2*pi/360];
0053 %     fmdl.nodes = fmdl.nodes*[cos(th),sin(th);-sin(th),cos(th)]; show_fem(fmdl);
0054 %
0055 % 2D cylinder. elecs at 0, 90 and 120 degrees
0056 %     fmdl= ng_mk_ellip_models([0,1.2,0.8],[0;90;120],[0.2,0,0.03]); show_fem(fmdl);
0057 %
0058 % 3D cylinder. Various elecs at 0, 30, 60, 90 in planes
0059 %     el_pos = [0,0.5;30,1;60,1.5;90,2.0];
0060 %     el_sz  = [0.2,0,0.1;0.1,0,0.05;0.2,0.2,0.02;0.2,0.4,0.5];
0061 %     fmdl= ng_mk_ellip_models([3,0.8,1.2],el_pos,el_sz); show_fem(fmdl);
0062 %
0063 % Simple 3D cylinder with a ball
0064 %     extra={'ball','solid ball = sphere(0.5,0.5,1;0.4);'};
0065 %     [fmdl,mat_idx]= ng_mk_ellip_models([2,1.2,0.8],[8,1],[0.1],extra);
0066 %     img= mk_image(fmdl, 1);
0067 %     img.elem_data(mat_idx{2}) = 2;   show_fem(img);
0068 
0069 
0070 
0071 % (C) Andy Adler, 2010. (C) Alistair Boyle 2013. Licenced under GPL v2 or v3
0072 % $Id: ng_mk_ellip_models.m 6970 2024-10-01 11:15:54Z aadler $
0073 
0074 if ischar(ellip_shape) && strcmp(ellip_shape,'UNIT_TEST'), do_unit_test, return, end
0075 if nargin < 4; extra_ng_code = {'',''}; end
0076 
0077 copt.cache_obj = { ellip_shape, elec_pos, elec_shape, extra_ng_code};
0078 copt.fstr = 'ng_mk_ellip_models';
0079 args = {ellip_shape, elec_pos, elec_shape, extra_ng_code};
0080 copt.cache_on_ng_opt = true;
0081 
0082 fmdl = eidors_cache(@mk_ellip_model,args,copt);
0083 
0084 mat_idx = fmdl.mat_idx;
0085 
0086 function fmdl = mk_ellip_model( ellip_shape, elec_pos, elec_shape, extra_ng_code );
0087 
0088    fnstem = tempname;
0089    geofn= [fnstem,'.geo'];
0090    ptsfn= [fnstem,'.msz'];
0091    meshfn= [fnstem,'.vol'];
0092 
0093    [tank_height, tank_radius, tank_maxh, is2D] = parse_shape(ellip_shape);
0094    [elecs, centres] = parse_elecs( elec_pos, elec_shape,  ...
0095                           tank_height, tank_radius, is2D );
0096 
0097    n_pts = write_geo_file(geofn, ptsfn, tank_height, tank_radius, ...
0098                   tank_maxh, elecs, extra_ng_code);
0099    if n_pts == 0 
0100       call_netgen( geofn, meshfn);
0101    else
0102       call_netgen( geofn, meshfn, ptsfn);
0103    end
0104 
0105    fmdl = ng_mk_fwd_model( meshfn, centres, 'ng', []);
0106 
0107    delete(geofn); delete(meshfn); delete(ptsfn); % remove temp files
0108    if is2D
0109       fmdl = mdl2d_from3d(fmdl);
0110    end
0111 
0112    % convert CEM to PEM if so configured
0113    % TODO shunt model is unsupported
0114    if isfield(fmdl,'electrode');
0115      fmdl.electrode = pem_from_cem(elecs, fmdl.electrode, fmdl.nodes);
0116    end
0117 
0118    % standard field order
0119    fmdl = eidors_obj('fwd_model',fmdl);
0120 
0121 % for the newest netgen, we can't call msz file unless there are actually points in  it
0122 function n_pts_elecs = write_geo_file(geofn, ptsfn, tank_height, tank_radius, ...
0123                         tank_maxh, elecs, extra_ng_code);
0124    fid=fopen(geofn,'w');
0125    write_header(fid,tank_height,tank_radius,tank_maxh,extra_ng_code);
0126 
0127    n_elecs = length(elecs);
0128    %  elecs(i).pos   = [x,y,z]
0129    %  elecs(i).shape = 'C' or 'R'
0130    %  elecs(i).dims  = [radius] or [width,height]
0131    %  elecs(i).maxh  = '-maxh=#' or '';
0132    pts_elecs_idx = []; 
0133 
0134    for i=1:n_elecs
0135       name = sprintf('elec%04d',i);
0136       pos = elecs(i).pos;
0137       % calculate the normal vector to the shape
0138       ab = tank_radius(1)/tank_radius(2);
0139       dirn= pos.*[inv(ab), ab, 0 ];
0140       switch elecs(i).shape
0141        case 'C'
0142          write_circ_elec(fid,name, pos, dirn,  ...
0143                elecs(i).dims, tank_radius, elecs(i).maxh);
0144        case 'R'
0145          write_rect_elec(fid,name, pos, dirn,  ...
0146                elecs(i).dims, tank_radius, elecs(i).maxh);
0147        case 'P'
0148          if 0 % Netgen doesn't put elecs where you ask
0149             pts_elecs_idx = [ pts_elecs_idx, i]; 
0150             continue; % DON'T print solid cyl
0151          end
0152          write_rect_elec(fid,name, pos, dirn,  ...
0153                elecs(i).dims, tank_radius, elecs(i).maxh);
0154 
0155        otherwise; error('huh? shouldnt get here');
0156       end
0157       fprintf(fid,'solid cyl%04d = %s and not bigcyl; \n',i,name);
0158    end
0159 
0160    % SHOULD tank_maxh go here?
0161    fprintf(fid,'tlo bigcyl;\n');
0162    for i=1:n_elecs
0163       if any(i == pts_elecs_idx); continue; end
0164       fprintf(fid,'tlo cyl%04d cyl -col=[1,0,0];\n ',i);
0165    end
0166 
0167    for i=1:length(extra_ng_code)-1
0168       if ~isempty(extra_ng_code{i})
0169          fprintf(fid,'tlo %s  -col=[0,1,0];\n',extra_ng_code{i});
0170       end
0171    end
0172 
0173    fclose(fid); % geofn
0174 % From Documentation: Syntax is
0175 % np
0176 % x1 y1 z1 h1
0177 % x2 y2 z2 h2
0178    n_pts_elecs= length(pts_elecs_idx);
0179    fid=fopen(ptsfn,'w');
0180    fprintf(fid,'%d\n',n_pts_elecs);
0181    for i = pts_elecs_idx;
0182       posxy = elecs(i).pos(1:2);
0183       fprintf(fid,'%10f %10f 0 %10f\n', posxy, elecs(i).dims(1) );
0184    end
0185    fclose(fid); % ptsfn
0186 
0187 function [tank_height, tank_radius, tank_maxh, is2D] = ...
0188               parse_shape(cyl_shape);
0189    tank_height = cyl_shape(1);
0190    tank_radius = [1,1];
0191    tank_maxh   = 0;
0192    is2D = 0;
0193    lcs = length(cyl_shape);
0194 
0195    if lcs == 2
0196       tank_radius(1)=cyl_shape(2);
0197    elseif lcs >= 3
0198       tank_radius=cyl_shape(2:3);
0199       if diff(tank_radius) == 0;
0200          warning(['Using ng_mk_ellip_models to create cylinder. This may fail for '...
0201                   'even electrode numbers. Recommend use ng_mk_cyl_models']);
0202       end
0203    end
0204    if length(cyl_shape)>=4; 
0205       tank_maxh  =cyl_shape(4);
0206    end
0207    if tank_height==0;
0208       is2D = 1;
0209 
0210       %Need some width to let netgen work, but not too much so
0211       % that it meshes the entire region
0212       tank_height = min(tank_radius)/5; % initial extimate
0213       if tank_maxh>0
0214          tank_height = min(tank_height,2*tank_maxh);
0215       end
0216    end
0217 
0218 % ELECTRODE POSITIONS:
0219 %  elec_pos = [n_elecs_per_plane,z_planes]
0220 %     OR
0221 %  elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)
0222 %
0223 % ELECTRODE SHAPES::
0224 %  elec_shape = [width,height, {maxsz}]  % Rectangular elecs
0225 %     OR
0226 %  elec_shape = [radius, {0, maxsz} ]  % Circular elecs
0227 %     maxsz  (OPT)  -> max size of mesh elems (default = courase mesh)
0228 %
0229 % OUTPUT:
0230 %  elecs(i).pos   = [x,y,z]
0231 %  elecs(i).shape = 'C' or 'R'
0232 %  elecs(i).dims  = [radius] or [width,height]
0233 %  elecs(i).maxh  = '-maxh=#' or '';
0234 function [elecs, centres] = parse_elecs(elec_pos, elec_shape, hig, rad, is2D );
0235    if isempty(elec_pos)
0236        elecs = [];
0237        centres = [];
0238        return
0239    end
0240 
0241    if is2D
0242       elec_pos(:,2) = hig/2;
0243    end
0244 
0245    % It never makes sense to specify only one elec
0246    % So elec_pos means the number of electrodes in this case
0247    if size(elec_pos,1) == 1
0248        % Parse elec_pos = [n_elecs_per_plane,z_planes]
0249       n_elecs= elec_pos(1); % per plane
0250       th = space_electrodes('ellipse', n_elecs, rad );
0251 
0252       on_elecs = ones(n_elecs, 1);
0253       el_th = []; 
0254       el_z  = []; 
0255       for i=2:length(elec_pos)
0256         el_th = [el_th; th];
0257         el_z  = [el_z ; on_elecs*elec_pos(i)];
0258       end
0259    else
0260       el_th = elec_pos(:,1)*2*pi/360;
0261       el_z  = elec_pos(:,2);
0262    end
0263       
0264    n_elecs= size(el_z,1); 
0265 
0266    if size(elec_shape,1) == 1
0267       elec_shape = ones(n_elecs,1) * elec_shape;
0268    end
0269 
0270    for i= 1:n_elecs
0271      row = elec_shape(i,:); 
0272      elecs(i) = elec_spec( row, is2D, hig, rad );
0273    end
0274    
0275 % Electrodes are numbered clockwise from top.
0276    centres = [rad(1)*sin(el_th),rad(2)*cos(el_th),el_z];   
0277    for i= 1:n_elecs; elecs(i).pos  = centres(i,:); end
0278 
0279    if n_elecs == 0
0280       elecs= struct([]); % empty
0281    end
0282 
0283 function elec = elec_spec( row, is2D, hig, rad )
0284   if     is2D
0285      if row(1) == 0;
0286         elec.shape = 'P';
0287 % To create a PEM, we make a square and take the corner. This isn't perfect, since
0288 % the elec isn't quite where we asked for it, but that's as good is I can do. I tried
0289 % asking for two rectangles to touch, but that freaks netgen out.
0290         elec.dims  =  [min(rad)/20, hig]; 
0291      else
0292         elec.shape = 'R';
0293         elec.dims  = [row(1),hig];
0294      end
0295   else
0296      if row(1) == 0
0297         elec.shape = 'P' 
0298         elec.dims  = [min(rad)/20, hig/10];
0299      elseif length(row)<2 || row(2) == 0 % Circular electrodes
0300         elec.shape = 'C';
0301         elec.dims  = row(1);
0302      elseif row(2)>0      % Rectangular electrodes
0303         elec.shape = 'R';
0304         elec.dims  = row(1:2);
0305      else
0306         error('negative electrode width');
0307      end
0308   end
0309 
0310   if length(row)>=3 && row(3) > 0
0311      elec.maxh = sprintf('-maxh=%f', row(3));
0312   else
0313      elec.maxh = '';
0314   end
0315 
0316 
0317 function write_header(fid,tank_height,tank_radius,maxsz,extra);
0318    if maxsz==0; 
0319       maxsz = '';
0320    else
0321       maxsz = sprintf('-maxh=%f',maxsz);
0322    end
0323 
0324    extra_ng= '';
0325    for i=1:length(extra)-1
0326       if ~isempty( extra{i} )
0327          extra_ng = sprintf(' %s and (not %s) ', ...
0328             extra_ng,extra{i});
0329       end
0330    end
0331 
0332    fprintf(fid,'#Automatically generated by ng_mk_ellip_models\n');
0333    fprintf(fid,'algebraic3d\n');
0334    fprintf(fid,['solid mainobj_bot= plane(0,0,0;0,0,-1);\n']);
0335    fprintf(fid,['solid mainobj_top= plane(0,0,%f;0,0,1);\n'], ...
0336                  tank_height);
0337    fprintf(fid,'%s\n',extra{end}); % Define extra stuff here
0338    fprintf(fid,'solid cyl=ellipticcylinder (0,0,0;%.4f,0,0;0,%.2f,0); \n', ...
0339             tank_radius);
0340    fprintf(fid,['solid bigcyl= mainobj_top and mainobj_bot and ' ...
0341                 'cyl %s %s;\n'],extra_ng,maxsz);  
0342 
0343 
0344 function write_rect_elec(fid,name,c, dirn,wh,d,maxh)
0345 % writes the specification for a netgen cuboid on fid, named name, centerd on c,
0346 % in the direction given by vector dirn,
0347 % hw = [height, width]  and depth d
0348 % direction is in the xy plane
0349    d= min(d);
0350    w = wh(1); h= wh(2);
0351    dirn(3) = 0; dirn = dirn/norm(dirn);
0352    dirnp = [-dirn(2),dirn(1),0];
0353    dirnp = dirnp/norm(dirnp);
0354 
0355    bl = c - (d/2)* dirn + (w/2)*dirnp - [0,0,h/2];
0356    tr = c + (d/2)* dirn - (w/2)*dirnp + [0,0,h/2];
0357    fprintf(fid,'solid %s  = ', name);
0358    fprintf(fid,' plane (%6.3f,%6.3f,%6.3f;0, 0, -1) and\n', ...
0359            bl(1),bl(2),bl(3));
0360    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0361            bl(1),bl(2),bl(3),-dirn(1),-dirn(2),0);
0362    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0363            bl(1),bl(2),bl(3),dirnp(1),dirnp(2),0);
0364    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;0, 0, 1) and\n', ...
0365            tr(1),tr(2),tr(3));
0366    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0367            tr(1),tr(2),tr(3),dirn(1),dirn(2),0);
0368    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f  )%s;\n', ...
0369            tr(1),tr(2),tr(3),-dirnp(1),-dirnp(2),0,maxh);
0370 
0371 function write_circ_elec(fid,name,c, dirn,rd,ln,maxh)
0372 % writes the specification for a netgen cylindrical rod on fid,
0373 %  named name, centerd on c,
0374 % in the direction given by vector d, radius rd  lenght ln
0375 % direction is in the xy plane
0376 % the direction vector
0377    dirn(3) = 0; dirn = dirn/norm(dirn);
0378 
0379    ln = min(ln);
0380  % This is hard to debug here, why does netgen sometime just fail
0381  % fails for 16 (but no 15 or 17) electrodes
0382  % The 'exact' fix seems to fix this, now. Leave comment above to test
0383    inpt = c - dirn.*(ln/6);
0384    outpt =c + dirn.*(ln/6);
0385 
0386    fprintf(fid,'solid %s  = ', name);
0387    fprintf(fid,'  plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0388          inpt(1),inpt(2),inpt(3),-dirn(1),-dirn(2),-dirn(3));
0389    fprintf(fid,'  plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0390          outpt(1),outpt(2),outpt(3),dirn(1),dirn(2),dirn(3));
0391    fprintf(fid,'  cylinder(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f;%6.3f) %s;\n', ...
0392          inpt(1),inpt(2),inpt(3),outpt(1),outpt(2),outpt(3), rd,maxh);
0393 
0394 
0395 function electrode = pem_from_cem(elecs, electrode, nodes)
0396 % elecs = electrode structure of model, from the parse_elecs function
0397 % electrode = the forward electrode model
0398 % nodes = the coordinates for the nodes
0399 % Can only have one node per electrode so we get a Point Electrode Model.
0400 % Choose the node with the greatest angle, so we atlest pick a consistent
0401 % side of the electrode: NetGen seems to give a random order to the nodes
0402 % in the electrode listing so we can't just pick the first one.
0403 % The nodes aside from those on the edges are not garanteed to be at any
0404 % particular location, so won't be consistent between meshes.
0405 % TODO should probably also adjust contact impedance too: its found later
0406 % by taking the average of the edges around the PEM's node, and those
0407 % will vary for each mesh -- should adjust so all electrodes get a
0408 % consistent effective impedance later.
0409   Ne = length(electrode);
0410   for i = 1:Ne
0411     if elecs(i).shape == 'P'
0412       % find the angles of the nodes for this electrode relative to (0,0)
0413       xy = nodes(electrode(i).nodes,:);
0414       ang = atan2(xy(:,2),xy(:,1));
0415       % if the angles cover more than 180 degrees, must be an angle
0416       % roll-over from -pi to +pi, so take all the negative angles
0417       % and move them up
0418       if (max(ang) - min(ang)) > pi
0419         ang = ang + (ang <0)*2*pi;
0420       end
0421       % choose the counter-clockwise most node only
0422       if size(xy,2) == 3 ; ang = ang - xy(:,3); end
0423       [jnk, ind] = max(ang);
0424       electrode(i).nodes = electrode(i).nodes(ind);
0425     end
0426   end
0427 
0428 
0429 function do_unit_test
0430    sp=1;     subplot(4,4,sp);
0431 % Simple 3D ellipse. Major, minor axes = [1.5, 0.8]. No electrodes
0432     fmdl= ng_mk_ellip_models([1,1.5,0.8],[0],[]);  show_fem(fmdl);
0433 
0434    sp=sp+1;  subplot(4,4,sp);
0435 % Simple 2D cylinder. Axes = [1.5,0.8]. Refined to 0.1
0436     fmdl= ng_mk_ellip_models([0,1.5,0.8,0.1],[],[]); show_fem(fmdl);
0437 
0438    sp=sp+1;  subplot(4,4,sp);
0439 % 3D cylinder. Axes = [1.5,0.8]. 2 planes of 8 elecs with radius 0.1
0440     fmdl= ng_mk_ellip_models([1,1.2,0.8],[8,0.3,0.7],[0.1]); show_fem(fmdl);
0441 
0442    sp=sp+1;  subplot(4,4,sp);
0443 % 3D cylinder. Axes= [1.3,1] = 1. 7 rect elecs with no refinement
0444     fmdl= ng_mk_ellip_models([3,1.3],[7,1],[0.2,0.3]); show_fem(fmdl);
0445 
0446    sp=sp+1;  subplot(4,4,sp);
0447 % 2D cylinder. Axes = [1.2,0.8]. 11 rect elecs with refinement. Rotated 45 degrees
0448     fmdl= ng_mk_ellip_models([0,1.2,0.8],[11],[0.2,0,0.05]); 
0449     th = 45* [2*pi/360];
0450     fmdl.nodes = fmdl.nodes*[cos(th),sin(th);-sin(th),cos(th)]; show_fem(fmdl);
0451 
0452    sp=sp+1;  subplot(4,4,sp);
0453 % 2D cylinder. elecs at 0, 90 and 120 degrees
0454     fmdl= ng_mk_ellip_models([0,1.2,0.8],[0;90;120],[0.2,0,0.03]); show_fem(fmdl);
0455 
0456    sp=sp+1;  subplot(4,4,sp);
0457 % 3D cylinder. Various elecs at 0, 30, 60, 90 in planes
0458     el_pos = [0,0.5;30,1;60,1.5;90,2.0];
0459     el_sz  = [0.2,0,0.1;0.1,0,0.05;0.2,0.2,0.02;0.2,0.4,0.5];
0460     fmdl= ng_mk_ellip_models([3,0.8,1.2],el_pos,el_sz); show_fem(fmdl);
0461 
0462    sp=sp+1;  subplot(4,4,sp);
0463 % Simple 3D cylinder with a ball
0464     extra={'ball','solid ball = sphere(0.3,0.3,1;0.4);'};
0465     [fmdl,mat_idx]= ng_mk_ellip_models([2,1.2,0.8],[8,1],[0.1],extra); 
0466     img= mk_image(fmdl, 1);
0467     img.elem_data(mat_idx{2}) = 2;   show_fem(img);
0468 
0469    sp=sp+1;  subplot(4,4,sp);
0470 % 3D cylinder with a two balls
0471     b1 = 'solid ball1= sphere(0.5, 0.5,1;0.2);';
0472     b2 = 'solid ball2= sphere(0.5,-0.5,1;0.2);';
0473     extra = {'ball1','ball2',[b1,b2]};
0474     [fmdl,mat_idx]= ng_mk_ellip_models([2,1.2,0.8],[8,1],[0.1],extra); 
0475     img= mk_image(fmdl, 1);
0476     img.elem_data(mat_idx{2}) = 2;
0477     img.elem_data(mat_idx{3}) = 0.5;
0478     show_fem(img);
0479      
0480    sp=sp+1;  subplot(4,4,sp);
0481 % Simple 3D cylinder with a ball
0482     extra={'ball','solid ball = sphere(0.3,0.3,1;0.4);'};
0483     [fmdl,mat_idx]= ng_mk_ellip_models([1.15,1.2,0.8],[8,1],[0.1],extra); 
0484     img= mk_image(fmdl, 1);
0485     img.elem_data(mat_idx{2}) = 2;   show_fem(img); view(-30,3);
0486 
0487    sp=sp+1;  subplot(4,4,sp);
0488 % Simple 3D cylinder with a ball
0489     extra={'ball',[ ...
0490        'solid topcut = plane(0,0,1.15;0,0,1);' ...
0491        'solid ball = sphere(0.3,0.3,1;0.4) and topcut;']};
0492     [fmdl,mat_idx]= ng_mk_ellip_models([1.15,1.2,0.8],[8,1],[0.1],extra); 
0493     img= mk_image(fmdl, 1);
0494     img.elem_data(mat_idx{2}) = 2;   show_fem(img); view(-30,3);