function [Electron_NEW_BS] = MC_sphere(Electron, TH_max, Probe_diam, TH_calc,Probe_calc,A,focus,R,dR,MultiCapsle)
% TH_max - apertura (max. dopadajici uhel theta)
% Probe_diam - spot size (prumer svazku)
% TH_calc - sposob vypoctu nah. cisla pre vyber dopadajuceho uhlu theta
% Probe_calc - sposob vypoctu nah. cisla pre dopadajuci uhol theta
% A - rozhonuti generujem nnove naklony a posuny alebo pri pohybe po kouli
%     nechavam tie sam
% dR = round(0.1*R);
Th1d = 0.034; % 34 mrad
Th2d = 0.039; % 39 mrad
Th3d = 0.054; % 54 mrad

R = R*10;
dR = dR*10;
focus = focus*10;

if dR > 0
    %    dR = round(Probe_diam/2);
    %dR = 0;   %bez posunu
elseif dR == 0
    %bez posunu
else
    errordlg('dR must be >0');
end


% dR_vec = 0:dR:8*R;
Electron_NEW = Electron;
[m,n] = size(Electron.Z);



if MultiCapsle == 1
    h=waitbar(0,'MC MultiCapsle preparation ...');
elseif MultiCapsle == 0
    h=waitbar(0,'MC Sphere ...');
elseif MultiCapsle > 2
    errordlg('MultiCapsle calculation must be 0/1  (on-1 / off-0) ')
end



for i = 1:m
    if MultiCapsle == 1 || MultiCapsle == 0
    waitbar(i/m);
    end
    %% Vypocet novych suradnic podla naklonu a posunu v ramci svazku
    ind = find(Electron.Z(i,:)~=0);  % len v miestach kolizie
    index = ind(end);
    beam_shift = 1;   %on/off
    j = 0;
    while beam_shift
        j = j+1;  % pocitadlo
        dR_vec(j) = -dR + j*dR ;
        Electron_NEW_BS(j).SHIFT = dR_vec(j);
        
        % ZASTAVENIE POSUNU KU HRANE
        if dR == 0  % ak nechceme posuvat beamom
            beam_shift = 0;
        elseif dR_vec(j) > R
            beam_shift = 0;
        end
        
        %GENERUJ NAKLON A POSUN
        % ak je A == 0 generuje sa len prvy krat, ak je rovno 1 generuje sa
        % stale iny naklon a posun
        
        if A == 0 && j == 1  % bolo i == 1
            [TH_nx,TH_ny,XE,YE] = generate(TH_max, Probe_diam, TH_calc,Probe_calc);
        elseif A == 1
            [TH_nx,TH_ny,XE,YE] = generate(TH_max, Probe_diam, TH_calc,Probe_calc);
        elseif A ~= 0 && A~=1
            error('A must be equal 0 or 1 (0- generates TH_n,XE,YE only once;  1 - generates in every step)')
        end
        
        % otocenie X,Z rotacnou maticou RMx
        RMx = [cos(TH_nx), -sin(TH_nx); sin(TH_nx),cos(TH_nx)];
        Mx = [Electron.X(i,:);Electron.Z(i,:)]';
        M_newx = Mx*RMx;
        
        Electron_NEW.Z(i,:) = M_newx(:,2)';
        Electron_NEW.X(i,:) = M_newx(:,1)';
        
        % otocenie Y,Z rotacnou maticou RMy
        RMy = [cos(TH_ny), -sin(TH_ny); sin(TH_ny),cos(TH_ny)];
        My = [Electron.Y(i,:); Electron.Z(i,:)]';
        M_newy = My*RMy;
        
        Electron_NEW.Z(i,:) = M_newy(:,2)';
        Electron_NEW.Y(i,:) = M_newy(:,1)';
        
        % posun zvazku kvoli nenulovemu focusu
        
        focus_shift_x = tan(TH_nx)*focus;
        focus_shift_y = tan(TH_ny)*focus;
        
        
        %% Vypocet prienikovych bodov (P0 a P1) paprsku a gule
        % VYpocet prienikovych bodov pre 3 typy fokusu
        %P0:
        if focus==R
            P0X = focus_shift_x + XE + dR_vec(j);
            P0Y = focus_shift_y + YE;
            P0Z = R;
            
            P1X = -focus_shift_x + XE + dR_vec(j);   % -P0X
            P1Y = -focus_shift_y + YE;             %- P0Y
            P1Z = -R;
            
        elseif focus == 2*R
            P0X = focus_shift_x + XE + dR_vec(j);
            P0Y = focus_shift_y + YE;
            P0Z = R;
            
            P1X = 0 + XE + dR_vec(j);
            P1Y = 0 + YE;
            P1Z = -R;
        elseif focus == 0
            P0X = 0 + XE + dR_vec(j);
            P0Y = 0 + YE;
            P0Z = R;
            
