clc
clear all
close all

Rh=20;
Th=20;
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=dz:1:dz_count;
ddz=ddz*dz;
Rhh=dz:dz:hmax;

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)

x=[0,100,200,300];
y=[0,0,20.001,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);
% %%
r=vp(x(1),x(2),Th,y(1),y(2),dz_count,dz,hmax);
ro=vo(x(1),x(2),dz,Th,y(1),y(2),r);
R=koef_od(x(2),x(1),y(2),y(1),ro,polarisation,Eps,Th,dz);
E_0=exp(-1j*k.*r)./r;
E_0od=R.*exp(-1j*k.*ro)./ro;
E_p=E_0+E_0od;
figure
plot(ddz,abs(E_p))
hold on
xlabel('vyska nad terenem[m]');
ylabel('E [V/m]');
%%
%dalsi vlnoplochy

for ind=3:1:length(x)
    ind
    E_n=E_n*0; E_n1=E_n1*0; E_n2=E_n2*0;
    K=sqrt(k*x(ind-1)/(2*pi*1j*x(ind)*(x(ind)-x(ind-1))));
    for indd=1:1:dz_count
        if y(ind)+dz*indd<=hmax
            r=vp(x(ind-1),x(ind),indd*dz,y(ind-1),y(ind),dz_count,dz,hmax);
            ro=vo(x(ind-1),x(ind),dz,indd*dz,y(ind-1),y(ind),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
    end
    E_p=E_n.*filter;
    plot(ddz,abs(E_p))
end
hold off

%%
% E_Th=abs(E_n(round(Rh/dz)))

%%
% %UTD
indd=1;
Rhhh=0:dz:hmax;
for Rh=Rhhh
    fi_edge=0;
    s_inc=hypot(x(3),y(3)-Th);
    s_diff=hypot(x(end)-x(3),y(3)-Rh);
    fi_inc_edge=asin(x(3)/s_inc);
    E_01=1*exp(-1j*k*s_inc)/s_inc;
    fi_inc_diff=2*pi-asin((x(end)-x(3))/s_diff);
    A=(s_inc/(s_diff*(s_inc+s_diff)))^0.5;
    D=diff_coeff_rrfc( fi_inc_edge,fi_inc_diff,fi_edge,f,polarisation,s_inc,s_diff);
    E(indd)=E_01*D*A*exp(-1j*k*s_diff);
    if Rh>Th+(x(end))/tan(fi_inc_edge)
        E(indd)=E(indd)-exp(-1j*k*hypot(Rh-Th,x(end)))/hypot(Rh-Th,x(end));
    end
    indd=indd+1;
end

figure('Name','Intenzita jako funkce vysky v poloze prijimace')
plot(ddz(1:round(length(filter)/2)+300),-20*log10(4*pi/lam)+20*log10(abs(E_n(1:round(length(filter)/2)+300))),'r')
hold on
plot(Rhhh(1:round(length(filter)/2)+300),20*log10(abs(E(1:round(length(filter)/2)+300)))-20*log10(4*pi/lam),'g')
title('Porovnani UTD a fyzikalni optiky')
ylabel('uroven vykonu v bode prijimace P [dB]');
xlabel('Vyska nad terenem [m]');
legend('Fyzikalni optika','UTD');
legend('Location','southeast')
hold off