function csvFile = model4_queue_aware_dcf_tgax(U, p, csvFile)
% Model 4:
% Queue-Aware DCF/TGax Delay Model
%
% Inputs  = p + 18
% Output  = 1 (targetDelay)

if nargin < 2, p = 6; end
if nargin < 3 || isempty(csvFile)
    csvFile = 'model4_queue_aware_dcf_tgax.csv';
end

N = numel(U.packetBytes);

slotMs = 0.009;
difsMs = 0.034;
sifsMs = 0.016;
ackMs = 0.044;
timeoutMs = 0.050;
cwMin = 15;
cwMax = 1023;

collisionProb = zeros(N,1);
packetDetected = zeros(N,1);
backoff0Slots = zeros(N,1);
retryCount = zeros(N,1);
queueLen = zeros(N,1);
waitMs = zeros(N,1);
serviceMs = zeros(N,1);
successFlag = zeros(N,1);
snrMarginDb = zeros(N,1);
rxPowerdB = zeros(N,1);
txTimeMs = zeros(N,1);
delayMs = zeros(N,1);

serverFreeMs = 0;
meanServiceMs = 0.8;

snrTable = [4 6 8 11 14 17 20 23 26 29 32 35];

for i = 1:N
    cfgHE = wlanHESUConfig;
    cfgHE.ChannelBandwidth = U.channelBandwidth;
    cfgHE.MCS = U.mcs(i);
    cfgHE.APEPLength = U.packetBytes(i);

    fs = wlanSampleRate(cfgHE);

    try
        psduBytes = getPSDULength(cfgHE);
    catch
        psduBytes = cfgHE.APEPLength;
    end

    txBits = randi([0 1], 8*psduBytes, 1);
    txWave = wlanWaveformGenerator(txBits, cfgHE);

    txTimeMs(i) = 1e3 * numel(txWave) / fs;
    snrMarginDb(i) = U.snrDb(i) - snrTable(U.mcs(i)+1);

    % Queueing from common arrival trace
    startServiceMs = max(U.arrivalTimeMs(i), serverFreeMs);
    waitMs(i) = startServiceMs - U.arrivalTimeMs(i);
    queueLen(i) = max(0, round(waitMs(i) / max(meanServiceMs, 0.05)));

    % Collision proxy from same load trace
    collisionProb(i) = min(max(1 - exp(-0.035*(U.nStations(i)-1)*U.busyRatio(i)), 0.001), 0.95);

    cw = cwMin;
    totalService = 0;
    ok = false;
    detAny = 0;
    rxPwrTmp = NaN;

    for k = 1:(U.maxRetries+1)
        if k == 1
            backoff0Slots(i) = randi([0 cwMin]);
            boSlots = backoff0Slots(i);
        else
            boSlots = randi([0 cw]);
        end

        tgax = wlanTGaxChannel;
        tgax.SampleRate = fs;
        tgax.ChannelBandwidth = U.channelBandwidth;
        tgax.DelayProfile = U.delayProfile;
        tgax.NumTransmitAntennas = 1;
        tgax.NumReceiveAntennas = 1;

        rxWave = tgax(txWave);
        rxPwrTmp = 10*log10(mean(abs(rxWave).^2) + eps);
        rxWave = awgn(rxWave, U.snrDb(i), 'measured');

        pktOffset = wlanPacketDetect(rxWave, cfgHE.ChannelBandwidth);
        detected = ~isempty(pktOffset) && all(isfinite(pktOffset)) && (pktOffset >= 0);
        detAny = max(detAny, double(detected));

        phyOk = detected && (U.snrDb(i) >= snrTable(U.mcs(i)+1));
        macOk = rand > collisionProb(i);

        totalService = totalService + difsMs + boSlots*slotMs + txTimeMs(i);

        if phyOk && macOk
            totalService = totalService + sifsMs + ackMs;
            ok = true;
            break;
        else
            totalService = totalService + timeoutMs;
            retryCount(i) = retryCount(i) + 1;
            cw = min(2*cw + 1, cwMax);
        end
    end

    packetDetected(i) = detAny;
    rxPowerdB(i) = rxPwrTmp;
    successFlag(i) = double(ok);
    serviceMs(i) = totalService;

    serverFreeMs = startServiceMs + serviceMs(i);
    meanServiceMs = 0.95*meanServiceMs + 0.05*serviceMs(i);

    delayMs(i) = waitMs(i) + serviceMs(i);
end

tbl = table( ...
    U.arrivalGapMs, ...
    U.arrivalTimeMs, ...
    U.packetBytes, ...
    U.nStations, ...
    U.busyRatio, ...
    U.snrDb, ...
    U.mcs, ...
    collisionProb, ...
    packetDetected, ...
    backoff0Slots, ...
    retryCount, ...
    queueLen, ...
    waitMs, ...
    serviceMs, ...
    txTimeMs, ...
    successFlag, ...
    rxPowerdB, ...
    snrMarginDb, ...
    delayMs, ...
    'VariableNames', { ...
    'arrivalGapMs','arrivalTimeMs','packetBytes','nStations','busyRatio', ...
    'snrDb','mcs','collisionProb','packetDetected','backoff0Slots', ...
    'retryCount','queueLen','waitMs','serviceMs','txTimeMs','successFlag', ...
    'rxPowerdB','snrMarginDb','delayMs'});

Treg = make_regression_table(tbl, 'delayMs', p, 1);
writetable(Treg, csvFile);

fprintf('Model 4 CSV saved: %s\n', csvFile);
fprintf('Model 4 inputs = %d, outputs = 1\n', width(Treg)-1);
end
