%% =========================
% Routine A: BER vs ROP plots
%% =========================
function S = plot_BER_vs_ROP(res, varargin)
% S = plot_BER_vs_ROP(res, 'fec', 3.8e-3, 'eval_list', [], 'colormap', 'Spectral')
%
% Creates one BER-vs-ROP figure per evaluated distance.
% Returns S struct with FEC crossings: S.Sffe, S.Sdfe, S.Svnle, S.Smlse, S.Sdbt
%
% res.* assumed: res.ffe{ch,rop,realiz,eval_distance} etc.

p = inputParser;
p.addParameter('fec', 3.8e-3, @(x)isnumeric(x)&&isscalar(x));
p.addParameter('eval_list', [], @(x)isnumeric(x));
p.addParameter('colormap', 'Spectral', @(x)ischar(x) || isstring(x));
p.parse(varargin{:});
fec = p.Results.fec;

% -------------------- Basic metadata --------------------
rop            = res.settings.rop(:);
wavelengthplan = res.settings.wavelengthplan(:);
distances      = res.eval_dist_km(:);

N_ch  = numel(wavelengthplan);
N_rop = numel(rop);

% -------------------- Determine dims robustly --------------------
dims = size(res.ffe);
if numel(dims) < 4, dims(end+1:4) = 1; end
N_distances = dims(4);

if numel(distances) ~= N_distances
    distances = (1:N_distances).';
end

% -------------------- Choose eval distances --------------------
eval_list = p.Results.eval_list;
if isempty(eval_list)
    eval_list = 1:N_distances;
end

% -------------------- Colors --------------------
try
    cols = cbrewer2(char(p.Results.colormap), N_ch);
catch
    cols = linspecer(N_ch);
end

% -------------------- Crossings containers --------------------
S = struct();
S.rop = rop;
S.wavelengthplan = wavelengthplan;
S.distances = distances;
S.fec = fec;

S.Sffe  = cell(N_distances, N_ch);
S.Sdfe  = cell(N_distances, N_ch);
S.Svnle = cell(N_distances, N_ch);
S.Smlse = cell(N_distances, N_ch);
S.Sdbt  = cell(N_distances, N_ch);

% -------------------- Plot per distance --------------------
for eval_ptr = eval_list

    figure('Name', sprintf('BER vs ROP @ %s', distLabel(distances, eval_ptr)));
    hold on;

    for ch = 1:N_ch

        % Extract cell slices
        ffe_cells  = sliceCells4D(res.ffe , ch, eval_ptr, N_rop);
        dfe_cells  = sliceCells4D(res.dfe , ch, eval_ptr, N_rop);
        vnle_cells = sliceCells4D(res.vnle, ch, eval_ptr, N_rop);
        mlse_cells = sliceCells4D(res.mlse, ch, eval_ptr, N_rop);
        dbt_cells  = sliceCells4D(res.dbt , ch, eval_ptr, N_rop);

        % Crossings (only depends on cells)
        [S.Sffe{eval_ptr,ch},  ~] = fecCrossings(rop, ffe_cells,  fec);
        [S.Sdfe{eval_ptr,ch},  ~] = fecCrossings(rop, dfe_cells,  fec);
        [S.Svnle{eval_ptr,ch}, ~] = fecCrossings(rop, vnle_cells, fec);
        [S.Smlse{eval_ptr,ch}, ~] = fecCrossings(rop, mlse_cells, fec);
        [S.Sdbt{eval_ptr,ch},  ~] = fecCrossings(rop, dbt_cells,  fec);

        % BER matrices (complete realizations only)
        ffe_mat = extractCompleteBER(ffe_cells);  % N_rop x K

        if ~isempty(ffe_mat)
            % Your current style: plot each realization curve (no boundedline)
            plot(rop, ffe_mat, 'Color', cols(ch,:), ...
                'DisplayName', sprintf('FFE @ %dnm', round(wavelengthplan(ch))));
        end
    end

    set(gca,'XScale','linear','YScale','log', ...
        'TickLabelInterpreter','latex','FontSize',11);

    yline([3.8e-3, 2.2e-4], 'HandleVisibility','off','LineWidth',1.5);

    ylabel('BER');
    xlabel('ROP [dB]');
    title(sprintf('BER vs. ROP @ %s', distLabel(distances, eval_ptr)));
    xlim([min(rop) max(rop)]);
    ylim([1e-5 0.3]);
    grid on;
    legend show;
    beautifyBERplot;
end
end



%% =========================
% Shared helpers (unchanged)
%% =========================

function cells2D = sliceCells4D(C4, ch, eval_ptr, N_rop)
dims = size(C4);
if numel(dims) < 4, dims(end+1:4) = 1; end
N_realiz = dims(3);

cells2D = cell(N_rop, N_realiz);

chMax   = dims(1);
ropMax  = dims(2);
rMax    = dims(3);
eMax    = dims(4);

if ch > chMax || eval_ptr > eMax
    return;
end

ropUse = min(N_rop, ropMax);
rUse   = min(N_realiz, rMax);

tmp = squeeze(C4(ch, 1:ropUse, 1:rUse, eval_ptr));
tmp = reshape(tmp, ropUse, rUse);

cells2D(1:ropUse, 1:rUse) = tmp;
end

function lbl = distLabel(distances, eval_ptr)
if eval_ptr <= numel(distances)
    d = distances(eval_ptr);
    if isfinite(d)
        lbl = sprintf('%.0f km', d);
        return;
    end
end
lbl = sprintf('eval\\_%d', eval_ptr);
end

function [S, noCrossingMask] = fecCrossings(rop, cellsNxR, fec)
Y = extractCompleteBER(cellsNxR);
if isempty(Y)
    S = [];
    noCrossingMask = [];
    return;
end

ok = mean(Y,1,'omitnan') <= 0.1;
Y  = Y(:, ok);
if isempty(Y)
    S = [];
    noCrossingMask = [];
    return;
end

nR = size(Y,2);
S  = nan(1,nR);
noCrossingMask = true(1,nR);

rop = rop(:);
for j = 1:nR
    y = Y(:,j);
    above = (y > fec);
    idx = find(above(1:end-1) & ~above(2:end), 1, 'first');
    if ~isempty(idx)
        x1 = rop(idx);   y1 = y(idx);
        x2 = rop(idx+1); y2 = y(idx+1);
        if isfinite(y1) && isfinite(y2) && y2 ~= y1
            t  = (fec - y1) / (y2 - y1);
            S(j) = x1 + t*(x2 - x1);
            noCrossingMask(j) = false;
        end
    end
end
end

function Y = extractCompleteBER(cellSlice)
if isempty(cellSlice), Y = []; return; end
nR = size(cellSlice,2);
keep = false(1,nR);
for r = 1:nR
    col = cellSlice(:,r);
    keep(r) = all(cellfun(@(c) ~isempty(c), col));
end
if ~any(keep), Y = []; return; end
Y = cellfun(@(c) c.metrics.BER, cellSlice(:,keep), 'UniformOutput', true);
end