matlab codes for rain attenuation for Ku
function [atten]=attfunc8(f,R01,th,hs,lat,tau,p);

%R01 : rain rate exceeded 0.01% of the time
%lat : station latitude (deg)
%hs  : ground station height (km)
%th  : path elevation angle (deg)
%f   : frequency (GHz)
%Re  : effective Earth radius
%tau : polarization tilt angle
%p   : percent of time that output rain rate is exceeded (can be an array
%      if desired)

Re=8500; % 4/3 Earth radius


% %Dar Es Salaam
 lat=6.52;   R01=99;     hs=0.053;     th=52;  ho=4.5;   f=11.356;   

% hfr=ho+0.36;

 hfr=5-0.075*(abs(lat)-23);

if (abs(lat)<23) hfr=5; end;
if (hfr<0) hfr=0; end;


Ls=(hfr-hs)/sin(th*pi/180);
if (th<5) Ls=2*(hfr-hs)/(sqrt(sin(th*pi/180)^2+2*(hfr-hs)/Re)+sin(th*pi/180)); end;
Lg=Ls*cos(th*pi/180);

% gam=raingam(f,th*pi/180,tau*pi/180,R01);

% gam=kei.*(R01).^alpha
gam=4.9643

rh=1/(1+0.78*sqrt(Lg*gam/f)-0.38*(1-exp(-2*Lg)));
zeta=atan2(hfr-hs,Lg*rh);
if (zeta>th*pi/180) 
    Lr=Lg*rh/cos((th*pi)/180);
else
    Lr=(hfr-hs)/sin(th*pi/180);
end
chi=0;

if (abs(lat)<36) chi=36-abs(lat); 
end


rv=1/( 1+ sqrt(sin(th*pi/180)) *((31*(1-exp(-th/(1+chi)))*sqrt(Lr*gam)/f^2)-0.45));

Le=Lr*rv;

A01=Le*gam
itu_p=[5 3 2 1 .5 .3 .2 .1 .05 .03 .02 .01 .005 .003 .002 .001];
     
  p=itu_p'; 
for i1=1:length(p)
    pp(i1)=p(i1)
    if (pp(i1)>1) 
        z=0;
    elseif (abs(lat)>36)
        z=0;
    elseif ((abs(lat)<36) & (th>25))
        z=-0.005*(abs(lat)-36);
    else
        z=-0.005*(abs(lat)-36)+1.8-4.25*sin(th*pi/180);
    end
    
    atten(i1)=A01*(pp(i1)./0.01)^-(0.655+0.033*log(pp(i1))-0.045*log(A01)-z*sin(th*pi/180)*(1-pp(i1)));
     
    
end

          
       semilogx(pp,atten,'rs-','LineWidth',1.5)
  
       
       hold on
% end
                
 %kibaha
 
  lat=6.50;   R01=94;     hs=0.167;     th=52;  ho=4.5;   f=11.356;   

% hfr=ho+0.36;

 hfr=5-0.075*(abs(lat)-23);

if (abs(lat)<23) hfr=5; end;
if (hfr<0) hfr=0; end;


Ls=(hfr-hs)/sin(th*pi/180);
if (th<5) Ls=2*(hfr-hs)/(sqrt(sin(th*pi/180)^2+2*(hfr-hs)/Re)+sin(th*pi/180)); end;
Lg=Ls*cos(th*pi/180);

% gam=raingam(f,th*pi/180,tau*pi/180,R01);

% gam=kei.*(R01).^alpha
gam=4.6645

rh=1/(1+0.78*sqrt(Lg*gam/f)-0.38*(1-exp(-2*Lg)));
zeta=atan2(hfr-hs,Lg*rh);
if (zeta>th*pi/180) 
    Lr=Lg*rh/cos((th*pi)/180);
else
    Lr=(hfr-hs)/sin(th*pi/180);
end
chi=0;

if (abs(lat)<36) chi=36-abs(lat); 
end


rv=1/( 1+ sqrt(sin(th*pi/180)) *((31*(1-exp(-th/(1+chi)))*sqrt(Lr*gam)/f^2)-0.45));

