minutesToSimulate = 10;
satConfig.stopTime = satConfig.startTime + minutes(minutesToSimulate);

snrs = 0:1:40;

if exist("satSim", "var")
    clear satSim;
end

satSim = create_satellite_sim(satConfig, groundStation);
packetsToTransmit = 1000;

sampleTimes = satConfig.startTime:seconds(satConfig.sampleTime):satConfig.stopTime;
packetsPerTimeSlot = packetsToTransmit / numel(sampleTimes);

packetTimes = satConfig.startTime + seconds(linspace(0, minutesToSimulate * 60, packetsToTransmit));

% viewer = satelliteScenarioViewer(satSim.sc);

preCompensations = {"none", "ideal", 10, 100, 1000, 5000};
preCompensationLabels = [
    "none"
    "ideal"
    "10Hz"
    "100Hz"
    "1000Hz"
    "5000Hz"
];

markers = 'ox*sd^v><ph+ox*sd^v><ph+';

if calc_task4
    distances = zeros(1, packetsToTransmit);
    elevationAngles = zeros(1, packetsToTransmit);
    dopplerShifts = zeros(1, packetsToTransmit);

    progressStartText = "Starting packet parameter calculation ...";
    f = waitbar(0, progressStartText, "Name", "Task 4 packet times");

    progressTotal = numel(packetTimes);
    progressQueue = parallel.pool.DataQueue;

    parfor_waitbar(f, progressTotal, progressStartText, true);
    progressListener = afterEach(progressQueue, @(info) parfor_waitbar(f, progressTotal, info, false));

    parfor timeIdx = 1 : numel(packetTimes)
        [~, elevationAngles(timeIdx), distances(timeIdx)] = aer(satSim.ue, satSim.sat, packetTimes(timeIdx));
        dopplerShifts(timeIdx) = dopplershift(satSim.sat, satSim.ue, packetTimes(timeIdx), Frequency=carrier_frequency);
        send(progressQueue, sprintf("Packet: %d", timeIdx));

        %{
        res = perform_access_analysis(satSim, packetTimes(timeIdx));
        
        distances(timeIdx) = res.distanceToSatellite;
        elevationAngles(timeIdx) = res.elevationAngleToSatellite;
        dopplerShifts(timeIdx) = res.fshift;
        send(progressQueue, sprintf("Packet: %d", timeIdx));
        %}
    end

    delete(progressListener);
    if ishandle(f)
        close(f);
    end

    accessAnalysisParameters.times = packetTimes;
    accessAnalysisParameters.distances = distances;
    accessAnalysisParameters.elevationAngles = elevationAngles;
    accessAnalysisParameters.dopplerShifts = dopplerShifts;
    accessAnalysisParameters.enablePilotTracking = false;

    simParameters.IncludeFreeSpacePathLoss = false;

    progressStartText = "Starting PER sweep ...";
    f = waitbar(0, progressStartText, "Name", "Task 4");

    progressTotal = numel(preCompensations) * numel(snrs);
    progressQueue = parallel.pool.DataQueue;

    parfor_waitbar(f, progressTotal, progressStartText, true);
    progressListener = afterEach(progressQueue, @(info) parfor_waitbar(f, progressTotal, info, false));

    errorCounts = zeros(numel(preCompensations), numel(snrs));
    packetCounts = zeros(numel(preCompensations), numel(snrs));
    errorRates = zeros(numel(preCompensations), numel(snrs));
    cfrMagnitudesBySnr = cell(numel(preCompensations), numel(snrs));

    for preCompensationIdx = 1 : numel(preCompensations)
        scenarioParameters = accessAnalysisParameters;
        scenarioParameters.preCompensation = preCompensations{preCompensationIdx};
        scenarioLabel = preCompensationLabels(preCompensationIdx);

        errorCountsForScenario = zeros(1, numel(snrs));
        packetCountsForScenario = zeros(1, numel(snrs));
        errorRatesForScenario = zeros(1, numel(snrs));
        cfrMagnitudesForScenario = cell(1, numel(snrs));

        parfor snrIdx = 1 : numel(snrs)
            localChannel = clone(channel);
            [numPacketErrors, numPacketsSimulated, cfrMagnitudes, ~, ~] = simulateTransmission(scenarioParameters, cfgEHT, localChannel, maxChDelay, snrs(snrIdx), inf, simParameters);
            errorCountsForScenario(snrIdx) = numPacketErrors;
            packetCountsForScenario(snrIdx) = numPacketsSimulated;
            errorRatesForScenario(snrIdx) = numPacketErrors / numPacketsSimulated;
            cfrMagnitudesForScenario{snrIdx} = cfrMagnitudes;
            send(progressQueue, sprintf("Pre-Compensation: %s, SNR: %d", char(scenarioLabel), snrs(snrIdx)));
        end

        errorCounts(preCompensationIdx, :) = errorCountsForScenario;
        packetCounts(preCompensationIdx, :) = packetCountsForScenario;
        errorRates(preCompensationIdx, :) = errorRatesForScenario;
        cfrMagnitudesBySnr(preCompensationIdx, :) = cfrMagnitudesForScenario;
    end

    delete(progressListener);
    if ishandle(f)
        close(f);
    end

    plotErrorRates = max(errorRates, 1 / packetsToTransmit);

    figure;
    hold on;
    for preCompensationIdx = 1 : numel(preCompensations)
        semilogy(snrs, plotErrorRates(preCompensationIdx, :), ['-' markers(preCompensationIdx)]);
    end

    grid on;
    xlabel("SNR [dB]");
    ylabel("PER");
    ylim([1 / packetsToTransmit 1]);
    legend(preCompensationLabels, "Location", "NorthEastOutside");
    title("PER over SNR with Doppler pre-compensation");

    % greatest difference between no pre compensation and ideal
    max(abs(errorRates(1,:) - errorRates(2,:)))
