%% ============================================================
%  IM/DD Power Fading Evolution vs Wavelength (L = 10 km)
% ============================================================

clear; close all; clc;

%% Fixed fiber parameters
lambda0 = 1310e-9;       % zero-dispersion wavelength [m]
S0      = 0.09;          % dispersion slope [ps/(nm^2·km)]
L       = 10e3;          % fiber length FIXED [m]
c       = physconst('lightspeed');

%% Frequency grid
f_max = 150e9;
f = linspace(0, f_max, 4000);  % [Hz]

%% Figure setup (stable axes for clean GIF)
fig = figure('Color','w');
ax = axes(fig);
hold(ax,'on'); grid(ax,'on'); box(ax,'on');
xlabel(ax,'Frequency [GHz]');
ylabel(ax,'Magnitude [dB]');
ylim(ax,[-30 0]);
xlim(ax,[0 f_max/1e9]);

%% GIF writer
g = GifWriter('Name','power_fading_vs_wavelength', ...
              'DelayTime',0.12, ...
              'Parallel',false);

%% Wavelength sweep (around ZDW)
lambda_vec = linspace(1260e-9, 1360e-9, 25);  % 1260–1360 nm

for k = 1:length(lambda_vec)

    lambda = lambda_vec(k);

    %% Dispersion for current wavelength
    D_lambda = (S0/4) * (lambda*1e9 - (lambda0*1e9)^4/(lambda*1e9)^3); % ps/(nm·km)
    D_si     = D_lambda * 1e-6;                                       % s/m^2
    b2       = -D_si * lambda^2 / (2*pi*c);                           % s^2/m

    %% Power fading transfer function
    phi = 2*pi^2 * b2 * f.^2 * L;
    H   = abs(cos(phi));
    HdB = 10*log10(max(H, 1e-12));

    cla(ax)
    plot(ax, f/1e9, HdB, 'LineWidth',1.8,'Color','black');

    %% First-null frequency
    if abs(D_si) > 0
        f_null = sqrt(c*(0.5)/(abs(D_si)*lambda^2*L));
        xline(ax, f_null/1e9, 'r--', 'LineWidth',1.2, ...
            'Label', sprintf('f_{null}=%.1f GHz', f_null/1e9), ...
            'LabelOrientation','horizontal', ...
            'LabelVerticalAlignment','bottom');
    end

    title(ax, sprintf('Power Fading for: 10 km @ %.0f nm', lambda*1e9));

    drawnow;
    g.addFrame(fig);
end

%% Compile GIF
g.compile(fig.Number);

disp(fullfile(g.OutputDir, sprintf('%s_fig_%d.gif', g.Name, fig.Number)));