difital modulation and BER output the code is normarl ?

Dear Friend i have a code thath ptoduce the BER for 16 qam moduation and send this message in a channell with noise m rician raeleigh finding . I use also a combination of COST 231 MODEL AND AL HURAMI for the attenuation

%% ===================== Παράμετροι =====================

clear; clc; close all;

M = 16;

k = log2(M);

Nsym = 1e5; % Simulation για χαμηλό/μεσαίο SNR

EbN0dB_sim = 0:2:12; % SNR για simulation

EbN0dB_th = 12:2:20; % SNR για θεωρητικό BER

K_Rician = 5;

fc = 2e9;

hBS = 50;

hUAV = 100;

d_mean = 500;

d_sigma = 200;

sigma_PL = 1;

%% ===================== Δημιουργία bits =====================

bits = randi([0 1], Nsym*k, 1);

tx_sym = qammod(bits, M, 'InputType','bit','UnitAveragePower',true);

%% ===================== Simulation =====================

BER_awgn_sim = zeros(size(EbN0dB_sim));

BER_rayleigh_sim = zeros(size(EbN0dB_sim));

BER_rician_sim = zeros(size(EbN0dB_sim));

for ii = 1:length(EbN0dB_sim)

snr_lin = 10^(EbN0dB_sim(ii)/10);

noise_var = 1/(2*k*snr_lin);

% Τυχαίες αποστάσεις UAV

d_horizontal = d_mean + d_sigma*randn(Nsym,1);

d_horizontal(d_horizontal<50) = 50;

d_km = d_horizontal/1000;

% COST-231 path loss

f_MHz = fc/1e6;

a_hUAV = (1.1*log10(f_MHz)-0.7)*hUAV - (1.56*log10(f_MHz)-0.8);

C = 3;

PL_COST231_dB = 46.3 + 33.9*log10(f_MHz) - 13.82*log10(hBS) - a_hUAV + ...

(44.9-6.55*log10(hBS))*log10(d_km) + C;

% L-Hourani correction

FSPL = 20*log10(4*pi*d_horizontal*fc/3e8);

p_LoS = 1./(1 + exp(-0.1*(hUAV - 20)));

eta = p_LoS*1 + (1-p_LoS)*20;

PL_combined_dB = PL_COST231_dB + (FSPL + eta - FSPL);

% Τυχαία path loss

PL_lin_random = 10.^(-PL_combined_dB/10) .* (1 + sigma_PL*randn(Nsym,1));

% --- AWGN ---

rx_awgn = sqrt(PL_lin_random).*tx_sym + sqrt(noise_var)*(randn(Nsym,1)+1i*randn(Nsym,1));

rx_bits_awgn = qamdemod(rx_awgn, M, 'OutputType','bit','UnitAveragePower',true);

BER_awgn_sim(ii) = sum(bits ~= rx_bits_awgn)/length(bits);

% --- Rayleigh ---

h_ray = (randn(Nsym,1)+1i*randn(Nsym,1))/sqrt(2);

rx_ray = sqrt(PL_lin_random).*(h_ray .* tx_sym) + sqrt(noise_var)*(randn(Nsym,1)+1i*randn(Nsym,1));

rx_ray_eq = rx_ray ./ h_ray;

rx_bits_ray = qamdemod(rx_ray_eq, M, 'OutputType','bit','UnitAveragePower',true);

BER_rayleigh_sim(ii) = sum(bits ~= rx_bits_ray)/length(bits);

% --- Rician ---

h_ric = sqrt(K_Rician/(K_Rician+1)) + sqrt(1/(K_Rician+1))*(randn(Nsym,1)+1i*randn(Nsym,1))/sqrt(2);

rx_ric = sqrt(PL_lin_random).*(h_ric .* tx_sym) + sqrt(noise_var)*(randn(Nsym,1)+1i*randn(Nsym,1));

rx_ric_eq = rx_ric ./ h_ric;

rx_bits_ric = qamdemod(rx_ric_eq, M, 'OutputType','bit','UnitAveragePower',true);

BER_rician_sim(ii) = sum(bits ~= rx_bits_ric)/length(bits);

end

%% ===================== Θεωρητικό BER =====================

snr_lin_th = 10.^(EbN0dB_th/10);

% AWGN

BER_awgn_th = (4/k)*(1-1/sqrt(M)) .* qfunc(sqrt(3*k/(M-1)*snr_lin_th));

% Rayleigh

gamma = 3*k/(M-1) * snr_lin_th;

BER_rayleigh_th = 0.5*(1 - sqrt(gamma./(1+gamma)));

% Rician (περίπου)

BER_rician_th = BER_rayleigh_th .* (0.5*(1+K_Rician/(K_Rician+1)));

%% ===================== Συνδυασμός Simulation + Θεωρία (χωρίς διπλό 12 dB) =====================

EbN0dB_th_unique = EbN0dB_th(EbN0dB_th > max(EbN0dB_sim));

EbN0dB_all = [EbN0dB_sim EbN0dB_th_unique];

BER_awgn_all = [BER_awgn_sim BER_awgn_th(EbN0dB_th > max(EbN0dB_sim))];

BER_rayleigh_all = [BER_rayleigh_sim BER_rayleigh_th(EbN0dB_th > max(EbN0dB_sim))];

BER_rician_all = [BER_rician_sim BER_rician_th(EbN0dB_th > max(EbN0dB_sim))];

%% ===================== Smooth interpolation =====================

EbN0dB_fine = 0:0.1:20;

BER_awgn_smooth = interp1(EbN0dB_all, BER_awgn_all, EbN0dB_fine, 'pchip');

BER_rayleigh_smooth = interp1(EbN0dB_all, BER_rayleigh_all, EbN0dB_fine, 'pchip');

BER_rician_smooth = interp1(EbN0dB_all, BER_rician_all, EbN0dB_fine, 'pchip');

%% ===================== Plot =====================

figure;

semilogy(EbN0dB_fine, BER_awgn_smooth, 'k-', 'LineWidth',2); hold on;

semilogy(EbN0dB_fine, BER_rayleigh_smooth, 'b-', 'LineWidth',2);

semilogy(EbN0dB_fine, BER_rician_smooth, 'r-', 'LineWidth',2);

grid on;

xlabel('Eb/N0 [dB]');

ylabel('Bit Error Rate (BER)');

title('16-QAM BER with UAV Motion & COST-231 + L-Hourani Path Loss');

legend('AWGN','Rayleigh','Rician (K=5)');

ylim([1e-6 1]);

this code is ok ??

thank you