From b5387f78c686571ad19f2587a660eee18acaacc3 Mon Sep 17 00:00:00 2001 From: "silas (home)" Date: Fri, 17 Oct 2025 08:02:30 +0200 Subject: [PATCH] Dispersion theory scripts --- Functions/Theory/dispersion_10km.m | 90 ++++++++++++++ .../Theory/dispersion_wavelength_notch.m | 32 +++++ Functions/Theory/dispersion_wdm.m | 111 ++++++++++++++++++ Functions/Theory/power_fading.m | 40 +++++++ Functions/convert_freq_lambda.m | 62 ++++++++++ 5 files changed, 335 insertions(+) create mode 100644 Functions/Theory/dispersion_10km.m create mode 100644 Functions/Theory/dispersion_wavelength_notch.m create mode 100644 Functions/Theory/dispersion_wdm.m create mode 100644 Functions/Theory/power_fading.m create mode 100644 Functions/convert_freq_lambda.m diff --git a/Functions/Theory/dispersion_10km.m b/Functions/Theory/dispersion_10km.m new file mode 100644 index 0000000..87a17a8 --- /dev/null +++ b/Functions/Theory/dispersion_10km.m @@ -0,0 +1,90 @@ +%% ============================================================ +% IM/DD Fading Notch – λ_null vs. Bandwidth (Fixed 10 km) +% ============================================================ + + + +%% Fiber and dispersion parameters +lambda0 = 1315e-9; % Zero-dispersion wavelength [m] +S0 = 0.08; % Dispersion slope at ZDW [ps/(nm²·km)] +L = 10e3; % Fiber length [m] +c = physconst('lightspeed'); + +%% Frequency sweep (defines the desired first-fading notch) +f_targets = linspace(40e9, 150e9, 200); % [Hz] +f_GHz = f_targets / 1e9; + +%% Compute wavelength λ_null for each target f_null +[lambda_vec, Dacc_vec] = lambda_for_first_null_full(f_targets, L, lambda0, S0); +lambda_nm = lambda_vec * 1e9; % Convert to nm + +%% ------------------------------------------------------------ +% Plot λ_null vs. f_null for 10 km fiber +% ------------------------------------------------------------ +% figure('Color','w'); +% plot(lambda_nm,f_GHz, 'LineWidth', 2); +% grid on; box on; + +cols = cbrewer2('Paired',10); +figure('Color','w');hold on + +plot(lambda_nm, f_GHz, 'LineWidth',2,'DisplayName',sprintf('%d km',L),'Color',cols(2,:)); + +yticks([56,75,90,112]) + +f_GHz = [56,75,90,112] * 1e9; +[lambda_vec, Dacc_vec] = lambda_for_first_null_full(f_GHz, L, lambda0, S0); +lambda_nm = lambda_vec * 1e9; % Convert to nm + +xticks(round(lambda_nm)) + +xlabel('$\Delta \lambda$ from ZDW [nm]'); +ylabel('$F_{null}$ [GHz]'); +grid on; box on; +lim=(lambda0.*1e9)-[8,40]; +xlim([lim(2) lim(1)]); +% ylim([40,130]) + +%% ------------------------------------------------------------ +% Helper function: lambda_for_first_null_full +% Stable, single-branch, clamped to O-band +% ------------------------------------------------------------ +function [lambda_vec, Dacc_vec] = lambda_for_first_null_full(f_target, L, lambda0, S0) +c = physconst('lightspeed'); +S0_si = S0 * 1e3; % ps/(nm²·km) -> s/(m³) + +% Define O-band boundaries [m] +lambda_min = 1260e-9; +lambda_max = 1360e-9; + +f_target = f_target(:); +N = numel(f_target); + +lambda_vec = zeros(N,1); +Dacc_vec = zeros(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; + + % Solve within normal-dispersion range + try + lambda_sol = fzero(fun, [lambda_min, lambda0 * 0.999]); + catch + lambda_sol = lambda_min; + end + + % Clamp to O-band + lambda_sol = min(max(lambda_sol, lambda_min), lambda_max); + lambda_vec(k) = lambda_sol; + + % 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 + Dacc_val = min(max(Dacc_val, -100), 100); + + Dacc_vec(k) = Dacc_val; +end +end diff --git a/Functions/Theory/dispersion_wavelength_notch.m b/Functions/Theory/dispersion_wavelength_notch.m new file mode 100644 index 0000000..0f4da7e --- /dev/null +++ b/Functions/Theory/dispersion_wavelength_notch.m @@ -0,0 +1,32 @@ +%% Dependency f_null vs Delta_lambda +lambda0 = 1310e-9; +S0 = 0.09; % ps/(nm²·km) +L = 10e3; % m +c = physconst('lightspeed'); + +% Convert slope to SI +S0_si = S0 * 1e3; % s/m³ + +Delta_lambda = linspace(5e-9, 80e-9, 300); % [m] detuning +f_null_2 = sqrt( c * 0.5 ./ (S0_si .* abs(Delta_lambda) .* lambda0.^2 .* L) ); +L = 2e3; % m +f_null_10 = sqrt( c * 0.5 ./ (S0_si .* abs(Delta_lambda) .* lambda0.^2 .* L) ); + +cols = cbrewer2('Paired',10); +figure('Color','w');hold on +cnt = 2; +for L = 10%[2,5,10] + f_null_10 = sqrt( c * 0.5 ./ (S0_si .* abs(Delta_lambda) .* lambda0.^2 .* L*1e3) ); + plot(1310-Delta_lambda*1e9, f_null_10/1e9, 'LineWidth',2,'DisplayName',sprintf('%d km',L),'Color',cols(cnt,:)); + cnt = cnt+2; +end +yticks([56,75,90,112]) +tickse = 1310-[7.5, 12, 17, 31.5]; +xticks(flip(tickse)); + +xlabel('$\Delta \lambda$ from ZDW [nm]'); +ylabel('$F_{null}$ [GHz]'); +grid on; box on; +lim=1310-[5,35]; +xlim([lim(2) lim(1)]); +ylim([40,130]) \ No newline at end of file diff --git a/Functions/Theory/dispersion_wdm.m b/Functions/Theory/dispersion_wdm.m new file mode 100644 index 0000000..f9ead0c --- /dev/null +++ b/Functions/Theory/dispersion_wdm.m @@ -0,0 +1,111 @@ +%% ============================================================ +% IM/DD Fading Notch Design Map +% Shows λ_null vs. bandwidth (f_target) and fiber length (L) +% ============================================================ + +clear; close all; clc; + +%% Parameters +lambda0 = 1310e-9; % Zero-dispersion wavelength [m] +S0 = 0.08; % Dispersion slope at ZDW [ps/(nm²·km)] +c = physconst('lightspeed'); + +% Frequency and length sweep +f_targets = linspace(20e9, 140e9, 80); % [Hz] → x-axis +L_values = linspace(0.5e3, 12e3, 80); % [m] → y-axis + +% Preallocate result matrices +lambda_surface = zeros(numel(L_values), numel(f_targets)); +Dacc_surface = zeros(numel(L_values), numel(f_targets)); + +%% Compute λ_null and Dacc for each (f_target, L) +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_surface(iL, :) = lambda_vec; % [m] + Dacc_surface(iL, :) = Dacc_vec; % [ps/nm] +end + +%% Convert to display units +lambda_surface_nm = lambda_surface * 1e9; % [nm] +L_km = L_values / 1000; % [km] +f_GHz = f_targets / 1e9; % [GHz] + +%% ------------------------------------------------------------ +% Contour plot (λ_null as function of f_null and L) +% ------------------------------------------------------------ +figure('Color','w'); + +% Define wavelength contour levels [nm] +lambda_levels = [1260:10:1290, 1290:5:1300, 1300:2:1310]; + +contourf(f_GHz, L_km, lambda_surface_nm, lambda_levels, ... + 'LineWidth', 1.5, ... + 'ShowText', 'on', ... + 'LabelFormat', '%1.1d nm'); + +% Colormap and colorbar +colormap(flip(cbrewer2('RdYlGn',100))); +clim([1260 1310]); +% c = colorbar; +% ylabel(c, 'λ_{null} [nm]', 'Rotation', 90); + +% Axis formatting +xlabel('Signal Bandwidth [GHz]'); +ylabel('Fiber length L [km]'); +% X-axis ticks (every 16 GHz starting at 56 GHz) +xticks(56:8:120); +xlim([56,120]) +grid on; box on; + +%% Optional overlay: accumulated dispersion contours +hold on; +[CS, h] = contour(f_GHz, L_km, Dacc_surface, 10, 'k--', 'LineWidth', 0.8); +clabel(CS, h, 'Color','k', 'FontSize',8); +legend('λ_{null} contours','|D_{acc}| [ps/nm]','Location','best'); + +%% ============================================================ +% Helper function: lambda_for_first_null_full +% Stable, single-branch, clamped to O-band +% ============================================================ +function [lambda_vec, Dacc_vec] = lambda_for_first_null_full(f_target, L, lambda0, S0) + c = physconst('lightspeed'); + S0_si = S0 * 1e3; % ps/(nm²·km) -> s/(m³) + + % Define O-band boundaries (in meters) + lambda_min = 1260e-9; + lambda_max = 1360e-9; + + % Force column vector + f_target = f_target(:); + N = numel(f_target); + + lambda_vec = zeros(N,1); + Dacc_vec = zeros(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; + + % Solve within the normal-dispersion range + try + lambda_sol = fzero(fun, [lambda_min, lambda0 * 0.999]); + catch + lambda_sol = lambda_min; + end + + % Clamp to O-band range + lambda_sol = min(max(lambda_sol, lambda_min), lambda_max); + lambda_vec(k) = lambda_sol; + + % 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 + + % Clamp to physical range + Dacc_val = min(max(Dacc_val, -100), 100); + Dacc_vec(k) = Dacc_val; + end +end diff --git a/Functions/Theory/power_fading.m b/Functions/Theory/power_fading.m new file mode 100644 index 0000000..0e14663 --- /dev/null +++ b/Functions/Theory/power_fading.m @@ -0,0 +1,40 @@ +%% ============================================================ +% 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'); diff --git a/Functions/convert_freq_lambda.m b/Functions/convert_freq_lambda.m new file mode 100644 index 0000000..e65f569 --- /dev/null +++ b/Functions/convert_freq_lambda.m @@ -0,0 +1,62 @@ +%% ============================================================ +% Wavelength–Frequency Conversion Utilities +% ============================================================ + +% Example usage: +% f = lambda2freq(1310e-9); % 1310 nm -> Hz +% lambda = freq2lambda(224e12); % 224 THz -> m +% delta_lambda_nm = df2dlambda(224e12, 400e9); % 400 GHz @ 224 THz -> nm +% delta_freq_GHz = dlambda2df(1310e-9, 3.45); % 3.45 nm @ 1310 nm -> GHz + +%% ---- Core conversion functions ---- +function f = lambda2freq(lambda) +% lambda2freq Convert wavelength [m] → frequency [Hz] + c = physconst('lightspeed'); + f = c ./ lambda; +end + +function lambda = freq2lambda(f) +% freq2lambda Convert frequency [Hz] → wavelength [m] + c = physconst('lightspeed'); + lambda = c ./ f; +end + +%% ---- Differential conversions ---- +function d_lambda = df2dlambda(f_center, d_f) +% df2dlambda Convert frequency spacing Δf [Hz] → wavelength spacing Δλ [m] +% around a given center frequency f_center [Hz]. +% Uses first-order differential: Δλ ≈ (c / f^2) * Δf + + c = physconst('lightspeed'); + d_lambda = (c ./ (f_center.^2)) .* d_f; +end + +function d_f = dlambda2df(lambda_center, d_lambda) +% dlambda2df Convert wavelength spacing Δλ [m] → frequency spacing Δf [Hz] +% around a given center wavelength λ_center [m]. +% Uses first-order differential: Δf ≈ (c / λ^2) * Δλ + + c = physconst('lightspeed'); + d_f = (c ./ (lambda_center.^2)) .* d_lambda; +end + +%% ============================================================ +% Example section (can be commented out) +% ============================================================ + +if ~isdeployed + fprintf('--- Example conversions ---\n'); + + lambda_nm = 1310; % nm + lambda = lambda_nm * 1e-9; % m + f = lambda2freq(lambda); % Hz + fprintf('λ = %.1f nm → f = %.3f THz\n', lambda_nm, f/1e12); + + d_f = 2000e9; % 400 GHz spacing + d_lambda = df2dlambda(f, d_f); % [m] + fprintf('Δf = %.0f GHz @ %.1f nm → Δλ = %.3f nm\n', d_f/1e9, lambda_nm, d_lambda*1e9); + + % Verify reverse direction + d_f_back = dlambda2df(lambda, d_lambda); + fprintf('Δλ = %.3f nm @ %.1f nm → Δf = %.0f GHz\n', d_lambda*1e9, lambda_nm, d_f_back/1e9); +end