Le=Lr*rv;

A01=Le*gam
itu_p=[5 3 2 1 .5 .3 .2 .1 .05 .03 .02 .01 .005 .003 .002 .001];
     
  p=itu_p'; 
for i1=1:length(p)
    pp(i1)=p(i1)
    if (pp(i1)>1) 
        z=0;
    elseif (abs(lat)>36)
        z=0;
    elseif ((abs(lat)<36) & (th>25))
        z=-0.005*(abs(lat)-36);
    else
        z=-0.005*(abs(lat)-36)+1.8-4.25*sin(th*pi/180);
    end
    
    atten(i1)=A01*(pp(i1)./0.01)^-(0.655+0.033*log(pp(i1))-0.045*log(A01)-z*sin(th*pi/180)*(1-pp(i1)));
     
    
end


          
       semilogx(pp,atten,'b*-','LineWidth',1.5)  

%Mtwara
lat=10.21;   R01=98;     hs=0.113;     th=50;  ho=4.5;   f=11.356;   

% hfr=ho+0.36;

 hfr=5-0.075*(abs(lat)-23);

if (abs(lat)<23) hfr=5; end;
if (hfr<0) hfr=0; end;


Ls=(hfr-hs)/sin(th*pi/180);
if (th<5) Ls=2*(hfr-hs)/(sqrt(sin(th*pi/180)^2+2*(hfr-hs)/Re)+sin(th*pi/180)); end;
Lg=Ls*cos(th*pi/180);

% gam=raingam(f,th*pi/180,tau*pi/180,R01);

% gam=kei.*(R01).^alpha
gam=4.9041

rh=1/(1+0.78*sqrt(Lg*gam/f)-0.38*(1-exp(-2*Lg)));
zeta=atan2(hfr-hs,Lg*rh);
if (zeta>th*pi/180) 
    Lr=Lg*rh/cos((th*pi)/180);
else
    Lr=(hfr-hs)/sin(th*pi/180);
end
chi=0;

if (abs(lat)<36) chi=36-abs(lat); 
end


rv=1/( 1+ sqrt(sin(th*pi/180)) *((31*(1-exp(-th/(1+chi)))*sqrt(Lr*gam)/f^2)-0.45));

Le=Lr*rv;

A01=Le*gam
itu_p=[5 3 2 1 .5 .3 .2 .1 .05 .03 .02 .01 .005 .003 .002 .001];
     
  p=itu_p'; 
for i1=1:length(p)
    pp(i1)=p(i1)
    if (pp(i1)>1) 
        z=0;
    elseif (abs(lat)>36)
        z=0;
    elseif ((abs(lat)<36) & (th>25))
        z=-0.005*(abs(lat)-36);
    else
        z=-0.005*(abs(lat)-36)+1.8-4.25*sin(th*pi/180);
    end
    
    atten(i1)=A01*(pp(i1)./0.01)^-(0.655+0.033*log(pp(i1))-0.045*log(A01)-z*sin(th*pi/180)*(1-pp(i1)));
     
    
end
          
       semilogx(pp,atten,'g+-','LineWidth',1.5)
       
 %Kilwa
 lat=8.55;   R01=95;     hs=0.014;     th=51;  ho=4.5;   f=11.356;   

% hfr=ho+0.36;

 hfr=5-0.075*(abs(lat)-23);

if (abs(lat)<23) hfr=5; end;
if (hfr<0) hfr=0; end;


Ls=(hfr-hs)/sin(th*pi/180);
if (th<5) Ls=2*(hfr-hs)/(sqrt(sin(th*pi/180)^2+2*(hfr-hs)/Re)+sin(th*pi/180)); end;
Lg=Ls*cos(th*pi/180);

% gam=raingam(f,th*pi/180,tau*pi/180,R01);

% gam=kei.*(R01).^alpha
gam=4.7242

rh=1/(1+0.78*sqrt(Lg*gam/f)-0.38*(1-exp(-2*Lg)));
zeta=atan2(hfr-hs,Lg*rh);
if (zeta>th*pi/180) 
    Lr=Lg*rh/cos((th*pi)/180);
