diff --git a/projects/FSO_transmission/first_analysis.m b/projects/FSO_transmission/first_analysis.m
index c535c19..c47a2b1 100644
--- a/projects/FSO_transmission/first_analysis.m
+++ b/projects/FSO_transmission/first_analysis.m
@@ -81,9 +81,23 @@ Pform = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"al
 Rx_matched = Pform.process(Scope_sig);
 Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1);
 
-% timing sync -> at this point we still have no symbol timing recovery, we
-% try to do this with 2sps EQ!
-[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1);
+
+%%
+
+sys = comm.SymbolSynchronizer('TimingErrorDetector', 'Gardner (non-data-aided)', ...
+                              'SamplesPerSymbol', 2, ...
+                              'DampingFactor', 0.7, ...
+                              'NormalizedLoopBandwidth', 0.01);
+Rx_symbolsync = Rx_matched;
+[Rx_symbolsync.signal, timing_error] = sys(Rx_matched.signal);
+
+plot(timing_error); % If this is a ramp, you have drift!
+
+%% timing sync -> at this point we still have no symbol timing recovery, we
+% % try to do this with 2sps EQ!
+
+[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_symbolsync.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1);
+
 
 %% not working.. 
 Rx_synced = Rx_synced_cell{1};
@@ -136,6 +150,22 @@ mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_styl
 [vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ...
     "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style);
 
-mlse_results.metrics.print
+mlse_results.metrics.print("description",'MLSE')
 fprintf('My EQ: %.1e \n',mlse_results.metrics.BER);
-fprintf('Paper: %.1e \n \n',ber_in_paper);
\ No newline at end of file
+fprintf('Paper: %.1e \n \n',ber_in_paper);
+
+%% -------------------- DB target --------------------
+mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels);
+
+eq_       = EQ("Ne",[50, 5, 5],"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5, ...
+    "K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
+
+dbt_results = duobinary_target(eq_,mlse_db_, M, Rx_synced, Symbols, Bits, ...
+    "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [],"eth_style_symbol_mapping",mapping_style);
+
+dbt_results.metrics.print("description",'Duobinary');
+mlse_results.metrics.print
+fprintf('My EQ: %.1e \n',dbt_results.metrics.BER);
+fprintf('Paper: %.1e \n \n',ber_in_paper);
+
+
diff --git a/projects/WDM/WDM_auswertung.m b/projects/WDM/WDM_auswertung.m
index 4c213fa..98b42b9 100644
--- a/projects/WDM/WDM_auswertung.m
+++ b/projects/WDM/WDM_auswertung.m
@@ -1,7 +1,58 @@
-filename = "C:\Users\Silas\Documents\MATLAB\Datensätze\FWM_2025\WDM_20260106_131152_n261_14408161_10km_16ch_400ghz_co_alpha0_5.mat";
-filename = "C:\Users\Silas\Documents\MATLAB\Datensätze\FWM_2025\WDM_20260106_131152_n262_14408165_10km_16ch_400ghz_co_alpha-0_5.mat";
-filename = "C:\Users\Silas\Documents\MATLAB\Datensätze\FWM_2025\WDM_20260106_134921_n223_14409338_10km_16ch_400ghz_co_alpha0_0.mat";
-res = load(filename,'res');
+
+base = "C:\Users\Silas\Nextcloud\Cluster";
+all_files = dir(fullfile(base, "**/*.mat"));
+
+schemes = ["co","pair","alt","seg"];
+
+% Preallocate as table (minimal + convenient)
+T = table('Size',[0 10], ...
+    'VariableTypes', ["string","string","string","datetime","double","double","double","double","double","double"], ...
+    'VariableNames', ["folder","file","scheme","date","node","jobid","L_km","Nch","df_GHz","alpha"]);
+
+% filename parser
+rx = "^WDM_(?<date>\d{8})_(?<time>\d{6})_n(?<node>\d+)_(?<jobid>\d+)_" + ...
+     "(?<L>\d+)km_(?<Nch>\d+)ch_(?<df>\d+)ghz_(?<scheme>[a-z]+)_alpha(?<alpha>\d+(?:_\d+)?)\.mat$";
+
+for k = 1:numel(all_files)
+    f = all_files(k);
+    folder = string(f.folder);
+    file   = string(f.name);
+
+    tok = regexp(file, rx, 'names');
+    if isempty(tok)
+        continue
+    end
+
+    % scheme from filename is the most reliable
+    scheme = string(tok.scheme);
+
+    % optional: ignore unexpected schemes
+    if ~any(strcmpi(scheme, schemes)), continue; end
+
+    node  = str2double(tok.node);
+    jobid = str2double(tok.jobid);
+    L_km  = str2double(tok.L);
+    Nch   = str2double(tok.Nch);
+    df_GHz = str2double(tok.df);
+    date = datetime(strcat(tok.date, tok.time), 'InputFormat','yyyyMMddHHmmss');
+
+
+    % alpha uses "_" as decimal separator in your filenames
+    alpha = str2double(strrep(tok.alpha, "_", "."));
+
+    T(end+1,:) = {folder, file,  scheme,date, node, jobid, L_km, Nch, df_GHz, alpha};
+end
+
+%%
+idx = strcmp(T.scheme,"co") & ...
+      T.alpha == 0.4 & ...
+      T.date >= datetime(2026,1,9) & ...
+      T.date <= datetime(2026,1,10,23,59,59);
+
+T_sel = T(idx,:);
+
+i = 2;
+res = load(fullfile(T_sel.folder(i),T_sel.file(i)),'res');
 res = res.res;
 
 %%
@@ -9,4 +60,4 @@ res = res.res;
 S = plot_BER_vs_ROP(res, 'fec', 3.8e-3);
 
 %% Routine B: violin plot (independent)
-plot_FEC_violin(res, 'tech','VNLE', 'fec',3.8e-3, 'eval_ptr',3, 'ylim',[-10 -3]);
+plot_FEC_violin(res, 'tech','VNLE', 'fec',3.8e-3, 'ylim',[-10 0]);
diff --git a/projects/WDM/WDM_model_10km_queue.m b/projects/WDM/WDM_model_10km_queue.m
index e81cbb2..a41cb02 100644
--- a/projects/WDM/WDM_model_10km_queue.m
+++ b/projects/WDM/WDM_model_10km_queue.m
@@ -6,7 +6,7 @@ function WDM_model_10km_queue(options)
 
 %%% Run parameters
 arguments
-    options.num_channels = 16;
+    options.num_channels = 8;
     options.channel_spacing = 400e9;
     options.fiber_length_km = 10;
     options.rand_key = 1;
@@ -26,7 +26,7 @@ warning('off','MATLAB:datetime:AmbiguousTimeZone');
 
 %% How many workers?
 cpus = str2double(getenv('SLURM_CPUS_PER_TASK'));
-if ~isfinite(cpus) || cpus < 1, cpus = max(1, feature('numcores')); end
+if ~isfinite(cpus) || cpus < 1, cpus = max(1, feature('numcores'))-1; end
 
 % Use a per-job, node-local JobStorageLocation to avoid stale locks on $HOME
 tmpbase = getenv('TMPDIR');
@@ -43,7 +43,7 @@ c.JobStorageLocation = jsl;
 p = gcp('nocreate');
 if isempty(p) || p.NumWorkers ~= cpus
     if ~isempty(p), delete(p); end
-    p = parpool(c, cpus);
+    p = parpool(c, cpus,"IdleTimeout",300);
 end
 fprintf('parpool up with %d workers; JobStorage=%s\n', p.NumWorkers, c.JobStorageLocation);
 
@@ -119,7 +119,7 @@ upsample_pow      = 2^nextpow2(upsample_required);
 s.f_opt     = fdac*kover*upsample_pow;
 s.f_opt_nyq = s.f_opt/2;
 
-s.rop = -10:0.5:0;
+s.rop = -10:1:0;
 
 %% ---------- Intermediate evaluation points (distance dimension) ----------
 segment_length = 1; % km (must match fiber loop below)
@@ -184,7 +184,7 @@ for realiz = 1:s.num_realiz
     for l = 1:N
 
         [Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ...
-            "fsym",fsym,"M",s.M,"order",17,"useprbs",0, ...
+            "fsym",fsym,"M",s.M,"order",15,"useprbs",0, ...
             "fs_out",fdac, ...
             "applyclipping",0,"clipfactor",1.5, ...
             "applypulseform",apply_pulsef,"pulseformer",Pform, ...
@@ -282,20 +282,21 @@ for realiz = 1:s.num_realiz
                 [F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt] = ...
                     throttle_inflight(F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, maxInFlight);
 
-                % Enqueue one job per (l,ri)
-                F(end+1,1) = parfeval(p, @rx_job, 5, ...
+                m = struct('l',l,'ri',ri,'realiz',realiz,'eval_ptr',eval_ptr);
+
+                [F, meta] = enqueue_atomic(F, meta, p, @rx_job, 5, ...
                     Opt_sig_wdm_demux{l}, ...
                     s.rop(ri), ...
                     Symbols{l}, Tx_bits{l}, ...
                     s, l, ri, realiz, eval_ptr, ...
                     fdac, kover, fadc, fsym, ...
                     len_tr, mu_dc, mu_ffe, mu_dfe, dfe_order, duob_mode, ...
-                    memlog_dir);
+                    memlog_dir, ...
+                    'META', m);
 
                 fprintf('Enqueued job %d/%d for realiz %d/%d, l=%d/%d, ri=%d/%d at 0 km (inflight=%d/%d)\n', ...
                     cnt, total, realiz, s.num_realiz, l, N, ri, length(s.rop), numel(F), maxInFlight);
 
-                meta(end+1) = struct('l',l,'ri',ri,'realiz',realiz,'eval_ptr',eval_ptr);
             end
         end
 
@@ -304,11 +305,6 @@ for realiz = 1:s.num_realiz
     end
     % -----------------------------------------------------
 
-    save_results_per_realization( ...
-    s, eval_dist_km, ...
-    output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, ...
-    output_root, fname);
-
     for seg = 1:nSegments
 
         fprintf('Realiz %d/%d: Segment %d/%d \n', realiz, s.num_realiz, seg, nSegments);
@@ -336,19 +332,21 @@ for realiz = 1:s.num_realiz
                     [F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt] = ...
                         throttle_inflight(F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, maxInFlight);
 
-                    F(end+1,1) = parfeval(p, @rx_job, 5, ...
+                    m = struct('l',l,'ri',ri,'realiz',realiz,'eval_ptr',eval_ptr);
+
+                    [F, meta] = enqueue_atomic(F, meta, p, @rx_job, 5, ...
                         Opt_sig_wdm_demux{l}, ...
                         s.rop(ri), ...
                         Symbols{l}, Tx_bits{l}, ...
                         s, l, ri, realiz, eval_ptr, ...
                         fdac, kover, fadc, fsym, ...
                         len_tr, mu_dc, mu_ffe, mu_dfe, dfe_order, duob_mode, ...
-                        memlog_dir);
+                        memlog_dir, ...
+                        'META', m);
 
                     fprintf('Seg: %d - Enqueued job %d/%d for realiz %d/%d, l=%d/%d, ri=%d/%d (inflight=%d/%d)\n', ...
                         seg, cnt, total, realiz, s.num_realiz, l, N, ri, length(s.rop), numel(F), maxInFlight);
 
-                    meta(end+1) = struct('l',l,'ri',ri,'realiz',realiz,'eval_ptr',eval_ptr);
                 end
             end
 
@@ -546,15 +544,21 @@ end
 first = true;
 while ~isempty(F)
 
+    assert(numel(F)==numel(meta))
     idxUn = strcmp({F.State}, 'unavailable');
     if any(idxUn)
+        for ii = find(idxUn)
+            m = meta(ii);
+            fprintf('[UNAVAILABLE] realiz=%d eval=%d l=%d ri=%d\n', ...
+                m.realiz, m.eval_ptr, m.l, m.ri);
+        end
         warning('collect_done: dropping %d unavailable futures.', sum(idxUn));
         F(idxUn) = [];
         meta(idxUn) = [];
+
         if isempty(F), break; end
     end
 
-
     if first
         timeout = timeout_first;
         first = false;
@@ -703,65 +707,146 @@ catch
 end
 end
 
+
 function [F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt] = ...
     throttle_inflight(F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, maxInFlight)
 %THROTTLE_INFLIGHT
-% Enforces a HARD cap on outstanding futures (running + queued + finished-not-yet-fetched).
-% Blocks until numel(F) < maxInFlight, then returns.
+% Hard cap on outstanding futures. Robust against:
+%   - 'unavailable' futures (ID=-1 or State unavailable)
+%   - finished-with-error futures (State may be 'finished' with non-empty Error)
+%   - rare F/meta desync
+%   - race-y property access on FevalFuture arrays (use snapshot Fs)
 
 while numel(F) >= maxInFlight
 
-    % Wait until at least one finishes OR fails
-    % (wait('finished',1) only wakes for finished, so we also handle failures via periodic checks)
+    % Wait for at least one to finish (or time out quickly)
     try
         wait(F, 'finished', 1);
     catch
-        % If wait fails for some reason, fall back to a short pause
         pause(0.1);
     end
 
-    % Harvest finished and failed futures (robust snapshot)
-    n = numel(F);
+    % ---- Ensure F/meta alignment (best-effort repair) ----
+    if numel(F) ~= numel(meta)
+        warning('throttle_inflight: F/meta desync (F=%d meta=%d). Attempting repair.', numel(F), numel(meta));
 
-    % Safety: meta must match F length
-    if numel(meta) ~= n
-        % Keep them aligned (rare; indicates earlier mismatch)
-        m = min(numel(meta), n);
-        F = F(1:m);
+        % Prefer dropping unusable futures first (likely the unmatched ones)
+        Fs = F; % snapshot
+        n  = numel(Fs);
+
+        st = cell(1,n);
+        idbad = false(1,n);
+        for k = 1:n
+            try, st{k} = Fs(k).State; catch, st{k} = ''; end
+            try, idbad(k) = (Fs(k).ID < 0); catch, idbad(k) = false; end
+        end
+        idxBad = find(strcmp(st,'unavailable') | idbad);
+
+        if numel(F) > numel(meta) && ~isempty(idxBad)
+            need = numel(F) - numel(meta);
+            idxBad = idxBad(1:min(numel(idxBad), need));
+            idxBad = idxBad(idxBad >= 1 & idxBad <= numel(F));
+            if ~isempty(idxBad)
+                warning('throttle_inflight: dropping %d bad futures to restore alignment.', numel(idxBad));
+                F(idxBad) = [];
+            end
+        end
+
+        % Final fallback: truncate both to min length
+        m = min(numel(F), numel(meta));
+        F    = F(1:m);
         meta = meta(1:m);
-        n = m;
     end
 
-    states = cell(1,n);
+    if isempty(F)
+        break;
+    end
+
+    % ===================== IMPORTANT: snapshot =====================
+    Fs = F;                 % stable snapshot for property access
+    n  = numel(Fs);
+    % ===============================================================
+
+    % ---- Compute unavailable indices from snapshot ONLY ----
+    st = cell(1,n);
+    idbad = false(1,n);
     for k = 1:n
-        states{k} = F(k).State;
+        try, st{k} = Fs(k).State; catch, st{k} = ''; end
+        try, idbad(k) = (Fs(k).ID < 0); catch, idbad(k) = false; end
     end
 
-    doneFinished = strcmp(states, 'finished');
-    doneFailed   = strcmp(states, 'failed');
-    done         = doneFinished | doneFailed;
+    idxUn = find(strcmp(st,'unavailable') | idbad);
 
+    % Sanitize indices (prevents out-of-range deletions even if something weird happens)
+    idxUn = idxUn(idxUn >= 1 & idxUn <= numel(F));
+
+    if ~isempty(idxUn)
+        for ii = idxUn(:).'
+            if ii <= numel(meta)
+                mm = meta(ii);
+                fprintf('[UNAVAILABLE@throttle] realiz=%d eval=%d l=%d ri=%d\n', ...
+                    mm.realiz, mm.eval_ptr, mm.l, mm.ri);
+            else
+                fprintf('[UNAVAILABLE@throttle] meta_missing for ii=%d\n', ii);
+            end
+        end
+        warning('throttle_inflight: dropping %d unavailable/bad futures.', numel(idxUn));
+
+        % Delete in descending order is safest for structs (not strictly necessary, but robust)
+        idxUn = sort(idxUn, 'descend');
+        F(idxUn) = [];
+        meta(idxUn) = [];
+        continue; % re-check cap after shrinking
+    end
+
+    % ---- Determine "done" (finished or finished-with-error) ----
+    isFinished = strcmp(st, 'finished');
+
+    hasErr = false(1,n);
+    for k = 1:n
+        try
+            hasErr(k) = ~isempty(Fs(k).Error);  % use snapshot object
+        catch
+            hasErr(k) = false;
+        end
+    end
+
+    done   = isFinished | hasErr;
     idxDone = find(done);
+    idxDone = idxDone(idxDone >= 1 & idxDone <= numel(F)); % sanitize
 
     if isempty(idxDone)
         pause(0.05);
         continue;
     end
 
-    % Fetch outputs for finished; mark failed as empty results
+    % ---- Harvest done futures ----
     for jj = 1:numel(idxDone)
         ii = idxDone(jj);
         m  = meta(ii);
 
-        if doneFinished(ii)
-            [ffe_r, dfe_r, vnle_r, mlse_r, dbt_r] = fetchOutputs(F(ii));
+        if isFinished(ii) && ~hasErr(ii)
+            try
+                [ffe_r, dfe_r, vnle_r, mlse_r, dbt_r] = fetchOutputs(F(ii));
+            catch ME
+                warning('throttle_inflight: fetchOutputs failed: realiz=%d eval=%d l=%d ri=%d | %s', ...
+                    m.realiz, m.eval_ptr, m.l, m.ri, ME.message);
+                ffe_r = []; dfe_r = []; vnle_r = []; mlse_r = []; dbt_r = [];
+            end
         else
-            % Failed
+            % finished-with-error
             try
                 err = F(ii).Error;
-                warning('rx_job failed: realiz=%d eval=%d l=%d ri=%d | %s', ...
-                    m.realiz, m.eval_ptr, m.l, m.ri, err.message);
+                if ~isempty(err)
+                    warning('rx_job error: realiz=%d eval=%d l=%d ri=%d | %s', ...
+                        m.realiz, m.eval_ptr, m.l, m.ri, err.message);
+                else
+                    warning('rx_job error: realiz=%d eval=%d l=%d ri=%d | (unknown error)', ...
+                        m.realiz, m.eval_ptr, m.l, m.ri);
+                end
             catch
+                warning('rx_job error: realiz=%d eval=%d l=%d ri=%d | (error unreadable)', ...
+                    m.realiz, m.eval_ptr, m.l, m.ri);
             end
             ffe_r = []; dfe_r = []; vnle_r = []; mlse_r = []; dbt_r = [];
         end
@@ -773,9 +858,43 @@ while numel(F) >= maxInFlight
         output_dbt{m.l, m.ri, m.realiz, m.eval_ptr}  = dbt_r;
     end
 
-    % Remove harvested futures/meta entries using numeric indices (safe)
+    % Remove harvested entries (descending for safety)
+    idxDone = sort(idxDone, 'descend');
     F(idxDone) = [];
     meta(idxDone) = [];
+
 end
 
 end
+
+
+
+
+function [F, meta] = enqueue_atomic(F, meta, p, fun, nOut, varargin)
+%ENQUEUE_ATOMIC Enqueue a parfeval job and append matching meta atomically.
+% Usage:
+%   m = struct('l',l,'ri',ri,'realiz',realiz,'eval_ptr',eval_ptr);
+%   [F, meta, f] = enqueue_atomic(F, meta, p, @rx_job, 5, args..., 'META', m);
+
+% Split varargin into args + meta struct
+idx = find(strcmp(varargin, 'META'), 1, 'last');
+if isempty(idx) || idx == numel(varargin)
+    error('enqueue_atomic: META struct must be provided as last named argument.');
+end
+args = varargin(1:idx-1);
+m    = varargin{idx+1};
+
+% Create future first (local variable), but don't mutate F/meta yet
+f = parfeval(p, fun, nOut, args{:});
+
+% If future is unusable, do not append (prevents meta shift)
+if f.ID < 0 || strcmp(f.State, 'unavailable')
+    warning('enqueue_atomic: got unusable future (ID=%d, State=%s). Dropping enqueue: realiz=%d eval=%d l=%d ri=%d', ...
+        f.ID, string(f.State), m.realiz, m.eval_ptr, m.l, m.ri);
+    return;
+end
+
+% Now append both together (atomic w.r.t. fprintf etc.)
+F(end+1,1)    = f;
+meta(end+1,1) = m;
+end
diff --git a/projects/WDM/plot_BER_vs_ROP.m b/projects/WDM/plot_BER_vs_ROP.m
index b8ac31f..86a1fdc 100644
--- a/projects/WDM/plot_BER_vs_ROP.m
+++ b/projects/WDM/plot_BER_vs_ROP.m
@@ -85,7 +85,6 @@ for eval_ptr = eval_list
         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
diff --git a/projects/WDM/plot_FEC_violin.m b/projects/WDM/plot_FEC_violin.m
index 5ce9217..7bf78e3 100644
--- a/projects/WDM/plot_FEC_violin.m
+++ b/projects/WDM/plot_FEC_violin.m
@@ -1,29 +1,38 @@
 function plot_FEC_violin(res, varargin)
-% plot_fec_violin_from_res  Self-contained violin plot (no dependency on other routines)
+% plot_FEC_violin  Violin plot of ROP-at-FEC crossings per wavelength/channel
+% Uses the Bechtold community violinplot (MATLAB Central 45134).
+%
+% Requirements:
+% - community violinplot must be on path, either as:
+%     (a) violinplot.m   (shadows MATLAB official), OR
+%     (b) violinplot_community.m (renamed + function name adjusted)
 %
 % Usage:
-%   plot_fec_violin_from_res('WDM_20260106_....mat', 'tech', 'FFE', 'fec', 3.8e-3, 'eval_ptr', 2, 'ylim', [-10 -6]);
-%   plot_fec_violin_from_res(res,                       'tech', 'FFE', 'fec', 3.8e-3, 'eval_ptr', 2);
-%
-% Notes:
-% - Works with partially-filled res.* cell arrays (missing eval distances / missing realizations).
-% - Uses only complete realizations per channel (all ROP points present).
-% - Requires violinplot.m on path.
+%   plot_FEC_violin(res, 'tech','VNLE', 'fec',3.8e-3, 'eval_ptr',[], 'ylim',[-10 -6], 'bandwidth',0.15);
 
+% -------- args --------
 p = inputParser;
-p.addParameter('tech', 'FFE', @(x)ischar(x) || isstring(x));
-p.addParameter('fec', 3.8e-3, @(x)isnumeric(x)&&isscalar(x));
-p.addParameter('eval_ptr', [], @(x)isnumeric(x)&&isscalar(x));
-p.addParameter('ylim', [], @(x)isnumeric(x)&&numel(x)==2);
-p.addParameter('title_prefix', 'ROP at FEC crossing', @(x)ischar(x)||isstring(x));
+p.addParameter('tech',      'VNLE', @(x)ischar(x)||isstring(x));
+p.addParameter('fec',       3.8e-3, @(x)isnumeric(x)&&isscalar(x));
+p.addParameter('eval_ptr',  [],     @(x)isnumeric(x)&&isscalar(x));
+p.addParameter('ylim',      [],     @(x)isnumeric(x)&&numel(x)==2);
+p.addParameter('bandwidth', [],     @(x)isnumeric(x)&&isscalar(x));
 p.parse(varargin{:});
+opt = p.Results;
 
-if exist('violinplot','file') ~= 2
-    error('violinplot.m not found on path. Add it to path first.');
+tech = upper(string(opt.tech));
+fec  = opt.fec;
+
+% -------- pick result field --------
+switch tech
+    case "FFE",  C4 = res.ffe;
+    case "DFE",  C4 = res.dfe;
+    case "VNLE", C4 = res.vnle;
+    case "MLSE", C4 = res.mlse;
+    case "DBT",  C4 = res.dbt;
+    otherwise, error('Unknown tech "%s". Use FFE/DFE/VNLE/MLSE/DBT.', tech);
 end
 
-
-% ---------- Metadata ----------
 rop            = res.settings.rop(:);
 wavelengthplan = res.settings.wavelengthplan(:);
 distances      = res.eval_dist_km(:);
@@ -31,103 +40,97 @@ distances      = res.eval_dist_km(:);
 N_ch  = numel(wavelengthplan);
 N_rop = numel(rop);
 
-% Determine distances dim robustly from the stored cell arrays (prefer ffe)
-dims = size(res.ffe);
+dims = size(C4);
 if numel(dims) < 4, dims(end+1:4) = 1; end
-N_distances = dims(4);
-if numel(distances) ~= N_distances
-    distances = (1:N_distances).';
+N_dist = dims(4);
+
+eval_ptr = opt.eval_ptr;
+if isempty(eval_ptr), eval_ptr = N_dist; end
+eval_ptr = max(1, min(N_dist, eval_ptr));
+
+% -------- compute crossings per channel + "completeness" --------
+S_cell   = cell(1, N_ch);  % crossings values (finite only)
+K_total  = zeros(1, N_ch); % number of complete realizations for that channel
+K_cross  = zeros(1, N_ch); % number of realizations that crossed
+
+for ch = 1:N_ch
+    cells = sliceCells4D(C4, ch, eval_ptr, N_rop); % N_rop x N_realiz (may contain [])
+    [S_cell{ch}, K_total(ch), K_cross(ch)] = crossings_and_counts(rop, cells, fec);
 end
 
-% Choose eval distance
-eval_ptr = p.Results.eval_ptr;
-if isempty(eval_ptr)
-    eval_ptr = N_distances; % default: last available
-end
-if eval_ptr < 1 || eval_ptr > N_distances
-    error('eval_ptr=%d out of range. Available: 1..%d', eval_ptr, N_distances);
+% -------- build vector+category representation (NO NaNs in xAll) --------
+catLabels = arrayfun(@(nm)sprintf('%d nm', round(nm)), wavelengthplan, 'UniformOutput', false);
+catsAll = categorical(strings(0,1), catLabels, 'Ordinal', false);
+xAll    = zeros(0,1);
+
+for ch = 1:N_ch
+    v = S_cell{ch}(:);
+    if ~isempty(v)
+        xAll    = [xAll; v];
+        catsAll = [catsAll; repmat( ...
+            categorical(string(catLabels{ch}), catLabels, 'Ordinal', false), ...
+            numel(v), 1)];
+    end
 end
 
-% Choose technique field
-tech = upper(string(p.Results.tech));
-switch tech
-    case "FFE",  C4 = res.ffe;
-    case "DFE",  C4 = res.dfe;
-    case "VNLE", C4 = res.vnle;
-    case "MLSE", C4 = res.mlse;
-    case "DBT",  C4 = res.dbt;
-    otherwise
-        error('Unknown tech "%s". Use one of: FFE, DFE, VNLE, MLSE, DBT.', tech);
-end
 
-fec = p.Results.fec;
-
-% ---------- Compute crossings per channel ----------
-
-% ---------- Plot ----------
-figure('Name', sprintf('%s violin @ %s', tech, distLabel(distances, eval_ptr)));
+% -------- plot --------
+figure('Name', sprintf('%s FEC violin @ %s', tech, distLabel(distances, eval_ptr)));
 hold on;
-c = linspecer(N_distances);
-for eval_ptr = 1:N_distances
-    S_cell = cell(1, N_ch);
-    
-    for ch = 1:N_ch
-        cells = sliceCells4D(C4, ch, eval_ptr, N_rop);    % N_rop x N_realiz cell
-        S_cell{ch} = fecCrossings_only(rop, cells, fec); % 1 x K crossings (may be empty)
-    end
-    
-    % Pad to matrix for violinplot: rows=realizations (max K), cols=channels
-    Kmax = max(cellfun(@numel, S_cell));
-    if isempty(Kmax) || Kmax == 0
-        warning('No FEC crossings found for %s at eval_ptr=%d (%s).', tech, eval_ptr, distLabel(distances, eval_ptr));
-        return;
-    end
-    
-    S_mat = NaN(Kmax, N_ch);
-    for ch = 1:N_ch
-        k = numel(S_cell{ch});
-        if k > 0
-            S_mat(1:k, ch) = S_cell{ch}(:);
-        end
-    end
-    
-    catLabels = arrayfun(@(nm) sprintf('%d nm', nm), round(wavelengthplan), 'UniformOutput', false);
-
-    if numel(S_mat)>numel(wavelengthplan)
-        S_mat(isnan(S_mat)) = 0;
-        v=violinplot(S_mat, catLabels, ...
-            'ViolinColor', c(eval_ptr,:), ...
-            'ViolinAlpha', 0.10, ...
-            'MarkerSize', 20, ...
-            'ShowMedian', true, ...
-            'EdgeColor', c(eval_ptr,:), ...
-            'ShowWhiskers', false, ...
-            'ShowData', true, ...
-            'ShowBox', false, ...
-            'Bandwidth', 0.05);
-    else
-        channs = 1:numel(wavelengthplan);
-        scatter( wavelengthplan ,S_mat,10,c(eval_ptr,:),'Marker','o','LineWidth',2,'DisplayName',distLabel(distances, eval_ptr));
-    
-        channs = 1:numel(wavelengthplan);
-        scatter( round(wavelengthplan(S_mat==0)),S_mat(S_mat==0),10,c(eval_ptr,:),'Marker','o','HandleVisibility','off');
-    end
 
+% Violin only if we have any data at all
+if ~isempty(xAll)
+    violinplot_community(xAll, catsAll,'Bandwidth',    0.15, ...                 % KDE bandwidth (critical!)
+    'ViolinColor',  [0.2 0.4 0.8], ...         % or Nx3 for per-channel colors
+    'ViolinAlpha',  0.15, ...
+    'EdgeColor',    [0.3 0.3 0.3], ...
+    'MarkerSize',   14, ...
+    'ShowData',     true, ...
+    'ShowMean',     false, ...
+    'ShowMedian',   true, ...
+    'ShowBox',      false, ...
+    'ShowWhiskers', false);
+else
+    warning('No FEC crossings found at eval_ptr=%d (%s). Plotting only markers.', eval_ptr, distLabel(distances, eval_ptr));
+    set(gca,'XTick',1:N_ch,'XTickLabel',catLabels);
 end
 
+% Marker X positions (violin groups are at 1..N_ch)
+xpos = 1:N_ch;
+
+% -------- overlay mean markers (black/red/blue) --------
+yBlue = -9;
+
+for ch = 1:N_ch
+    if K_total(ch) == 0 || K_cross(ch) == 0
+        % no complete realizations OR none crossed -> blue X at -9
+        scatter(xpos(ch), yBlue, 80, 'x', 'LineWidth', 2, 'MarkerEdgeColor', [0 0 1], 'HandleVisibility','off');
+        continue;
+    end
+
+    mu = mean(S_cell{ch}, 'omitnan');
+
+    if K_cross(ch) < K_total(ch)
+        % some complete realizations did not cross -> red X
+        scatter(xpos(ch), mu, 80, 'x', 'LineWidth', 2, 'MarkerEdgeColor', [1 0 0], 'HandleVisibility','off');
+    else
+        % all complete realizations crossed -> black X
+        scatter(xpos(ch), mu, 80, 'x', 'LineWidth', 2, 'MarkerEdgeColor', [0 0 0], 'HandleVisibility','off');
+    end
+end
+
+title(sprintf('%s | BER %.2e | %s', tech, fec, distLabel(distances, eval_ptr)));
 ylabel('ROP at FEC crossing');
-title(sprintf('%s: %s | BER %.2e | %s', ...
-    p.Results.title_prefix, tech, fec, distLabel(distances, eval_ptr)));
 grid on; box on;
 
-if ~isempty(p.Results.ylim)
-    ylim(p.Results.ylim);
+if ~isempty(opt.ylim)
+    ylim(opt.ylim);
 end
 
 end
 
 
-%% ================= helpers (self-contained) =================
+% ================= helpers =================
 
 function cells2D = sliceCells4D(C4, ch, eval_ptr, N_rop)
 dims = size(C4);
@@ -136,68 +139,64 @@ 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
+if ch > dims(1) || eval_ptr > dims(4)
     return;
 end
 
-ropUse = min(N_rop, ropMax);
-rUse   = min(N_realiz, rMax);
+ropUse = min(N_rop, dims(2));
+rUse   = min(N_realiz, dims(3));
 
 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);
+function [S, K_total, K_cross] = crossings_and_counts(rop, cellsNxR, fec)
+% counts "complete" realizations and how many of those crossed
+% S returns only finite crossings (no NaN).
+
+Y = extractCompleteBER(cellsNxR);   % N_rop x K_total
+K_total = size(Y,2);
+
+if K_total == 0
+    S = [];
+    K_cross = 0;
+    return;
 end
 
-function S = fecCrossings_only(rop, cellsNxR, fec)
-% Returns 1xK crossing positions for all complete realizations (drops bad ones).
-Y = extractCompleteBER(cellsNxR);  % N_rop x K
-if isempty(Y), S = []; return; end
-
-% Optional quality gate (same as your previous scripts)
+% optional quality gate (keep your style)
 ok = mean(Y,1,'omitnan') <= 0.1;
-Y  = Y(:, ok);
-if isempty(Y), S = []; return; end
+Y  = Y(:,ok);
+K_total = size(Y,2);
 
-nR = size(Y,2);
-S  = nan(1,nR);
+if K_total == 0
+    S = [];
+    K_cross = 0;
+    return;
+end
 
 rop = rop(:);
-for j = 1:nR
+S = nan(1, K_total);
+
+for j = 1:K_total
     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);
+        if isfinite(y1) && isfinite(y2) && (y2 ~= y1)
+            t    = (fec - y1) / (y2 - y1);
             S(j) = x1 + t*(x2 - x1);
         end
     end
 end
 
-S = S(isfinite(S)); % keep only successful crossings
+K_cross = sum(isfinite(S));
+S = S(isfinite(S));
 end
 
 function Y = extractCompleteBER(cellSlice)
-% Keep only those realization columns where ALL N_rop entries are non-empty
 if isempty(cellSlice), Y = []; return; end
 nR = size(cellSlice,2);
 keep = false(1,nR);
@@ -208,6 +207,13 @@ for r = 1:nR
 end
 
 if ~any(keep), Y = []; return; end
-
 Y = cellfun(@(c) c.metrics.BER, cellSlice(:,keep), 'UniformOutput', true);
 end
+
+function lbl = distLabel(distances, eval_ptr)
+if eval_ptr <= numel(distances) && isfinite(distances(eval_ptr))
+    lbl = sprintf('%.0f km', distances(eval_ptr));
+else
+    lbl = sprintf('eval_%d', eval_ptr);
+end
+end