set_3d_currents

function [I,Ib]=set_3d_currents(protocol,elec,vtx,gnd_ind,no_pl);

This function sets current patterns in a system with (no_pl) planes of 
equal number of electrodes according to "opposite" or "adjacent" protocols, 
or their 3D similar.



protocol= The selected protocol '{op}' or '{ad}'
elec    = The electrodes (only the number of electrodes is used)
vtx     = The vertices
gnd_ind = the index of the ground node
no_pl   = The number of planes
Ib      = The current patterns
I       = The RHS vectors, i.e., the current patterns padded with zeroes

0001 function [I,Ib] = set_3d_currents(protocol,elec,vtx,gnd_ind,no_pl);
0002 %function [I,Ib]=set_3d_currents(protocol,elec,vtx,gnd_ind,no_pl);
0003 %
0004 %This function sets current patterns in a system with (no_pl) planes of
0005 %equal number of electrodes according to "opposite" or "adjacent" protocols,
0006 %or their 3D similar.
0007 %
0008 %
0009 %
0010 %protocol= The selected protocol '{op}' or '{ad}'
0011 %elec    = The electrodes (only the number of electrodes is used)
0012 %vtx     = The vertices
0013 %gnd_ind = the index of the ground node
0014 %no_pl   = The number of planes
0015 %Ib      = The current patterns
0016 %I       = The RHS vectors, i.e., the current patterns padded with zeroes
0017 
0018 
0019 [vr,vc] = size(vtx);
0020    
0021 [el_no,q] = size(elec); 
0022 
0023 el_pp = el_no/no_pl;
0024 
0025 a=1:el_no;
0026 
0027 X = reshape(a,el_pp,no_pl)';
0028 
0029 if protocol == '{op}'
0030    
0031    Ib = [];
0032    
0033     for i=1:no_pl
0034    
0035        this_plane = X(i,:);
0036       
0037       for j=this_plane(1):this_plane(8)
0038          
0039          Ip = zeros(el_no,1);
0040          Ip(j) = 1;
0041          Ip(j+ el_pp/2) = -1;
0042          Ib = [Ib,Ip];
0043       end
0044    
0045    end 
0046    
0047     
0048    Is_supl = zeros(vr,size(Ib,2));
0049 
0050    I = [Is_supl;Ib];
0051 
0052    I(gnd_ind,:) = 0;
0053 
0054 elseif protocol == '{ad}'
0055    
0056    Ib = [];
0057    
0058    for i=1:no_pl
0059    
0060      this_plane = X(i,:); 
0061    
0062         for j=this_plane(1):this_plane(el_pp-1)
0063            
0064            Ip = zeros(el_no,1);
0065            Ip(j) = 1;
0066            Ip(j+1) = -1;
0067            Ib =[Ib,Ip];
0068            
0069            if j==this_plane(el_pp-1) %the ring pattern
0070               
0071               Ip = zeros(el_no,1);
0072               
0073               Ip(j+1) = 1;
0074               Ip(this_plane(1)) = -1;
0075               Ib = [Ib,Ip];
0076            end
0077            
0078            
0079         end
0080         
0081      end
0082      
0083     Is_supl = zeros(vr,size(Ib,2));
0084 
0085     I = [Is_supl;Ib];
0086 
0087     I(gnd_ind,:) = 0;
0088 
0089 else
0090    error(['protocol ',protocol,' is not recognized']);
0091 end %protocol
0092          
0093 
0094 
0095 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0096 % This is part of the EIDORS suite.
0097 % Copyright (c) N. Polydorides 2003
0098 % Copying permitted under terms of GNU GPL
0099 % See enclosed file gpl.html for details.
0100 % EIDORS 3D version 2.0
0101 % MATLAB version 5.3 R11
0102 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%