here i need to calculate the BER of UAV based fso system using lognormal model.but after writing the code, i got BER is increasing wrt pt. can anyone suggest me where i got the mistake ?

clear; close all; clc;

%% Define constants and variables
% FSO Parameters
q = 1.602e-19; % Charge of electron
F = 2.75; % Excess noise factor of APD
ra = 10; % Receiver lens radius in cm
wz = 400; % Receiver beam width in cm
sigma_R = 0.2; % Rytov variance for weak turbulence modeled by LN
lambda = 1550e-9; % Optical wavelength
h_planck = 6.626e-34; % Planck’s constant in m^2-kg/s
k_boltz = 1.38064852e-23; % Boltzmann constant in W/K
c = 3e8; % Speed of light
O_f = c / lambda; % Optical frequency
Mew = 0.8; % APD quantum efficiency
G = 30; % Average APD Gain
Be = 1e9; % Electrical bandwidth
Bo = 1e-9; % Optical filter bandwidth
Nb = 1e-3; % Spectral radiance of background radiation
Rl = 1e3; % Receiver Load Resistance
Tr = 300; % Receiver Temperature
Rb = 1; % Data rate in Gbps
Z = 500; % Link length in meters

theta_FOV = 40e-3; % Field of view variations in radians
Pt_dBm = -10:2:50; % Transmit power in dBm
pt = (10^-3) * 10.^(Pt_dBm / 10); % Transmit power in linear scale
eta = (q * G * Mew) / (h_planck * O_f);
Pb = 0.25 * ra^2 * Bo * Nb * theta_FOV^2; % Background power
sigma_th = (4 * k_boltz * Tr * Be) / Rl; % variance of thermal noise
sigma_b = 2 * q * G * F * Be * eta * Pb; % noise variance due to background radiation

% UAV Parameters
sigma_to = 1e-3; % Standard deviation of orientation fluctuation of Tx in radians
sigma_ro = 1e-3; % Standard deviation of orientation fluctuation of Rx in radians
sigma_tp = 40; % Standard deviation of position fluctuation of Tx in centimeters
sigma_rp = 40; % Standard deviation of position fluctuation of Rx in centimeters

% Channel Coefficient Formulas Calculation
hl_dB = 10; % Channel loss in dB
hl = 10^(hl_dB / 10);
v = sqrt(pi / 2) * (ra / wz);
wz_eq = wz^2 * ((sqrt(pi) * erf(v)) / (2 * v * exp(-v^2)));
Ao = (erf(v))^2;
T_base = (wz_eq) / (4 * ((Z * sigma_to)^2 + sigma_tp^2 + sigma_rp^2));
C1 = (T_base / ((Ao * hl)^T_base)) * exp(0.5 * sigma_R * T_base * (1 + T_base));
C2 = T_base – 1;
C3 = 0.5 * sigma_R * (1 + 2 * T_base);

a = [5/24, 4/24, 1/24];
ai_prime = [2, 11/20, 1/2];
M = 2;

% Preallocate Pe array
Pe = zeros(size(Pt_dBm));

% Calculate h1 and b1
h1 = Ao * hl * exp(-C3);

% Loop over different transmit powers
for idx = 1:length(pt)
current_pt = pt(idx);
b1 = eta * current_pt / sqrt(8 * (sigma_th^2 + sigma_b^2));

% Initialize Pe for current pt
Pe(idx) = 0;

% First part of the equation
Pe(idx) = Pe(idx) + 0.5 * exp(-theta_FOV^2 / (2 * (sigma_to^2 + sigma_ro^2)));

% Second part of the equation
sum_term1 = 0;
for i = 1:3
term1 = (a(i) * h1^(C2 + 1)) / 2;
term2 = (1 / (C2 + 1))+sqrt((sigma_R * pi) / ai_prime(i)) * exp((sigma_R * (C2 + 1)^2) / (4 * ai_prime(i))) * …
erf((sqrt(sigma_R) * (C2 + 1)) / (2 * sqrt(ai_prime(i))));
sum_term2 = 0; % Initialize sum_term2 for each i
for m = 0:M
term3 = (2 / sqrt(pi)) * ((-1)^m * (b1 * h1)^(2 * m + 1) / (factorial(m) * (2 * m + 1)));
term4 = (-1) / (C2 + 2 * m + 2);
term5 = sqrt((sigma_R*pi) / ai_prime(i)) * exp((sigma_R * (C2 + 2 * m + 2)^2 )/ (4 * ai_prime(i))) * …
erf((sqrt(sigma_R) * (C2 + 2 * m + 2)) / (2 * sqrt(ai_prime(i))));
sum_term2 = sum_term2 + term3 * (term4 + term5);
end
sum_term1 = sum_term1 + term1 * (term2 + sum_term2);
end
Pe(idx) = Pe(idx) + C1 * (1 – exp(-theta_FOV^2 / (2 * (sigma_to^2 + sigma_ro^2)))) * sum_term1;
end

% Normalize the BER values to decrease from 10^0
Pe = Pe / max(Pe);

% Plot Pt_dBm versus normalized BER
figure;
semilogy(Pt_dBm, Pe, ‘o-‘, ‘LineWidth’, 2);
xlabel(‘Transmit Power (dBm)’);
ylabel(‘Normalized Bit Error Rate (BER)’);
title(‘Transmit Power vs. Normalized BER’);
grid on;clear; close all; clc;

a = [5/24, 4/24, 1/24];
ai_prime = [2, 11/20, 1/2];
M = 2;

% Preallocate Pe array
Pe = zeros(size(Pt_dBm));

% Calculate h1 and b1
h1 = Ao * hl * exp(-C3);

% Loop over different transmit powers
for idx = 1:length(pt)
current_pt = pt(idx);
b1 = eta * current_pt / sqrt(8 * (sigma_th^2 + sigma_b^2));

% Initialize Pe for current pt
Pe(idx) = 0;

% First part of the equation
Pe(idx) = Pe(idx) + 0.5 * exp(-theta_FOV^2 / (2 * (sigma_to^2 + sigma_ro^2)));

% Normalize the BER values to decrease from 10^0
Pe = Pe / max(Pe);

a = [5/24, 4/24, 1/24];
ai_prime = [2, 11/20, 1/2];
M = 2;

% Preallocate Pe array
Pe = zeros(size(Pt_dBm));

% Calculate h1 and b1
h1 = Ao * hl * exp(-C3);

% Loop over different transmit powers
for idx = 1:length(pt)
current_pt = pt(idx);
b1 = eta * current_pt / sqrt(8 * (sigma_th^2 + sigma_b^2));

% Initialize Pe for current pt
Pe(idx) = 0;

% First part of the equation
Pe(idx) = Pe(idx) + 0.5 * exp(-theta_FOV^2 / (2 * (sigma_to^2 + sigma_ro^2)));

% Normalize the BER values to decrease from 10^0
Pe = Pe / max(Pe);

Cart

Cart