stim_pattern_geophys

 STIM_PATTERN_GEOPHYS: Create Geophysical Stimulation Patterns

 stim= stim_pattern_geophys( n_elec, type, options_list )
  n_elec = Number of electrodes (in single line)
  type   = name of stimulation pattern (string)
  options_list = options for the pattern

  options = {'current', 0.1}  % 0.1 Amps current (default 1)
  options = {'gain', 2.0}     % 2.0 gain         (default 1)
  options = {'spacings',[1,2,3,4]} % list of spacings 'a'
  options = {'multiples',[1,2,3,4]} % list of multiples 'n'
  options = {'reciprocal_meas',1} % do reciprocal (default 0);
  options = {'circumferential_meas',1} % do a circumferential protocol (default 0 = in line protocol);

 Supported stimulation pattern types are:
  'Wenner'       - a Wenner array
   [0,3,1,2]*a + 1; % as [s+ s- m+ m-] -->  s+  m+ m-  s-
   same as Schlumberger with n=1
 'Schlumberger' - a Schlumberger array
   [0,(2*n+1)*a,n*a,(n+1)*a]+1;
 'DipoleDipole' - a dipole-dipole array
   [0,a,a+n*a,2*a+n*a]+1; % -->  s+ s-  m+     m-

0001 function [stim,S]= stim_pattern_geophys( n_elec, pat_type,  options )
0002 % STIM_PATTERN_GEOPHYS: Create Geophysical Stimulation Patterns
0003 %
0004 % stim= stim_pattern_geophys( n_elec, type, options_list )
0005 %  n_elec = Number of electrodes (in single line)
0006 %  type   = name of stimulation pattern (string)
0007 %  options_list = options for the pattern
0008 %
0009 %  options = {'current', 0.1}  % 0.1 Amps current (default 1)
0010 %  options = {'gain', 2.0}     % 2.0 gain         (default 1)
0011 %  options = {'spacings',[1,2,3,4]} % list of spacings 'a'
0012 %  options = {'multiples',[1,2,3,4]} % list of multiples 'n'
0013 %  options = {'reciprocal_meas',1} % do reciprocal (default 0);
0014 %  options = {'circumferential_meas',1} % do a circumferential protocol (default 0 = in line protocol);
0015 %
0016 % Supported stimulation pattern types are:
0017 %  'Wenner'       - a Wenner array
0018 %   [0,3,1,2]*a + 1; % as [s+ s- m+ m-] -->  s+  m+ m-  s-
0019 %   same as Schlumberger with n=1
0020 % 'Schlumberger' - a Schlumberger array
0021 %   [0,(2*n+1)*a,n*a,(n+1)*a]+1;
0022 % 'DipoleDipole' - a dipole-dipole array
0023 %   [0,a,a+n*a,2*a+n*a]+1; % -->  s+ s-  m+     m-
0024 
0025 % (C) 2011 Andy Adler. License: GPL version 2 or version 3
0026 % $Id: stim_pattern_geophys.m 5225 2016-05-08 06:35:38Z alistair_boyle $
0027 
0028 if ischar(n_elec) && strcmp(n_elec,'UNIT_TEST'); do_unit_test; return; end
0029 if nargin<3; options = {}; end
0030 options = parse_options(options);
0031 
0032 switch(upper(pat_type))
0033     case 'WENNER';          S = stim_pattern_wenner(n_elec, options);
0034     case 'SCHLUMBERGER';    S = stim_pattern_schlumberger(n_elec, options);
0035     case 'DIPOLEDIPOLE';    S = stim_pattern_dipoledipole(n_elec, options);
0036   otherwise;
0037     error('No pattern of type "%s" available', upper(pat_type));
0038 end
0039 stim = stim_meas_list( S, n_elec, options.current, options.gain); 
0040 
0041 function pp= parse_options(opts)
0042 % Set defaults here
0043 %   pp.spacings = [1,2,4,8,16,32,64,128,256,1024,2048];
0044   pp.spacings= 1:256;
0045   pp.multiples= ones(1,256);
0046   pp.reciprocal_meas = 0;
0047   pp.circumferential_meas = 0;
0048   pp.current = 1;
0049   pp.gain = 1;
0050 
0051   for i= 1:2:length(opts)-1;
0052     pp.( lower( opts{i} ) ) = opts{i+1};
0053   end
0054 
0055 function stim = stim_pattern_wenner(n_elec, options)
0056   stim=[];
0057   for a= options.spacings;
0058     block =  [0,3,1,2]*a + 1;
0059     if options.circumferential_meas 
0060         if max(block(1))>n_elec
0061            break;
0062         end
0063         n_meas= n_elec-1;
0064     else
0065         if max(block)>n_elec
0066            break;
0067         end
0068         n_meas= n_elec-max(block);
0069     end 
0070     for i=(0:n_meas)
0071        stim = [stim; i+block];
0072        if options.reciprocal_meas
0073           stim = [stim; i+block([3,4,1,2])];
0074        end
0075     end
0076     if options.circumferential_meas 
0077        stim(stim>n_elec)= stim(stim>n_elec)-n_elec;
0078        stim= reshape(stim,[],4);
0079     end
0080 
0081   end
0082 
0083 function stim = stim_pattern_schlumberger(n_elec, options)
0084   stim=[];
0085   a= options.spacings;
0086   n= options.multiples;
0087   if length(a)~=length(n)
0088    error('spacings and multiples must be vector with a same size');
0089   end
0090   for k= 1:length(a);
0091       block =  [0,round((2*n(k)+1)*a(k)),round(n(k)*a(k)),round((n(k)+1)*a(k))]+1;
0092     if options.circumferential_meas 
0093         if max(block(1))>n_elec
0094            break;
0095         end
0096         n_meas= n_elec-1;
0097     else
0098         if max(block)>n_elec
0099            break;
0100         end
0101         n_meas= n_elec-max(block);
0102     end 
0103     for i=(0:n_meas)
0104        stim = [stim; i+block];
0105        if options.reciprocal_meas
0106           stim = [stim; i+block([3,4,1,2])];
0107        end
0108     end
0109     if options.circumferential_meas 
0110        stim(stim>n_elec)= stim(stim>n_elec)-n_elec;
0111        stim= reshape(stim,[],4);
0112     end
0113   end  
0114   
0115   
0116 function stim = stim_pattern_dipoledipole(n_elec, options)
0117   stim=[];
0118   a= options.spacings;
0119   n= options.multiples;
0120   if length(a)~=length(n)
0121    error('spacings and multiples must be vector with a same size');
0122   end
0123   for k= 1:length(a);
0124       block =  [0,a(k),a(k)+n(k)*a(k),2*a(k)+n(k)*a(k)]+1;
0125     if options.circumferential_meas 
0126         if max(block(1))>n_elec
0127            break;
0128         end
0129         n_meas= n_elec-1;
0130     else
0131         if max(block)>n_elec
0132            break;
0133         end
0134         n_meas= n_elec-max(block);
0135     end 
0136     for i=(0:n_meas)
0137        stim = [stim; i+block];
0138        if options.reciprocal_meas
0139           stim = [stim; i+block([3,4,1,2])];
0140        end
0141     end
0142     if options.circumferential_meas
0143        while any(any(stim>n_elec))
0144           stim(stim>n_elec)= stim(stim>n_elec)-n_elec;
0145        end
0146        stim= reshape(stim,[],4);
0147     end
0148   end  
0149   
0150 function do_unit_test
0151   unit_test_wenner;
0152   unit_test_schlumberger;
0153   unit_test_dipoledipole;
0154 
0155 function unit_test_wenner
0156   stim= stim_pattern_wenner( 13, parse_options({}) ); 
0157   test1= [1,4,2,3; 2,5,3,4;
0158 3,6,4,5; 4,7,5,6; 5,8,6,7; 6,9,7,8; 7,10,8,9; 8,11,9,10;
0159 9,12,10,11; 10,13,11,12; 1,7,3,5; 2,8,4,6; 3,9,5,7;
0160 4,10,6,8; 5,11,7,9; 6,12,8,10; 7,13,9,11; 1,10,4,7; 2 11 5 8; 3 12 6 9;
0161 4 13 7 10; 1,13,5,9];
0162   unit_test_cmp('WENNER #1', stim, test1);
0163   stim= stim_pattern_geophys( 13,'wenner', {});
0164   unit_test_cmp('WENNER #2', stim, stim_meas_list(test1));
0165 
0166 
0167   stim= stim_pattern_geophys( 5,'wenner', {'reciprocal_meas',0});
0168   unit_test_cmp('WENNER #3a',stim, stim_meas_list([1,4,2,3;2,5,3,4]));
0169   stim= stim_pattern_geophys( 5,'wenner', {'reciprocal_meas',1});
0170   unit_test_cmp('WENNER #3b',stim, stim_meas_list([1,4,2,3;2,3,1,4;2,5,3,4;3,4,2,5]));
0171 
0172   stim= stim_pattern_geophys( 5,'wenner', {'current',0.1});
0173   unit_test_cmp('WENNER #4a',stim, stim_meas_list([1,4,2,3;2,5,3,4],5,0.1));
0174   stim= stim_pattern_geophys( 5,'wenner', {'gain',0.1});
0175   unit_test_cmp('WENNER #4b',stim, stim_meas_list([1,4,2,3;2,5,3,4],5,1,0.1));
0176 
0177   stim= stim_pattern_geophys( 5,'wenner', {'spacings',[1,2,3]});
0178   unit_test_cmp('WENNER #5a',stim, stim_meas_list([1,4,2,3;2,5,3,4]));
0179   stim= stim_pattern_geophys( 7,'wenner', {'spacings',[2,1,3]});
0180   unit_test_cmp('WENNER #5b',stim, stim_meas_list(...
0181        [1,7,3,5; 1,4,2,3; 2,5,3,4; 3,6,4,5; 4,7,5,6]));
0182   stim= stim_pattern_geophys( 8,'wenner', {'spacings',[1],'circumferential_meas',1});
0183   unit_test_cmp('WENNER #5c',stim, stim_meas_list(...
0184        [1,4,2,3; 2,5,3,4; 3,6,4,5; 4,7,5,6;5,8,6,7;6,1,7,8;7,2,8,1;8,3,1,2]));
0185   stim= stim_pattern_geophys( 8,'wenner', {'spacings',[1,2],'circumferential_meas',1});
0186   unit_test_cmp('WENNER #5d',stim, stim_meas_list(...
0187        [1,4,2,3; 2,5,3,4; 3,6,4,5; 4,7,5,6;5,8,6,7;6,1,7,8;7,2,8,1;8,3,1,2;1,7,3,5;2,8,4,6;3,1,5,7;4,2,6,8;...
0188        5,3,7,1;6,4,8,2;7,5,1,3;8,6,2,4]));
0189 
0190 
0191 function unit_test_schlumberger
0192     stim= stim_pattern_schlumberger( 13, parse_options({}) );
0193     test1= [1,4,2,3; 2,5,3,4;
0194         3,6,4,5; 4,7,5,6; 5,8,6,7; 6,9,7,8; 7,10,8,9; 8,11,9,10;
0195         9,12,10,11; 10,13,11,12; 1,7,3,5; 2,8,4,6; 3,9,5,7;
0196         4,10,6,8; 5,11,7,9; 6,12,8,10; 7,13,9,11; 1,10,4,7; 2 11 5 8; 3 12 6 9;
0197         4 13 7 10; 1,13,5,9];
0198     unit_test_cmp('SCHLUMBERGER #1', stim, test1);
0199     stim= stim_pattern_geophys( 13,'schlumberger', {});
0200     unit_test_cmp('SCHLUMBERGER #2', stim, stim_meas_list(test1));
0201 
0202 
0203     stim= stim_pattern_geophys( 5,'schlumberger', {'reciprocal_meas',0});
0204     unit_test_cmp('SCHLUMBERGER #3a',stim, stim_meas_list([1,4,2,3;2,5,3,4]));
0205     stim= stim_pattern_geophys( 5,'schlumberger', {'reciprocal_meas',1});
0206     unit_test_cmp('SCHLUMBERGER #3b',stim, stim_meas_list([1,4,2,3;2,3,1,4;2,5,3,4;3,4,2,5]));
0207 
0208     stim= stim_pattern_geophys( 5,'schlumberger', {'current',0.1});
0209     unit_test_cmp('SCHLUMBERGER #4a',stim, stim_meas_list([1,4,2,3;2,5,3,4],5,0.1));
0210     stim= stim_pattern_geophys( 5,'schlumberger', {'gain',0.1});
0211     unit_test_cmp('SCHLUMBERGER #4b',stim, stim_meas_list([1,4,2,3;2,5,3,4],5,1,0.1));
0212 
0213     stim= stim_pattern_geophys( 6,'schlumberger', {'spacings',[1,1],'multiples',[1,2]});
0214     unit_test_cmp('SCHLUMBERGER #5a',stim, stim_meas_list([1,4,2,3;2,5,3,4;3,6,4,5;1,6,3,4]));
0215     stim= stim_pattern_geophys( 11,'schlumberger', {'spacings',[1,2],'multiples',[3,2]});
0216     unit_test_cmp('SCHLUMBERGER #5b',stim, stim_meas_list(...
0217         [1,8,4,5;2,9,5,6;3,10,6,7;4,11,7,8;1,11,5,7]));
0218    
0219 % options.spacing= [1 1 1 2 3 4 6 8 8 11 12 14 17];
0220 % options.multiples= [1 2 3 2 5/3 6/4 7/6 1 10/8 1 13/12 15/14 1];
0221 
0222 function unit_test_dipoledipole
0223     stim= stim_pattern_dipoledipole( 13, parse_options({}) );
0224     test1= [1,2,3,4; 2,3,4,5;
0225         3,4,5,6; 4,5,6,7; 5,6,7,8; 6,7,8,9; 7,8,9,10; 8,9,10,11;
0226         9,10,11,12; 10,11,12,13; 1,3,5,7; 2,4,6,8; 3,5,7,9;
0227         4,6,8,10; 5,7,9,11; 6,8,10,12; 7,9,11,13; 1,4,7,10; 2 5 8 11; 3 6 9 12;
0228         4 7 10 13; 1,5,9,13];
0229     unit_test_cmp('DIPOLEDIPOLE #1', stim, test1);
0230     stim= stim_pattern_geophys( 13,'dipoledipole', {});
0231     unit_test_cmp('DIPOLEDIPOLE #2', stim, stim_meas_list(test1));
0232 
0233 
0234     stim= stim_pattern_geophys( 5,'dipoledipole', {'reciprocal_meas',0});
0235     unit_test_cmp('DIPOLEDIPOLE #3a',stim, stim_meas_list([1,2,3,4;2,3,4,5]));
0236     stim= stim_pattern_geophys( 5,'dipoledipole', {'reciprocal_meas',1});
0237     unit_test_cmp('DIPOLEDIPOLE #3b',stim, stim_meas_list([1,2,3,4;3,4,1,2;2,3,4,5;4,5,2,3]));
0238 
0239     stim= stim_pattern_geophys( 5,'dipoledipole', {'current',0.1});
0240     unit_test_cmp('DIPOLEDIPOLE #4a',stim, stim_meas_list([1,2,3,4;2,3,4,5],5,0.1));
0241     stim= stim_pattern_geophys( 5,'dipoledipole', {'gain',0.1});
0242     unit_test_cmp('DIPOLEDIPOLE #4b',stim, stim_meas_list([1,2,3,4;2,3,4,5],5,1,0.1));
0243 
0244     stim= stim_pattern_geophys( 6,'dipoledipole', {'spacings',[1,1],'multiples',[1,2]});
0245     unit_test_cmp('DIPOLEDIPOLE #5a',stim, stim_meas_list([1,2,3,4;2,3,4,5;3,4,5,6;1,2,4,5;2,3,5,6]));
0246     stim= stim_pattern_geophys( 9,'dipoledipole', {'spacings',[1,2],'multiples',[3,2]});
0247     unit_test_cmp('DIPOLEDIPOLE #5b',stim, stim_meas_list(...
0248         [1,2,5,6;2,3,6,7;3,4,7,8;4,5,8,9;1,3,7,9]));
0249     
0250     
0251