%             P1X = tan(TH_nx)*2*R + XE + dR_vec(j);
%             P1Y = tan(TH_ny)*2*R + YE;
%             P1Z = -R;
            P1X = tan(TH_nx)*2*R + XE + dR_vec(j);
            P1Y = tan(TH_ny)*2*R + YE;
            P1Z = -R;
        else
            error('Focus must be equal 0(top of sphere) or R(middle) or 2R(bottom)');
        end
        
        % priesecniky
        P0 = [P0X,P0Y,P0Z];
        P1 = [P1X,P1Y,P1Z];
        
        S = [0,0,0];      % sphere center
        
        %Vypocet priesecnikov:
        dx = P1(1)-P0(1); dy = P1(2)-P0(2);  dz = P1(3)-P0(3);
        a = dx^2 + dy^2 + dz^2;
        b = 2*dx*(P0(1)-S(1)) + 2*dy*(P0(2)-S(2)) + 2*dz*(P0(3)-S(3));
        c = S(1)^2 + S(2)^2 + S(3)^2 + P0(1)^2 + P0(2)^2 + P0(3)^2 - 2*(S(1)*P0(1) + S(2)*P0(2) + S(3)*P0(3)) - R^2;
        D = b^2 - 4*a*c;
        t1 = (-b -sqrt(D))/(2*a);
        t2 = (-b +sqrt(D))/(2*a);
        
        % SU DVA PRIENIKY ?
        if t1 ~= t2 && t1 ~= 0 && t2 ~= 0 && isreal(t1)~=0 && isreal(t2)~=0 % dva priesecniky
            %% su dva priesecniky 
            % PR1:
            Xn = P0(1)+t1*dx;   Yn = P0(2)+t1*dy;   Zn = P0(3)+t1*dz;
            % PR2:
            Xn2 = P0(1)+t2*dx;   Yn2 = P0(2)+t2*dy;   Zn2 = P0(3)+t2*dz;
            
            % PRIPOCITAJ PRIENIK KU SURADNICIAM
            Electron_NEW_BS(j).X(i,1:index) = Electron_NEW.X(i,1:index) + Xn;
            Electron_NEW_BS(j).Y(i,1:index) = Electron_NEW.Y(i,1:index) + Yn;
            Electron_NEW_BS(j).Zcalc(i,1:index) = Zn - Electron_NEW.Z(i,1:index) ;  % pre overenie polohy bodu -> vo vnutri / vonku
            Electron_NEW_BS(j).Z(i,1:index) = Electron_NEW.Z(i,1:index) + (R-Zn) ;   % ulozenie tloustky
            
            Electron_NEW_BS(j).Th(i,1:index)= Electron_NEW.Th(i,1:index) + TH_nx; % naklon od osi z v smere x
            Electron_NEW_BS(j).Ps(i,1:index)= Electron_NEW.Ps(i,1:index) + TH_ny;
            Electron_NEW_BS(j).Scat = 1;
            Electron_NEW_BS(j).dE(i,1:index) = Electron_NEW.dE(i,1:index);
            Electron_NEW_BS(j).SE(i,1:index) = Electron_NEW.SE(i,1:index);
            
            Electron_NEW_BS(j).Xn(i,1) = Xn; Electron_NEW_BS(j).Xn2(i,1) = Xn2;
            Electron_NEW_BS(j).Yn(i,1) = Yn; Electron_NEW_BS(j).Yn2(i,1) = Yn2;
            Electron_NEW_BS(j).Zn(i,1) = Zn; Electron_NEW_BS(j).Zn2(i,1) = Zn2;
            
            
            % KTORY ROZPTYL JE UZ MIMO GULU?
            Limit = sqrt(Electron_NEW_BS(j).X(i,1:index).^2 + Electron_NEW_BS(j).Y(i,1:index).^2 + Electron_NEW_BS(j).Zcalc(i,1:index).^2);
            Limit(Limit <=R) = 0;
            Limit(Limit > R) = 1;
            [maxind, LimIndex] = max(Limit,[],2);
            
            %% ULOZENIE THETA, vypocet signalu
            if maxind == 0  % BSE alebo absorb
                Electron_NEW_BS(j).THETA(i,1) = Electron_NEW_BS(j).Th(i,index);
            elseif maxind == 1 && LimIndex == 1 % hned prvy rozptyl je mimo gulu
                Electron_NEW_BS(j).THETA(i,1) = TH_nx;
            elseif maxind == 1 && LimIndex > 1 % nastal rozptyl, elektron vychadza vonku az ponom
                Electron_NEW_BS(j).THETA(i,1) = Electron_NEW_BS(j).Th(i,LimIndex);
            else
                disp('nezaradeny rozptyl');
            end
            
            % VYPOCET SIGNALU
            Electron_NEW_BS(j).iSE(i,1) = 1;
            Electron_NEW_BS(j).iTH_X(i,1) = TH_nx;
            
            
            %% TYP ELEKTRONU PODLA JEHO POLOHY
            Electron_NEW_BS(j).IndexDraw(i,1) = LimIndex+1;  % pre vykreslenie
            Electron_NEW_BS(j).Index(i,1) = LimIndex;    % pre vypocty
            %                  if length(LimInd)>=1 && Electron_NEW_BS(j).Zcalc(i,LimIndex)>0
            if LimIndex==1 && Electron_NEW_BS(j).Zcalc(i,LimIndex)>0
                Electron_NEW_BS(j).BSE(i,1) = 1;
                Electron_NEW_BS(j).SCE(i,1) = 0;
                Electron_NEW_BS(j).USE(i,1) = 0;
            elseif LimIndex>1 && Electron_NEW_BS(j).Zcalc(i,LimIndex)> Electron_NEW_BS(j).Zcalc(i,LimIndex-1)
                Electron_NEW_BS(j).BSE(i,1) = 1;
                Electron_NEW_BS(j).SCE(i,1) = 0;
                Electron_NEW_BS(j).USE(i,1) = 0;
            elseif LimIndex>1 && Electron_NEW_BS(j).Zcalc(i,LimIndex)<0
                Electron_NEW_BS(j).SCE(i,1) = 1;
                Electron_NEW_BS(j).USE(i,1) = 0;
                Electron_NEW_BS(j).BSE(i,1) = 0;
                Electron_NEW_BS(j).iTH_X(i,1) = Electron_NEW_BS(j).Th(i,LimIndex); % bolo LimIndex-1
            elseif LimIndex == 1 && Electron_NEW_BS(j).Zcalc(i,LimIndex)<0
                Electron_NEW_BS(j).USE(i,1) = 1;
                Electron_NEW_BS(j).BSE(i,1) = 0;
                Electron_NEW_BS(j).SCE(i,1) = 0;
                
            end
            
        else  % nieje priesecnik s gulou, ukladaju sa vysoke hodnoty kvoli analyze
            %             beam_shift = 0;
            Electron_NEW_BS(j).X(i,1:index) = 3*R; %NaN;
            Electron_NEW_BS(j).Y(i,1:index) = 3*R; %NaN;
            Electron_NEW_BS(j).Zcalc(i,1:index) = -3*R; %NaN ;  % bolo 3*R .. pre overenie polohy bodu -> vo vnutri / vonku
            Electron_NEW_BS(j).Z(i,1:index) = 3*R; %NaN;   % ulozenie tloustky
            
            Electron_NEW_BS(j).Th(i,1:index)= TH_nx; % naklon od osi z v smere x
            Electron_NEW_BS(j).ThY(i,1:index)= TH_ny;
            Electron_NEW_BS(j).Ps(i,1:index)= 0;
            Electron_NEW_BS(j).Scat = 0;
            Electron_NEW_BS(j).dE = Electron_NEW.dE;
            Electron_NEW_BS(j).SE = Electron_NEW.SE;
            
            Electron_NEW_BS(j).Xn(i,1) = NaN; Electron_NEW_BS(j).Xn2(i,1) = NaN;
            Electron_NEW_BS(j).Yn(i,1) = NaN; Electron_NEW_BS(j).Yn2(i,1) = NaN;
            Electron_NEW_BS(j).Zn(i,1) = NaN; Electron_NEW_BS(j).Zn2(i,1) = NaN;
            Electron_NEW_BS(j).Index(i,1) = NaN;
            Electron_NEW_BS(j).IndexDraw(i,1) = NaN;
            Electron_NEW_BS(j).USE(i,1) = 1;
            Electron_NEW_BS(j).BSE(i,1) = 0;
            Electron_NEW_BS(j).SCE(i,1) = 0;
            
            % Ulozenie SIGNALU
            Electron_NEW_BS(j).iSE(i,1) = 1;
            Electron_NEW_BS(j).iTH_X(i,1) = TH_nx;
            Electron_NEW_BS(j).THETA(i,1) = TH_nx;
            
            Electron_NEW_BS(j).BF(i,1) = 1;
            Electron_NEW_BS(j).DF(i,1) = 0;
            Electron_NEW_BS(j).HAADF(i,1) = 0;
        end
        
        
        
    end
