Strengthen WDM and IM/DD optical tests

Functions/Theory/Dispersion/dispersion_power_fading.m

@@ -0,0 +1,165 @@
+%% Chromatic Dispersion Power Fading Demonstration
+% ------------------------------------------------------------
+% This script computes and visualizes power fading after
+% photodiode detection caused by chromatic dispersion in IM/DD links.
+%
+% It also determines the wavelength λ that produces the first
+% fading null at a specified RF frequency f_target using the
+% full physical dispersion model:
+%
+%   D(λ) = (S0/4) * (λ - λ0^4 / λ^3)
+%
+% and compares the analytic null frequency with simulation.
+% ------------------------------------------------------------
+
+% clear; close all; clc;
+
+%% Fiber and wavelength parameters
+lambda0  = 1310e-9;   % Zero-dispersion wavelength (ZDW) [m]
+S0       = 0.08;      % Dispersion slope at ZDW [ps/(nm^2·km)]
+L        = 10000;      % Fiber length [m]
+alpha_dB = 0;         % Attenuation [dB/m] (ignored here)
+
+%% Target null frequency
+f_targets = linspace(55e9,58e9,10);
+f_targets = 56e9;
+% f_targets = 80e9;
+% Compute wavelength that gives the first null at f_target
+[lambda_vec, Dacc_vec] = lambda_for_first_null_full(f_targets, L, lambda0, S0);
+% lambda_vec = 1293e-9;
+
+
+fprintf('\n----------------------------------------------\n');
+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');
+
+%% Frequency grid
+f_simu = 500e9;                % Simulation bandwidth [Hz]
+N_freq = 500000;
+faxis  = linspace(-f_simu/2, f_simu/2, N_freq);
+
+%% Derived fiber parameters
+c  = physconst('lightspeed');
+S0_si = S0 * 1e3;  % ps/(nm²·km) -> s/m³
+
+% Convert wavelengths to nm for the D(lambda) model
+lambda_nm  = lambda_vec(end) * 1e9;
+lambda0_nm = lambda0 * 1e9;
+
+% Dispersion parameter [ps/(nm·km)]
+D_lambda = (S0/4) * (lambda_nm - (lambda0_nm^4)/(lambda_nm^3));
+
+% Convert to [s/m²]
+D_si = D_lambda * 1e-6;
+
+% β2 in [s²/m]
+b2 = -D_si * lambda_vec(end)^2 / (2*pi*c);
+
+%% IM/DD intensity response (simulation)
+phi = 2*pi^2*b2*faxis.^2*L;
+H_field_pos = exp(-1j*phi);  % +f sideband
+H_field_neg = exp(+1j*phi);  % -f sideband
+H_intensity = 0.5 * (H_field_pos + H_field_neg);  % PD beating term
+H_sim = abs(H_intensity);
+
+%% Theoretical analytical IM/DD response
+phi = 2*pi^2 * abs(b2) * faxis.^2 * L;
+H_theoretical = abs(cos(phi));
+
+%% Analytic first null (for verification)
+f_null_analytic = sqrt(c*(0.5)/(abs(D_si)*lambda_vec(end)^2*L));
+fprintf('Analytic first null from D,λ,L: %.2f GHz\n\n', f_null_analytic/1e9);
+
+%% Plot
+cols = linspecer(5);
+figure('Color','w'); hold on; grid on; box on;
+plot(faxis*1e-9, 10*log10(H_sim), 'DisplayName','$|H_{sim}|$ (IM/DD simulation)','Color',cols(1,:));
+plot(faxis*1e-9, 10*log10(H_theoretical), 'DisplayName','|cos($\phi$)| (theory)','Color',cols(2,:),'LineStyle','--');
+xline(f_targets(end)/1e9,'k:','LineWidth',1.2,'DisplayName','Target null (56 GHz)');
+xline(f_null_analytic/1e9,'Color',[0.2 0.6 0.2],'LineStyle','-.','LineWidth',1.2,'DisplayName','Analytic null');
+xlabel('Frequency [GHz]');
+ylabel('Magnitude [dB]');
+title(sprintf('Power Fading for %.2f nm,  L = %.1f km',lambda_nm,L/1000));
+legend('Location','best'); ylim([-30 0]);
+
+%% Plot Bandwidth vs Lambda max
+
+figure();
+hold on;
+plot(lambda_vec.*1e6,f_targets.*1e-9)
+xlabel('wavelength');
+ylabel('max. Bandwidth')
+
+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