# Error in the visualizing Capacity bounds of FSO Links over Gamma-Gamma Atmospheric Turbulence Channels Using OOK Signaling?

I simply want to visualize 【Capacity bounds and mutual information】 versus 【SNR $$\gamma$$(dB) or probability $$p$$ (The input to the channel is a binomial distribution with parameter $$p$$)】

But I cannot reproduce the results in a paper (On the Capacity of FSO Links over Gamma-Gamma Atmospheric Turbulence Channels Using OOK Signaling) for some reason.

Is there some mistake in my Mathematica code?

\[Gamma] = 10;(*\[Gamma]=P_{opt}/\sqrt{N_o W} is the SNR \
definition*)Clear["Global*"];
p = 1/2;
\[Xi] = 1;
(*\[Xi] represents the square of the increment in Euclidean distance \
due to the use of a pulse shape of high PAOPR(peak-to-average optical \
power ratio),alternative to the classical rectangular pulse.*)
\[Theta] = 0.9028;
\[Kappa] = 2;(*with 0<\[Theta]<1,representing the fact that the \
channel under study is constrained to \[Kappa]=2 W T_b degrees of \
freedom*)\[Alpha] = 4;
\[CapitalTheta] = \[Alpha]^(3 \[Alpha]/
2)*(1 - \[Alpha])^((1 - \[Alpha])/
2)*(1 - \[Alpha]/2)^(1 - \[Alpha]/2);
\[Beta] = 4;
(*\[Alpha] and \[Beta] are related to the atmospheric conditions*)
SI = 1/\[Alpha] + 1/\[Beta] +
1/(\[Alpha]*\[Beta]);(*scintillation index,a parameter of interest \
used to describe the strength of atmospheric fading*)Capacity =
MeijerG[{{1,
1, (1 - \[Alpha])/2, (2 - \[Alpha])/2, (1 - \[Beta])/
2, (2 - \[Beta])/2}, {}}, {{1}, {0}},
16*(1/p -
1)*\[Xi]*\[Theta]*\[Kappa]*\[Gamma]^2/(\[Alpha]^2*\[Beta]^2)]*(\
2^(\[Alpha] + \[Beta])/
Log[E, 2]/(2*Pi*Gamma[\[Alpha]]*Gamma[\[Beta]])) /. {\[Gamma] ->
10*10^\[Gamma]};
Subscript[C, "nt"] =
1/2*Log[2,
1 + (1/p - 1)*\[Kappa]*\[Xi]*\[Theta]*\[Gamma]^2] /. {\[Gamma] ->
10*10^\[Gamma]};
(*nt means:no atmospheric turbulence*)
Subscript[C, "AWGN"] =
1/2*Log[2, 1 + 2*\[Xi]*\[Theta]*\[Kappa]*\[Gamma]^2] /. {\[Gamma] ->
10*10^\[Gamma]};
Subscript[C, "H"] =
Log[2, (Sqrt[\[Xi]*\[Theta]*\[Kappa]]*\[Gamma] + 2)*
Sqrt[E/(2 Pi)]] /. {\[Gamma] -> 10*10^\[Gamma]};
Subscript[C, "F"] =
Log[2, Sqrt[
E^2/(4 Pi)]^\[Alpha]*(Sqrt[\[Xi]*\[Theta]*\[Kappa]]*\[Gamma])^\
\[Alpha]*1/\[CapitalTheta]] /. {\[Gamma] -> 10*10^\[Gamma]};
Subsuperscript[C, "AWGN", "turb"] =
MeijerG[{{1,
1, (1 - \[Alpha])/2, (2 - \[Alpha])/2, (1 - \[Beta])/
2, (2 - \[Beta])/2}, {}}, {{1}, {0}},
32*\[Xi]*\[Theta]*\[Kappa]*\[Gamma]^2/(\[Alpha]^2*\[Beta]^2)]*(2^(\
\[Alpha] + \[Beta] - 2)/
Log[E, 2]/(2*Pi*Gamma[\[Alpha]]*Gamma[\[Beta]])) /. {\[Gamma] ->
10*10^\[Gamma]}; plot1 =
Plot[{Capacity, Subscript[C, "nt"], Subscript[C, "AWGN"],
Subscript[C, "H"], Subscript[C, "F"],
Subsuperscript[C, "AWGN", "turb"]}, {\[Gamma], 0.01, 50},
PlotLegends -> "Expressions", PlotRange -> Full,
AxesLabel -> {"SNR \[Gamma](in dB)", "Capacity(bits/channel use)"}];
Subscript[f, "I"] =
2*(\[Alpha]*\[Beta])^((\[Alpha] + \[Beta])/2)/(Gamma[\[Alpha]]*
Gamma[\[Beta]])*i^((\[Alpha] + \[Beta])/2 - 1)*
BesselK[\[Alpha] - \[Beta], 2 Sqrt[\[Alpha]*\[Beta]*i]];
CGammap[\[Gamma]_] :=
NIntegrate[
1/2*Log[2, (1 + (1/p - 1)*\[Kappa]*\[Xi]*\[Theta]*\[Gamma]^2*i^2)]*
Subscript[f, "I"] /. {\[Gamma] -> 10*10^\[Gamma]}, {i, 0,
Infinity}];
plot2 = ListPlot[{#, CGammap@#} & /@ Range[5, 50, 0.1],
PlotLegends -> "anotherBound"];
Show[plot1, plot2]


Clear["Global*"];
\[Xi] = 1;
\[Theta] = 0.9028;
\[Kappa] = 2;

fYwhenxequals1[y_, p_, i_, \[Gamma]_] :=
1/Sqrt[2*Pi]*
Exp[-(y - 1/p*\[Gamma]*Sqrt[\[Xi]*\[Theta]*\[Kappa]]*i^2)^2/2];
fYwhenxequals0[y_] := 1/Sqrt[2*Pi]*Exp[-y^2/2];

MutualInformation[i_, p_, \[Gamma]_] :=
Module[{denominator},
denominator[y_] := (1 - p)*fYwhenxequals0[y] +
p*fYwhenxequals1[y, p, i, \[Gamma]];
(1 - p)*
NIntegrate[
fYwhenxequals0[y]*
Log[2, fYwhenxequals0[y]/denominator[y]], {y, -65, 109}] +
p*NIntegrate[
fYwhenxequals1[y, p, i, \[Gamma]]*
Log[2, fYwhenxequals1[y, p, i, \[Gamma]]/
denominator[y]], {y, -65, 109}]];
listplotWhenSNRequalsMinus5dB =
ListPlot[{#, MutualInformation[1, #, 10^(-5/10)]} & /@
Range[0, 1, 0.00625]];
listplotWhenSNRequalsMinus10dB =
ListPlot[{#, MutualInformation[1, #, 10^(-10/10)]} & /@
Range[0, 1, 0.00625]];
listplotWhenSNRequalsMinus1dB =
ListPlot[{#, MutualInformation[1, #, 10^(-1/10)]} & /@
Range[0, 1, 0.00625]];
Show[listplotWhenSNRequalsMinus5dB, listplotWhenSNRequalsMinus10dB, \
listplotWhenSNRequalsMinus15dB,
PlotLegends -> {"SNR is -5 dB", "SNR is -10 dB", "SNR is -1 dB"}]


Clear["Global*"];
2*(\[Alpha]*\[Beta])^((\[Alpha] + \[Beta])/2)/(Gamma[\[Alpha]]*
Gamma[\[Beta]])*i^((\[Alpha] + \[Beta])/2 - 1)*
BesselK[\[Alpha] - \[Beta], 2 Sqrt[\[Alpha]*\[Beta]*i]];
GammaGammaDistributionPDF[1, 1]}, {i, 0.001, 10}, PlotRange -> Full,
PlotLabel ->
Row[{"Gamma-Gamma distribution for different \[Alpha] and \[Beta]",
""}], AxesLabel -> {"i", Row[{Subscript["f", "I"], "(i)"}]},
PlotLegends -> "Expressions"]

(*For plane wave,here sigmaR^2 is the variance of laser \
scintillation (Rytov variance?),and its magnitude determines the strength of \
atmospheric turbulence.The larger its value,the stronger the \
turbulence intensity.*)
1/(Exp[(0.49 sigmaR^2)/((1 + 1.11 sigmaR^(12/5))^(7/6))] - 1),
1/(Exp[(0.51 sigmaR^2)/((1 + 0.69 sigmaR^(12/5))^(5/6))] - 1)]
PlotRange -> Full,
PlotLabel ->
Row[{"Gamma-Gamma distribution for different ",
Subsuperscript[\[Sigma], "R", ""], ""}],
AxesLabel -> {"i", Row[{Subscript["f", "I"], "(i)"}]},
PlotLegends -> "Expressions"]


By the way,

LogNormal model for weak turbulence

LogNormalDistributionForWeakTurbulence[MeanOfLogorithmLightIntensity_,
StandardDeviationOfLogorithmLightIntensity_] :=
Module[{MeanOfChi, VarianceOfChi, StandardDeviationOfChi},
MeanOfChi = MeanOfLogorithmLightIntensity/2;
VarianceOfChi = StandardDeviationOfLogorithmLightIntensity^2/4;
StandardDeviationOfChi = Sqrt[VarianceOfChi];
PDF[LogNormalDistribution[MeanOfChi, StandardDeviationOfChi], i]]

Plot[Table[
LogNormalDistributionForWeakTurbulence[0.05,
StandardDeviationOfLogorithmLightIntensity], \
{StandardDeviationOfLogorithmLightIntensity, {0.2, 0.4, 0.6, 0.8} //
Sqrt}] // Evaluate, {i, 0, 4}, Filling -> Axis,
PlotRange -> Full, PlotLegends -> {0.2, 0.4, 0.6, 0.8},
PlotLabel ->
Row[{"LogNormal distribution for different ",
Subsuperscript[\[Sigma], "l", "2"], ""}]]
`