close all;
clear all;

h1 = 0.8;
h2 = 1;
r = 10;
P = 100;

% Izotropní zářič
D = ones(1,360);



%% Prijmaci antena - charky
 
T_both = xlsread('3142_BiConiLog_97.xls'); 
 T_ver = T_both(:,1:8);
 T_hor = T_both(:,9:16);


%% r = 10 m

r = 10;

for h2 = 1:0.05:4
    i = 1;
   
for freq = 30e6:10e6:1000e6
    
    if freq <= 400e6
         Tx_hor = db2mag(T_hor(:,1));
         Tx_ver = db2mag(T_ver(:,1));
    elseif freq <= 500e6
        Tx_hor = db2mag(T_hor(:,2));
        Tx_ver = db2mag(T_ver(:,2));
    elseif freq <= 600e6
        Tx_hor = db2mag(T_hor(:,3));
        Tx_ver = db2mag(T_ver(:,3));
    elseif freq <= 700e6
        Tx_hor = db2mag(T_hor(:,4));
        Tx_ver = db2mag(T_ver(:,4));
    elseif freq <= 800e6
        Tx_hor = db2mag(T_hor(:,5));
        Tx_ver = db2mag(T_ver(:,5));
    elseif freq <= 900e6
        Tx_hor = db2mag(T_hor(:,6));
        Tx_ver = db2mag(T_ver(:,6));
    elseif freq <= 1000e6
        Tx_hor = db2mag(T_hor(:,7));
        Tx_ver = db2mag(T_ver(:,7));
    end
    
r1 = sqrt((h2-h1)^2 + r^2);
r2 = sqrt((h1+h2)^2 + r^2);

del_r = r2 - r1;

lam = 3e8/freq;
k = 2*pi/lam;

alfaD = (180/pi)*asin((h2-h1)/r1);
betaD = (180/pi)*asin((h2+h1)/r2);

Up = round(alfaD);
Uo = round(360-betaD);

Upt = round(alfaD);
Uot = round(betaD);


 E1(i) = (sqrt(30*P*D(Up))/r1);
 E2(i) = (sqrt(30*P*D(Uo))/r2);

 Ev(i) = abs( Tx_ver(1)*E1(i) + Tx_ver(1)*E2(i)*exp(j*k*del_r));
 Eh(i) = abs( Tx_hor(1)*E1(i) - Tx_hor(1)*E2(i)*exp(j*k*del_r));
  
 Evt(i) = abs( Tx_ver(Upt)*E1(i) + Tx_ver(Uot)*E2(i)*exp(j*k*del_r));
 Eht(i) = abs( Tx_hor(Upt)*E1(i) - Tx_hor(Uot)*E2(i)*exp(j*k*del_r));
 
 if h2 == 1
      Evmax(i) = Ev(i);
      Ehmax(i) = Eh(i);
      
      
      Evmaxt(i) = Evt(i);
      Ehmaxt(i) = Eht(i);   
 
 else
    if Evmax(i) < Ev(i)
       Evmax(i) = Ev(i);
       
    end
    
    if Evmaxt(i) < Evt(i)
       Evmaxt(i) = Evt(i);
       
    end

    if Ehmax(i) < Eh(i)
       Ehmax(i) = Eh(i)
    end
    
    if Ehmaxt(i) < Eht(i)
       Ehmaxt(i) = Eht(i)
    end
    
 end
 
 
i = i+1;
end

i = 1;
end

Evmax_dB = 20*LOG10(Evmax*1e6);
Ehmax_dB = 20*LOG10(Ehmax*1e6);
Evmaxt_dB = 20*LOG10(Evmaxt*1e6);
Ehmaxt_dB = 20*LOG10(Ehmaxt*1e6);


%% r = 3 m

r = 3;

for h2 = 1:0.05:4
    i = 1;
   
for freq = 30e6:10e6:1000e6
    
    if freq <= 400e6
        Tx_hor = db2mag(T_hor(:,1));
        Tx_ver = db2mag(T_ver(:,1));
    elseif freq <= 500e6
        Tx_hor = db2mag(T_hor(:,2));
        Tx_ver = db2mag(T_ver(:,2));
    elseif freq <= 600e6
        Tx_hor = db2mag(T_hor(:,3));
        Tx_ver = db2mag(T_ver(:,3));
    elseif freq <= 700e6
        Tx_hor = db2mag(T_hor(:,4));
        Tx_ver = db2mag(T_ver(:,4));
    elseif freq <= 800e6
        Tx_hor = db2mag(T_hor(:,5));
        Tx_ver = db2mag(T_ver(:,5));
    elseif freq <= 900e6
        Tx_hor = db2mag(T_hor(:,6));
        Tx_ver = db2mag(T_ver(:,6));
    elseif freq <= 1000e6
        Tx_hor = db2mag(T_hor(:,7));
        Tx_ver = db2mag(T_ver(:,7));
    end

 
    
