I'am looking for demonstration of ITU P.618-9(Propagation data and prediction methods required for the design of Earth-space telecommunication systems) P.618-9 (Year 2007 version), in MATLAB or Mathematica.

I have written some Mathematica code, any help would be appreciated.

(*ITU-R rain attenuation prediction model, ITU P.618-9 version, year 2007*)

\[CurlyPhi] = 56.5;(* the lattitute of the Earth station \[CurlyPhi],in degree*)
Subscript[h, "R"] = 
   If[0 <= \[CurlyPhi] && \[CurlyPhi] < 36, 3.0 + 0.028 \[CurlyPhi] , 
    If[\[CurlyPhi] >= 36, 
     4.0 - 0.075 (\[CurlyPhi] - 36)]]];(*Calculate equivalent rain roof heights using empirical formulas h_R*)
Subscript[h, "S"] = 2000;(*Average elevation of Earth stations h_S,in km*)
Subscript[R, "e"] = 8500;(*Equivalent Earth radius taking into account refraction, R_e,in km*)
\[Theta] = 4.5;(*Earth station antenna elevation angle \[Theta],in degree*)
Subscript[L, "S"] = 
   If[\[Theta] >= 5, (Subscript[h, "R"] - Subscript[h, "S"])/
    If[\[Theta] < 5, 
     2 (Subscript[h, "R"] - 
         Subscript[h, "S"])/(Sqrt[
          Sin[\[Theta]]^2 + 
           2 (Subscript[h, "R"] - Subscript[h, "S"])/
             Subscript[R, "e"]] + 
         Sin[\[Theta]])]]];(*Length of inclined path below rainfall height L_S,in km*)
Subscript[L, "G"] = 
  Subscript[L, "S"] Cos[\[Theta]]; (*Horizontal projection of an inclined path L_G,in km*)

\[Tau] = 32; (*Polarization angle of the received point wave*)
Subscript[\[Kappa], "H"] = 31231;
Subscript[\[Kappa], "V"] = 31546;
Subscript[\[Alpha], "H"] = 86935;
Subscript[\[Alpha], "V"] = 464534;
\[Kappa] = (Subscript[\[Kappa], "H"] + 
     "V"] + (Subscript[\[Kappa], "H"] - Subscript[\[Kappa], "V"])*
     Cos[\[Theta]]^2*Cos[2 \[Tau]])/2
\[Alpha] = (Subscript[\[Kappa], "H"] Subscript[\[Alpha], "H"] + 
      Subscript[\[Kappa], "V"] Subscript[\[Alpha], 
        "V"] + (Subscript[\[Kappa], "H"] Subscript[\[Alpha], "H"] - 
         Subscript[\[Kappa], "V"] Subscript[\[Alpha], "V"])*
       Cos[\[Theta]]^2*Cos[2 \[Tau]])/2 \[Kappa];

are coefficients related to frequency, polarization mode, raindrop size, rain dielectric constant, etc., \[Kappa]_H, \[Alpha]_H and \[Kappa]_V, \[Alpha]_\
v correspond to the horizontally polarized and vertically polarized components of the wave, respectively*)

Subscript[R, "0.01"] = 
  1234.5;(*Local rainfall rate in mm/h at the point of time probability of exceeding 0.01% of the annual average, in 1-minute integration time*)
Subscript[Gamma, "R"] = \[Kappa]*
     "0.01"]^\[Alpha];(*Characteristic attenuation of rainfall rate at 0.01% annual average time probability point\gamma_R in dB/km*)

Subscript[r, "0.01"] = 
   Subscript[L, "0"] = 35*Exp[-0.015*Subscript[R, "0.01"]](*km*); 
   1/(1 + Subscript[L, "G"]/Subscript[L, "0"])];
(*Path shortening factor for annual average 0.01% time probability points r_{0.01}*)

Subscript[A, "0.01"] = 
  Subscript[\[Gamma], "R"]*Subscript[L, "S"]*Subscript[r, "0.01"];
(* Rainfall attenuation value A_{0.01} at the point of 0.01% annual average time probability, in dB  *)
ApFunction[p_] := 
  Subscript[A, "0.01"]*0.12*p^(-(0.546 + 0.043 Log[10, p]));
ApFunction[A, "p"][0.01];  
(*Calculate the attenuation value A_p(dB) from A_{0.01} for any time probability p in the annual averaging time range of 0.001%~1%, here we take 
0.01% to see the result.*)
ApFunction[A, "p"][0.02];  (*Take 0.02% and look at the results*)

