0001 function [elec_nodes, refine_nodes] = dm_mk_elec_nodes( elec_posn, ...
0002 elec_width, refine_level);
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021 [ne,nd]= size(elec_posn);
0022 elec_width= ones(ne,1).*elec_width;
0023 refine_level= ones(ne,1).*refine_level;
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
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
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
0097
0098
0099
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