%%FWM analysis from "Analytical Calculation of the Number of
%%Four-Wave-Mixing Products in Optical Multichannel Communication Systems"

N_ = [4,8,16];

df_hz = 400e9;
center_nm = 1310;

figure()
for i = 1:length(N_)

    vec = (2*N_(i)-1:-1:2-N_(i)) -(N_(i)/2+0.5);
    channelplan_hz = nm2hz(center_nm) + (vec * df_hz) ;
    channelplan_nm = hz2nm(channelplan_hz);

    [Mndg,Mdg] = getProducts(N_(i));
    total(i) = sum(Mndg) + sum(Mdg);
    subplot(1,length(N_),i)
    xline(calcWavelengthPlan(N_(i), df_hz, center_nm));
    hold on
    stem(channelplan_nm,(Mndg+Mdg),'filled','LineWidth',1,'Marker','o','MarkerSize',2)
    stem(channelplan_nm,(Mdg),'filled','LineWidth',1,'Marker','none');
    ylim([0,100])
    xlim([1260, 1365]);
    grid off
    xlabel('O-band wavelength region in nm');
    ylabel('Number of FWM products');
    title([num2str(N_(i)),' ch.'])

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