else
    Lr=(hfr-hs)/sin(th*pi/180);
end
chi=0;

if (abs(lat)<36) chi=36-abs(lat); 
end


rv=1/( 1+ sqrt(sin(th*pi/180)) *((31*(1-exp(-th/(1+chi)))*sqrt(Lr*gam)/f^2)-0.45));

Le=Lr*rv;

A01=Le*gam
itu_p=[5 3 2 1 .5 .3 .2 .1 .05 .03 .02 .01 .005 .003 .002 .001];
     
  p=itu_p'; 
for i1=1:length(p)
    pp(i1)=p(i1)
    if (pp(i1)>1) 
        z=0;
    elseif (abs(lat)>36)
        z=0;
    elseif ((abs(lat)<36) & (th>25))
        z=-0.005*(abs(lat)-36);
    else
        z=-0.005*(abs(lat)-36)+1.8-4.25*sin(th*pi/180);
    end
    
    atten(i1)=A01*(pp(i1)./0.01)^-(0.655+0.033*log(pp(i1))-0.045*log(A01)-z*sin(th*pi/180)*(1-pp(i1)));
     
    
end
          
       semilogx(pp,atten,'k*-','LineWidth',1.5)
       
 %Pemba 
 lat=6.25;   R01=109;     hs=0.046;     th=52;  ho=4.5;   f=11.356;   

% hfr=ho+0.36;

 hfr=5-0.075*(abs(lat)-23);

if (abs(lat)<23) hfr=5; end;
if (hfr<0) hfr=0; end;


Ls=(hfr-hs)/sin(th*pi/180);
if (th<5) Ls=2*(hfr-hs)/(sqrt(sin(th*pi/180)^2+2*(hfr-hs)/Re)+sin(th*pi/180)); end;
Lg=Ls*cos(th*pi/180);

% gam=raingam(f,th*pi/180,tau*pi/180,R01);

% gam=kei.*(R01).^alpha
gam=5.5728

rh=1/(1+0.78*sqrt(Lg*gam/f)-0.38*(1-exp(-2*Lg)));
zeta=atan2(hfr-hs,Lg*rh);
if (zeta>th*pi/180) 
    Lr=Lg*rh/cos((th*pi)/180);
else
    Lr=(hfr-hs)/sin(th*pi/180);
end
chi=0;

if (abs(lat)<36) chi=36-abs(lat); 
end


rv=1/( 1+ sqrt(sin(th*pi/180)) *((31*(1-exp(-th/(1+chi)))*sqrt(Lr*gam)/f^2)-0.45));

Le=Lr*rv;

A01=Le*gam
itu_p=[5 3 2 1 .5 .3 .2 .1 .05 .03 .02 .01 .005 .003 .002 .001];
     
  p=itu_p'; 
for i1=1:length(p)
    pp(i1)=p(i1)
    if (pp(i1)>1) 
        z=0;
    elseif (abs(lat)>36)
        z=0;
    elseif ((abs(lat)<36) & (th>25))
        z=-0.005*(abs(lat)-36);
    else
        z=-0.005*(abs(lat)-36)+1.8-4.25*sin(th*pi/180);
    end
    
    atten(i1)=A01*(pp(i1)./0.01)^-(0.655+0.033*log(pp(i1))-0.045*log(A01)-z*sin(th*pi/180)*(1-pp(i1)));
     
    
end

          
       semilogx(pp,atten,'yo-','LineWidth',1.5)
       
  %Tanga
       
 lat=5.5;   R01=103;     hs=0.049;     th=52;  ho=4.5;   f=11.356;   

% hfr=ho+0.36;

 hfr=5-0.075*(abs(lat)-23);

if (abs(lat)<23) hfr=5; end;
if (hfr<0) hfr=0; end;


Ls=(hfr-hs)/sin(th*pi/180);
if (th<5) Ls=2*(hfr-hs)/(sqrt(sin(th*pi/180)^2+2*(hfr-hs)/Re)+sin(th*pi/180)); end;
Lg=Ls*cos(th*pi/180);