r1 = sqrt((h2-h1)^2 + r^2);
r2 = sqrt((h1+h2)^2 + r^2);


del_r = r2 - r1;

lam = 3e8/freq;
k = 2*pi/lam;

alfaD = (180/pi)*asin((h2-h1)/r1);
betaD = (180/pi)*asin((h2+h1)/r2);

Up = round(alfaD);
Uo = round(360-betaD);

Upt = round(alfaD);
Uot = round(betaD);


 E1(i) = (sqrt(30*P*D(Up))/r1);
 E2(i) = (sqrt(30*P*D(Uo))/r2);

 Ev3(i) = abs( Tx_ver(1)*E1(i) + Tx_ver(1)*E2(i)*exp(j*k*del_r));
 Eh3(i) = abs( Tx_hor(1)*E1(i) - Tx_hor(1)*E2(i)*exp(j*k*del_r));
  
 Evt3(i) = abs( Tx_ver(Upt)*E1(i) + Tx_ver(Uot)*E2(i)*exp(j*k*del_r));
 Eht3(i) = abs( Tx_hor(Upt)*E1(i) - Tx_hor(Uot)*E2(i)*exp(j*k*del_r));
 
 if h2 == 1
      Evmax3(i) = Ev3(i);
      Ehmax3(i) = Eh3(i);
      
      
      Evmaxt3(i) = Evt3(i);
      Ehmaxt3(i) = Eht3(i);   
 
 else
    if Evmax3(i) < Ev3(i)
       Evmax3(i) = Ev3(i);
       
    end
    
    if Evmaxt3(i) < Evt3(i)
       Evmaxt3(i) = Evt3(i);
       
    end

    if Ehmax3(i) < Eh3(i)
       Ehmax3(i) = Eh3(i)
    end
    
    if Ehmaxt3(i) < Eht3(i)
       Ehmaxt3(i) = Eht3(i)
    end
    
 end
 
 
i = i+1;
end

Xplot = 30:10:1000


i = 1;
end

Evmax_dB3 = 20*LOG10(Evmax3*1e6);
Ehmax_dB3 = 20*LOG10(Ehmax3*1e6);
Evmaxt_dB3 = 20*LOG10(Evmaxt3*1e6);
Ehmaxt_dB3 = 20*LOG10(Ehmaxt3*1e6);


%%
figure(1)
set(gcf, 'color', [1 1 1])
subplot(2,1,1)
semilogx(Xplot, Evmax_dB3 - Evmaxt_dB3, 'r','LineWidth',2)
hold on;
semilogx(Xplot, Evmax_dB - Evmaxt_dB, 'g','LineWidth',2)
set(gca,'FontSize',25)
legend('3 m', '10 m ', 'Location','northwest','Orientation','vertical')
title('Vertikální polarizace');
xlabel('Frekvence [MHz]');
ylabel('Rozdíl [dB]');
set(gca,'FontSize',20)

subplot(2,1,2)
semilogx(Xplot, Ehmax_dB3 - Ehmaxt_dB3, 'r','LineWidth',2)
hold on;
semilogx(Xplot, Ehmax_dB - Ehmaxt_dB, 'g','LineWidth',2)
set(gca,'FontSize',25)
legend('3 m', '10 m ', 'Location','northwest','Orientation','vertical')
%ylim([5 15])
title('Horizontální polarizace');
xlabel('Frekvence [MHz]');
ylabel('Rozdíl [dB]');
set(gca,'FontSize',20)


%%

    
   freq_400 = 400e6
         Tx_hor_400 = db2mag(T_hor(:,1));
         Tx_ver_400 = db2mag(T_ver(:,1));
  
    freq_1000 = 1000e6
        Tx_hor_1000 = db2mag(T_hor(:,7));
        Tx_ver_1000 = db2mag(T_ver(:,7));
    
        
     i = 1;
   
for h2 = 1:0.005:4
         
r1_10 = sqrt((h2-h1)^2 + 100);
r2_10 = sqrt((h1+h2)^2 + 100);
r1_3 = sqrt((h2-h1)^2 + 9);
r2_3 = sqrt((h1+h2)^2 + 9);

del_r_10 = r2_10 - r1_10;
del_r_3 = r2_3 - r1_3;

lam_400 = 3e8/freq_400;
k_400 = 2*pi/lam_400;
lam_1000 = 3e8/freq_1000;
k_1000 = 2*pi/lam_1000;