end % calc_task4


% 4.4
carrierFrequenciesF = [2.4e9, 5e9, 6e9, 10e9, 20e9, 40e9, 45e9, 50e9, 55e9 60e9];
preCompensationsF = {"none", "ideal"};
preCompensationLabelsF = [
    "none"
    "ideal"
];

if calc_task4_4
    errorCountsF = zeros(numel(carrierFrequenciesF), numel(preCompensationsF), numel(snrs));
    packetCountsF = zeros(numel(carrierFrequenciesF), numel(preCompensationsF), numel(snrs));
    errorRatesF = zeros(numel(carrierFrequenciesF), numel(preCompensationsF), numel(snrs));
    cfrMagnitudesBySnrF = cell(numel(carrierFrequenciesF), numel(preCompensationsF), numel(snrs));

    for carrierFrequencyIdxF = 1 : numel(carrierFrequenciesF)
        carrierFrequencyF = carrierFrequenciesF(carrierFrequencyIdxF);

        [cfgEHTF, simParametersF, channelF, ~, maxChDelayF, ~, ~] = create_baseline_configuration(profile, bandwidth, false, delay_model, carrierFrequencyF, mobile_speed, mobile_altitude, satellite_altitude, random_stream, seed, tx_antenna_count, rx_antenna_count, apep);
        cfgEHTF.User{1}.MCS = 3;
        simParametersF.TransmitPower = txPower;
        simParametersF.IncludeFreeSpacePathLoss = false;

        distancesF = zeros(1, packetsToTransmit);
        elevationAnglesF = zeros(1, packetsToTransmit);
        dopplerShiftsF = zeros(1, packetsToTransmit);

        progressStartTextF = sprintf("Starting packet parameter calculation for %.1f GHz ...", carrierFrequencyF / 1e9);
        waitbarF = waitbar(0, progressStartTextF, "Name", "Task 4.4 packet times");

        progressTotalF = numel(packetTimes);
        progressQueueF = parallel.pool.DataQueue;

        parfor_waitbar(waitbarF, progressTotalF, progressStartTextF, true);
        progressListenerF = afterEach(progressQueueF, @(info) parfor_waitbar(waitbarF, progressTotalF, info, false));

        parfor timeIdxF = 1 : numel(packetTimes)
            [~, elevationAnglesF(timeIdxF), distancesF(timeIdxF)] = aer(satSim.ue, satSim.sat, packetTimes(timeIdxF));
            dopplerShiftsF(timeIdxF) = dopplershift(satSim.sat, satSim.ue, packetTimes(timeIdxF), Frequency=carrierFrequencyF);
            send(progressQueueF, sprintf("Packet: %d", timeIdxF));
        end

        delete(progressListenerF);
        if ishandle(waitbarF)
            close(waitbarF);
        end

        accessAnalysisParametersF.times = packetTimes;
        accessAnalysisParametersF.distances = distancesF;
        accessAnalysisParametersF.elevationAngles = elevationAnglesF;
        accessAnalysisParametersF.dopplerShifts = dopplerShiftsF;
        accessAnalysisParametersF.enablePilotTracking = false;

        progressStartTextF = sprintf("Starting PER sweep for %.1f GHz ...", carrierFrequencyF / 1e9);
        waitbarF = waitbar(0, progressStartTextF, "Name", "Task 4.4");

        progressTotalF = numel(preCompensationsF) * numel(snrs);
        progressQueueF = parallel.pool.DataQueue;

        parfor_waitbar(waitbarF, progressTotalF, progressStartTextF, true);
        progressListenerF = afterEach(progressQueueF, @(info) parfor_waitbar(waitbarF, progressTotalF, info, false));

        for preCompensationIdxF = 1 : numel(preCompensationsF)
            scenarioParametersF = accessAnalysisParametersF;
            scenarioParametersF.preCompensation = preCompensationsF{preCompensationIdxF};
            scenarioLabelF = preCompensationLabelsF(preCompensationIdxF);

            errorCountsForScenarioF = zeros(1, numel(snrs));
            packetCountsForScenarioF = zeros(1, numel(snrs));
            errorRatesForScenarioF = zeros(1, numel(snrs));
            cfrMagnitudesForScenarioF = cell(1, numel(snrs));

            parfor snrIdxF = 1 : numel(snrs)
                localChannelF = clone(channelF);
                [numPacketErrorsF, numPacketsSimulatedF, cfrMagnitudesF, ~, ~] = simulateTransmission(scenarioParametersF, cfgEHTF, localChannelF, maxChDelayF, snrs(snrIdxF), inf, simParametersF);
                
                errorCountsForScenarioF(snrIdxF) = numPacketErrorsF;
                packetCountsForScenarioF(snrIdxF) = numPacketsSimulatedF;
                errorRatesForScenarioF(snrIdxF) = numPacketErrorsF / numPacketsSimulatedF;
                cfrMagnitudesForScenarioF{snrIdxF} = cfrMagnitudesF;
                
                send(progressQueueF, sprintf("%.1f GHz | Pre-Compensation: %s | SNR: %d", carrierFrequencyF / 1e9, char(scenarioLabelF), snrs(snrIdxF)));
            end

            errorCountsF(carrierFrequencyIdxF, preCompensationIdxF, :) = reshape(errorCountsForScenarioF, 1, 1, []);
            packetCountsF(carrierFrequencyIdxF, preCompensationIdxF, :) = reshape(packetCountsForScenarioF, 1, 1, []);
            errorRatesF(carrierFrequencyIdxF, preCompensationIdxF, :) = reshape(errorRatesForScenarioF, 1, 1, []);
            cfrMagnitudesBySnrF(carrierFrequencyIdxF, preCompensationIdxF, :) = reshape(cfrMagnitudesForScenarioF, 1, 1, []);
        end

        delete(progressListenerF);
        if ishandle(waitbarF)
            close(waitbarF);
        end
    end

    plotErrorRatesF = max(errorRatesF, 1 / packetsToTransmit);
    carrierFrequencyLabelsF = string(carrierFrequenciesF / 1e9) + " GHz";

    figure;
    tiledlayout(1, numel(preCompensationsF));
    for preCompensationIdxF = 1 : numel(preCompensationsF)
        nexttile;
        hold on;
        for carrierFrequencyIdxF = 1 : numel(carrierFrequenciesF)
            semilogy(snrs, squeeze(plotErrorRatesF(carrierFrequencyIdxF, preCompensationIdxF, :)), ['-' markers(carrierFrequencyIdxF)]);
        end

        grid on;
        xlabel("SNR [dB]");
        ylabel("PER");
        ylim([1 / packetsToTransmit 1]);
        legend(carrierFrequencyLabelsF, "Location", "NorthEastOutside");
        title("Pre-compensation: " + preCompensationLabelsF(preCompensationIdxF));
    end

    targetPERF = 1e-1;
    minErrorRatesF = min(errorRatesF, [], 3);
    usableCarrierFrequenciesF = nan(1, numel(preCompensationsF));

    for preCompensationIdxF = 1 : numel(preCompensationsF)
        usableCarrierFrequencyIdxF = find(minErrorRatesF(:, preCompensationIdxF) <= targetPERF, 1, "last");
        if ~isempty(usableCarrierFrequencyIdxF)
            usableCarrierFrequenciesF(preCompensationIdxF) = carrierFrequenciesF(usableCarrierFrequencyIdxF);
        end
    end

    fprintf("Maximum usable carrier frequency without pre-compensation: %.1f GHz\n", usableCarrierFrequenciesF(1) / 1e9);
    fprintf("Maximum usable carrier frequency with ideal pre-compensation: %.1f GHz\n", usableCarrierFrequenciesF(2) / 1e9);
end % calc_task4_4