%Copyright WaVe 2012
f= 86; % vband 71-76 GHz
p = 1013; %hPa
t = 0; %degrees celsius
pho = 7.5; %water vapor density in g./m3
El = 49.8; %elevation angle in degrees
h1 = 0;   %altitude 1 Earth Station Height
h2 = 53;  %altitude 2 Height of atmosphere
El = El*pi/180;

%calculation of the dry air attenuation
rp = (p./1013);
rt = 288./(273+t);
eta1 = 6.957566.*(rp.^-0.3461).*(rt.^0.2535).*exp(1.3766.*(1-rt))-1;
eta2 = 42.1309.*(rp.^-0.3068).*(rt.^1.2023).*exp(2.5147.*(1-rt))-1;
c = log(eta2./eta1)./log(3.5);
d = (4.^c)./eta1;
gamma_o_dash = 1.935.*(rp.^1.6657).*(rt.^-3.3714).*exp(-4.1643.*(1-rt));

%dry air attenuation from the ITU model in dB./km
gamma_o = (((0.2296.*d.*gamma_o_dash) ./ (((f-66).^c) + d)) +...
    ((0.286.*(rp.^2).*(rt.^3.8))./((f-118.75).^2+2.97.*rp.^2.*rt.^1.6))).*(f.^2).*10.^-3;

%calculation of attenuation due to water vapor
etaw1 = 0.9544 .* rp.* (rt.^0.69) + 0.0061.*pho;
etaw2 = 0.95.*rp.*(rt.^0.64) + 0.0067.*pho;
etaw3 = 0.9561.*rp.*(rt.^0.67) + 0.0059.*pho;
etaw4 = 0.9543.*rp.*(rt.^0.68) + 0.0061.*pho;
etaw5 = 0.955.*rp.*(rt.^0.68) + 0.006.*pho;

g22 = 1+(((f-22.235).^2)./((f+22.235).^2));
g557 = 1+(((f-557).^2)./((f+557).^2));
g752 = 1+(((f-752).^2)./((f+752).^2));

%water vapor attenuation in dB./km
gamma_w = ((3.13.*10.^-2).*rp.*rt.^2 + 1.76 .* (10.^-3) .* pho .* (rt.^8.5) + ...
(rt.^2.5).*((3.84.*etaw1.*g22.*exp(2.23.*(1-rt))./((f-22.235).^2 + 9.42.*etaw1.^2))+...
(10.48.*etaw2.*exp(0.7.*(1-rt))./((f-183.31).^2+9.48.*etaw2.^2)) + ...
(0.078.*etaw3.*exp(6.4385.*(1-rt))./((f - 321.226).^2+6.29.*etaw3.^2)) +...
(3.76.*etaw4.*exp(1.6.*(1-rt))./((f-325.153).^2 + 9.22.*etaw4.^2)) +...
(26.36.*etaw5.*exp(1.09.*(1-rt))./((f-380).^2)) +...
(17.87.*etaw5.*exp(1.46.*(1-rt))./((f-448).^2)) +...
(883.7.*etaw5.*g557.*exp(0.17.*(1-rt))./((f-557).^2)) +...
(302.6.*etaw5.*g752.*exp(0.41.*(1-rt))./((f-752).^2)))) .* (f.^2) .* pho .* (10.^-4);

%slant path attenuation
ho = f.*((0.039581-1.19751.*(10.^-3).*f +...
    9.14810.*(10.^-6).*(f.^2))./(1-0.028687.*f+2.07858.*(10.^-4).*(f.^2))) +...
    (90.6./((f-60).^2));
hw = 1.65 .* (1 + (1.61./((f-22.23).^2 +2.91)) + (3.33./((f-183.3).^2 + 4.58))+...
    (1.90./((f-325.1).^2)+3.34));

ho_dash = ho.*(exp(-h1./ho) - exp(-h2./ho));
hw_dash = hw.*(exp(-h1./hw) - exp(-h2./hw));

%angle of elevation between 10 degrees and 90 degrees
%attenuation is given by
A_dryair_watervapor = (ho_dash .* gamma_o + gamma_w .* hw_dash)./sin(El)











