function [hdm,sim_tl,fit] = CzSimulation(cfg)
% CZSIMULATION Creates the extended segmentation, then computes
% the headmodel and simulates some dipoles, all data will be saved to given
% folder 
%
%Call:
%   [hdm,sim_tl,fit]    = CzSimulation(cfg) Create extended segmentation,
%                                           generate BEM-Model simulate
%                                           data and apply fit
%   [hdm,sim_tl]        = CzSimulation(cfg) Create extended segmentation,
%                                           generate BEM-Model amd simulate
%                                           data
%   [hdm]               = CzSimulation(cfg) Create extended segmentation
%                                           and headmodel only
%                         CzSimulation(cfg) Same as first, save all data to
%                                           file
%
%Returns:
%   hdm         -   BEM headmodel
%   sim_tl      -   Simulated data after applying timelock analysis
%   fit         -   Fitted data
%
%Parameters:
%   cfg         -   Configuration structure with fields:
%       folder      -   Folder containing the segmented structures
%       file        -   File suffix
%           or
%       seg         -   Structure with segmentation from SPM8
%       mri         -   MRI data structure


if(nargout == 0)
    disp([cfg.file,': Keine Ausgabeparameter, speichere Ausgabe']);
end

if(isfield(cfg,'folder') && isfield(cfg, 'file'))
    seg_cfg.segmented = true;
    seg_cfg.file = cfg.file;
    seg_cfg.folder = cfg.folder;
    if(isfield(cfg, 'anatomy_threshold')); seg_cfg.anatomy_threshold = cfg.anatomy_threshold; end
elseif(isfield(cfg, 'seg') && isfield(cfg, 'mri'))
    seg_cfg.segmented = true;
    seg_cfg.seg = cfg.seg;
    seg_cfg.mri = cfg.mri;
elseif(isfield(cfg, 'ftseg'))
    disp('Segmentation already done. Skipping.');
else
    error('Invalid configuration structure');
end

if(isfield(cfg, 'ftseg') == false)
    seg = CreateExtendedSegmentation(seg_cfg);
else
    seg = cfg.ftseg;
end

bem_cfg.mri = seg;

model = CreateBEMModel(bem_cfg);

try 
    nipa = open([cfg.folder, 'nipa', cfg.file, '.mat']);
    nipa = nipa.nipa;
catch
    disp('Failed to open nasion-inion-pre-auricular data');
end

[s,t] = FiducialSystem(nipa.nasion, nipa.left, nipa.right);
electrodes = ElectrodesToFT(Get1020(seg.seg > 0, nipa.nasion, nipa.inion, nipa.left, nipa.right));

if(nargout > 0)
    hdm = model;
end

if(nargout ~= 1)
    [time,sig] = GetP1N1(1,200,100,10,20,120,10,15);

    sim_cfg.elec = electrodes;
    sim_cfg.ntrials = 5;
    sim_cfg.triallength = 1;
    sim_cfg.dip.signal = sig;
    sim_cfg.absnoise = 3;
    sim_cfg.relnoise = 0.1;
    
    sim_cfg.vol = model.vol;
    
    sim_tl = cell(0);
    
    timelock_cfg.covariance = 'yes';
    timelock_cfg.covariancewindow = [0 200];
    
    timelock_cfg.blcovariance = 'yes';
    timelock_cfg.blcovariancewindow = [0 200];
    
    for i=1:length(cfg.dipoles)
        sim_cfg.dip.pos = s*cfg.dipoles(i).pos'+t;      %convert position
        sim_cfg.dip.mom = cfg.dipoles(i).mom;        
        
        sim_raw = ft_dipolesimulation(sim_cfg);
        sim_preproc = ft_preprocessing([], sim_raw);
        sim_timelock = ft_timelockanalysis(timelock_cfg, sim_preproc);

        sim_tl(i) = {sim_timelock};
        
        if(nargout == 0)
            fileout.sim(i).raw = sim_raw;
            fileout.sim(i).preproc = sim_preproc;
        else
            sim(i) = sim_raw;
            
        end
    end
end

if(nargout==0 || nargout==3)
    fit_cfg.numdipoles = 1;
    fit_cfg.vol = model.vol;
    fit_cfg.elec = electrodes;
    fit_cfg.grid.resolution = 1;
    
    for i=1:length(sim_tl)
        dipole_fit = ft_dipolefitting(fit_cfg, sim_tl{i});
        
        if(nargout == 0)
            fileout.fit(i) = dipole_fit;
        else
            fit(i) = dipole_fit;
        end
    end
end


if(nargout == 0)
    fileout.hdm = model;
    fileout.seg = seg;
    
    filename = [cfg.folder, 'msf', cfg.file, '.mat'];
    
    save(filename, 'fileout');
end

end