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


%% Fiber and system parameters
lambda0 = 1309e-9;       % zero-dispersion wavelength [m]
lambda  = 1370e-9;       % operating wavelength [m]
S0      = 0.092;          % dispersion slope [ps/(nm²·km)]
L       = 30e3;          % 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

Dacc = D_lambda * L;
fprintf('Accumulated Dispersion: %.2f ps/nm \n', Dacc / 1e3);

%% Frequency grid
f_max = 200e9;
f = linspace(0, f_max, 5000);                        % [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,'Color','black');
grid on; box on;
xlabel('Frequency [GHz]');
ylabel('Magnitude [dB]');
% title(sprintf('IM/DD Power Fading: 10 km; 1275nm', lambda*1e9, L/1000),"Interpreter","latex");
ylim([-20 0]);

%% Mark analytic null frequencies up to order 5
max_order = 3;
for n = 0:max_order
    f_null_n = sqrt( c*(2*n + 1)/(2*abs(D_si)*lambda^2*L) );
    xline(f_null_n/1e9, '--', 'LineWidth', 1.2, ...
        'Color', [0.1216, 0.4706, 0.7059]);
    text(f_null_n/1e9, -15, sprintf('$f_{\\mathrm{null}, %d}=%.1f$ GHz', n, f_null_n/1e9), ...
        'BackgroundColor', 'w', 'EdgeColor', 'k', 'Interpreter', 'latex', ...
        'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle');
end

% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\dispersion\power_fading.tikz')