clear all
clc
close all

%%
%nacteni hodnot a print terenu
vys='vysledky.xlsx';
fid=fopen('Zidovice_Nabdin_Zabaged_100m_Vinc.txt','r');
sizeA=[2 Inf];
Data=fscanf(fid,'%f',sizeA);
x=Data(1,:)*1E3;
y=Data(2,:);
xx=x(1:7:length(x));
yy=y(1:7:length(y));
yy(2)=y(2);
yy(end)=y(end);
figure(1)
% plot(x,y)

plot(x,y)
hold on
xlabel('x [m]');
ylabel('y nadmorska vyska [m]');
title('Tvar terenu');

ref_grad=-6;
ke=1/(1+6378*ref_grad*(10e-6));
Re=ke*6371;
for index=1:1:length(xx)
bb(index)=xx(index)*(xx(end)-xx(index))/(2*Re*10e3)
yyr(index)=yy(index)+bb(index);
end
plot(xx,yy)
legend('Puvodni','Aproximace')
legend('Location','south')
hold off
qwe
%%
Th=12.9;
Rh=5.2;
f=2e9;
c=3e8;
lam=c/f;
k=2*pi/lam;
dz=lam/8;
hmax=350;
dz_count=round((hmax-min(y))/dz);
Eps=10;
polarisation='vertical';
ddz=1:1:dz_count;
ddz=ddz*dz;
Rhh=0:dz:500;

filter=(1+cos(-pi+2*pi/(dz_count):2*pi/(dz_count):pi))/2;
filter(1:round(length(filter)/2))=1;

r=zeros(1,dz_count); E_0=zeros(1,dz_count); E_n=zeros(1,dz_count); r_0=zeros(1,dz_count); E_n1=zeros(1,dz_count);E_n2=zeros(1,dz_count);
Mat_mesh=zeros(length(xx),dz_count);

%%
r=vp(xx(1),xx(2),Th,yy(1),yy(2),dz_count,dz,hmax);
ro=vo(xx(1),xx(2),dz,Th,yy(1),yy(2),r);
R=koef_od(xx(2),xx(1),yy(2),yy(1),ro,polarisation,Eps,Th,dz);
E_0(1:round((hmax-yy(2))/dz))=exp(-1j*k.*r(1:round((hmax-yy(2))/dz)))./r(1:round((hmax-yy(2))/dz));

% E_0(~isfinite(E_0))=0;
E_0od=R.*exp(-1j*k.*ro)./ro;
E_p=E_0+E_0od;

Mat_mesh(2,dz_count-round((hmax-yy(2))/dz)+1:dz_count)=E_p(1:round((hmax-yy(2))/dz));
figure(2)
plot(20*log10(abs(E_p)))

for ind=3:1:length(xx)
    ind
    E_n=E_n*0; E_n1=E_n1*0; E_n2=E_n2*0;
    K=sqrt(k*xx(ind-1)/(2*pi*1j*xx(ind)*(xx(ind)-xx(ind-1))));
    for indd=1:1:dz_count
        if yy(ind)+indd*dz<=hmax
            r=vp(xx(ind-1),xx(ind),indd*dz,yy(ind-1),yy(ind),dz_count,dz,hmax);
            ro=vo(xx(ind-1),xx(ind),dz,indd*dz,yy(ind-1),yy(ind),r);
            R=koef_od(xx(ind),xx(ind-1),yy(ind),yy(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
    figure
    plot(20*log10(abs(E_n)))
    E_p=E_n.*filter;
    Mat_mesh(ind,dz_count-round((hmax-yy(ind))/dz)+1:dz_count)=E_p(1:round((hmax-yy(ind))/dz));
    
end

E_Th=abs(E_n(round(Rh/dz)));
L_noref=20*log10(1/E_Th)
mesh(20*log10(abs(Mat_mesh))-20*log10(4*pi/lam)+20*log10(4*pi*x(end)/lam))
%%
% % ddz=dz:1:dz_count;
% % ddz=ddz*dz;
% % dddz=0:50;
% % figure
% % plot(ddz,20*log10(abs(E_p))-20*log10(4*pi/lam)+20*log10(4*pi*x(end)/lam))
% % E_k=Mat_mesh(end,:);
% % figure
% % plot(ddz,20*log10(abs(E_k))-20*log10(4*pi/lam)+20*log10(4*pi*x(end)/lam))
% %     for i=1:1:length(Rhh)
% %         utmer2FO(i)=20*log10(abs(E_k(round(Rhh/dz))));
% %     end
% %     
% %     figure
% %     
% %     for ind=1:1:length(Rhhh)
% % ut2gg(ind)=SS(round(Rhhh(ind)/dz+50/dz));
% % end
% % figure
% % hold on
% % plot(Rhhh,Ltot(:,1),'-.r*')
% % plot(Rhhh,utmer2,'-.b+')
% % plot(Rhhh,deygout2,'-.go')
% % plot(Rhh,-ut2gg,'-.o')
% % hold off
% % xlabel('vyska prijimace Rxh [m]'); 
% % ylabel('L difracni utlum [dB]'); 
% % title('Porovnani s merenim f=2GHz');
% % legend('UTD','Mereni','Deygout','Fyzikalni optika') 