carrier_frequency = 2.4e9;
wifi_be_subcarrier_spacing = 78.125e3;
leo_altitude = 600e3;
bandwidth = 20; 
mcs_values = [2];
apep = 1000;
tx_antenna_count = 1;
rx_antenna_count = 1;

snrs = (0:1:40);
packets_per_snr = 250;
maxNumErrors = 50;

% elevation_angle = 50;
elevationAngles = [ ...
    linspace(10, 90, ceil(packets_per_snr / 2)), ...
    linspace(90, 10, floor(packets_per_snr / 2)) ...
];

meanElevationAngle = mean(elevationAngles);
slantRanges = arrayfun(@(el) slantRangeCircularOrbit(el, leo_altitude, 0), elevationAngles);
dopplerShifts = arrayfun(@(el) dopplerShiftCircularOrbit(el, leo_altitude, 0, carrier_frequency), elevationAngles);

random_stream = "mt19937ar with seed";
seed = 666;

ntn_tdl_delay_profiles = ["A" "B" "C" "D"];

[cfgEHT, simParameters, channel, chInfo, maxChDelay, satelliteDopplerShift, chanBW] = ...
    create_baseline_configuration(ntn_tdl_delay_profiles(3), bandwidth, false, "None", ...
                       carrier_frequency, 0, ...
                       0, leo_altitude, random_stream, ...
                       seed, tx_antenna_count, ...
                       rx_antenna_count, apep);

meanSlantRange = mean(slantRanges);
meanSatelliteDopplerShift = mean(dopplerShifts);

meanNormalizedCFO = meanSatelliteDopplerShift / wifi_be_subcarrier_spacing;

disp("Mean Elevation Angle: " + meanElevationAngle + "°")
disp("Mean Slant Range: " + meanSlantRange + "m")
disp("Mean Satellite Doppler Shift: " + meanSatelliteDopplerShift + "Hz")
disp("Normalized CFO: " + meanNormalizedCFO)

errorRates = zeros(numel(mcs_values), numel(snrs));
for imcs = 1 : numel(mcs_values)
    cfgEHT.User{1}.MCS = mcs_values(imcs);
    parfor snr_idx = 1 : length(snrs)
        local_channel = clone(channel)
        [errorCount, packetCount] = simulateTransmission(cfgEHT, local_channel, maxChDelay, snrs(snr_idx), packets_per_snr, maxNumErrors, simParameters, chanBW, snr_idx, elevationAngles);
        
        errorCounts(snr_idx) = errorCount;
        packetCounts(snr_idx) = packetCount;
        errorRates(imcs, snr_idx) = errorCount / packetCount;
    end
end

markers = 'ox*sd^v><ph+ox*sd^v><ph+';
figure;
for imcs = 1 : numel(mcs_values) 
    semilogy(snrs, errorRates(imcs, :).',['-' markers(imcs)]);
    hold on;
end
grid on;
xlabel('SNR (dB)');
ylabel('PER');
% task 1
dataStr = arrayfun(@(x)sprintf('MCS %d',x), 2, 'UniformOutput', false);
% task 4
% dataStr = arrayfun(@(x)'MCS Sweep Values: [0 2 4 8 10 12 13]', 'UniformOutput', false);
legend(dataStr, 'Location', 'NorthEastOutside');
title(['PER (EHT MU), ' num2str(cfgEHT.ChannelBandwidth) ', NTN-TDL-C, ' num2str(tx_antenna_count) '-by-' num2str(rx_antenna_count)]);