                                                                     
clear all;
close all;
clc; 
 

%% nahrání naměřených dat 
load L1_Gs_d.mat;                        % data Gaussův svazek
data_G = L1_Gs;                          % převedení do proměnné


     load L1_BI1s_d.mat;                % data Top hat svazek
data_Th = L1_BI1s;                      % převedení do proměnné

%% zadání vstupních konstant
Pskutec_G = 2.6e-3;                   % výkon optického svazku změřený měřičem optického výkonu pro Gaussův svazek
proc_sumu_G = 0.08;                    % velikost odfiltrovaného šumu z Gaussova svazku
konst_wx = 0.9255;                     % 2Wx/2Wy - konstanta určující velikost násobku průměru Gaussova svazku v ose x oproti ose y- 
                                       % pro kruhový Gaussův svazek je velikost konstanty = 1
                                        
Pskutec_Th = 1.62e-3;                  % výkon optického svazku změřený měřičem optického výkonu pro Top hat svazek
 proc_sumu_Th = 0.1;                   %nastavení velikosti filtru šumu Top hat svazku 
           fd = 0.20;                   %nastavení dolní rozhodovací úrovně pro metodu delta/2  (oříznutí šumu pro detekci náběžných a sestupných hran)
           fh = 0.8;                   %nastavení horní rozhodovací úrovně pro metodu delta/2  (oříznutí vrchu Top hat svazku pro detekci náběžných a sestupných hran) 
        
      

%% definování prostoru pro vykreslení

x = 9.0515e-6;                           % velikost jednoho pixelu v horizontálním směru v mm
y = 8.3e-6;                              % velikost jednoho pixelu ve vertikálním směru v mm
podstava = 8.3e-6*9.0515e-6;             % plocha jednoho pixelu v mm

xa = x:x:640*x;                         
ya = y:y:480*y;

% matice pro vykreslení
[Xm,Ym] = meshgrid(xa,ya);

%% Přepočet velikosti I naměřených Beam Profilerem na I [W/m^2]_Gaussův svazek

% uložení naměřených dat do matice z (640x480)
k = 1;
for i = 1:length(ya),
  for j = 1:length(xa),
    z_tmp = data_G((k),3);
    z(i,j) = z_tmp;                     
    k = k+1;
  end
end

Pmer = sum(sum(z.*podstava));            %celkový výkon obsažený v naměřeném optickém svazku

konstP = Pskutec_G/Pmer;                %konstanta pro přenásobení matice z 
z_G = z.*konstP;                        %matice (640x480) optické intenzity I [W/m^2]

Imax_G = (max(max(z_G)))                %maximální optická intezita měřeného Gaussova svazku Imax [W/m^2]

% vykreslení Gaussova svazku 
figure;                                 
mesh(Xm,Ym,z_G)                          
title('Naměřený Gaussův svazek')
xlabel('x [m]');ylabel('y [m]');zlabel('I [W / m^2]')

z_Gn= z_G/max(max(z_G));                %normování optické intenzity I [W/m^2] na I/I0 [-]        

% vykreslení normovaného Gaussova svazku 
figure;
mesh(Xm,Ym,z_Gn)
title('Naměřený Gaussův svazek')
xlabel('x [m]');ylabel('y [m]');zlabel('I / I0 [-]')

%% Zanedbání šumu u Gaussova svazku

noise_G = Imax_G * proc_sumu_G;         
k = 1;
for i = 1:length(ya),
  for j = 1:length(xa),
    if z_G(i,j)< noise_G
        z_G_filt(i,j) = 0;
    else
        z_G_filt(i,j) = z_G(i,j);
    end
  end
end

z_G_filt;                                %matice (640x480) optické intenzity I [W/m^2] s odfiltrovaným šumem
P_G_filt = sum(sum(z_G_filt.*podstava)); %optický výkon Gaussova svazku s odfiltrovaným šumem

%% Výpočet ideálního Gaussova svazku

w0 = sqrt( P_G_filt/(0.5*pi*Imax_G));   %výpočet pološířky ideálního Gaussova svazku

w_x = w0/(sqrt(konst_wx));              %přepočet poměru šířky Gaussova svazku v ose x oproti šířce Gausova svazku v ose y
w_y = w0*(sqrt(konst_wx));


%% Vykreslení ideálního Gaussova svazku

sigma_x = w_x/2;                
sigma_y = w_y/2;

% Vstupní data
Lx = x*480;                         % délka osy x
Ly = y*640;                         % délka osy y
amp = Imax_G;                       % maximální amplituda ideálního Gaussova svazku
stred=[0  0];                       % souřadnice vrcholu ideálního Gaussova svazku

