function grad = yokogawa2grad(hdr)

% YOKOGAWA2GRAD converts the position and weights of all coils that
% compromise a gradiometer system into a structure that can be used
% by FieldTrip. This implementation uses the old "yokogawa" toolbox.
%
% See also CTF2GRAD, BTI2GRAD, FIF2GRAD, MNE2GRAD, ITAB2GRAD,
% FT_READ_SENS, FT_READ_HEADER

% Copyright (C) 2005-2008, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if ~ft_hastoolbox('yokogawa')
    ft_error('cannot determine whether Yokogawa toolbox is present');
end

if isfield(hdr, 'label')
  label = hdr.label; % keep for later use
end

if isfield(hdr, 'orig')
  hdr = hdr.orig; % use the original header, not the FieldTrip header
end

% The "channel_info" contains
% 1  channel number, zero offset
% 2  channel type, type of gradiometer
% 3  position x (in m)
% 4  position y (in m)
% 5  position z (in m)
% 6  orientation of first coil (theta in deg)
% 7  orientation of first coil (phi in deg)
% 8  orientation from the 1st to 2nd coil for gradiometer (theta in deg)
% 9  orientation from the 1st to 2nd coil for gradiometer (phi in deg)
% 10 coil size (in m)
% 11 baseline (in m)

handles    = definehandles;
isgrad     = (hdr.channel_info(:,2)==handles.AxialGradioMeter | ...
              hdr.channel_info(:,2)==handles.PlannerGradioMeter | ...
              hdr.channel_info(:,2)==handles.MagnetoMeter | ...
              hdr.channel_info(:,2)==handles.RefferenceAxialGradioMeter);
% reference channels are excluded because the positions are not specified
%              hdr.channel_info(:,2)==handles.RefferencePlannerGradioMeter
%              hdr.channel_info(:,2)==handles.RefferenceMagnetoMeter
isgrad_handles = hdr.channel_info(isgrad,2);
ismag      = (isgrad_handles(:)==handles.MagnetoMeter | isgrad_handles(:)==handles.RefferenceMagnetoMeter);
grad.coilpos   = hdr.channel_info(isgrad,3:5)*100;    % cm
grad.unit  = 'cm';

% Get orientation of the 1st coil
ori_1st   = hdr.channel_info(find(isgrad),[6 7]);
% polar to x,y,z coordinates
ori_1st = ...
  [sin(ori_1st(:,1)/180*pi).*cos(ori_1st(:,2)/180*pi) ...
  sin(ori_1st(:,1)/180*pi).*sin(ori_1st(:,2)/180*pi) ...
  cos(ori_1st(:,1)/180*pi)];
grad.coilori = ori_1st;

% Get orientation from the 1st to 2nd coil for gradiometer
ori_1st_to_2nd   = hdr.channel_info(find(isgrad),[8 9]);
% polar to x,y,z coordinates
ori_1st_to_2nd = ...
  [sin(ori_1st_to_2nd(:,1)/180*pi).*cos(ori_1st_to_2nd(:,2)/180*pi) ...
  sin(ori_1st_to_2nd(:,1)/180*pi).*sin(ori_1st_to_2nd(:,2)/180*pi) ...
  cos(ori_1st_to_2nd(:,1)/180*pi)];
% Get baseline
baseline = hdr.channel_info(isgrad,size(hdr.channel_info,2));

% Define the location and orientation of 2nd coil
info = hdr.channel_info(isgrad,2); 
for i=1:sum(isgrad)
  if (info(i) == handles.AxialGradioMeter || info(i) == handles.RefferenceAxialGradioMeter )
    grad.coilpos(i+sum(isgrad),:) = [grad.coilpos(i,:)+ori_1st(i,:)*baseline(i)*100];
    grad.coilori(i+sum(isgrad),:) = -ori_1st(i,:);
  elseif (info(i) == handles.PlannerGradioMeter || info(i) == handles.RefferencePlannerGradioMeter)
    grad.coilpos(i+sum(isgrad),:) = [grad.coilpos(i,:)+ori_1st_to_2nd(i,:)*baseline(i)*100];
    grad.coilori(i+sum(isgrad),:) = -ori_1st(i,:);
  else
    grad.coilpos(i+sum(isgrad),:) = [0 0 0];
    grad.coilori(i+sum(isgrad),:) = [0 0 0];  
  end
