function [symbols_for_lvl,avg_for_lvl] = show2Dconstellation(eq_signal,ref_symbols,options)
arguments
    eq_signal
    ref_symbols
    options.fignum (1,1) double = NaN % Default to NaN if not provided
    options.displayname (1,:) char = '' % Default to an empty string if not provided
end

plot_shit = 1;

if isa(eq_signal,'Signal')
    eq_signal = eq_signal.signal;
end
if isa(ref_symbols,'Signal')
    ref_symbols = ref_symbols.signal;
end


if plot_shit
    % Determine the figure number to use or create a new figure
    if isnan(options.fignum)
        fig = figure; % Create a new figure and get its handle
    else
        fig = figure(options.fignum); % Use the specified figure number
        clf;
    end
end


rx_symbols = eq_signal; %./ rms(eq_signal);
correct_symbols = ref_symbols;

col = cbrewer2('Paired',numel(unique(correct_symbols))*2);
ccnt = -1;

levels = unique(correct_symbols);
M = numel(levels);
symbols_for_lvl = NaN(numel(levels),length(correct_symbols));
start = 1;
ende = length(correct_symbols);




for l = 1:numel(levels)
    ccnt = ccnt+2;

    level_amplitude = levels(l);

    symbols_for_lvl(l,correct_symbols==level_amplitude) = rx_symbols(correct_symbols==level_amplitude);

    statistical_mean(l) = mean(symbols_for_lvl(l,:),'omitnan');

end

tx_even = correct_symbols(1:2:end);
tx_odd = correct_symbols(2:2:end);

rx_even = rx_symbols(1:2:end);
rx_odd = rx_symbols(2:2:end);

D_even = tx_even - rx_even;
D_odd  = tx_odd - rx_odd;

D = sqrt(D_even.^2 + D_odd.^2);

[X,Y] = meshgrid(levels, levels);
[X_,Y_] = meshgrid(statistical_mean, statistical_mean);

hold on;
scatter(rx_even,rx_odd,5*ones(1,length(D)),D,'.','DisplayName',['Even/Odd PAM-',num2str(M)],'MarkerEdgeColor',col(2,:));
% colormap(gca,flip(cbrewer2('Spectral',100)))
colormap(gca,'hsv');
scatter(X_(:), Y_(:), 2, 'x', 'LineWidth', 10, 'MarkerEdgeColor', col(6,:),'DisplayName','Statistical Rx Levels');
scatter(X(:), Y(:), 2, 'x', 'LineWidth', 10, 'MarkerEdgeColor', col(4,:),'DisplayName','Tx Levels');
xlim([floor(min(X_(:)))-1, ceil(max(X_(:)))+1]);
ylim([floor(min(Y_(:)))-1, ceil(max(Y_(:)))+1]);

yticks((levels(1:end-1) + levels(2:end)) / 2);
xticks((levels(1:end-1) + levels(2:end)) / 2);
legend 
xlabel('even symbols');
ylabel('odd symbols');
axis equal; grid on;
end