%% Validate the analytical FWM product count against a brute-force reference
% The paper counts channel combinations, not only unique output frequencies:
%   non-degenerate: i < j, k ~= i, k ~= j, n = i + j - k
%   degenerate:     i == j, k ~= i,       n = 2*i - k

clear;
clc;

N_values = [4, 8, 16];
plot_N = 8;

fprintf('Validating analytical FWM product count from Goebel and Hanik (2008)\n');

for idxN = 1:numel(N_values)
    N = N_values(idxN);
    fprintf('\nN = %d\n', N);

    [Mndg_ana, Mdg_ana] = getProducts(N);
    [Mndg_brute, Mdg_brute, n_values] = getProductsBruteForce(N);

    diff_ndg = Mndg_ana - Mndg_brute;
    diff_dg = Mdg_ana - Mdg_brute;

    fprintf('  Analytical total   : %d\n', sum(Mndg_ana) + sum(Mdg_ana));
    fprintf('  Brute-force total  : %d\n', sum(Mndg_brute) + sum(Mdg_brute));

    if all(diff_ndg == 0) && all(diff_dg == 0)
        fprintf('  Match              : yes\n');
    else
        fprintf('  Match              : no\n');
        fprintf('  Non-degenerate diff: %s\n', mat2str(diff_ndg));
        fprintf('  Degenerate diff    : %s\n', mat2str(diff_dg));
    end

    if N == plot_N
        plotComparison(n_values, Mndg_ana, Mdg_ana, Mndg_brute, Mdg_brute, N);
    end
end

function [Mndg, Mdg, n_values] = getProductsBruteForce(N)
    n_values = (2 - N):(2*N - 1);
    Mndg = zeros(size(n_values));
    Mdg = zeros(size(n_values));

    for idx = 1:numel(n_values)
        n = n_values(idx);

        % Degenerate products: two identical pumps and one different channel.
        for i = 1:N
            k = 2*i - n;
            if isValidChannel(k, N) && (k ~= i)
                Mdg(idx) = Mdg(idx) + 1;
            end
        end

        % Non-degenerate products: unordered pump pair plus one third channel.
        for i = 1:N
            for j = (i + 1):N
                k = i + j - n;
                if isValidChannel(k, N) && (k ~= i) && (k ~= j)
                    Mndg(idx) = Mndg(idx) + 1;
                end
            end
        end
    end
end

function tf = isValidChannel(channel_idx, N)
    tf = (channel_idx >= 1) && (channel_idx <= N) && (channel_idx == round(channel_idx));
end

function plotComparison(n_values, Mndg_ana, Mdg_ana, Mndg_brute, Mdg_brute, N)
    figure;

    subplot(1,2,1);
    stem(n_values, Mndg_ana + Mdg_ana, 'filled', 'LineWidth', 1, 'Marker', 'o', 'MarkerSize', 2);
    hold on;
    stem(n_values, Mdg_ana, 'filled', 'LineWidth', 1, 'Marker', 'none');
    title(['Analytical (N=', num2str(N), ')']);
    xlabel('Product index n');
    ylabel('Number of FWM products');
    legend('Total', 'Degenerate');
    grid on;

    subplot(1,2,2);
    stem(n_values, Mndg_brute + Mdg_brute, 'filled', 'LineWidth', 1, 'Marker', 'o', 'MarkerSize', 2);
    hold on;
    stem(n_values, Mdg_brute, 'filled', 'LineWidth', 1, 'Marker', 'none');
    title(['Brute force (N=', num2str(N), ')']);
    xlabel('Product index n');
    ylabel('Number of FWM products');
    legend('Total', 'Degenerate');
    grid on;
end

function [Mndg, Mdg] = getProducts(N)
    s = abs(2 - N - 1);

    for n = 2 - N:2*N - 1
        if n < -N
            Mdg(n + s) = 0;
        elseif (-N <= n) && (n <= 0)
            Mdg(n + s) = N - ceil((N - n)/2);
        elseif (1 <= n) && (n <= N)
            Mdg(n + s) = N - 1 - floor(n/2) - ceil((N - n)/2);
        elseif (N < n) && (n <= 2*N)
            Mdg(n + s) = N - floor(n/2);
        elseif n > 2*N
            Mdg(n + s) = 0;
        end

        if n < -N
            Mndg(n + s) = 0;
        elseif (-N <= n) && (n < 1)
            Mndg(n + s) = ceil((N^2 + n^2 - 2*N - 2*n + 2*N*n)/4);
        elseif (1 <= n) && (n <= N)
            Mndg(n + s) = ceil(((N^2 - 6*N - 2*n^2 + 2*n + 4)/4) + floor((N*n)/2));
        elseif (N < n) && (n <= 2*N)
            Mndg(n + s) = floor(N^2 + n^2/4 - N*n);
        elseif n > 2*N
            Mndg(n + s) = 0;
        end
    end
end