clear all 
close all
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% load electrodes positions
elec = [];
load electrodes1020
% positions of 33 electrodes
i_elec=[4 6 18 19 21 23 25 27 28 31 33 35 37 40 41 43 45 47 49 50 53 55 57 59 62 63 65 67 69 71 72 84 86];
% load concentric spheres or bem 
vol = [];

% the values of this vectors correspond
% to the number of spheres
vol.r = [88 92 100]; % radii of spheres
vol.c = [1 1/80 1];       % conductivity

vol.o = [0 -2 0];          % center of sphere
% Boundary Element Methods
% Real
% load standard_vol.mat
% Spheric
% vol.cond = [1 1/80 1];
% load spheric_vol_162.mat
% load spheric_vol_642.mat
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
cfg      = [];
cfg.vol  = vol;            
% cfg.elec = elec;
cfg.elec.pnt = elec.pnt(i_elec,1:3);
cfg.elec.label = elec.label(i_elec);
% initialization simulation parameters
fs=250;
ns=250;
cfg.fsample = fs;            
time = (1:ns)/fs;           
signal1 = 1*sin(2*time*6*pi);
signal2 = 1*cos(2*time*12*pi);
signal3 = 1*cos(2*time*8*pi);
% head model
[vols, sens, cfg] = prepare_headmodel(cfg, []);
% source positions
dipin=[];
[x,y,z]=meshgrid(linspace(-40,40,4), linspace(-40,40,4), linspace(0,30,3) );
dipin.pos=[x(:) y(:) z(:)];
mom_t=zeros(3*length(dipin.pos),ns);
% inicialization
x_orig=mom_t;
n_orig=mom_t;
n_orig(20,:)=signal1;
n_orig(80,:)=signal2;
n_orig(40,:)=signal3;
% dynamic model parameters
a11=1.2;%1.2
b11=0.05;%0.05
M11=3*length(dipin.pos);
A = a11*diag(ones(M11,1)) + b11*(diag(ones(M11-1,1),1) + diag(ones(M11-1,1),-1));
% b must be zero for a linear model
b12=0;%0.01
b13=0;%-0.001
c11=-0.9;%-0.9
for b=1:ns,
    if b==1,
        x_orig(:,b)=n_orig(:,b);
    end
    if b==2,
        x_orig(:,b)=A*x_orig(:,b-1)+b12*x_orig(:,b-1).^2-b13*x_orig(:,b-1).^3+n_orig(:,b);
    end
    if b>2,
        x_orig(:,b)=A*x_orig(:,b-1)+b12*x_orig(:,b-1).^2+b13*x_orig(:,b-1).^3+c11*x_orig(:,b-2)+n_orig(:,b);
    end
end
% leadfield computing
C=compute_leadfield(dipin.pos, sens, vols);
% EEG simulation
dat=C*x_orig;
figure(1)
plot(dat')