alfaD_10 = (180/pi)*asin((h2-h1)/r1_10);
betaD_10 = (180/pi)*asin((h2+h1)/r2_10);
alfaD_3 = (180/pi)*asin((h2-h1)/r1_3);
betaD_3 = (180/pi)*asin((h2+h1)/r2_3);


U_prm_10 = round(alfaD_10);
U_odr_10 = round(betaD_10);
U_prm_3 = round(alfaD_3);
U_odr_3 = round(betaD_3);

V3_400(i) =  Tx_ver_400(U_prm_3);
V3_1000(i) =  Tx_ver_1000(U_prm_3);

H3_400(i) =  Tx_hor_400(U_prm_3);
H3_1000(i) =  Tx_hor_1000(U_prm_3);

V10_400(i) =  Tx_ver_400(U_prm_10);
V10_1000(i) =  Tx_ver_1000(U_prm_10);

H10_400(i) =  Tx_hor_400(U_prm_10);
H10_1000(i) =  Tx_hor_1000(U_prm_10);

i = i+1;
end

V3_400db =   20*LOG10(V3_400);
V3_1000db =   20*LOG10(V3_1000);

H3_400db =   20*LOG10(H3_400);
H3_1000db =   20*LOG10(H3_1000);

V10_400db =   20*LOG10(V10_400);
V10_1000db =  20*LOG10(V10_1000);

H10_400db =   20*LOG10(H10_400);
H10_1000db =   20*LOG10(H10_1000);

Xplot = 1:0.005:4

O = 1;
for Q = 1:1:length(Xplot)
      for i = 1:1:1
        Xplots(O) = Xplot(Q);
        
        V3_400dbs(O) = V3_400db(Q);
        V3_1000dbs(O) = V3_1000db(Q);
        
        H3_400dbs(O) = H3_400db(Q);
        H3_1000dbs(O) = H3_1000db(Q);
        
        V10_400dbs(O) = V10_400db(Q);
        V10_1000dbs(O) = V10_1000db(Q);
        
        H10_400dbs(O) = H10_400db(Q);
        H10_1000dbs(O) = H10_1000db(Q);
                         
            
       O = O+1;
    end
end
   W = window(@hann, 10);
   Xplotn = filtfilt(W, 4.5, Xplots);
   
   V3_400dbn = filtfilt(W, 4.5, V3_400dbs);
     V3_1000dbn = filtfilt(W, 4.5, V3_1000dbs);
       H3_400dbn = filtfilt(W, 4.5, H3_400dbs);
         H3_1000dbn = filtfilt(W, 4.5, H3_1000dbs); 
         
           V10_400dbn = filtfilt(W, 4.5, V10_400dbs);
             V10_1000dbn = filtfilt(W, 4.5, V10_1000dbs);
               H10_400dbn = filtfilt(W, 4.5, H10_400dbs);
                 H10_1000dbn = filtfilt(W, 4.5, H10_1000dbs);
   
                 p = 5;
                 
  pV3_400db = polyfit(Xplot,V3_400db,p);
  pV3_1000db = polyfit(Xplot,V3_1000db,p);
  pH3_400db = polyfit(Xplot,H3_400db,p);
  pH3_1000db = polyfit(Xplot,H3_1000db,p);

  V3_400dbnn = polyval(pV3_400db,Xplot);
    V3_1000dbnn = polyval(pV3_1000db,Xplot);
      H3_400dbnn = polyval(pH3_400db,Xplot);
      H3_1000dbnn = polyval(pH3_1000db,Xplot);


figure(2)
set(gcf, 'color', [1 1 1])
subplot(2,1,1)


plot(Xplotn, V3_400dbnn, 'r','LineWidth',2)
hold on;
plot(Xplotn, V3_1000dbnn, 'g','LineWidth',2)

set(gca,'FontSize',25)
legend('400 MHz', '1000 MHz', 'Location','northwest','Orientation','vertical')
title('Vertikální polarizace');
xlabel('Výška antény [m]');
ylabel('Rozdíl [dB]');
set(gca,'FontSize',20)

subplot(2,1,2)

plot(Xplotn, H3_400dbnn, 'r','LineWidth',2)
hold on;
plot(Xplotn, H3_1000dbnn, 'g','LineWidth',2)

set(gca,'FontSize',25)
legend('400 MHz', '1000 MHz', 'Location','northwest','Orientation','vertical')
%ylim([5 15])
title('Horizontální polarizace');
xlabel('Výška antény [m]');
ylabel('Rozdíl [dB]');
set(gca,'FontSize',20)


