dm_mk_elec_nodes

 DM_MK_ELEC_NODES: create node points for dm_mk_fwd_model
 [elec_nodes, refine_nodes] = dm_mk_elec_nodes( n_elecs, elec_width);
 INPUT:
  elec_posn:        vector N x [x,y,{z}] of centre of each electrode
  elec_width:        vector N x 1 of width of each electrode
  refine_level:      vector N x 1 of refine_level for each electrode
      refine_level = 0 -> no refinement
      refine_level = 1 -> more refinement etc

 elec_width and refine_level may be a scalar

 OUTPUT:
  elec_nodes:        cell of matrix N x [x,y,{z}] for each electrode
  refine_nodes:      vector of fixed nodes to add to model to refine it

0001 function [elec_nodes, refine_nodes] = dm_mk_elec_nodes( elec_posn, ...
0002                     elec_width, refine_level);
0003 % DM_MK_ELEC_NODES: create node points for dm_mk_fwd_model
0004 % [elec_nodes, refine_nodes] = dm_mk_elec_nodes( n_elecs, elec_width);
0005 % INPUT:
0006 %  elec_posn:        vector N x [x,y,{z}] of centre of each electrode
0007 %  elec_width:        vector N x 1 of width of each electrode
0008 %  refine_level:      vector N x 1 of refine_level for each electrode
0009 %      refine_level = 0 -> no refinement
0010 %      refine_level = 1 -> more refinement etc
0011 %
0012 % elec_width and refine_level may be a scalar
0013 %
0014 % OUTPUT:
0015 %  elec_nodes:        cell of matrix N x [x,y,{z}] for each electrode
0016 %  refine_nodes:      vector of fixed nodes to add to model to refine it
0017 
0018 % (C) 2008 Andy Adler. License: GPL version 2 or version 3
0019 % $Id: dm_mk_elec_nodes.m 1758 2009-03-26 00:10:34Z aadler $
0020 
0021 [ne,nd]= size(elec_posn);
0022 elec_width=   ones(ne,1).*elec_width;   % expand scalar if required
0023 refine_level= ones(ne,1).*refine_level; % expand scalar if required
0024 
0025 if nd==2
0026    [elec_nodes, refine_nodes]= mk_elec_nodes_2d( ...
0027              elec_posn, elec_width, refine_level, ne); 
0028 elseif nd==3
0029    error('can`t solve 3D problem, yet');
0030 else
0031    error('elec_posn isn`t 2 or 3D');
0032 end
0033 
0034 function [elec_nodes, refine_nodes]= mk_elec_nodes_2d( ...
0035              elec_posn, elec_width, refine_level, ne); 
0036 
0037 % electrodes start top and go clockwise
0038    refine_nodes= [];
0039    for i=1:ne
0040       [ctr, radius] = find_ctr_rad( i, elec_posn);
0041       th= (i-1)*2*pi/ne;
0042       th_delta = elec_width(i)/2/pi/radius;
0043       switch refine_level(i)
0044         case 0,
0045            the= th;
0046            thr= th;
0047         case 1,
0048            the= th + th_delta*[-1;0;1];
0049            thr= th + th_delta*[-3;-1;0;1;3];
0050         case 2,
0051            the= th + th_delta*[-1;0;1];
0052            thr= th + th_delta*[-5;-2;-1;0;1;2;5];
0053         case 3,
0054            the= th + th_delta*[-1;-.5;0;.5;1];
0055            thr= th + th_delta*[-5;-2;-1;-.5;0;.5;1;2;5];
0056         case 4,
0057            the= th + th_delta*[-1;-.5;0;.5;1];
0058            thr= th + th_delta*[-5;-3;-2;-1.5;-1;-.5;0;.5;1;1.5;2;3;5];
0059         otherwise
0060            error('refine level = %d not understood', refine_level(i));
0061       end
0062       this_e_node= radius*[sin(the),cos(the)];
0063       this_e_node= this_e_node + ones(size(this_e_node,1),1)*ctr;
0064       elec_nodes{i} = this_e_node;
0065 
0066       csthr= [sin(thr),cos(thr)];
0067       switch refine_level(i)      
0068         case 0,
0069            % no refine_nodes
0070            refine_new= zeros(0,2);
0071         case 1,
0072            refine_new= [radius*csthr([1,5],:); ...
0073                     .97*radius*csthr([3],:)];
0074         case 2,
0075            refine_new= [radius*csthr([1:2,6:7],:); ...
0076                     .98*radius*csthr([1:7],:); ...
0077                     .99*radius*csthr([1:3:7],:)];
0078         case 3,
0079            refine_new= [radius*csthr([1:2,8:9],:); ...
0080                     .98*radius*csthr([1:9],:); ...
0081                     .95*radius*csthr([1:2:9],:); ...
0082                     .90*radius*csthr([1:4:9],:)];
0083         case 4,
0084            refine_new= [radius*csthr([1:4,10:13],:); ...
0085                     .98*radius*csthr([1:13],:); ...
0086                     .96*radius*csthr([1:2:13],:); ...
0087                     .93*radius*csthr([1,4,6,8,10,13],:); ...
0088                     .90*radius*csthr([2:5:12],:)];
0089         otherwise
0090            error('refine level = %d not understood', refine_level(i));
0091       end
0092       refine_new = refine_new + ones(size(refine_new,1),1)*ctr;
0093       refine_nodes= [refine_nodes; refine_new];
0094    end
0095 
0096 % Find the ctr and radius of this_elec and the two closest
0097 %  ones. This will allow fitting the electrode to the curvature
0098 %  locally.
0099 % Alg: http://www.geocities.com/kiranisingh/center.html
0100 function [ctr, rad] = find_ctr_rad( idx, elec_posn);
0101    nelec= size(elec_posn,1);
0102    idx = rem(idx+[-1,0,1]+nelec-1,nelec)+1;
0103    x = elec_posn(idx,1);
0104    y = elec_posn(idx,2);
0105    s21= x(2)^2 + y(2)^2 - x(1)^2 - y(1)^2;
0106    s31= x(3)^2 + y(3)^2 - x(1)^2 - y(1)^2;
0107    x21 = x(2) - x(1);
0108    x31 = x(3) - x(1);
0109    y21 = y(2) - y(1);
0110    y31 = y(3) - y(1);
0111    n1 = det([s21,y21;s31,y31]);
0112    n2 = det([x21,s21;x31,s31]);
0113    D  = det([x21,y21;x31,y31]);
0114    ctr= [n1,n2]/2/D;
0115    rad= sqrt((x-ctr(1)).^2 + (y-ctr(2)).^2);
0116    if std(rad)/mean(rad)>.001; error('PROBLEM WITH ALGORITHM');end
0117    rad=mean(rad);
0118 %  disp([idx(2),ctr,rad]);