end

% Define the pair of 1st and 2nd coils for each gradiometer
grad.tra = repmat(diag(ones(1,size(grad.coilpos,1)/2),0),1,2);

% for mangetometers change tra as there is no second coil
if any(ismag)
    sz_pnt = size(grad.coilpos,1)/2;
    % create logical variable
    not_2nd_coil = ([diag(zeros(sz_pnt),0)' ismag']~=0);
    grad.tra(ismag,not_2nd_coil) = 0;
end

% the gradiometer labels should be consistent with the channel labels in
% read_yokogawa_header, the predefined list of channel names in ft_senslabel
% and with ft_channelselection
% ONLY consistent with read_yokogawa_header as NO FIXED relation between
% channel index and type of channel exists for Yokogawa systems. Therefore
% all have individual label sequences: No useful support in ft_senslabel possible
if ~isempty(label)
    grad.label = label(isgrad);
else
    % this is only backup, if something goes wrong above.
    label = cell(size(isgrad));
    for i=1:length(label)
    label{i} = sprintf('AG%03d', i);
    end
    grad.label = label(isgrad);    
end
grad.unit = 'cm';


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this defines some usefull constants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function handles = definehandles
handles.output = [];
handles.sqd_load_flag = false;
handles.mri_load_flag = false;
handles.NullChannel         = 0;
handles.MagnetoMeter        = 1;
handles.AxialGradioMeter    = 2;
handles.PlannerGradioMeter  = 3;
handles.RefferenceChannelMark = hex2dec('0100');
handles.RefferenceMagnetoMeter       = bitor( handles.RefferenceChannelMark, handles.MagnetoMeter );
handles.RefferenceAxialGradioMeter   = bitor( handles.RefferenceChannelMark, handles.AxialGradioMeter );
handles.RefferencePlannerGradioMeter = bitor( handles.RefferenceChannelMark, handles.PlannerGradioMeter );
handles.TriggerChannel      = -1;
handles.EegChannel          = -2;
handles.EcgChannel          = -3;
handles.EtcChannel          = -4;
handles.NonMegChannelNameLength = 32;
handles.DefaultMagnetometerSize       = (4.0/1000.0);       % ????4.0mm???????`
handles.DefaultAxialGradioMeterSize   = (15.5/1000.0);      % ???a15.5mm???~??
handles.DefaultPlannerGradioMeterSize = (12.0/1000.0);      % ????12.0mm???????`
handles.AcqTypeContinuousRaw = 1;
handles.AcqTypeEvokedAve     = 2;
handles.AcqTypeEvokedRaw     = 3;
handles.sqd = [];
handles.sqd.selected_start  = [];
handles.sqd.selected_end    = [];
handles.sqd.axialgradiometer_ch_no      = [];
handles.sqd.axialgradiometer_ch_info    = [];
handles.sqd.axialgradiometer_data       = [];
handles.sqd.plannergradiometer_ch_no    = [];
handles.sqd.plannergradiometer_ch_info  = [];
handles.sqd.plannergradiometer_data     = [];
handles.sqd.nullchannel_ch_no   = [];
handles.sqd.nullchannel_data    = [];
handles.sqd.selected_time       = [];
handles.sqd.sample_rate         = [];
handles.sqd.sample_count        = [];
handles.sqd.pretrigger_length   = [];
handles.sqd.matching_info   = [];
handles.sqd.source_info     = [];
handles.sqd.mri_info        = [];
handles.mri                 = [];

