function [R,s,sensnum,th,rings] = sensor_ring(r,th,rings,version);
%SENSOR_RING - Generate rings of sensors about a sphere roughly uniformly spaced.
% function [R,s,sensnum,th,rings] = sensor_ring(r,th,rings,version);
% For radius 'r', angle separation 'theta' (radians), and 
%  'rings' number  of sensor rings:
% Generate R (mx3) location and s (mx3) orientation of 
%  radially oriented sensors.
% If rings is given as a string, 'h' or 'f', then theta is adjusted to give
%   exactly half or full head coverage. Adjusted values are returned in 
%   the output as th and rings.
%
% The first sensor is on the z-axis, at [0 0 r];  Each ring is incremented
%  theta radians down from the z-axis, and the first sensor of each ring is
%  positioned on the x-axis (in spherical coords, the first sensor is at
%  phi=0). 
% The desired sensor spacing is set by the distance from sensor 1 to 2.
% Sensnum is a vector of the number of sensors per ring.  The sensors are
%  spaced about each ring such that the distance between sensors in a given
%  ring are as close as possible to the desired spacing.
% Optionally: if theta has two components, then the first component is the
%   initial offset from the z-axis for the first ring, and the second
%   component is increment between rings. (Option has no effect for 'old'
%   mode, but the second component will be properly interpreted as the
%   increment.) 
%
% Example:
%   [R,s,sensnum] = sensor_ring(0.12,12*(pi/180),4); produces a simulated
%   version of the BTi configuration, which has four rings of sensors,
%   including the center ring of 1 sensor, then 6, then 12, then 18, for a
%   total of 37 sensors.  Each ring is separated from the adjacent ring by
%   about 12 degrees, and the entire array sits on a virtual sphere of radius
%   12 cm (0.12 meters).
%
% OTHER ANGLES AND MORE RINGS MAY GENERATE DIFFERENT PATTERNS THAN THE
%  OLD SENSOR_RING PROGRAM.  The old routine may be activated by the 
%  string 'old' in version, which builds rings of 1, 6, 12, 18, 24, ...
%  Default is the new routine, which builds each ring with a suitable number
%  of sensors to maintain the desired sensor spacing.  For the BTi pattern
%  example above, there is virtually no difference between old and new modes.
%
% See also CONN_DIST for the distance between sensor points, 
%  and SENSOR_SPACING for a convenient method of calling this routine.

%<autobegin> ---------------------- 26-May-2004 11:34:20 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Forward Modeling
%
% At Check-in: $Author: Mosher $  $Revision: 14 $  $Date: 5/26/04 10:02a $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 24-May-2004
% 
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
% 
% See BrainStorm website at http://neuroimage.usc.edu for further information.
% 
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 26-May-2004 11:34:20 -----------------------

% ----------------------------- Script History ---------------------------------
% Author John C. Mosher, Ph.D. ~1992
% June 8, 1994 JCM: modified from old sensor_ring 6/8/94 to give more uniform
%  spacing over hemisphere .
% July 28, 1994 JCM: added two element theta option
% January 18, 1995 JCM: reversed ordering around each ring to correspond to
%  HRL numbering scheme.
% March 2, 1995 JCM: adjusted 'new' routine to interleave subsequent rings
%  and therefore make the true distances somewhat more comparable.
%  Also added 'new' rings option 'half' and 'full', for upperhemisphere
%  coverage and full sphere coverage.
% August 31, 1995, bug fix with 'full' coverage, was only giving 'half'.
% 19-May-2004 JCM Comments Cleaning
% ----------------------------- Script History ---------------------------------


if(exist('version')~=1),	% user gave no version
  version = 'new';		% default is not old
end

if(strcmp(lower(version),'old')), % user wants old version
  OLD_MODE = 1;
else
  OLD_MODE = 0;			% default mode
end

th_initial = 0;			% initial angle from z-axis
if(length(th)==2),		% user using option
  th_initial = th(1);
  th = th(2);
end

% design first for unit length, then scale to radius

if(OLD_MODE),
  
  R(1,:) = [0 0 1];		% first sensor, "ring" one
  degrees = 180/pi;		% conversion

  deg = 60;			% phi increment in second ring of six sensors

  for j = 2:rings
    phi = -[0:deg:359]'./degrees;% phi index this ring, 1/18/95 added '-'sign
    thj = (j-1)*th;		% next theta increment

    % append next ring of sensors
    R = [R;...
	[cos(phi)*sin(thj) ...
	sin(phi)*sin(thj) ones(length(phi),1)*cos(thj)]];
    deg = 60/j;			% phi increment next ring of sensors
  end
  
  sensnum = [1 [1:(rings-1)]*6]; % old pattern of rings

else
  
  if(isstr(rings)),		% user wants us to calculate parameters
    % first figure for half hemisphere coverage, but not exceeding.
    tmp1 = round((pi/2 - th_initial)/th); % number of rings, except top
    th = ((pi/2)-th_initial)/tmp1;		% new value
    nrings = tmp1;		% number of rings for half coverage
    if(th_initial==0),		% need a top single point also
      nrings = nrings + 1;
    end
    if(rings(1)=='f'), 		% want full coverage
      nrings = tmp1 + nrings; 	% including bottom set ring of 1 sensor
    end
    % now map number of rings over the string rings
    rings = nrings;
  end				% if rings given as string
      
  [x1,y1,z1] = sph2cart(0,pi/2 - th_initial,1); % first sensor, first ring
  [x2,y2,z2] = sph2cart(0,pi/2 - th_initial - th,1); %first sensor, second ring
  d = norm([x2,y2,z2] - [x1,y1,z1]); % distance between, desired distance

  sensnum = zeros(1,rings);	% number of sensors per ring
  R = [];
  for j = 1:rings,
    thj = (j-1)*th + th_initial;% theta this ring
    circum = 2*pi*sin(thj);	% circumference of circle
    sens = max(round(circum/d),1);% number of sensors this ring, minimum 1
    sensnum(j) = sens;
    ph = -[0:(sens-1)]'/sens*2*pi;	% phi increments this ring,
    % 1/18/95 added '-' sign to ph
    % 3/2/95 added ph/2 factor for odd numbered rings
    if(rem(j,2) & (sens > 1)),		% odd ring of more than one sensor
      ph = ph + ph(2)/2;	% bump half an increment
    end

    % append next ring of sensors
    R = [R;...
	[cos(ph)*sin(thj) ...
	sin(ph)*sin(thj) ones(length(ph),1)*cos(thj)]];
  end

end

s = R;				% radial orientation of unit length

R = R*r;			% scale to radius

return
