clc
clear
tic
%brainstorm %norgui
%%
WGN = 1; % Using white Gaussian Nosie (WGN = 1) or Human EEG noise (WGN  = 0);
LFNormlization = 0; % Whether normalize the LeadField Matrix
Uniform = 1; % Uniform/NonUniform Sources
% tic
%% Export the Channel,Cortex and HeadModel to workspace
if WGN
    channelselect=[1:32,34:42,44:64]; %Select EEG data
%     channelselect=[1:32 34:42 44:59 61:63];
else
    channelselect=[1:32 34:42 44:59 61:63]; % for Real noise simulations
end
% [sStudy, iStudy] = bst_get('StudyWithCondition', 'LiDaoli/Stim_2');

[sStudy, iStudy] = bst_get('StudyWithCondition','Subject01/Simulation');
% [sStudy, iStudy] = bst_get('StudyWithCondition','Subject01/Simulation');%GaussianSources');%Extents');%SNRs');
% [sStudy, iStudy] =bst_get('StudyWithCondition','Subject01/mind004_050924_median01_raw_clean_notch');

index = 1;
bst_call(@export_matlab, {char(sStudy.Data(1,index).FileName)},'data');

%=======================I021mport the LFM and Cortex==================================%
[sSurface, iSurface] = bst_get('SurfaceFileByType',[],'Cortex');
% [sSurface, iSurface] = bst_get('SurfaceFileByType',2,'Cortex');
bst_call(@export_matlab, {char(sSurface.FileName)},'Cortex');
Atlas = Cortex.Atlas(2);

[sHeadModel] = bst_get('HeadModelForStudy', iStudy);
bst_call(@export_matlab, {char(sHeadModel.FileName)},'model');
Gain=model.Gain(channelselect,:);
GridLoc=model.GridLoc;
GridOrient=model.GridOrient;

Gain = bst_gain_orient(Gain,GridOrient);
% load('MEGLeadfield.mat')
clear GridOrient
% load 'MEGLeadfield.mat';
% load 'CortexMEG.mat';
[nSensor,nSource] = size(Gain);

%% Output path
outpath = 'D:\Data\1';


%% Iteration
Miter = 1000;
Eccentricity = sqrt(sum(GridLoc.^2,2));
Ec = find(Eccentricity > 70*1e-3);

se = 1;

for iter = 1 : Miter
    ind = randperm(numel(Ec));
    for iteration = 1 : 10
    seedvox = Ec(ind(iteration));
    tau = [0.1 0.35 0.5 0.6];omega = [0.1 0.15 0.15 0.15];%[0.07 0.05 0.10 0.10];%[0.035 0.035 0.035 0.035];%*max(Time(Activetime));
    %tau = [0.1 0.2 0.5 0.6];omega = [0.1 0.15 0.15 0.15];%[0.07 0.05 0.10 0.10];%[0.035 0.035 0.035 0.035];%*max(Time(Activetime));
    f = [10 11 8 9];%10*ones(1,4);%5 (Hz);
    Amp = 1e-8;
    StimTime = find(abs(data.Time) == min(abs(data.Time)));
    TimeLen = 100;
    Time = data.Time(StimTime-0.5*TimeLen:StimTime+0.5*TimeLen-1);

    OPTIONS.DefinedArea    = 18*1e-4;
    OPTIONS.seedvox        = seedvox;
    OPTIONS.frequency      = f;
    OPTIONS.tau            = tau;
    OPTIONS.omega          = omega;
    OPTIONS.Amp            = Amp;
    OPTIONS.GridLoc        = GridLoc;
    OPTIONS.uniform       = Uniform;
    OPTIONS.WGN           = WGN;
    OPTIONS.SNR           = 5;
    OPTIONS.ar            = 0;
    OPTIONS.params(:,:,1) = [ 0.8    0    0 ;
                                0  0.9  0.5 ;
                              0.4    0  0.5];

    OPTIONS.params(:,:,2) = [-0.5    0    0 ;
                                0 -0.8    0 ;
                                0    0 -0.2];

    OPTIONS.noisecov      = [ 0.3    0    0 ;
                                0    1    0 ;
                                0    0  0.2];
    OPTIONS.SinglePoint   = 0;

    [B,s_real,Result] = Simulation_Data_Generate(LFM,Cortex,Time,OPTIONS);
    ActiveVoxSeed = Result.ActiveVoxSeed;
    % %=======================================================================%
    fprintf('Actual SNR is %g\n',20*log10(norm(Gain*s_real,'fro')/norm(B-Gain*s_real,'fro')));
    
    data_name = sprintf('data_%d.mat', se);
    sreal_name = sprintf('sreal_%d.mat', se);
    
    save([outpath '\' data_name], 'B', 's_real')
    
    se = se + 1;
    end
end

save([outpath '\' 'GridLoc.mat'], 'GridLoc')
save([outpath '\' 'Cortex.mat'], 'Cortex')
save([outpath '\' 'L.mat'], 'Gain')
    
    
    
    
    