close all
clear all
clc

Rh=10;
Th=10;
f=1e9;
c=3e8;
lam=c/f;
k=2*pi/lam;
dz=lam/8;
hmax=60;
dz_count=round(hmax/dz);
polarisation='vertical';
Eps=10;

ddz=1:1:dz_count;
ddz=ddz*dz;

Rhh=dz:dz:60;
% filter=ones(1,dz_count);
filter=(1+cos(-pi+2*pi/(dz_count):2*pi/(dz_count):pi))/2;
filter(1:round(length(filter)/2))=1;
% figure
% plot(Rhh, filter)
% title('Filtracni funkce')
% ylabel('Amplituda [-]');
% xlabel('vyska nad terenem [m]');
% qwe
D=500;
x=[0,100,200,300,400,500];
y=[0,0,0,0,0,0];
r=zeros(1,length(x)); E_0=zeros(1,dz_count); E_n1=zeros(1,dz_count);E_n2=zeros(1,dz_count);
E_n=zeros(1,dz_count);E_0od=zeros(1,dz_count);
%%
%prvni vlnoplocha
r=vp(x(1),x(2),Th,0,0,dz_count,dz,hmax);
ro=vo(x(1),x(2),dz,Th,0,0,r);
R=koef_od(x(2),x(1),y(2),y(1),ro,polarisation,Eps,Rh,dz);
E_0=exp(-1j*k.*r)./r;
E_0od=R.*exp(-1j*k.*ro)./ro;
E_p=E_0+E_0od; %pomocny vektor do ktereho se ukladaji mezivysledky
plot(ddz,abs(E_p))
hold on
xlabel('vyska nad terenem[m]');
ylabel('E [V/m]');
%%

%dalsi vlnoplochy
for ind=3:1:length(x)
    
    K=sqrt(k*x(ind-1)/(2*pi*1j*x(ind)*(x(ind)-x(ind-1))));
    for indd=1:1:dz_count
        
        r=vp(x(ind-1),x(ind),indd*dz,0,0,dz_count,dz,hmax);
        ro=vo(x(ind-1),x(ind),dz,indd*dz,0,0,r);
        R=koef_od(x(ind),x(ind-1),y(ind),y(ind-1),ro,polarisation,Eps,indd*dz,dz);
        E_n1(indd)=K*(sum(E_p.*exp(-1j*k.*r))*dz);
        E_n2(indd)=K*(sum(E_p.*exp(-1j*k.*ro).*R)*dz);
        E_n(indd)=E_n1(indd)+E_n2(indd);
    end
    E_p=E_n.*filter;
    plot(ddz,abs(E_p))
end

hold off
%%
E_Th=abs(E_n(round(Rh/dz)))
E_Th_v=abs(exp(-1j*k*D)/D)
%%
r=vp(x(1),D,Th,0,0,dz_count,dz,hmax);
ro=vo(x(1),D,dz,Th,0,0,r);R=zeros(1,length(ro));
R=koef_od(D,x(1),y(2),y(1),ro,polarisation,Eps,Th,dz);


E_2=exp(-1j*k.*r)./r+R.*exp(-1j*k.*ro)./ro;
figure('Name','Intenzita jako funkce vysky v poloze prijimace')
plot(ddz(1:round(length(filter)/2)),abs(E_n(1:round(length(filter)/2))),'r')
hold on
plot(ddz(1:round(length(filter)/2)),abs(E_2(1:round(length(filter)/2))),'g')
title('Obecny koeficient odrazu')
ylabel('E [V/m]');
xlabel('Vyska nad terenem [m]');
legend('Fyzikalni optika','Dvoupaprskovy model');
legend('Location','southeast')
hold off

%%
%prvni vlnoplocha
r=vp(x(1),x(2),Th,0,0,dz_count,dz,hmax);
ro=vo(x(1),x(2),dz,Th,0,0,r);
R=1;
E_0=exp(-1j*k.*r)./r;
E_0od=R.*exp(-1j*k.*ro)./ro;
E_p=E_0+E_0od; %pomocny vektor do ktereho se ukladaji mezivysledky
%%
%dalsi vlnoplochy
for ind=3:1:length(x)
    K=sqrt(k*x(ind-1)/(2*pi*1j*x(ind)*(x(ind)-x(ind-1))));
    for indd=1:1:dz_count
        r=vp(x(ind-1),x(ind),indd*dz,0,0,dz_count,dz,hmax);
        ro=vo(x(ind-1),x(ind),dz,indd*dz,0,0,r);
        E_n1(indd)=K*(sum(E_p.*exp(-1j*k.*r))*dz);
        E_n2(indd)=K*(sum(E_p.*exp(-1j*k.*ro).*R)*dz);
        E_n(indd)=E_n1(indd)+E_n2(indd);
    end
    E_p=E_n.*filter;
    
end
r=vp(x(1),D,Th,0,0,dz_count,dz,hmax);
ro=vo(x(1),D,dz,Th,0,0,r);
E_2=exp(-1j*k.*r)./r+R.*exp(-1j*k.*ro)./ro;
figure('Name','Intenzita jako funkce vysky v poloze prijimace')
plot(ddz(1:round(length(filter)/2)),abs(E_n(1:round(length(filter)/2))),'r')
hold on
plot(ddz(1:round(length(filter)/2)),abs(E_2(1:round(length(filter)/2))),'g')
title('R=1')
ylabel('E [V/m]');
xlabel('vyska nad terenem [m]');
legend('Fyzikalni optika','Dvoupaprskovy model');
legend('Location','southeast')
hold off


