%% ------------------------------------------------------------
% Contour plot: λ_null as function of bandwidth (f_target) and reach (L)
% ------------------------------------------------------------

% Parameters
lambda0 = 1310e-9;   % [m]
S0      = 0.08;      % [ps/(nm²·km)]
c       = physconst('lightspeed');

% Sweep dimensions
f_targets = linspace(50e9, 120e9, 100);   % [Hz]  (x-axis)
L_values  = linspace(0.5e3, 10e3, 100);   % [m]   (y-axis)

lambda_surface = zeros(numel(L_values), numel(f_targets));
Dacc_surface   = zeros(numel(L_values), numel(f_targets));

% Outer loop over fiber length (since L must be scalar)
for iL = 1:numel(L_values)
    L = L_values(iL);
    [lambda_vec, Dacc_vec] = lambda_for_first_null_full(f_targets, L, lambda0, S0);
    lambda_vec = 2*abs(lambda0 - lambda_vec);

    if 0
    fprintf('\n- %d km ------------------------------------\n',L);
    fprintf('   f_null [GHz]     lambda [nm]     Dacc [ps/nm]\n');
    fprintf('----------------------------------------------\n');
    fprintf('%10.1f    %8.2f       %+8.3f\n',[f_targets(:)/1e9, lambda_vec(:)*1e9, Dacc_vec(:)].');
    fprintf('----------------------------------------------\n\n');
    end

    lambda_surface(iL, :) = lambda_vec;   % λ for each f_target
    Dacc_surface(iL, :)   = Dacc_vec;     % corresponding accumulated dispersion
end

% Convert for plotting
lambda_surface_nm = lambda_surface * 1e9;  % [nm]
L_km              = L_values / 1000;       % [km]
f_GHz             = f_targets / 1e9;       % [GHz]

%% Contour plot
figure('Color','w');

% Define wavelength contour levels [nm]
lambda_levels = [1260:10:1290, 1290:5:1300, 1300:2.5:1310];
lambda_levels = [100:-20:50, 50:-10:30,30:-5:0];

% Contour plot
contour(f_GHz, L_km, lambda_surface_nm, lambda_levels, ...
    'LineWidth', 1.5, ...
    'ShowText', 'on', ...
    'LabelFormat', '%.0f nm');

% Colormap and colorbar
colormap((cbrewer2('RdYlGn',100)));
colorbar;
clim([0 100]);

% Axis formatting
xlabel('Signal Bandwidth [GHz]');
ylabel('Fiber length [km]');
legend('$\Delta \lambda$')

% X-axis ticks at 56 : 16 : 150 GHz
xticks(56:8:150);

grid on; box on;


%% Optional: overlay accumulated-dispersion contours
if 0 
    hold on;
    [CS, h] = contour(f_GHz, L_km, Dacc_surface, 10, 'k--', 'LineWidth', 0.8);
    clabel(CS, h, 'Color','k', 'FontSize',8);
end

function [lambda_vec, Dacc_vec] = lambda_for_first_null_full(f_target, L, lambda0, S0)
% lambda_for_first_null_full  (stable, single-branch + validity checks)
% --------------------------------------------------------------------
% Computes the wavelength(s) at which the first IM/DD fading null
% occurs at frequency/ies f_target using the full dispersion model:
%
%   D(lambda) = (S0/4)*(lambda - lambda0^4 / lambda^3)
%
% Restricted to the NORMAL-dispersion branch (λ < λ0),
% and valid only in the O-band (1260–1360 nm).
%
% Inputs:
%   f_target - scalar or vector of target null frequencies [Hz]
%   L        - fiber length [m]
%   lambda0  - zero-dispersion wavelength (ZDW) [m]
%   S0       - dispersion slope at ZDW [ps/(nm²·km)]
%
% Outputs:
%   lambda_vec - wavelength(s) [m] where first null occurs (clamped to O-band)
%   Dacc_vec   - accumulated dispersion(s) [ps/nm] (NaN if out of valid range)
% --------------------------------------------------------------------

    c = physconst('lightspeed');
    S0_si = S0 * 1e3;   % ps/(nm²·km) -> s/(m³)

    % Define O-band boundaries (in meters)
    lambda_min = 1255e-9;
    lambda_max = 1361e-9;

    % Force column vector
    f_target = f_target(:);
    N = numel(f_target);

    lambda_vec = NaN(N,1);
    Dacc_vec   = NaN(N,1);

    for k = 1:N
        RHS = c * 0.5 / (f_target(k)^2 * L);

        % Normal-dispersion branch (λ < λ0)
        fun = @(lambda) -(S0_si/4).*(lambda - (lambda0^4)./(lambda.^3)).*lambda.^2 - RHS;

        % Limit the search to [λ_min, λ0)
        try
            lambda_sol = fzero(fun, [lambda_min, lambda0 * 0.999]);
        catch
            % If the zero is not within bounds, skip this point
            lambda_sol = NaN;
        end

        % Validate solution
        if isnan(lambda_sol) || lambda_sol < lambda_min || lambda_sol > lambda_max
            lambda_vec(k) = NaN;
            Dacc_vec(k)   = NaN;
            continue
        end

        % Compute D(lambda) and accumulated dispersion
        D_lambda = (S0_si/4) * (lambda_sol - (lambda0^4)/(lambda_sol^3)) / 1e-6; % ps/(nm·km)
        Dacc_val = D_lambda * (L/1000);                                         % ps/nm

        % Sanity bound on dispersion (avoid unphysical > ±100 ps/nm)
        if abs(Dacc_val) > 100
            lambda_vec(k) = NaN;
            Dacc_vec(k)   = NaN;
        else
            lambda_vec(k) = lambda_sol;
            Dacc_vec(k)   = Dacc_val;
        end
    end
end