re = (4/3)*6371000;
hr = 5:5:100;
ht = 500
lambda = 0.03
f = 10000
r = 10000
fi = pi
p = (2/sqrt(3))*sqrt(re.*(ht+hr)+((r.^2)/4))
ksi = asin((2*re.*r.*(ht-hr))/p.^3)
r1 = (r/2) - (p.*sin(ksi./3))
r2 = r - r1
phi1 = r1/re
phi2 = r2/re
R1 = sqrt((hr.^2)+4*re.*(re+hr).*(sin(phi1./2)).^2)
R2 = sqrt((ht.^2)+4*re.*(re+ht).*(sin(phi2./2)).^2)
Rd = sqrt(((ht-hr).^2)+4*(re+ht).*(re+hr).*(sin((phi1+phi2)/2)).^2)
psi = asin((ht./R1)-(R1./(2*re)))
deltaR = (4*R1.*R2.*(sin(psi)).^2)/(R1+R2+Rd)
eps = 15 + ((i * 0.15)/(2*pi.*f))
arg1 = eps-((cos(psi)).^2);
arg2 = sqrt(arg1);
arg3 = sin(psi);
arg4 = eps.*arg3;
rv = (arg4-arg2)./(arg4+arg2)
rh = (arg3-arg2)./(arg3+arg2);
ro = rv/(exp(i*fi))
deltaPhi = ((2*pi)./lambda).*deltaR
alpha = deltaPhi + fi
F = abs(1+ro.*exp(i*alpha))
%figure(1)
plot(hr, F)
axis( [0 400 0 4]) 
grid on;
xlabel('Wysokość nadajnika');
ylabel('Współczynnik F');