A = 123.2;(*Rainfall attenuation, in dB*)
Subscript[T, "r"] = 273.15 + 23;(* Rainfall temperature in Kelvin *) 
Subscript[DeltaT, "r"] = (1 - 10^(-A/10))*Subscript[T, "r"];
(* Equivalent noise temperature caused by rainfall\Delta T_r (in Kelvin) *)

(* Cross-polarization discrimination XPD is defined as the ratio of the main polarization component of this channel to the cross-polarization component generated by another signal within this channel, and it exists only in dual-polarization systems. Since current satellite communications are generally \
As the current satellite communications are generally dual-polarization system, and XPD can directly reflect the degree of interference of other channels to this channel, so generally use XPD to measure polarization interference. *)
Subscript[A, "p"] = 0.12;(*attenuation of rainfall A_p*)
Subscript[C, "f"] = 0.34;(*frequency factor C_f*)
Subscript[C, \[Tau]] = 0.56; (*Line polarization improvement factor C_\tau*)
Subscript[C, \[Theta]] = 0.78;(*Geographic gain factor C_\theta*)
Subscript[C, \[Delta]] = 0.89;(*raindrop inclination factor C_\delta*)
f = 9;(*carrier frequency f,in GHz*)
XPD = Module[{U, V}, 
   U = Subscript[C, "f"] + Subscript[C, \[Tau]] + 
     Subscript[C, \[Theta]] + Subscript[C, \[Delta]]; 
   V = If[8 <= f && f < 20, 12.8*f^0.19, If[20 <= f && f < 35, 22.6]];
   U - V*Log[10, Subscript[A, "p"]]
   ];(*Cross-polarization discrimination XPD*)

Subscript[\[Gamma], "w"] = 1.23;(*Water vapor loss factor \gamma_w,in dB/km*)
Subscript[h, "w"] = 4.56;(*Effective height of water vapor h_w,in km*)

Subscript[A, "w"] = 
  Module[{F}, F[x_] = 1/(0.661 x + 0.339 Sqrt[x^2 + 5.31]); 
   If[\[Theta] > 10, 
    Subscript[h, "w"]*Subscript[\[Gamma], "w"]/Sin[\[Theta]], 
    If[\[Theta] <= 
      10, (Subscript[\[Gamma], "w"]*
        Sqrt[Subscript[R, "e"]*Subscript[h, "w"]]/Cos[\[Theta]])*
      F[Tan[\[Theta]]*Sqrt[Subscript[R, "e"]/Subscript[h, "w"]]]]]];
(*Absorption loss of water vapor A_w,in dB*)

Subscript[\[Gamma], "o"] = 3.21;(*Oxygen loss factor \gamma_o,in dB/km*)
Subscript[h, "o"] = 6.54;(*Effective height of oxygen h_o,in km*)

Subscript[A, "o"] = 
 Module[{F}, F[x_] = 1/(0.661 x + 0.339 Sqrt[x^2 + 5.31]); 
  If[\[Theta] > 10, 
   Subscript[h, "o"]*Subscript[\[Gamma], "o"]/Sin[\[Theta]], 
   If[\[Theta] <= 
     10, (Subscript[\[Gamma], "o"]*
       Sqrt[Subscript[R, "e"]*Subscript[h, "o"]]/Cos[\[Theta]])*
     F[Tan[\[Theta]]*Sqrt[Subscript[R, "e"]/Subscript[h, "o"]]]]]];
(*Oxygen uptake losses A_o,in dB*)

L = 1.2;(*Cloud thickness L,in km*)

(*use this https://www.random-science-tools.com/electronics/water_dielectric.htm to calculate the DIELECTRIC CONSTANT OF WATER*)

Subscript[\[CurlyEpsilon], "'"] = 46.43;(*The real part of the dielectric constant of water*)
Subscript[\[CurlyEpsilon], "''"] = 41.15;(*The imaginary part of the dielectric constant of water*)
\[Eta] = (2 + 
    Subscript[\[CurlyEpsilon], "'"])/(Subscript[\[CurlyEpsilon], 
Subscript[A, "c"] = 
      "''"]*(1 + \[Eta]^2) Sin[\[Theta]]); (*cloud caused attenuation A_c,in dB*)
Subscript[A, "t"] = 
 Subscript[A, "o"] + 
  Sqrt[(Subscript[A, "w"] + Subscript[A, "c"] + Subscript[A, "s"])^2 +
     Subscript[A, "r"]^2];(*total decay value A_t*)



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.