% gam=raingam(f,th*pi/180,tau*pi/180,R01);

% gam=kei.*(R01).^alpha
gam=5.2063

rh=1/(1+0.78*sqrt(Lg*gam/f)-0.38*(1-exp(-2*Lg)));
zeta=atan2(hfr-hs,Lg*rh);
if (zeta>th*pi/180) 
    Lr=Lg*rh/cos((th*pi)/180);
else
    Lr=(hfr-hs)/sin(th*pi/180);
end
chi=0;

if (abs(lat)<36) chi=36-abs(lat); 
end


rv=1/( 1+ sqrt(sin(th*pi/180)) *((31*(1-exp(-th/(1+chi)))*sqrt(Lr*gam)/f^2)-0.45));

Le=Lr*rv;

A01=Le*gam
itu_p=[5 3 2 1 .5 .3 .2 .1 .05 .03 .02 .01 .005 .003 .002 .001];
     
  p=itu_p'; 
for i1=1:length(p)
    pp(i1)=p(i1)
    if (pp(i1)>1) 
        z=0;
    elseif (abs(lat)>36)
        z=0;
    elseif ((abs(lat)<36) & (th>25))
        z=-0.005*(abs(lat)-36);
    else
        z=-0.005*(abs(lat)-36)+1.8-4.25*sin(th*pi/180);
    end
    
    atten(i1)=A01*(pp(i1)./0.01)^-(0.655+0.033*log(pp(i1))-0.045*log(A01)-z*sin(th*pi/180)*(1-pp(i1)));
     
    
end

          
       semilogx(pp,atten,'cd-','LineWidth',1.5)
       
  %Unguja
  lat=6.13;   R01=112;     hs=0.018;     th=52;  ho=4.5;   f=11.356;   

% hfr=ho+0.36;

 hfr=5-0.075*(abs(lat)-23);

if (abs(lat)<23) hfr=5; end;
if (hfr<0) hfr=0; end;


Ls=(hfr-hs)/sin(th*pi/180);
if (th<5) Ls=2*(hfr-hs)/(sqrt(sin(th*pi/180)^2+2*(hfr-hs)/Re)+sin(th*pi/180)); end;
Lg=Ls*cos(th*pi/180);

% gam=raingam(f,th*pi/180,tau*pi/180,R01);

% gam=kei.*(R01).^alpha
gam=5.7576

rh=1/(1+0.78*sqrt(Lg*gam/f)-0.38*(1-exp(-2*Lg)));
zeta=atan2(hfr-hs,Lg*rh);
if (zeta>th*pi/180) 
    Lr=Lg*rh/cos((th*pi)/180);
else
    Lr=(hfr-hs)/sin(th*pi/180);
end
chi=0;

if (abs(lat)<36) chi=36-abs(lat); 
end


rv=1/( 1+ sqrt(sin(th*pi/180)) *((31*(1-exp(-th/(1+chi)))*sqrt(Lr*gam)/f^2)-0.45));

Le=Lr*rv;

A01=Le*gam
itu_p=[5 3 2 1 .5 .3 .2 .1 .05 .03 .02 .01 .005 .003 .002 .001];
     
  p=itu_p'; 
for i1=1:length(p)
    pp(i1)=p(i1)
    if (pp(i1)>1) 
        z=0;
    elseif (abs(lat)>36)
        z=0;
    elseif ((abs(lat)<36) & (th>25))
        z=-0.005*(abs(lat)-36);
    else
        z=-0.005*(abs(lat)-36)+1.8-4.25*sin(th*pi/180);
    end
    
    atten(i1)=A01*(pp(i1)./0.01)^-(0.655+0.033*log(pp(i1))-0.045*log(A01)-z*sin(th*pi/180)*(1-pp(i1)));
     
    
end

          
       semilogx(pp,atten,'k<-','LineWidth',1.5)
        grid on
        title('Rain Attenuation Cumulative Distribution')
        xlabel('Percentage of time exceeded (%)','FontSize',14)
        ylabel('Rain attenuation(dB)','FontSize',14)band

