%% ============================================================
%  Minimal IM/DD Power Fading Plot
% ============================================================


%% Fiber and system parameters
lambda0 = 1310e-9;       % zero-dispersion wavelength [m]
lambda  = 1275e-9;       % operating wavelength [m]
S0      = 0.08;          % dispersion slope [ps/(nm²·km)]
L       = 10e3;          % fiber length [m]
c       = physconst('lightspeed');

%% Derived quantities
S0_si = S0 * 1e3;                                    % → s/m³
D_lambda = (S0/4) * (lambda*1e9 - (lambda0*1e9)^4/(lambda*1e9)^3); % ps/(nm·km)
D_si = D_lambda * 1e-6;                              % → s/m²
b2   = -D_si * lambda^2 / (2*pi*c);                  % s²/m

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

%% IM/DD transfer function (power fading)
phi = 2*pi^2 * b2 * f.^2 * L;
H = abs(cos(phi));

%% Plot
figure('Color','w');
plot(f/1e9, 10*log10(H), 'LineWidth', 1.8);
grid on; box on;
xlabel('Frequency [GHz]');
ylabel('Magnitude [dB]');
title(sprintf('IM/DD Power Fading |H| for λ = %.1f nm, L = %.1f km', lambda*1e9, L/1000));
ylim([-30 0]);

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