space_electrodes

 space_electrodes: equal spacing of electrodes around body
 [thrad,thdeg] = space_electrodes(spacing_type, n_elecs, params);

 OUTPUT: thrad,thdeg = spacing angle (rad,degrees) from centre

 INPUT: n_elecs = number of electrodes
  spacing_type = 'ellipse':
    params = eliptical radii [r_x r_y]

0001 function [thrad,thdeg] = space_electrodes(spacing_type, n_elecs, params);
0002 % space_electrodes: equal spacing of electrodes around body
0003 % [thrad,thdeg] = space_electrodes(spacing_type, n_elecs, params);
0004 %
0005 % OUTPUT: thrad,thdeg = spacing angle (rad,degrees) from centre
0006 %
0007 % INPUT: n_elecs = number of electrodes
0008 %  spacing_type = 'ellipse':
0009 %    params = eliptical radii [r_x r_y]
0010 
0011 % (C) 2024 Andy Adler . License: GPL v2 or v3.
0012 % $Id: space_electrodes.m 7002 2024-11-24 13:11:35Z aadler $
0013 
0014 if ischar(spacing_type) && strcmp(spacing_type,'UNIT_TEST'); do_unit_test; return; end
0015 
0016 switch spacing_type
0017    case 'ellipse';
0018       thrad = ellip_space_elecs( n_elecs, params );
0019    otherwise
0020       error('spacing_type not recognized')
0021 end
0022 
0023 thdeg = 180/pi*thrad;
0024 
0025 % equally space n_elecs around an ellipse of outer radius rad(1),rad(2)
0026 function th = ellip_space_elecs( n_elecs, rad )
0027    % The radius is the integral of sqrt((r1*sin(th))^2 + (r2*cos(th))^2)
0028    %  This integral is the incomplete_elliptic_integral(th, 1-r2/r1)*sqrt(r1)
0029    %  Unfortunately, STUPID MATLAB, doesn't have incomplete elliptic integrals
0030    %  by default. So, rather than install a toolkit for it, we integrate numerically.
0031    if n_elecs==0; th=[]; return; end
0032    
0033    th = linspace(0,2*pi, 100*(n_elecs)); th(1)=[]; % Accuracy to 100x spacing
0034    len = cumsum( sqrt( rad(1)*cos(th).^2 + rad(2)*sin(th).^2 ) );
0035    len = len/max(len);
0036    xi = linspace(0,1,n_elecs+1); xi(1)= []; xi(end)=[];
0037    yi = interp1(len,th,xi);
0038 
0039    th= [0;yi(:)];
0040    for exact = 0:3;
0041       eth = exact/2*pi;
0042       ff = abs(th-eth)<1e-3;
0043       th(ff) = eth;
0044    end
0045 
0046 
0047 function do_unit_test
0048   th = space_electrodes('ellipse',16,[1,1]); 
0049   unit_test_cmp('circle', th(1), 0, 1e-14);
0050   unit_test_cmp('circle', diff(th), 2*pi/16, 1e-14);
0051 
0052   th = space_electrodes('ellipse',16,[1,2]); 
0053   unit_test_cmp('ellipse [1,2]', th(1), 0, 1e-14);
0054   unit_test_cmp('ellipse [1,2]', [th(1), min(diff(th)) max(diff(th))], ...
0055         [ 0   0.340379981348740   0.462371789340561 ], 1e-14);