# Any demonstration of ITU P.618-9(Propagation data and prediction methods required for the design of Earth-space telecommunication systems)?

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*)

Clear["Global*"];
\[CurlyPhi] = 56.5;(* the lattitute of the Earth station \[CurlyPhi]，in degree*)
Subscript[h, "R"] =
Module[{},
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"] =
Module[{},
If[\[Theta] >= 5, (Subscript[h, "R"] - Subscript[h, "S"])/
Sin[\[Theta]],
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"] +
Subscript[\[Kappa],
"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];

(*\[Kappa]\[Alpha]\
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]*
Subscript[R,
"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"] =
Module[{},
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"] =
0.4095*f*L/(Subscript[\[CurlyEpsilon],
"''"]*(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*)

`