%**********************************
%*********READING dataFIC**********
%**********************************
cfg = [];                                           % empty configuration
cfg.dataset                 = 'MAX_29_es.fif';       % name of CTF dataset  
cfg.trialdef.eventtype      = 'STI101';
cfg.trialdef.prestim        = 1;
cfg.trialdef.poststim       = 2;
cfg.trialdef.eventvalue     = 1;                    % event value of FIC
cfg = ft_definetrial(cfg);   
cfg.channel   = {'MEG','-MEG2442','-MEG0143','-MEG2543'};   % read all MEG channels except...
cfg.baseline = [-1 0];
dataTransparent = ft_preprocessing(cfg);
save dataTransparent dataTransparent


%**********************************
%*********WINDOW OF INTEREST*******
%**********************************

cfg = [];                                           
cfg.toilim = [-0.5 0];                       
dataPre = ft_redefinetrial(cfg, dataTransparent); 
cfg.toilim = [0.8 1.3];                       
dataPost = ft_redefinetrial(cfg, dataTransparent);


%**********************************
%********CROSS SPECTRAL MATRIX*****
%**********************************
cfg = [];
cfg.method    = 'mtmfft';
cfg.output    = 'powandcsd';
cfg.tapsmofrq = 4;
cfg.foilim    = [18 18];
freqPre = ft_freqanalysis(cfg, dataPre);
cfg = [];
cfg.method    = 'mtmfft';
cfg.output    = 'powandcsd';
cfg.tapsmofrq = 4;
cfg.foilim    = [18 18];
freqPost = ft_freqanalysis(cfg, dataPost);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%******************************************
%*****FORDWARD MODEL, LEADFIELD MATRIX*****
%******************************************
 mri = ft_read_mri('C:\Users\ilavado\Desktop\31_CW1AFMRI_ESG_ENDIKA_SANCHEZ_GOMEZ\1_028_t1_mpr_tra_p2_iso_20110314\31_CW1AFMRI_ESG_ENDIKA_SANCHEZ_GOMEZ_20110314_001_028_t1_mpr_tra_p2_iso.img');
%******* SHELL MODEL *******
cfg                = [];
cfg.coordsys       = 'spm'; %We are working with img data (SPM), axes seem to be ok
[segmentedmri]     = ft_volumesegment(cfg, mri)

%The segmentation fits the mri images perfectly. 
cfg                = [];
test               = segmentedmri;
test.avg.pow       = test.gray+test.white;
test.anatomy       = mri_real.anatomy;
cfg.funparameter   = 'avg.pow';
cfg.interactive    = 'yes';
ft_sourceplot(cfg,test);

%With the code that you include (volumerealign), the segmentation does not
%fit and the vol can not be calculated properly.

%cfg = []; 
%cfg.coordsys       = 'neuromag';
% cfg.interactive    = 'yes';
% mri_real           = ft_volumerealign(cfg,mri);
% [segmentedmri]     = ft_volumesegment(cfg, mri_real)

%************** CREATE THE HEAD MODEL ***************
cfg = [];
cfg.sourceunits='mm';
vol = ft_prepare_singleshell(cfg, segmentedmri);
ft_plot_vol(vol);  %At this point, the segmentation is well done, we have a 3D brain.

%************************** MAKE THE LEADFIELD *************

vol=ft_convert_units(vol,'cm')
cfg                 = [];
cfg.grad            = freqPre.grad;  %The sensor position of the fif file, from neuromag
cfg.vol             = vol;
cfg.reducerank      = 2;
cfg.grid.resolution = 1;   
[grid] = ft_prepare_leadfield(cfg);

%figure of the headmodel, sensor positions and grid positions
cfg                = [];
cfg.vol            = vol;
cfg.grid           = grid;
cfg.grad     =  freqPre.grad;
figure;
ft_headmodelplot(cfg);  %Here is the error, the brain does not match the coordinates of the helmet. It should be 
%higher and more frontal.