Matlab codes for rain attenuation at Ka band

function [atten]=attfunc9(f,R01,th,hs,lat,tau,p);

%R01 : rain rate exceeded 0.01% of the time
%lat : station latitude (deg)
%hs  : ground station height (km)
%th  : path elevation angle (deg)
%f   : frequency (GHz)
%Re  : effective Earth radius
%tau : polarization tilt angle
%p   : percent of time that output rain rate is exceeded (can be an array
%      if desired)

Re=8500; % 4/3 Earth radius


% %Dar Es Salaam
 lat=6.52;   R01=99;     hs=0.053;     th=52;  ho=4.5;   f=21.749;   

% This is from Allnutt paper, not sure if same
% as p.839

% hfr=ho+0.36;

 hfr=5-0.075*(abs(lat)-23);

if (abs(lat)<23) hfr=5; end;
if (hfr<0) hfr=0; end;


Ls=(hfr-hs)/sin(th*pi/180);
if (th<5) Ls=2*(hfr-hs)/(sqrt(sin(th*pi/180)^2+2*(hfr-hs)/Re)+sin(th*pi/180)); end;
Lg=Ls*cos(th*pi/180);

% gam=raingam(f,th*pi/180,tau*pi/180,R01);

% gam=kei.*(R01).^alpha
gam=13.1169

rh=1/(1+0.78*sqrt(Lg*gam/f)-0.38*(1-exp(-2*Lg)));
zeta=atan2(hfr-hs,Lg*rh);
if (zeta>th*pi/180) 
    Lr=Lg*rh/cos((th*pi)/180);
else
    Lr=(hfr-hs)/sin(th*pi/180);
end
chi=0;

if (abs(lat)<36) chi=36-abs(lat); 
end


rv=1/( 1+ sqrt(sin(th*pi/180)) *((31*(1-exp(-th/(1+chi)))*sqrt(Lr*gam)/f^2)-0.45));

Le=Lr*rv;

A01=Le*gam
itu_p=[5 3 2 1 .5 .3 .2 .1 .05 .03 .02 .01 .005 .003 .002 .001];
     
  p=itu_p'; 
for i1=1:length(p)
    pp(i1)=p(i1)
    if (pp(i1)>1) 
        z=0;
    elseif (abs(lat)>36)
        z=0;
    elseif ((abs(lat)<36) & (th>25))
        z=-0.005*(abs(lat)-36);
    else
        z=-0.005*(abs(lat)-36)+1.8-4.25*sin(th*pi/180);
    end
    
    atten(i1)=A01*(pp(i1)./0.01)^-(0.655+0.033*log(pp(i1))-0.045*log(A01)-z*sin(th*pi/180)*(1-pp(i1)));
     
    
end

          
       semilogx(pp,atten,'rs-','LineWidth',1)
  
       
       hold on
% end
                
 %kibaha
 
  lat=6.50;   R01=94;     hs=0.167;     th=52;  ho=4.5;   f=21.749;   


% hfr=ho+0.36;

 hfr=5-0.075*(abs(lat)-23);

if (abs(lat)<23) hfr=5; end;
if (hfr<0) hfr=0; end;


Ls=(hfr-hs)/sin(th*pi/180);
if (th<5) Ls=2*(hfr-hs)/(sqrt(sin(th*pi/180)^2+2*(hfr-hs)/Re)+sin(th*pi/180)); end;
Lg=Ls*cos(th*pi/180);

% gam=raingam(f,th*pi/180,tau*pi/180,R01);

% gam=kei.*(R01).^alpha
gam=12.4313

rh=1/(1+0.78*sqrt(Lg*gam/f)-0.38*(1-exp(-2*Lg)));
zeta=atan2(hfr-hs,Lg*rh);
if (zeta>th*pi/180) 
    Lr=Lg*rh/cos((th*pi)/180);
else
    Lr=(hfr-hs)/sin(th*pi/180);
end
chi=0;

if (abs(lat)<36) chi=36-abs(lat); 
end


rv=1/( 1+ sqrt(sin(th*pi/180)) *((31*(1-exp(-th/(1+chi)))*sqrt(Lr*gam)/f^2)-0.45));

