fc = 300e9; % Carrier frequency (300 GHz for terahertz)

% Parameters

fc = 300e9; % Carrier frequency (300 GHz for terahertz)

c = 3e8; % Speed of light

lambda = c / fc; % Wavelength

Pt_mean = 1; % Mean transmit power (W)

Pt_std = 0.1; % Standard deviation of transmit power

Gt = 10^(10/10); % Transmit antenna gain (linear scale)

Gr = 10^(10/10); % Receive antenna gain (linear scale)

numElements_mean = 100; % Mean number of elements in RIS

numElements_std = 10; % Standard deviation of number of elements in RIS

RIS_gain = 10^(20/10); % Gain of RIS (linear scale)

SNR_threshold = 10; % SNR threshold for acceptable BER (dB)

latency_per_meter = 5e-9; % Latency per meter (s)

% Constants

k = 1.38e-23; % Boltzmann constant

T = 290; % Noise temperature (K)

B = 1e9; % Bandwidth (1 GHz)

N0 = k * T * B; % Noise power

% Functions for path loss

path_loss_los = @(d) (lambda / (4 * pi * d))^2;

path_loss_ris = @(d1, d2, numElements) (numElements * lambda^2) / (16 * pi^2 * d1 * d2);

% Custom Q-function approximation

Q = @(x) 0.5 * erfc(x / sqrt(2));

% Distances

distances = linspace(10, 500, 50);

% Number of Monte Carlo iterations

num_iterations = 1000;

% Preallocate arrays for results

received_power_los_mc = zeros(length(distances), num_iterations);

received_power_ris_mc = zeros(length(distances), num_iterations);

received_power_los_ris_mc = zeros(length(distances), num_iterations);

data_rate_los_mc = zeros(length(distances), num_iterations);

data_rate_ris_mc = zeros(length(distances), num_iterations);

data_rate_los_ris_mc = zeros(length(distances), num_iterations);

ber_los_mc = zeros(length(distances), num_iterations);

ber_ris_mc = zeros(length(distances), num_iterations);

ber_los_ris_mc = zeros(length(distances), num_iterations);

latency_los_mc = zeros(length(distances), num_iterations);

latency_ris_mc = zeros(length(distances), num_iterations);

latency_los_ris_mc = zeros(length(distances), num_iterations);

for iter = 1:num_iterations

% Randomly vary the transmit power and number of RIS elements

Pt = normrnd(Pt_mean, Pt_std);

numElements = round(normrnd(numElements_mean, numElements_std));

for i = 1:length(distances)

d = distances(i);

d1 = d / 2; % Distance from Tx to RIS

d2 = d / 2; % Distance from RIS to Rx

% LoS only

Pr_los = Pt * Gt * Gr * path_loss_los(d);

SNR_los = Pr_los / N0;

data_rate_los_mc(i, iter) = B * log2(1 + SNR_los); % Shannon capacity

ber_los_mc(i, iter) = Q(sqrt(2 * SNR_los)); % BER for BPSK using custom Q-function

received_power_los_mc(i, iter) = Pr_los;

latency_los_mc(i, iter) = d / c + latency_per_meter * d;

% RIS only

Pr_ris = Pt * Gt * Gr * RIS_gain * path_loss_ris(d1, d2, numElements);

SNR_ris = Pr_ris / N0;

data_rate_ris_mc(i, iter) = B * log2(1 + SNR_ris); % Shannon capacity

ber_ris_mc(i, iter) = Q(sqrt(2 * SNR_ris)); % BER for BPSK using custom Q-function

received_power_ris_mc(i, iter) = Pr_ris;

latency_ris_mc(i, iter) = (d1 + d2) / c + latency_per_meter * (d1 + d2);

% LoS + RIS

Pr_los_ris = Pr_los + Pr_ris; % Sum of received powers

SNR_los_ris = Pr_los_ris / N0;

data_rate_los_ris_mc(i, iter) = B * log2(1 + SNR_los_ris); % Shannon capacity

ber_los_ris_mc(i, iter) = Q(sqrt(2 * SNR_los_ris)); % BER for BPSK using custom Q-function

received_power_los_ris_mc(i, iter) = Pr_los_ris;

latency_los_ris_mc(i, iter) = d / c + latency_per_meter * d;

end

end

% Calculate mean and standard deviation for each parameter

mean_received_power_los = mean(received_power_los_mc, 2);

std_received_power_los = std(received_power_los_mc, 0, 2);

mean_received_power_ris = mean(received_power_ris_mc, 2);

std_received_power_ris = std(received_power_ris_mc, 0, 2);

mean_received_power_los_ris = mean(received_power_los_ris_mc, 2);

std_received_power_los_ris = std(received_power_los_ris_mc, 0, 2);

mean_data_rate_los = mean(data_rate_los_mc, 2);

std_data_rate_los = std(data_rate_los_mc, 0, 2);

mean_data_rate_ris = mean(data_rate_ris_mc, 2);

std_data_rate_ris = std(data_rate_ris_mc, 0, 2);

mean_data_rate_los_ris = mean(data_rate_los_ris_mc, 2);

std_data_rate_los_ris = std(data_rate_los_ris_mc, 0, 2);

mean_ber_los = mean(ber_los_mc, 2);

std_ber_los = std(ber_los_mc, 0, 2);

mean_ber_ris = mean(ber_ris_mc, 2);

std_ber_ris = std(ber_ris_mc, 0, 2);

mean_ber_los_ris = mean(ber_los_ris_mc, 2);

std_ber_los_ris = std(ber_los_ris_mc, 0, 2);

mean_latency_los = mean(latency_los_mc, 2);

std_latency_los = std(latency_los_mc, 0, 2);

mean_latency_ris = mean(latency_ris_mc, 2);

std_latency_ris = std(latency_ris_mc, 0, 2);

mean_latency_los_ris = mean(latency_los_ris_mc, 2);

std_latency_los_ris = std(latency_los_ris_mc, 0, 2);

% Plotting

figure;

subplot(3, 1, 1);

errorbar(distances, mean_received_power_los, std_received_power_los, 'r', 'LineWidth', 2);

hold on;

errorbar(distances, mean_received_power_ris, std_received_power_ris, 'g', 'LineWidth', 2);

errorbar(distances, mean_received_power_los_ris, std_received_power_los_ris, 'b', 'LineWidth', 2);

xlabel('Distance (m)');

ylabel('Received Power (dBm)');

legend('LoS Only', 'RIS Only', 'LoS + RIS');

title('Received Power vs Distance with Monte Carlo Simulations');

subplot(3, 1, 2);

errorbar(distances, mean_data_rate_los / 1e6, std_data_rate_los / 1e6, 'r', 'LineWidth', 2);

hold on;

errorbar(distances, mean_data_rate_ris / 1e6, std_data_rate_ris / 1e6, 'g', 'LineWidth', 2);

errorbar(distances, mean_data_rate_los_ris / 1e6, std_data_rate_los_ris / 1e6, 'b', 'LineWidth', 2);

xlabel('Distance (m)');

ylabel('Data Rate (Mbps)');

legend('LoS Only', 'RIS Only', 'LoS + RIS');

title('Data Rate vs Distance with