function [R,s,u1,u2] = sensor_planar(R1,baseline);
%SENSOR_PLANAR - Generate planar gradiometers
% function [R,s,u1,u2] = sensor_planar(R1,baseline);
% Given sensor sites in R1, generate planar point gradiometers, 
%  oriented in directions orthogonal to R1 with a separation of baseline.
%  Optionally, if R1 has six columns, the last three are presumed to be the
%  desired planar normal.
% 
% For example, the Neuromag 122 system can be roughly approximated by
% spacing the sensors approximately 3 cm apart on a 10.5 cm shell, with
% a 1.5 cm baseline, for a total of 78 sensors about the entire hemisphere.
% R1 = sensor_ring(0.105,16.6*pi/180,6);  % generate sensor sites
% [R,s] = sensor_planar(R1,[0.015]);
%
% The first three columns of R and s will be the positions of the coils
%   to one side of R1.  The next three columns correspond to the second coil in
%   each sensor, a distance of baseline away.  The combined set for a channel.
%
% If R1 has mR rows, then R will have mR*2 rows.  Each pair of rows
%  corresponds to the primary and secondary channels at each site in R1.
%  where the secondary channel is oriented orthogonal
%  to the primary channel.
% s will give the unit radial direction of R1/norm(R1).
%
% U1 and U2 are calculated from spherical coordinates.
% Each row of u1 will give the theta direction of the primary channel.
% Each row of u2 will give the phi direction of the secondary channel.
% Each row of u1 and u2 is orthogonal to each other and 
%  the corresponding direction in R1.
%
% See also SENSOR_RING, SENSOR_AXIAL, SENSOR_122

%<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 ---------------------------------
% 1994 by John C. Mosher, Ph.D.
% 19-May-2004 JCM Comments Cleaning
% ----------------------------- Script History ---------------------------------

mR = size(R1,1);		% number of sensors
nCoils = 2;			% two per channel in planar config


% generate primary and secondary directions
if(size(R1,2) == 6),		% user has supplied desired orientation
  [ph,th,ra] = cart2sph(R1(:,4),R1(:,5),R1(:,6));
  s1 = R1(:,4:6);
else
  % get spherical coordinates for rho, theta, phi
  [ph,th,ra] = cart2sph(R1(:,1),R1(:,2),R1(:,3));
  % generate radial vectors
  s1 = R1 ./ (ra(:,[1 1 1])); 	% normalized radial orientation
end

% th is returned as elevation from xy plane.  Convert to angle from z-axis.
th = pi/2 - th;

% generate primary in theta direction
u1 = [cos(th).*cos(ph) cos(th).*sin(ph) -sin(th)];

% generate secondary in the phi direction
u2 = [-sin(ph) cos(ph) zeros(mR,1)];

% primary channels
pR = [(R1(:,1:3) + baseline*u1/2) (R1(:,1:3) - baseline*u1/2)];
% secondary channels
sR = [(R1(:,1:3) + baseline*u2/2) (R1(:,1:3) - baseline*u2/2)];

% interleave
R([[1:2:(mR*2)] [2:2:(mR*2)]],:) = [pR;sR];
s([[1:2:(mR*2)] [2:2:(mR*2)]],:) = [s1 s1;s1 s1];

return