Le=Lr*rv;

A01=Le*gam
itu_p=[5 3 2 1 .5 .3 .2 .1 .05 .03 .02 .01 .005 .003 .002 .001];
     
  p=itu_p'; 
for i1=1:length(p)
    pp(i1)=p(i1)
    if (pp(i1)>1) 
        z=0;
    elseif (abs(lat)>36)
        z=0;
    elseif ((abs(lat)<36) & (th>25))
        z=-0.005*(abs(lat)-36);
    else
        z=-0.005*(abs(lat)-36)+1.8-4.25*sin(th*pi/180);
    end
    
    atten(i1)=A01*(pp(i1)./0.01)^-(0.655+0.033*log(pp(i1))-0.045*log(A01)-z*sin(th*pi/180)*(1-pp(i1)));
     
    
end


          
       semilogx(pp,atten,'b*-','LineWidth',1)
     
%Mtwara
lat=10.21;   R01=98;     hs=0.113;     th=50;  ho=4.5;   f=21.749;   

% hfr=ho+0.36;

 hfr=5-0.075*(abs(lat)-23);

if (abs(lat)<23) hfr=5; end;
if (hfr<0) hfr=0; end;


Ls=(hfr-hs)/sin(th*pi/180);
if (th<5) Ls=2*(hfr-hs)/(sqrt(sin(th*pi/180)^2+2*(hfr-hs)/Re)+sin(th*pi/180)); end;
Lg=Ls*cos(th*pi/180);

% gam=raingam(f,th*pi/180,tau*pi/180,R01);

% gam=kei.*(R01).^alpha
gam=12.9797

rh=1/(1+0.78*sqrt(Lg*gam/f)-0.38*(1-exp(-2*Lg)));
zeta=atan2(hfr-hs,Lg*rh);
if (zeta>th*pi/180) 
    Lr=Lg*rh/cos((th*pi)/180);
else
    Lr=(hfr-hs)/sin(th*pi/180);
end
chi=0;

if (abs(lat)<36) chi=36-abs(lat); 
end


rv=1/( 1+ sqrt(sin(th*pi/180)) *((31*(1-exp(-th/(1+chi)))*sqrt(Lr*gam)/f^2)-0.45));

Le=Lr*rv;

A01=Le*gam
itu_p=[5 3 2 1 .5 .3 .2 .1 .05 .03 .02 .01 .005 .003 .002 .001];
     
  p=itu_p'; 
for i1=1:length(p)
    pp(i1)=p(i1)
    if (pp(i1)>1) 
        z=0;
    elseif (abs(lat)>36)
        z=0;
    elseif ((abs(lat)<36) & (th>25))
        z=-0.005*(abs(lat)-36);
    else
        z=-0.005*(abs(lat)-36)+1.8-4.25*sin(th*pi/180);
    end
    
    atten(i1)=A01*(pp(i1)./0.01)^-(0.655+0.033*log(pp(i1))-0.045*log(A01)-z*sin(th*pi/180)*(1-pp(i1)));
     
    
end
          
       semilogx(pp,atten,'g+-','LineWidth',1)
       
 %Kilwa
 lat=8.55;   R01=95;     hs=0.014;     th=51;  ho=4.5;   f=21.749;   

% hfr=ho+0.36;

 hfr=5-0.075*(abs(lat)-23);

if (abs(lat)<23) hfr=5; end;
if (hfr<0) hfr=0; end;


Ls=(hfr-hs)/sin(th*pi/180);
if (th<5) Ls=2*(hfr-hs)/(sqrt(sin(th*pi/180)^2+2*(hfr-hs)/Re)+sin(th*pi/180)); end;
Lg=Ls*cos(th*pi/180);

% gam=raingam(f,th*pi/180,tau*pi/180,R01);

% gam=kei.*(R01).^alpha
gam=12.5683

rh=1/(1+0.78*sqrt(Lg*gam/f)-0.38*(1-exp(-2*Lg)));
zeta=atan2(hfr-hs,Lg*rh);
if (zeta>th*pi/180) 
    Lr=Lg*rh/cos((th*pi)/180);
