function [simulated] = dipolesimulation(cfg)

% DIPOLESIMULATION computes the field or potential of a simulated dipole
% and returns a datastructure identical to the PREPROCESSING function.
%
% Use as
%   data = dipolesimulation(cfg)
%
% Where the configuration should contain information about the sensors
% and the forward model
%   cfg.channel    = cell-array (optional, see CHANNELSELECTION)
%   cfg.hdmfile    = string (required)
%   cfg.gradfile   = string (optional, can already be present in data)
%   cfg.elecfile   = string (optional, can already be present in data)
% or alternatively
%   cfg.vol        = structure with volume conduction model
%   cfg.grad       = structure with gradiometer definition
%   cfg.elec       = structure with electrode definition
%
% The dipoles position and orientation have to be specified with
%   cfg.dip.pos     = [Rx Ry Rz]
%   cfg.dip.mom     = [Qx Qy Qz]
%
% The timecourse of the dipole activity is given as a single vector or as a
% cell-array with one vectors per trial
%   cfg.dip.signal
% or by specifying a sine-wave signal
%   cfg.dip.frequency    in Hz
%   cfg.dip.phase        in radians
%   cfg.ntrials          number of trials
%   cfg.triallength      time in seconds
%   cfg.fsample          sampling frequency in Hz
%
% Random white noise can be added to the data in each trial, either by 
% specifying an absolute or a relative noise level
%   cfg.relnoise    = add noise with level relative to simulated signal
%   cfg.absnoise    = add noise with absolute level

% Copyright (C) 2004, Robert Oostenveld
%
% $Log: dipolesimulation.m,v $
% Revision 1.12  2005/11/10 10:02:14  roboos
% minor change in whitespace and help
%

% set the defaults
if ~isfield(cfg, 'dip'),        cfg.dip = [];             end
if ~isfield(cfg.dip, 'pos'),    cfg.dip.pos = [-5 0 15];  end
if ~isfield(cfg.dip, 'mom'),    cfg.dip.mom = [1 0 0]';   end
if ~isfield(cfg, 'fsample'),    cfg.fsample = 250;        end
if ~isfield(cfg, 'relnoise'),   cfg.relnoise = 0;         end
if ~isfield(cfg, 'absnoise'),   cfg.absnoise = 0;         end
if ~isfield(cfg, 'feedback'),   cfg.feedback = 'text';    end

Ndipoles = size(cfg.dip.pos,1);
if size(cfg.dip.pos,2)~=3
  error('dipole positions should be specified as [x, y, z]');
end

% prepare the volume conductor and the sensor array
[vol, sens, cfg] = prepare_vol_sens(cfg, []);

if ~isfield(cfg, 'ntrials') 
  if isfield(cfg.dip, 'signal')
    cfg.ntrials = length(cfg.dip.signal);
  else
    cfg.ntrials = 20;
  end
end
Ntrials  = cfg.ntrials;

if all(size(cfg.dip.mom)==[1 3])
  cfg.dip.mom = cfg.dip.mom';
end

try
  % this should be a column vector
  cfg.dip.frequency = cfg.dip.frequency(:);
end

try
  % this should be a column vector
  cfg.dip.phase = cfg.dip.phase(:);
end

% no signal was given, compute a cosine-wave signal as timcourse for the dipole
if ~isfield(cfg.dip, 'signal')
  % set some additional defaults if neccessary
  if ~isfield(cfg.dip, 'frequency')
    cfg.dip.frequency = ones(Ndipoles,1)*10;
  end
  if ~isfield(cfg.dip, 'phase')
    cfg.dip.phase = zeros(Ndipoles,1);
  end
  if ~isfield(cfg, 'triallength')
    cfg.triallength = 1;
  end
  % compute a cosine-wave signal wit the desired phase and frequency
  nsamples = round(cfg.triallength*cfg.fsample);
  time = (0:(nsamples-1))/cfg.fsample;
  cfg.dip.signal = cos(cfg.dip.frequency*time*2*pi + repmat(cfg.dip.phase,1,nsamples));
end

% construct the timecourse of the dipole activity for each individual trial
if ~iscell(cfg.dip.signal)
  dipsignal = {};
  time      = {};
  nsamples  = length(cfg.dip.signal);
  for trial=1:Ntrials
    % each trial has the same dipole signal 
    dipsignal{trial} = cfg.dip.signal;
    time{trial} = (0:(nsamples-1))/cfg.fsample;
  end
else
  dipsignal = {};
  time      = {};
  for trial=1:Ntrials
    % each trial has a different dipole signal 
    dipsignal{trial} = cfg.dip.signal{trial};
    time{trial} = (0:(length(dipsignal{trial})-1))/cfg.fsample;
  end
end

dippos    = cfg.dip.pos;
dipmom    = cfg.dip.mom;

if ~iscell(dipmom)
  dipmom = {dipmom};
end

if ~iscell(dippos)
  dippos = {dippos};
end

if length(dippos)==1
  dippos = repmat(dippos, 1, Ntrials);
elseif length(dippos)~=Ntrials
  error('incorrect number of trials specified in the dipole position');
end

if length(dipmom)==1
  dipmom = repmat(dipmom, 1, Ntrials);
elseif length(dipmom)~=Ntrials
  error('incorrect number of trials specified in the dipole moment');
end

simulated.trial  = {};
simulated.time   = {};
simulated.offset = [];
progress('init', cfg.feedback, 'computing simulated data');
for trial=1:Ntrials
  progress(trial/Ntrials, 'computing simulated data for trial %d\n', trial);
  lf = compute_leadfield(dippos{trial}, sens, vol);
  simulated.trial{trial}  = lf * dipmom{trial} * dipsignal{trial};
  simulated.time{trial}   = time{trial};
  simulated.offset(trial) = 0;
end
progress('close');

if strcmp(sens.type, 'gradiometers')
  simulated.grad = sens;
elseif strcmp(sens.type, 'electrodes')
  simulated.elec = sens;
end

% determine RMS value of simulated data
ss = 0;
sc = 0;
for trial=1:Ntrials
  ss = ss + sum(simulated.trial{trial}(:).^2);
  sc = sc + length(simulated.trial{trial}(:));
end
rms = sqrt(ss/sc);
fprintf('RMS value of simulated data is %g\n', rms);

% add noise to the simulated data
for trial=1:Ntrials
  relnoise = randn(size(simulated.trial{trial})) * cfg.relnoise * rms;
  absnoise = randn(size(simulated.trial{trial})) * cfg.absnoise;
  simulated.trial{trial} = simulated.trial{trial} + relnoise + absnoise;
end

simulated.fsample = cfg.fsample;
simulated.label   = sens.label;

% add version details to the configuration
try
  % get the full name of the function
  cfg.version.name = mfilename('fullpath');
catch
  % required for compatibility with Matlab versions prior to release 13 (6.5)
  [st, i] = dbstack;
  cfg.version.name = st(i);
end
cfg.version.id   = '$Id: dipolesimulation.m,v 1.12 2005/11/10 10:02:14 roboos Exp $';
% remember the configuration details of the input data
try, cfg.previous = data.cfg; end
% remember the exact configuration details in the output 
simulated.cfg = cfg;

