function [Electron_MULTILAYER] = MC_multilayer(Thicknesses,varargin)
l = length(Thicknesses);
Thicknesses = Thicknesses.*10;
if nargin == 1 || nargin==2 || l~=(nargin-1) || nargin>5
    error('Wrong input! 1. Vector of thicknesses, 2. - 5. Structure with layer simulated')
else
    switch nargin
        case 3
            Electron(1).Layer = varargin{1};
            Electron(2).Layer = varargin{2};
            [m(1),n(1)] = size(Electron(1).Layer.X);
            [m(2),n(2)] = size(Electron(2).Layer.X);
        case 4
            Electron(1).Layer = varargin{1};
            Electron(2).Layer = varargin{2};
            Electron(3).Layer = varargin{3};
            [m(1),n(1)] = size(Electron(1).Layer.X);
            [m(2),n(2)] = size(Electron(2).Layer.X);
            [m(3),n(3)] = size(Electron(3).Layer.X);
        case 5
            Electron(1).Layer = varargin{1};
            Electron(2).Layer = varargin{2};
            Electron(3).Layer = varargin{3};
            Electron(4).Layer = varargin{4};
            [m(1),n(1)] = size(Electron(1).Layer.X);
            [m(2),n(2)] = size(Electron(2).Layer.X);
            [m(3),n(3)] = size(Electron(3).Layer.X);
            [m(4),n(4)] = size(Electron(4).Layer.X);
        case 6
            Electron(1).Layer = varargin{1};
            Electron(2).Layer = varargin{2};
            Electron(3).Layer = varargin{3};
            Electron(4).Layer = varargin{4};
            Electron(5).Layer = varargin{5};
            [m(1),n(1)] = size(Electron(1).Layer.X);
            [m(2),n(2)] = size(Electron(2).Layer.X);
            [m(3),n(3)] = size(Electron(3).Layer.X);
            [m(4),n(4)] = size(Electron(4).Layer.X);
            [m(5),n(5)] = size(Electron(5).Layer.X);
    end
end
m = min(m);
ncols = sum(n);

Electron_MULTILAYER.X  = zeros(m,ncols);
Electron_MULTILAYER.Y  = zeros(m,ncols);
Electron_MULTILAYER.Z  = zeros(m,ncols);
Electron_MULTILAYER.Th = zeros(m,ncols);
Electron_MULTILAYER.Ps = zeros(m,ncols);
Electron_MULTILAYER.dE = zeros(m,ncols);
Electron_MULTILAYER.SE = zeros(m,ncols);
Electron_MULTILAYER.PX = zeros(m,1);
Electron_MULTILAYER.PY = zeros(m,1);
Electron_MULTILAYER.PZ = zeros(m,1);
[m1,n] = size(Electron(1).Layer(1).X);
Electron_MULTILAYER.X(1:m,1:n)  = Electron_MULTILAYER.X(1:m,1:n) + Electron(1).Layer(1).X;
Electron_MULTILAYER.Y(1:m,1:n)  = Electron_MULTILAYER.Y(1:m,1:n) + Electron(1).Layer(1).Y;
Electron_MULTILAYER.Z(1:m,1:n)  = Electron_MULTILAYER.Z(1:m,1:n) + Electron(1).Layer(1).Z;
Electron_MULTILAYER.Th(1:m,1:n) = Electron_MULTILAYER.Th(1:m,1:n) + Electron(1).Layer(1).Th;
Electron_MULTILAYER.Ps(1:m,1:n) = Electron_MULTILAYER.Ps(1:m,1:n) + Electron(1).Layer(1).Ps;
Electron_MULTILAYER.dE(1:m,1:n) = Electron_MULTILAYER.dE(1:m,1:n) + Electron(1).Layer(1).dE;
Electron_MULTILAYER.SE(1:m,1:n) = Electron_MULTILAYER.SE(1:m,1:n) + Electron(1).Layer(1).SE;