else
    Lr=(hfr-hs)/sin(th*pi/180);
end
chi=0;

if (abs(lat)<36) chi=36-abs(lat); 
end


rv=1/( 1+ sqrt(sin(th*pi/180)) *((31*(1-exp(-th/(1+chi)))*sqrt(Lr*gam)/f^2)-0.45));

Le=Lr*rv;

A01=Le*gam
itu_p=[5 3 2 1 .5 .3 .2 .1 .05 .03 .02 .01 .005 .003 .002 .001];
     
  p=itu_p'; 
for i1=1:length(p)
    pp(i1)=p(i1)
    if (pp(i1)>1) 
        z=0;
    elseif (abs(lat)>36)
        z=0;
    elseif ((abs(lat)<36) & (th>25))
        z=-0.005*(abs(lat)-36);
    else
        z=-0.005*(abs(lat)-36)+1.8-4.25*sin(th*pi/180);
    end
    
    atten(i1)=A01*(pp(i1)./0.01)^-(0.655+0.033*log(pp(i1))-0.045*log(A01)-z*sin(th*pi/180)*(1-pp(i1)));
     
    
end
          
       semilogx(pp,atten,'k*-','LineWidth',1)

       
      %Pemba 
       lat=6.25;   R01=109;     hs=0.046;     th=52;  ho=4.5;   f=21.749;   

% hfr=ho+0.36;

 hfr=5-0.075*(abs(lat)-23);

if (abs(lat)<23) hfr=5; end;
if (hfr<0) hfr=0; end;


Ls=(hfr-hs)/sin(th*pi/180);
if (th<5) Ls=2*(hfr-hs)/(sqrt(sin(th*pi/180)^2+2*(hfr-hs)/Re)+sin(th*pi/180)); end;
Lg=Ls*cos(th*pi/180);

% gam=raingam(f,th*pi/180,tau*pi/180,R01);

% gam=kei.*(R01).^alpha
gam=14.4917

rh=1/(1+0.78*sqrt(Lg*gam/f)-0.38*(1-exp(-2*Lg)));
zeta=atan2(hfr-hs,Lg*rh);
if (zeta>th*pi/180) 
    Lr=Lg*rh/cos((th*pi)/180);
else
    Lr=(hfr-hs)/sin(th*pi/180);
end
chi=0;

if (abs(lat)<36) chi=36-abs(lat); 
end


rv=1/( 1+ sqrt(sin(th*pi/180)) *((31*(1-exp(-th/(1+chi)))*sqrt(Lr*gam)/f^2)-0.45));

Le=Lr*rv;

A01=Le*gam
itu_p=[5 3 2 1 .5 .3 .2 .1 .05 .03 .02 .01 .005 .003 .002 .001];
     
  p=itu_p'; 
for i1=1:length(p)
    pp(i1)=p(i1)
    if (pp(i1)>1) 
        z=0;
    elseif (abs(lat)>36)
        z=0;
    elseif ((abs(lat)<36) & (th>25))
        z=-0.005*(abs(lat)-36);
    else
        z=-0.005*(abs(lat)-36)+1.8-4.25*sin(th*pi/180);
    end
    
    atten(i1)=A01*(pp(i1)./0.01)^-(0.655+0.033*log(pp(i1))-0.045*log(A01)-z*sin(th*pi/180)*(1-pp(i1)));
     
    
end

          
       semilogx(pp,atten,'yo-','LineWidth',1)
       
       %Tanga
       
 lat=5.5;   R01=103;     hs=0.049;     th=52;  ho=4.5;   f=21.749;   


% hfr=ho+0.36;

 hfr=5-0.075*(abs(lat)-23);

if (abs(lat)<23) hfr=5; end;
if (hfr<0) hfr=0; end;


Ls=(hfr-hs)/sin(th*pi/180);
if (th<5) Ls=2*(hfr-hs)/(sqrt(sin(th*pi/180)^2+2*(hfr-hs)/Re)+sin(th*pi/180)); end;
Lg=Ls*cos(th*pi/180);

