%% clear
clear all global
close all

%% load data....
data_folder = 'beamformer';
load(fullfile(data_folder, 'dataFIC.mat'));
load(fullfile(data_folder, 'vol.mat'));
mri = ft_read_mri(fullfile(data_folder, 'Subject01.mri'));

sens = dataFIC.grad;

%% make grid...
cfg = [];
cfg.grid.resolution = 10;
cfg.headmodel = vol;

grid = ft_prepare_sourcemodel(cfg);

%% set vars....
fsample = 200;
ntrials = 100;
dip = 850;
freq = 10;

%% simulate a dipole.....
cfg = [];
cfg.grad = sens;
cfg.headmodel = vol;
cfg.channel = 'MEG';
cfg.dip.pos = grid.pos(min(find(cumsum(grid.inside) == dip)), :);
cfg.dip.mom = [1 0 0]';
cfg.dip.frequency = freq;
cfg.dip.phase = 0;
cfg.dip.amplitude = 1;
cfg.ntrials = ntrials;
cfg.triallength = 1;
cfg.fsample = fsample;
cfg.absnoise = .1;

data = ft_dipolesimulation(cfg);

cfg.dip.amplitude = 0;

data_bl = ft_dipolesimulation(cfg);

%% make leadfield...
cfg = [];
cfg.grid = grid;
cfg.grad = sens;
cfg.vol = vol;

lf = ft_prepare_leadfield(cfg, data);

%% do mtmfft for beamformers....
cfg = [];
cfg.method = 'mtmfft';
cfg.output = 'fourier';
cfg.foi = [freq-1 freq+1];
cfg.taper = 'hanning';

data_fft = ft_freqanalysis(cfg, data);
data_fft_bl = ft_freqanalysis(cfg, data_bl);

%% do source analysis....
cfg = [];
cfg.method = 'dics';
cfg.grid = lf;
cfg.frequency = freq;
cfg.headmodel = vol;
cfg.dics.projectnoise = 'yes';
cfg.dics.lambda = '5%';

data_dics = ft_sourceanalysis(cfg, data_fft);
data_dics_bl = ft_sourceanalysis(cfg, data_fft_bl);

data_dics_norm_bl = data_dics;
data_dics_norm_bl.avg.pow = data_dics.avg.pow ./ data_dics_bl.avg.pow;

data_dics_norm_nai = data_dics;
data_dics_norm_nai.avg.pow = data_dics.avg.pow ./ data_dics.avg.noise;


dics_desc = ft_sourcedescriptives([], data_dics);
%% plot...
toplot = 'pow';

cfg = [];
cfg.parameter = toplot;

source = ft_sourceinterpolate(cfg, data_dics_norm_bl, mri);

cfg = [];
cfg.funparameter = toplot;
cfg.location = grid.pos(min(find(cumsum(grid.inside) == dip)), :);

ft_sourceplot(cfg, source);

%% check distance between max activity in source space and chosen dipole
[~, i] = max(data_dics_norm_bl.avg.(toplot))

rms(grid.pos(i, :) - grid.pos(min(find(cumsum(grid.inside) == dip)), :))

%% try dipolefit....
cfg = [];
cfg.grid = lf;
cfg.headmodel = vol;
cfg.frequency = [freq freq];

dipoles = ft_dipolefitting(cfg, data_fft);