h = waitbar(0, 'Mc MultiLayer')
L = length(Electron);
for i = 1:L-1
    T = cumsum(Thicknesses(1:i));
    
    Zi  = Electron_MULTILAYER.Z;
    Zii = Electron_MULTILAYER.Z;
    Xi  = Electron_MULTILAYER.X;
    Yi  = Electron_MULTILAYER.Y;
    Thi = Electron_MULTILAYER.Th;
    dEi = Electron_MULTILAYER.dE;
    
   
    Zi(Zi>T(i)) = 0;
    Zi(Zi>0)= 1;
    Zi = ~Zi;
    [Z1,IndZ1] = max(Zi,[],2);
    IndZ1 = IndZ1-1;
    inds0 = find(IndZ1 == 0);
    %            IndZ1 = IndZ1+1;
    IndZ1(IndZ1==0) = 1;
    
    [s1,s2] = size(Zi);
    v = (1:s1)';
    inds = (IndZ1-1).*s1 + v;
    X1 = Xi(inds);
    Y1 = Yi(inds);
    Th1 = Thi(inds);
    dE1 = dEi(inds);
    Z1 = Zii(inds);
    IndZ2 = IndZ1 + 1;
    inds2 = (IndZ2-1).*s1 + v;
    
    if ~isempty(inds0)
        X1(inds0) = 0;
        Y1(inds0) = 0;
        Th1(inds0) = 0;
        dE1(inds0) = 0;
        Z1(inds0) = 0;
        IndZ1(inds0)=1;
        inds2(inds0) = inds2(inds0)-s1;
    end
    
    %     inds2(Zi(:,1)==0) = inds2(Zi(:,1)==0)-s1;
    X2 = Xi(inds2);
    Y2 = Yi(inds2);
    Th2 = Thi(inds2);
    dE2 = dEi(inds2);
    Z2 = Zii(inds2);
        
    % body svazku pred  P4 a za P5 prechodem materialu
    P4 = [X1, Y1, Z1];
    P5 = [X2, Y2, Z2];
        
    % Body roviny
    P1 = [10 10 T(i)];
    P2 = [20 30 T(i)];
    P3 = [30 20 T(i)];
       
    for j = 1:length(X1)
        waitbar(j/length(X1)*(L-1));
        % Priesecnik
        L1 = P4(j,:);
        L2 = P5(j,:);
        
        if L2(3)<L1(3)  % BSE urcite nepretina T
            % ulozenie
            Electron_MULTILAYER.Z(j,1:IndZ2(j))  = Electron_MULTILAYER.Z(j,1:IndZ2(j)); % tu bolo +PZ
            Electron_MULTILAYER.X(j,1:IndZ2(j))  = Electron_MULTILAYER.X(j,1:IndZ2(j));
            Electron_MULTILAYER.Y(j,1:IndZ2(j))  = Electron_MULTILAYER.Y(j,1:IndZ2(j));
            Electron_MULTILAYER.Th(j,1:IndZ2(j)) = Electron_MULTILAYER.Th(j,1:IndZ2(j));
            Electron_MULTILAYER.Ps(j,1:IndZ2(j)) = Electron_MULTILAYER.Ps(j,1:IndZ2(j));
            Electron_MULTILAYER.dE(j,1:IndZ2(j)) = Electron_MULTILAYER.dE(j,1:IndZ2(j));
            Electron_MULTILAYER.PX(j,i+1) = 0;
            Electron_MULTILAYER.PY(j,i+1) = 0;
            Electron_MULTILAYER.PZ(j,i+1) = 0;
            
        else % front scattered
            
            % vypocet prusecniku:
            N = cross(P2-P1, P3 - P1);
            Intersec = L1 + (dot(N, P1 - L1) / dot(N, L2 - L1)) * (L2 - L1);
            % Body priesecniku
            PX = Intersec(1); PY = Intersec(2); PZ = Intersec(3);
            
            
            [Min,IndMin] = min(Electron(i+1).Layer.Z,[],2);
            IndMin = IndMin - 1;
            % naklon podla thety poslednej
            %                                  THn = Electron(i).Layer.Th(j,IndZ1(j));
            THn = Th1(j);
            if THn ~= 0
                RMx = [cos(THn), -sin(THn); sin(THn),cos(THn)];
                Mx = [Electron(i+1).Layer.X(j,:) ;  Electron(i+1).Layer.Z(j,:)]';
                M_newx = Mx*RMx;
                Electron(i+1).Layer.Z(j,:) = M_newx(:,2)';
                Electron(i+1).Layer.X(j,:) = M_newx(:,1)';
            else
            end
            
            % prepocet suradnic
%             vec = find(Electron_MULTILAYER.Z(j,:)>T(i)); % i+1
            vec = IndZ2(j);
            if isempty(vec)==0
                poz = vec(1);
                if poz==0
                    poz=1;
                end
                %         Electron_MULTILAYER.Z(j,poz:poz+poz2) = Electron(i).Layer.Z(j,IndZ2(j):end);
                %             dEn = Electron(i).Layer.dE(j,IndZ1(j));
                dEn = dE1(j);
                indn0 = find(Electron(i+1).Layer.Z(j,:)~=0);  % len v miestach kolizie
                index0 = indn0(end);
                Electron(i+1).Layer.Z(j,1:index0) = Electron(i+1).Layer.Z(j,1:index0)+PZ;
                Electron(i+1).Layer.X(j,1:index0) = Electron(i+1).Layer.X(j,1:index0)+PX;
                Electron(i+1).Layer.Y(j,1:index0) = Electron(i+1).Layer.Y(j,1:index0)+PY;
                Electron(i+1).Layer.Th(j,1:index0) = Electron(i+1).Layer.Th(j,1:index0)+THn;
                Electron(i+1).Layer.dE(j,1:index0) = Electron(i+1).Layer.dE(j,1:index0)+dEn;
                
                [rows,cols] = size(Electron(i+1).Layer.Z(j,:));
                Electron_MULTILAYER.Z(j,poz:poz+cols-1)  = Electron(i+1).Layer.Z(j,:); % tu bolo +PZ
                Electron_MULTILAYER.X(j,poz:poz+cols-1)  = Electron(i+1).Layer.X(j,:);
                Electron_MULTILAYER.Y(j,poz:poz+cols-1)  = Electron(i+1).Layer.Y(j,:);
                Electron_MULTILAYER.Th(j,poz:poz+cols-1) = Electron(i+1).Layer.Th(j,:);
                Electron_MULTILAYER.Ps(j,poz:poz+cols-1) = Electron(i+1).Layer.Ps(j,:);
                Electron_MULTILAYER.dE(j,poz:poz+cols-1) = Electron(i+1).Layer.dE(j,:);
                Electron_MULTILAYER.PX(j,i+1) = PX;
                Electron_MULTILAYER.PY(j,i+1) = PY;
                Electron_MULTILAYER.PZ(j,i+1) = PZ;
            else
                
            end
            
            
        end
    end
end
close(h) 
%ulozenie
% dt = string(datetime);
% dt(dt==':') = '_';
% NAME = strcat('Electron_MULTILAYER_',dt,'.mat');
% save(NAME,'Electron_MULTILAYER')
% 
% figure
% plot(SHIFT,FrUSE,'b')
% hold on
% plot(SHIFT,FrBSE,'g')
% plot(SHIFT,FrSCE,'r')
% hold off
% legend('Unscattered','BackScattered','FrontScattered')
% title('MultiLayer');