% gam=raingam(f,th*pi/180,tau*pi/180,R01);

% gam=kei.*(R01).^alpha
gam=13.6662

rh=1/(1+0.78*sqrt(Lg*gam/f)-0.38*(1-exp(-2*Lg)));
zeta=atan2(hfr-hs,Lg*rh);
if (zeta>th*pi/180) 
    Lr=Lg*rh/cos((th*pi)/180);
else
    Lr=(hfr-hs)/sin(th*pi/180);
end
chi=0;

if (abs(lat)<36) chi=36-abs(lat); 
end


rv=1/( 1+ sqrt(sin(th*pi/180)) *((31*(1-exp(-th/(1+chi)))*sqrt(Lr*gam)/f^2)-0.45));

Le=Lr*rv;

A01=Le*gam
itu_p=[5 3 2 1 .5 .3 .2 .1 .05 .03 .02 .01 .005 .003 .002 .001];
     
  p=itu_p'; 
for i1=1:length(p)
    pp(i1)=p(i1)
    if (pp(i1)>1) 
        z=0;
    elseif (abs(lat)>36)
        z=0;
    elseif ((abs(lat)<36) & (th>25))
        z=-0.005*(abs(lat)-36);
    else
        z=-0.005*(abs(lat)-36)+1.8-4.25*sin(th*pi/180);
    end
    
    atten(i1)=A01*(pp(i1)./0.01)^-(0.655+0.033*log(pp(i1))-0.045*log(A01)-z*sin(th*pi/180)*(1-pp(i1)));
     
    
end

          
       semilogx(pp,atten,'cd-','LineWidth',1)
       
  %Unguja
  lat=6.13;   R01=112;     hs=0.018;     th=52;  ho=4.5;   f=21.749;   


% hfr=ho+0.36;

 hfr=5-0.075*(abs(lat)-23);

if (abs(lat)<23) hfr=5; end;
if (hfr<0) hfr=0; end;


Ls=(hfr-hs)/sin(th*pi/180);
if (th<5) Ls=2*(hfr-hs)/(sqrt(sin(th*pi/180)^2+2*(hfr-hs)/Re)+sin(th*pi/180)); end;
Lg=Ls*cos(th*pi/180);

% gam=raingam(f,th*pi/180,tau*pi/180,R01);

% gam=kei.*(R01).^alpha
gam=14.9051

rh=1/(1+0.78*sqrt(Lg*gam/f)-0.38*(1-exp(-2*Lg)));
zeta=atan2(hfr-hs,Lg*rh);
if (zeta>th*pi/180) 
    Lr=Lg*rh/cos((th*pi)/180);
else
    Lr=(hfr-hs)/sin(th*pi/180);
end
chi=0;

if (abs(lat)<36) chi=36-abs(lat); 
end


rv=1/( 1+ sqrt(sin(th*pi/180)) *((31*(1-exp(-th/(1+chi)))*sqrt(Lr*gam)/f^2)-0.45));

Le=Lr*rv;

A01=Le*gam
itu_p=[5 3 2 1 .5 .3 .2 .1 .05 .03 .02 .01 .005 .003 .002 .001];
     
  p=itu_p'; 
for i1=1:length(p)
    pp(i1)=p(i1)
    if (pp(i1)>1) 
        z=0;
    elseif (abs(lat)>36)
        z=0;
    elseif ((abs(lat)<36) & (th>25))
        z=-0.005*(abs(lat)-36);
    else
        z=-0.005*(abs(lat)-36)+1.8-4.25*sin(th*pi/180);
    end
    
    atten(i1)=A01*(pp(i1)./0.01)^-(0.655+0.033*log(pp(i1))-0.045*log(A01)-z*sin(th*pi/180)*(1-pp(i1)));
     
    
end

          
       semilogx(pp,atten,'k<-','LineWidth',1)
        grid on
        title('Rain Attenuation Cumulative Distribution','fontsize',16)
        xlabel('Percentage of time exceeded (%)','fontsize',16)
        ylabel('Rain attenuation(dB)','fontsize',16)