end
Electron_NEW_BS(1).R = R;
disp('MC Sphere done');



    
    [m1, m2]=size(Electron_NEW_BS);
    
    for k = 1:m2
        [m3, m4]  =size(Electron_NEW_BS(k).X);
        FrUSE(k) = sum(Electron_NEW_BS(k).USE)/m3;
        FrBSE(k) = sum(Electron_NEW_BS(k).BSE)/m3;
        FrSCE(k) = sum(Electron_NEW_BS(k).SCE)/m3;
        SHIFT(k) = sum(Electron_NEW_BS(k).SHIFT);
        for l = 1:m3
            THETA = abs(Electron_NEW_BS(k).THETA(l,1));
            if THETA <= Th1d
                Electron_NEW_BS(k).BF(l,1) = 1;
                Electron_NEW_BS(k).DF(l,1) = 0;
                Electron_NEW_BS(k).HAADF(l,1) = 0;
            elseif THETA > Th1d && THETA <= Th3d
                Electron_NEW_BS(k).DF(l,1) = 1;
                Electron_NEW_BS(k).BF(l,1) = 0;
                Electron_NEW_BS(k).HAADF(l,1) = 0;
            elseif THETA > Th3d && THETA <= pi/2
                Electron_NEW_BS(k).HAADF(l,1) = 1;
                Electron_NEW_BS(k).DF(l,1) = 0;
                Electron_NEW_BS(k).BF(l,1) = 0;
            end
        end
        
        BF(k) = sum(Electron_NEW_BS(k).BF(:,1))/m3;
        DF(k) = sum(Electron_NEW_BS(k).DF(:,1))/m3;
        HAADF(k) = sum(Electron_NEW_BS(k).HAADF(:,1))/m3;
        SHIFT(k) = Electron_NEW_BS(k).SHIFT;
        
        
    end
 if MultiCapsle == 0
 % ulozenie vysledku do *.MAT file   
% dt = string(datetime);
% dt(dt==':') = '_';
% NAME = strcat('Electron_Sphere_',dt,'.mat');
% save(NAME,'Electron_NEW_BS');
    
%     figure
%     plot(SHIFT,FrUSE,'b')
%     hold on
%     plot(SHIFT,FrBSE,'g')
%     plot(SHIFT,FrSCE,'r')
%     hold off
%     legend('Unscattered','BackScattered','FrontScattered')
%     title('Sphere');
%     xlabel('dR [nm]'); ylabel('Types of electrons')
    
    
    figure
    plot(SHIFT,BF,'b')
    hold on
    plot(SHIFT,DF,'g')
    plot(SHIFT,HAADF,'r')
    hold off
    legend('BF','DF','HAADF');
    title('Sphere');
    
    ylim([0,1])
    xlabel('dR [A]'); ylabel('Signal')
    close(h)
    
elseif MultiCapsle == 1
%     NAME = strcat('Electron_FilledCapsle_Result',dt,'.mat');
    close(h)
else
    
end

% save ('Electron_NEW_BS.mat', 'Electron_NEW_BS')