% Osy, souřadnice
Nx = 480;
Ny = 640;
X = Lx*((0:Nx-1)-Nx/2+.5)/Nx;       % středy intervalů na ose x od -Lx/2+dx/2 do +Lx/2-dx/2 - neobsahují x = 0
Y = Ly*((0:Ny-1)-Ny/2+.5)/Ny;       % středy intervalů na ose y od -Ly/2+dy/2 do +Ly/2-dy/2 - neobsahují y = 0

% Gaussovské rozložení
 wx = amp*exp(((X-stred(1)).^2)/(-2*sigma_x^2));
 wy = amp*exp(((Y-stred(2)).^2)/(-2*sigma_y^2));
  
 wxy = wy'*wx;                      
 wxyta = (wxy')./(amp*amp);         
 z_G_id = wxyta.*amp;                % matice (640x480) optické intenzity I [W/m^2] ideálního Gaussova svazku            

P_G_id = sum(sum(z_G_id.*podstava)); %optický výkon ideálního Gaussova svazku 

% vykreslení ideálního Gaussova svazku 3D
 figure;
 mesh(Xm,Ym,z_G_id);
  
title('Ideální Gaussův svazek')
xlabel('x [m]');ylabel('y [m]');zlabel('I [W / m^2]')
 

% data v ose X a ose Y (úprava proměnných pro vykreslení)
za = z_G(:,320);
za_n = za/(max(max(za)));
zb = z_G(240,:);
zb_n = zb/(max(max(zb)));
zao = z_G_filt(:,320);
zao_n = zao/(max(max(zao)));
zbo = z_G_filt(240,:);
zbo_n = zbo/(max(max(zbo)));
wxz = z_G_id(:,320);
wxz_n = wxz/(max(max(wxz)));
wyz = z_G_id(240,:);
wyz_n = wyz/(max(max(wyz)));


% vykreslení ideálního Gaussova svazku 2D
figure;
subplot(2,1,2),plot(ya,za_n,'r','LineWidth',2); title('Intenzitní profil id. a měř. Gaussova svazku ve vertikální rovinně');xlabel('y [m]');ylabel('I / I0 (-)'); 
hold on
plot(ya,wxz_n,'LineWidth',2); 
subplot(2,1,1),plot(xa,zb_n,'r','LineWidth',2); title('Intenzitní profil id. a měř. Gaussova svazku ve horizontální rovinně');xlabel('x [m]');ylabel('I / I0 (-)'); 
hold on
plot(xa,wyz_n,'LineWidth',2);

figure;
subplot(2,1,2),plot(ya,zao_n,'r','LineWidth',2); title('Intenzitní profil id. a měř. Gaussova svazku s filt. šumem ve vertikální rovinně');xlabel('y [m]');ylabel('I / I0 (-)'); 
hold on
plot(ya,wxz_n,'LineWidth',2); 
subplot(2,1,1),plot(xa,zbo_n,'r','LineWidth',2); title('Intenzitní profil id. a měř. Gaussova svazku s filt. šumem v horizontáln rovinně');xlabel('x [m]');ylabel('I / I0 (-)'); 
hold on
plot(xa,wyz_n,'LineWidth',2);

%% Výpočet parametru kvality

sm_odch = zeros(480,640);
z_G_filt_n = z_G_filt/max(max(z_G_filt));                   %normování matice dat (640x480) Gaussova svazku s odfiltrovaným šumem 
z_G_id_n = z_G_id/max(max(z_G_id));                         %normování matice dat (640x480) ideálního Gaussova svazku 

poc_nenul_prvku_G = 0;                                      %výpočet počtu nenulových prvků v matici filtrovaného Gaussova svazku
for i = 1:length(ya),
  for j = 1:length(xa),
      if z_G_filt_n(i,j) ~= 0,                                      
            poc_nenul_prvku_G = poc_nenul_prvku_G + 1;
      end
  end
end

 for i = 1:length(ya),
  for j = 1:length(xa),
    z_tmp = z_G_filt_n(i,j);
    w_tmp = z_G_id_n(i,j);
    sm_odch(i,j) = sqrt((z_tmp-w_tmp)^2/(poc_nenul_prvku_G));         % výpočet směrodatné odchylky
    
  end
 end

 var = sum(sum(sm_odch));                                           % součet všech směrodatných odchilek
 smerodatna_odchylka_G = var 
 figure
 mesh(Xm,Ym,sm_odch)
 title('Směrodatné odchylky Gaussova svazku')
 shading interp;xlabel('x [m]');ylabel('y [m]');zlabel('sm. odchylka (x,y)');
 
% %% výpočet procentuelního parametru zvlnění optické intenzity Gaussova svazku
% % z důvodu velké doby výpočtu je tento výpočet defaultně zakomentován
% kk = 0;
% konst_zv = 1.5;                                                     %konstanta pro porovnání velikosti každeého bodu v matici
% konst_posun = 8;                                                    %konstanta pro posun v matici 640x480
% mnoz_zvl = 0;                                                       %matice pro ukládání počtu bodů, které jsou větší než nasobek konst_zv 
% for i = 1:length(ya),       
%   for j = 1:length(xa),
%       kk = kk + 1;
%       prvek_zv(kk) = z_G_filt(i,j);
%       if kk > konst_posun
%         if (prvek_zv(kk-konst_posun)> (prvek_zv(kk)*konst_zv));
%             mnoz_zvl = mnoz_zvl+1;
%         end
%       end
%         
%   end
% end
% param_zvl_G = (mnoz_zvl / 307200)*100                               %výpočet parametru zvlnění v %
% 
% 
 
%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TOP HAT SVAZEK  TOP HAT SVAZEK    TOP HAT SVAZEK    TOP HAT SVAZEK
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Přepočet velikosti I naměřených Beam Profilerem na I [W/m^2]_Top hat svazek

% uložení naměřených dat do matice data_th (640x480)
k = 1;
for i = 1:length(ya),
  for j = 1:length(xa),
    z_tmp_th = data_Th((k),3);
    data_th(i,j) = z_tmp_th;
    k = k+1;
  end
end

Pmer_th = sum(sum(data_th.*podstava));          %celkový výkon obsažený v naměřeném optickém svazku

konstP_th = Pskutec_Th/Pmer_th;                 %konstanta pro přenásobení matice data_th 
data_th_k = data_th.*(konstP_th);                 %matice (640x480) optické intenzity I [W/m^2]

Imax_Th = (max(max(data_th_k)))                 %maximální optická intezita měřeného Top hat svazku Imax [W/m^2]

% vykreslení Top hat svazku
figure;
mesh(Xm,Ym,data_th_k)
title('Naměřený Top hat svazek')
xlabel('x [m]');ylabel('y [m]');zlabel('I [W / m^2]')

data_th_kn = data_th_k/max(max(data_th_k));      %normování optické intenzity I [W/m^2] na I/I0 [-]        

% vykreslení normovaného Top hat svazku 
figure;
mesh(Xm,Ym,data_th_kn)
title('Naměřený Top hat svazek')
xlabel('x [m]');ylabel('y [m]');zlabel('I / I0 [-]')

%% Zanedbání šumu u Top hat svazku
  
noise_Th = Imax_Th *  proc_sumu_Th;
k = 1;
for i = 1:length(ya),
  for j = 1:length(xa),
    if data_th_k(i,j)< noise_Th
        data_th_k_filt(i,j) = 0;
    else
        data_th_k_filt(i,j) = data_th_k(i,j);
    end
  end
end

data_th_k_filt;                                 %matice (640x480) optické intenzity I [W/m^2] s odfiltrovaným šumem
P_Th_filt = sum(sum(data_th_k_filt.*podstava)); %optický výkon Top hat svazku s odfiltrovaným šumem

data_th_k_filtn = data_th_k_filt/max(max(data_th_k_filt));      %normování optické intenzity I [W/m^2] na I/I0 [-]

figure;
mesh(Xm,Ym,data_th_k_filtn)
title('Naměřený Top hat svazek-filt. šum')
xlabel('x [m]');ylabel('y [m]');zlabel('I / I0 [-]')

%% průměrování Top hat svazku i+4 / 5 
for i = 1:(640-5),
    prumer_th_X1(i)=data_th_k(240,i)+data_th_k(240,i+1)+data_th_k(240,i+2)+data_th_k(240,i+3)+data_th_k(240,i+4);%+data_th(320,i+5)+data_th(320,i+6)+data_th(320,i+7)+data_th(320,i+8)+data_th(320,i+9);
    prumer_th_X1(i)=prumer_th_X1(i)/5;
end

for i = 1:(480-5),
    prumer_th_Y1(i)=data_th_k(i,320)+data_th_k(i+1,320)+data_th_k(i+2,320)+data_th_k(i+3,320)+data_th_k(i+4,320);%+data_th(i+5,240)+data_th(i+6,240)+data_th(i+7,240)+data_th(i+8,240)+data_th(i+9,240);
    prumer_th_Y1(i) = prumer_th_Y1(i)/5;
end

xa_2 = x:x:635*x;                         
ya_2 = y:y:475*y;

% vykreslení průmervaného svazku
figure;
subplot(2,1,1),plot(xa_2,prumer_th_X1,'LineWidth',2); title('Porovnání měřeného a vyhlazeného Top hat svazku - horizontální rovinna');xlabel('x [m]');;ylabel('I (W/m^2)'); 
hold on
plot(xa,data_th_k(240,:),'r','LineWidth',1); 
subplot(2,1,2),plot(ya_2,prumer_th_Y1,'LineWidth',2); title('Porovnání měřeného a vyhlazeného Top hat svazku - vertikální rovinna');xlabel('y [m]');ylabel('I (W/m^2)');
hold on
plot(ya,data_th_k(:,320),'r','LineWidth',1);
%% průměrování Top hat svazku i+4
i = 0;
j = 0;
while i< (635),
    i = i+5;
    j = j+1;
    prumer_th_X2(j) = data_th_k(240,i)+data_th_k(240,i+1)+data_th_k(240,i+2)+data_th_k(240,i+3)+data_th_k(240,i+4);%+data_th(i+5,240)+data_th(i+6,240)+data_th(i+7,240)+data_th(i+8,240)+data_th(i+9,240);
    
end

i = 0;
j = 0;
while i< (475),
    i = i+5;
    j = j+1;
    prumer_th_Y2(j) = data_th_k(i,320)+data_th_k(i+1,320)+data_th_k(i+2,320)+data_th_k(i+3,320)+data_th_k(i+4,320);%+data_th(320,i+5)+data_th(320,i+6)+data_th(320,i+7)+data_th(320,i+8)+data_th(320,i+9);
end

%% Výpočet náběžné a sestupné hrany ideálního Top hat svazku metodou I0/2

% výpočet náběžné hrany pro osu x id. Top hat metodou I0/2
uroven1 = 2;                                       % hodnota pro dělení Imax
AB_ch1 = prumer_th_X2;                             % mocnina pro zlepšení detekce
d1_ch1 = AB_ch1 - max(prumer_th_X2)/uroven1;       % nastavení rozhodovací úrovně
d2_ch1 = sign(d1_ch1);                             % kde průchodu nulou
%rídící proměnné
m = 0;
i = 1;
while i < length(d2_ch1)    
    if (d2_ch1(i)==1)
        index2_ch1(i) = i;
        i = i + 100;                                % posuv z důvodu dvou průchodu těsně za sebou 
        m = m + 1;                                  % počet průchodů
    else
        i = i + 1;                                  % inkrementace z důvodu posuvu bez záchytu signalu
    end
end

index2_ch1 = sort(index2_ch1,'descend');            % řazení od největšího
index2a_ch1 = index2_ch1(1:m);                      % zmenšení pole              
hrana = sort(index2a_ch1);                          % seřazení pole 

% výpočet sestupné hrany pro osu x id. Top hat metodou I0/2
m = 0;
i = length(d2_ch1);
n = 0;
while (i > 0)    
    n = n + 1;
    if (d2_ch1(length(d2_ch1)-n)==1)
        index2_ch11(i) = i;
        i = i - 100;                                % posuv z důvodu dvou průchodu těsně za sebou 
        m = m + 1;                                  % počet průchodů
    else
        i = i-1;                                    % inkrementace z důvodu posuvu bez záchytu signalu
    end
end

index2_ch11 = sort(index2_ch11,'descend');        	% řazení od největšího
index2a_ch11 = index2_ch11(1:m);                  	% zmenšení pole              
hrana1 = sort(index2a_ch11);                        % seřazení pole 

[cc dd] = max(hrana1);
prumer_svazku_x = ((cc-hrana(1)).*5)*x;             % průměr svazku v ose x

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% výpočet náběžné hrany pro osu y id. Top hat metodou I0/2
uroven1 = 2;                                        % hodnota pro dělení Imax
AB_ch1_2 = prumer_th_Y2;                            % mocnina pro zlepšení detekce
d1_ch1_2 = AB_ch1_2 - max(prumer_th_Y2)/uroven1;    % nastavení rozhodovací úrovně
d2_ch1_2 = sign(d1_ch1_2);                          % kde průchodu nulou

%rídící proměnné
m = 0;
i = 1;
while i < length(d2_ch1_2)    
    if (d2_ch1_2(i)==1)
        index2_ch1_2(i) = i;
        i = i + 100;                               % posuv z důvodu dvou průchodu těsně za sebou 
        m = m + 1;                                 % počet průchodů
    else
        i = i + 1;                                 % inkrementace z důvodu posuvu bez záchytu signalu
    end
end

index2_ch1_2 = sort(index2_ch1_2,'descend');       % řazení od největšího
index2a_ch1_2 = index2_ch1_2(1:m);                 % zmenšení pole             
hrana = sort(index2a_ch1_2);                       % seřazení pole 

% výpočet sestupné hrany pro osu y id. Top hat metodou I0/2
m = 0;
i = length(d2_ch1_2);
n = 0;
while (i > 0)    
    n = n + 1;
    if (d2_ch1_2(length(d2_ch1_2)-n)==1)
        index2_ch11_2(i) = i;
        i = i - 100;                                % posuv z důvodu dvou průchodu těsně za sebou
        m = m + 1;                                  % počet průchodů
    else
        i = i-1;                                    % inkrementace z důvodu posuvu bez záchytu signalu        
    end
end

index2_ch11_2 = sort(index2_ch11_2,'descend');       % řazení od největšího
index2a_ch11_2 = index2_ch11_2(1:m);                 % zmenšení pole           
hrana1_2 = sort(index2a_ch11_2);                     % seřazení pole  

[cc dd] = max(hrana1_2);
prumer_svazku_y = ((cc-hrana(1)).*5)*y;              % průměr svazku v ose y



dve_a = prumer_svazku_x;                            % šířka elipsy 2a
dve_b = prumer_svazku_y;                            % šířka elipsy 2b

%% Výpočet náběžné a sestupné hrany ideálního Top hat svazku metodou D/2

%uložení filtrovaných svazků o fh a fd pro osu X
%hranice detekce pro metodu D/2--definováno ve vstupních proměnných

maxAmp = (max(prumer_th_X1));                       %maximální I0 z vyhlazeného svazku
A90 = maxAmp * fh;                                  %nastavení horní hranice detekce z max I0
k = 1;

  for j = 1:length(prumer_th_X1),
    if prumer_th_X1(j) > A90
        data_th_90X(j) = A90;                       
    else
        data_th_90X(j) = prumer_th_X1(j);           %matice vyhlazeného Top hat svazku v ose x oříznutá o fh
    end
  end

A10 = maxAmp * fd;
k = 1;
   for j = 1:length(prumer_th_X1),
    if prumer_th_X1(j)< A10
        data_th_10X(j) = A10;
    else
         data_th_10X(j) = prumer_th_X1(j);          %matice vyhlazeného Top hat svazku v ose x oříznutá o fd
    end
   end

% výpočet náběžné hrany v ose x 
skok = 50;                                          %skok zabezpečující přeskočení náhodných špiček šumu
                                                    %nebo propadů v horní části Top hat svazku

j = 0;
while j < length(data_th_10X),
     j = j + 1;    
    if data_th_10X(j+1)~= data_th_10X(j)            %porovnávání kdy se změní hodnota z 0 na kladnou hodnotu
        if (data_th_10X(j+skok)== data_th_10X(j))   %oštření aby cyklus přeskočil náhodné špičky šumu
            j= j + skok;
        else
            x1d = j;                                %dolní bod náběžné hrany v ose x
            break;
        end
    end
   
end

j = round(length(data_th_90X)/2);                   %nastaví j do poloviny pole
while j > 0
    
    if data_th_90X(j-1)~= data_th_90X(j)            %porovnávání kdy se změní honota z maxima na menší
        if (data_th_90X(j-skok)== data_th_90X(j))   %oštření aby cyklus přeskočil náhodné propady v horní části filt. Top hat
            j = j - skok;
        else
            x1h = j;                                 %horní bod náběžné hrany v ose x                        
            break;
        end;
    end
    j = j - 1;
end

% výpočet sestupné hrany v ose x 
skok = 50;                                          %skok zabezpečující přeskočení náhodných špiček šumu
                                                    %nebo propadů v horní části Top hat svazku

jj = length(data_th_10X);                           %nastaví jj na konec pole
while jj > 0                                      
         if (data_th_10X(jj-1)~= data_th_10X(jj))        %porovnávání kdy se změní hodnota z 0 na kladnou hodnotu
            if (data_th_10X(jj-skok)== data_th_10X(jj))  %oštření aby cyklus přeskočil náhodné špičky šumu
                jj = jj - skok;
            else
                x2d = jj;                                %dolní bod sestupné hrany v ose x
                break
            end;
        end
    jj = jj - 1;
end
                                                    
j = 0;
while j < length(data_th_90X)
    if j > round(length(data_th_90X)/2)                 %nastaví j do poloviny pole
        if data_th_90X(j-1)~= data_th_90X(j)            %porovnávání kdy se změní honota z maxima na menší
            if(data_th_90X(j+skok)== data_th_90X(j))    %oštření aby cyklus přeskočil náhodné propady v horní části filt. Top hat
                j = j + skok;
            else
                x2h = j;                                %horní bod sestupné hrany v ose x                     
                break;
            end
        end
    end
    j = j + 1;
end

% výpočet D/2 v ose x pro náběžnou a sestupnou hranu
hranaX1 = (x1h + x1d)/2;
hranaX2 = (x2h + x2d)/2; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%uložení filtrovaných svazků o fh a fd pro osu Y
%hranice detekce pro metodu D/2--definováno ve vstupních proměnných

maxAmp = (max(prumer_th_Y1));                           %maximální I0 z vyhlazeného svazku
A90 = maxAmp * fh;                                      %nastavení horní hranice detekce z max I0
k = 1;

  for j = 1:length(prumer_th_Y1),
    if prumer_th_Y1(j) > A90
        data_th_90Y(j) = A90;
    else
        data_th_90Y(j) = prumer_th_Y1(j);                 %matice vyhlazeného Top hat svazku v ose y oříznutá o fh
    end
  end

A10 = maxAmp * fd;
k = 1;
   for j = 1:length(prumer_th_Y1),
    if prumer_th_Y1(j)< A10
        data_th_10Y(j) = A10;
    else
         data_th_10Y(j) = prumer_th_Y1(j);                %matice vyhlazeného Top hat svazku v ose y oříznutá o fd
    end
    end
  
  
% výpočet náběžné hrany v ose y 
skok = 50;                                              %skok zabezpečující přeskočení náhodných špiček šumu   
                                                        %nebo propadů v horní části Top hat svazku            
j = 0;
while j < length(data_th_10Y),
     j = j + 1;    
    if data_th_10Y(j+1)~= data_th_10Y(j)                %porovnávání kdy se změní hodnota z 0 na kladnou hodnotu
        if (data_th_10Y(j+skok)== data_th_10Y(j))       %oštření aby cyklus přeskočil náhodné špičky šumu
            j= j + skok;
        else
            y1d = j;                                    %dolní bod náběžné hrany v ose y
            break;
        end
    end
   
end

j = round(length(data_th_90Y)/2);                       %nastaví j do poloviny pole
while j > 0
    
    if data_th_90Y(j-1)~= data_th_90Y(j)                %porovnávání kdy se změní honota z maxima na menší
        if (data_th_90Y(j-skok)== data_th_90Y(j))       %oštření aby cyklus přeskočil náhodné propady v horní části filt. Top hat
            j = j - skok;
        else
            y1h = j;                                    %horní bod náběžné hrany v ose y
            break;
        end;
    end
    j = j - 1;
end

% výpočet sestupné hrany v ose y 
skok = 50;                                          %skok zabezpečující přeskočení náhodných špiček šumu
                                                    %nebo propadů v horní                                                   
                                                    %části Top hat svazku
                                                    
jj = length(data_th_10Y);                           %nastaví jj na konec pole
while jj > 0
      
        if (data_th_10Y(jj-1)~= data_th_10Y(jj))        %porovnávání kdy se změní hodnota z 0 na kladnou hodnotu
            if (data_th_10Y(jj-skok)== data_th_10Y(jj)) %oštření aby cyklus přeskočil náhodné špičky šumu
                jj = jj - skok;
            else
                y2d = jj;                               %dolní bod sestupné hrany v ose y
                break
            end;
        end
    jj = jj - 1;
end

j = 0;
while j < length(data_th_90Y)
    if j > round(length(data_th_90Y)/2)                 %nastaví j do poloviny pole
        if data_th_90Y(j-1)~= data_th_90Y(j)            %porovnávání kdy se změní honota z maxima na menší
            if(data_th_90Y(j+skok)== data_th_90Y(j))    %oštření aby cyklus přeskočil náhodné propady v horní části filt. Top hat
                j = j + skok;
            else
                y2h = j;                                %horní bod sestupné hrany v ose y
                break;
            end
        end
    end
    j = j + 1;
end


% výpočet D/2 v ose y pro náběžnou a sestupnou hranu
hranaY1 = (y1h + y1d)/2;
hranaY2 = (y2h + y2d)/2; 

dve_a_2 = (hranaX2 - hranaX1)*x;                    % šířka elipsy 2a
dve_b_2 = (hranaY2 - hranaY1)*y;                    % šířka elipsy 2b

%% Výpočet a vykreslení id. Top hat svazku metodou I0/2
a = dve_a/2;                                        %délka hlavní poloosy elipsy (podstava id.Top hat je elipsa)
b = dve_b/2;                                        %délka vedlejší poloosy elipsy

r = sqrt(b*a);                                      %poloměr id Top hat 
Imax_Th_id1 = P_Th_filt / ((pi*r^2));               %výpočet Imax id. Top hat

stred=[320*x  240*y];                               % souřadnice vrcholu id. Top hatu
for i = 1:length(xa)
     for j = 1:length(ya)
           if sqrt((b^2*(xa(i)-stred(1))^2)+(a^2*(ya(j)-stred(2))^2))>(a*b)     %výpočet id. Top hatu
                wxy(i,j)=0;
           else
               wxy(i,j)=Imax_Th_id1;
           end
     end
end

wxyt = wxy';

%Vykreslení ideálního Top hat svazku metodou Io/2
 figure;                                                
 mesh(Xm,Ym,wxyt);                                                                                        
 xlabel('x [m]');ylabel('y [m]');zlabel('I [W / m^2]')
 title('Ideální Top hat svazek-metoda Io/2');
 
 %Normování dat
 maxAmp= max(max(data_th_k_filt));
 if maxAmp < Imax_Th_id1
    norma_1 = Imax_Th_id1;
else
    norma_1 = maxAmp;
 end
wxytn = wxyt/max(max(norma_1));
data_th_k_filt_n1 = data_th_k_filt /max(max(norma_1));

%vykreslení měř. a id. Top hat svazku metoda I0/2 
figure;
subplot(2,1,2),plot(ya,wxytn(:,320),'LineWidth',2); title('Intenzitní profil id. a měř. Top hat s filt. šumem-metoda Io/2-vert. rovinna');xlabel('y [m]'); 
hold on
plot(ya,data_th_k_filt_n1(:,320),'r','LineWidth',1); xlabel('y [m]');ylabel('I / I0 [-]');
subplot(2,1,1),plot(xa,wxytn(240,:),'LineWidth',2); title('Intenzitní profil id. a měř. Top hat s filt. šumem-metoda Io/2-horiz. rovinna');xlabel('x [m]');ylabel('I / I0 [-]');
hold on
plot(xa,data_th_k_filt_n1(240,:),'r')

%% Výpočet id. Top hat svazku metodou D/2

a_2 = dve_a_2/2;                                        %délka hlavní poloosy elipsy (podstava id.Top hat je elipsa)
b_2 = dve_b_2/2;                                        %délka vedlejší poloosy elipsy

r_2 = sqrt(b_2*a_2);                                    %poloměr id Top hat 
Imax_Th_id2 = P_Th_filt / ((pi*r_2^2));                 %výpočet Imax id. Top hat

stred=[320*x  240*y];                                   % souřadnice vrcholu id. Top hatu
for i = 1:length(xa)
     for j = 1:length(ya)
           if sqrt((b_2^2*(xa(i)-stred(1))^2)+(a_2^2*(ya(j)-stred(2))^2))>(a_2*b_2) %výpočet id. Top hatu
                wxy_2(i,j)=0;
           else
               wxy_2(i,j)=Imax_Th_id2;
           end
     end
end

wxyt_2 = wxy_2';

%Vykreslení ideálního Top hat svazku metodou D/2
 figure;                                                
 mesh(Xm,Ym,wxyt_2);                                                                                       
 xlabel('x [m]');ylabel('y [m]');zlabel('I [W / m^2]')
 title('Ideální Top hat svazek-metoda D/2');

%Normování dat
maxAmp= max(max(data_th_k_filt));
if maxAmp < Imax_Th_id2
    norma_2 = Imax_Th_id2;
else
    norma_2 = maxAmp;
end
wxyt_2 = wxyt_2/max(max(norma_2));
data_th_k_filt_n2 = data_th_k_filt /max(max(norma_2));

%vykreslení měř. a id. Top hat svazku metoda D/2 
figure;
subplot(2,1,2),plot(ya,wxyt_2(:,320),'LineWidth',2); title('Intenzitní profil id. a měř. Top hat s filt. šumem-metoda D/2-vert. rovinna');xlabel('y [m]');  
hold on
plot(ya,data_th_k_filt_n2(:,320),'r','LineWidth',1); xlabel('y [m]');ylabel('I / I0 [-]');
subplot(2,1,1),plot(xa,wxyt_2(240,:),'LineWidth',2);title('Intenzitní profil id. a měř. Top hat s filt. šumem-metoda D/2-horiz. rovinna');xlabel('x [m]');ylabel('I / I0 [-]');
hold on
plot(xa,data_th_k_filt_n2(240,:),'r')




%% Výpočet parametru kvality Top hat svazku pro metodu Io/2

sm_odch_Th1 = zeros(480,640);
wxytn=wxyt/max(max(wxyt));                                    %normování ideálního Top hat svazku - metoda I0/2

poc_nenul_prvku_TH1 = 0;                                      %výpočet počtu nenulových prvků v matici filtrovaného Top hat svazku
for i = 1:length(ya),
  for j = 1:length(xa),
      if data_th_k_filt(i,j) ~= 0,                                      
            poc_nenul_prvku_TH1 = poc_nenul_prvku_TH1 + 1;
      end
  end
end

 for i = 1:length(ya),
  for j = 1:length(xa),
    z_tmp = data_th_k_filtn(i,j);
    w_tmp = wxytn(i,j);
    sm_odch_Th1(i,j) = sqrt((z_tmp-w_tmp)^2/(poc_nenul_prvku_TH1));  % výpočet směrodatné odchylky
  end
 end

 var = sum(sum(sm_odch_Th1));                                        % součet všech směrodatných odchilek
 param_kvality_Th_I0_2 = var 
 figure;
 mesh(Xm,Ym,sm_odch_Th1)
 title('Směrodatné odchylky Top hat svazku-I_0/2')
shading interp;xlabel('x [m]');ylabel('y [m]');zlabel('sm. odchylka (x,y)');

%% Výpočet parametru kvality Top hat svazku pro metodu D/2

sm_odch_Th2 = zeros(480,640);
wxytn=0;
wxytn=wxyt_2/max(max(wxyt_2));                                 %normování ideálního Top hat svazku - metoda D/2

poc_nenul_prvku_TH2 = 0;                                      %výpočet počtu nenulových prvků v matici filtrovaného Top hat svazku
for i = 1:length(ya),
  for j = 1:length(xa),
      if data_th_k_filt(i,j) ~= 0,                                      
            poc_nenul_prvku_TH2 = poc_nenul_prvku_TH2 + 1;
      end
  end
end

 for i = 1:length(ya),
  for j = 1:length(xa),
    z_tmp = data_th_k_filtn(i,j);
    w_tmp = wxytn(i,j);
    sm_odch_Th2(i,j) = sqrt((z_tmp-w_tmp)^2/(poc_nenul_prvku_TH2)); % výpočet směrodatné odchylky
  end
 end

 var_D = sum(sum(sm_odch_Th2));                                 % součet všech směrodatných odchilek
 param_kvality_Th_D_2 = var_D
 figure;
 mesh(Xm,Ym,sm_odch_Th2)
 title('Směrodatné odchylky Top hat svazku-D/2')
 shading interp;xlabel('x [m]');ylabel('y [m]');zlabel('sm. odchylka (x,y)');

%  %% výpočet procentuelního parametru zvlnění optické intenzity Top hat svazku
%  % z důvodu velké doby výpočtu je tento výpočet defaultně zakomentován
% kk = 0;
% konst_zv = 1.5;                                                         %konstanta pro porovnání velikosti každeého bodu v matici
% konst_posun = 8;                                                        %konstanta pro posun v matici 640x480
% mnoz_zvl = 0;                                                           %matice pro ukládání počtu bodů, které jsou větší než nasobek konst_zv 
% for i = 1:length(ya),
%   for j = 1:length(xa),
%       kk = kk + 1;
%       prvek_zv(kk) = data_th_k_filt(i,j);
%       if kk > konst_posun
%         if (prvek_zv(kk-konst_posun)> (prvek_zv(kk)*konst_zv));
%             mnoz_zvl = mnoz_zvl+1;
%         end
%       end
%         if data_th_k_filt(i,j) ~= 0,
%             pocetnulovychprvku = pocetnulovychprvku + 1;
%         end
%   end
% end
% param_zvl_Th = (mnoz_zvl / 307200)*100                      %výpočet parametru zvlnění v %
%  
%  
% 
% 
