silas/imdd_silas
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

CLEANUP - changes to folder structure

Browse Source
This commit is contained in:
Silas Oettinghaus
2026-03-25 10:57:48 +01:00
parent 0c5ad28f0a
commit 0ae846d3c3
351 changed files with 405 additions and 1294 deletions
Download Patch File Download Diff File Expand all files Collapse all files

89
Theory/Optical/Dispersion/dispersion_10km.m Normal file
View File

@@ -0,0 +1,89 @@
%% ============================================================
% IM/DD Fading Notch – λ_null vs. Bandwidth (Fixed 10 km)
% ============================================================
clear; clc;
%% Fiber and dispersion parameters
lambda0 = 1310e-9; % Zero-dispersion wavelength [m]
S0 = 0.09; % 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
Dacc = Dacc_vec; % [ps/nm]
%% ------------------------------------------------------------
% Plot λ_null vs. f_null for 10 km fiber
% ------------------------------------------------------------
cols = cbrewer2('Paired',10);
figure('Color','w'); hold on;
hLine = plot(lambda_nm, f_GHz, ...
'LineWidth', 2, ...
'DisplayName', sprintf('L = %.1f km', L/1000), ...
'Color', cols(2,:));
xlabel('Wavelength λ [nm]');
ylabel('First fading notch f_{null} [GHz]');
title('IM/DD Fading Notch Position vs. Wavelength');
grid on; box on;
lim = (lambda0.*1e9) - [8, 40];
xlim([lim(2) lim(1)]);
yticks([56,75,90,112]);
%% ------------------------------------------------------------
% Custom DataTip Template
% ------------------------------------------------------------
% Add accumulated dispersion value to the DataTip
hLine.DataTipTemplate.DataTipRows(1).Label = 'λ [nm]';
hLine.DataTipTemplate.DataTipRows(2).Label = 'f_{null} [GHz]';
% Create a new row for Dacc
dRow = dataTipTextRow('D_{acc} [ps/nm]', Dacc);
hLine.DataTipTemplate.DataTipRows(end+1) = dRow;
%% ------------------------------------------------------------
% Helper function: lambda_for_first_null_full
% ------------------------------------------------------------
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³)
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);
fun = @(lambda) -(S0_si/4).*(lambda - (lambda0^4)./(lambda.^3)).*lambda.^2 - RHS;
try
lambda_sol = fzero(fun, [lambda_min, lambda0 * 0.999]);
catch
lambda_sol = lambda_min;
end
lambda_sol = min(max(lambda_sol, lambda_min), lambda_max);
lambda_vec(k) = lambda_sol;
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

77
Theory/Optical/Dispersion/dispersion_contour.m Normal file
View File

@@ -0,0 +1,77 @@
% Festen Betriebsparameter
lambda = 1290; % nm
L = 1; % km
mu_zwd = 1317; % nm
sigma_zwd = 2; % nm
mu_s0 = 0.0872; % ps / nm2 km
sigma_s0 = 0.0012; % ps / nm2 km
rho = -0.5; % Korrelation
% Gitter für lambda0 und S0
lambda0_vec = linspace(mu_zwd-10, mu_zwd+10, 100);
S0_vec = linspace(mu_s0-0.01, mu_s0+0.01, 100);
% Korrigierte meshgrid Reihenfolge
[S0, Lambda0] = meshgrid(S0_vec, lambda0_vec);
% Dispersion berechnen
D = (S0./4) .* ( lambda - (Lambda0.^4)./(lambda^3) ) * L;
%% 2D-Konturplot
figure('Color','w');
hold on
% --- 1. Hintergrund: Dispersions-Konturlinien (Gerundet für TikZ) ---
numLevels_D = 10;
% Erzeuge glatte, auf 1 Nachkommastelle gerundete Werte
levels_D = round(linspace(min(D(:)), max(D(:)), numLevels_D), 1);
levels_D = unique(levels_D); % Falls durch Rundung doppelte Werte entstehen
% Colormap in der Länge der verbliebenen Level erstellen
cmap_bg = cbrewer2('Blues', length(levels_D));
for i = 1:length(levels_D)
% WICHTIG: Das Level als [Wert, Wert] übergeben!
[C,h] = contour(S0, Lambda0, D, [levels_D(i), levels_D(i)], ...
'Color', cmap_bg(i,:),...
'LineWidth', 1.5, ...
'ShowText', 'on', ...
'LabelFormat', '%0.1f');
h.LabelColor = [0,0,0];
end
% --- NEU: Parameter für die Verteilungen ---
mu_zwd = 1317; % nm
sigma_zwd = 2; % nm
mu_s0 = 0.0872; % ps / nm2 km
sigma_s0 = 0.0012; % ps / nm2 km
rho = -0.5; % Korrelation
% --- 3. Randverteilungen (1D Gauss) an den Achsen ---
s0_vals = linspace(min(S0_vec), max(S0_vec), 500);
gauss_s0 = exp(-0.5*((s0_vals - mu_s0)/sigma_s0).^2);
scale_s0 = 4; % Skalierung für die Höhe in der Ansicht
plot(s0_vals, min(lambda0_vec) + gauss_s0 * scale_s0, 'k', 'LineWidth', 2);
zwd_vals = linspace(min(lambda0_vec), max(lambda0_vec), 500);
gauss_zwd = exp(-0.5*((zwd_vals - mu_zwd)/sigma_zwd).^2);
scale_zwd = 0.005; % Skalierung für die Auslenkung in der Ansicht
plot(min(S0_vec) + gauss_zwd * scale_zwd, zwd_vals, 'k', 'LineWidth', 2);
% Hilfslinien für die Mittelwerte
xline(mu_s0, '--k', 'Alpha', 0.4);
yline(mu_zwd, '--k', 'Alpha', 0.4);
% --- Achsenbeschriftung, Titel & Formatierung ---
xlabel('S0 ', 'FontSize', 12);
ylabel('ZDW [nm]', 'FontSize', 12);
% title(sprintf('Dispersion: %d km; %d nm', L, lambda), 'FontSize', 14);
axis([min(S0_vec) max(S0_vec) min(lambda0_vec) max(lambda0_vec)]);
grid on
hold off
%% Für den LaTeX Export
mat2tikz_improved("C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/02_IMDD_System/tikz/dispersion/dispersion_contour.tikz")

38
Theory/Optical/Dispersion/dispersion_first_notch_10km.m Normal file
View File

@@ -0,0 +1,38 @@
%% ------------------------------------------------------------
% Plot: Maximum usable IM/DD bandwidth vs wavelength
% ------------------------------------------------------------
% Fiber and dispersion parameters
lambda0 = 1310e-9; % [m]
S0 = 0.08; % [ps/(nm²·km)]
L = 10000; % [m]
c = physconst('lightspeed');
% Wavelength range around ZDW
lambda_vec = linspace(1250e-9, 1350e-9, 200); % [m]
% Compute D(lambda) using full model
lambda_nm = lambda_vec * 1e9;
lambda0_nm = lambda0 * 1e9;
D_lambda = (S0/4) .* (lambda_nm - (lambda0_nm.^4) ./ (lambda_nm.^3)); % [ps/(nm·km)]
% Convert D to [s/m²]
D_si = D_lambda * 1e-6;
% Compute first null frequency (f₀) for each wavelength
f_null = sqrt(c*(0.5) ./ (abs(D_si).*lambda_vec.^2*L)); % [Hz]
% Plot
figure('Color','w');
plot(lambda_vec*1e9, f_null/1e9, 'LineWidth', 1.6);
grid on; box on;
xlabel('Wavelength [nm]');
ylabel('First Fading Null Frequency [GHz]');
title(sprintf('IM/DD Bandwidth Limit vs. Wavelength (L = %.1f km)', L/1000));
% Highlight useful bandwidth thresholds
yline(25, '--', '25 GHz','Color',[0.4 0.4 0.4],'LabelHorizontalAlignment','left');
yline(50, '--', '50 GHz','Color',[0.2 0.6 0.2],'LabelHorizontalAlignment','left');
yline(100,'--', '100 GHz','Color',[0.6 0.2 0.2],'LabelHorizontalAlignment','left');
legend('First fading notch (f_{null})','Location','best');

165
Theory/Optical/Dispersion/dispersion_power_fading.m Normal file
View File

@@ -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

32
Theory/Optical/Dispersion/dispersion_wavelength_notch.m Normal file
View File

@@ -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
cols = [0.3467 0.5360 0.6907;...
0.9153 0.2816 0.2878;...
0.4416 0.7490 0.4322];
figure('Color','w');hold on
cnt = 1;
for L = [2,10,40]
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+1;
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,60];
xlim([lim(2) lim(1)]);
ylim([10,130])
legend

111
Theory/Optical/Dispersion/dispersion_wdm.m Normal file
View File

@@ -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

68
Theory/Optical/Dispersion/power_fading_gif.m Normal file
View File

@@ -0,0 +1,68 @@
%% ============================================================
% IM/DD Power Fading Evolution GIF (1 km -> 20 km)
% Uses the provided GifWriter (serial mode)
% ============================================================
clear; close all; clc;
%% Fiber and system parameters
lambda0 = 1310e-9; % zero-dispersion wavelength [m]
lambda = 1275e-9; % operating wavelength [m]
S0 = 0.09; % dispersion slope [ps/(nm^2·km)]
c = physconst('lightspeed');
%% Derived quantities (length-independent)
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
%% Frequency grid
f_max = 150e9;
f = linspace(0, f_max, 4000); % [Hz]
%% Figure setup (keep it stable for nicer GIFs)
fig = figure('Color','w');
ax = axes(fig); %#ok<LAXES>
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 (serial mode; simplest)
g = GifWriter('Name','power_fading_evolution', 'DelayTime',0.12, 'Parallel',false);
%% Loop: 1 km to 20 km
L = [1:20,19:-1:1];
for L_km = L
L_meter = L_km * 1e3; % [m]
% IM/DD transfer function (power fading)
phi = 2*pi^2 * b2 * f.^2 * L_meter;
H = abs(cos(phi));
HdB = 10*log10(max(H, 1e-12)); % avoid -Inf for deep notches
% Clear and redraw (stable axes)
cla(ax);
plot(ax, f/1e9, HdB, 'LineWidth', 1.8, 'Color','black');
% Analytic first-null frequency marker
f_null = sqrt(c*(0.5)/(abs(D_si)*lambda^2*L_meter));
xline(ax, f_null/1e9, 'r--', 'LineWidth', 1.2, ...
'Label', sprintf('f_{null}=%.1f GHz', f_null/1e9), ...
'LabelOrientation','horizontal', ...
'LabelVerticalAlignment','bottom');
title(ax, sprintf('Power Fading for: %.0f km @ 1275 nm', L_km));
drawnow;
% Add frame to GIF
g.addFrame(fig);
end
%% Done
g.compile(fig.Number);
disp(fullfile(g.OutputDir, sprintf('%s_fig_%d.gif', g.Name, fig.Number)));

69
Theory/Optical/Dispersion/power_fading_gif_lambda_variation.m Normal file
View File

@@ -0,0 +1,69 @@
%% ============================================================
% 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)));

116
Theory/Optical/PMD_maxwellshape/dgd_considerations.m Normal file
View File

@@ -0,0 +1,116 @@
if 0
pmd = 0.2 * ( 1e-12 / sqrt(1e3) ) ; % 0.1 ps/sqrt(km) -> 1e-9 -> s/sqrt(m)
L = 1000; % m
corr_len = 100;
num_wave_plates = 100;
wp_len = L / num_wave_plates; %waveplate length
dgd = pmd * sqrt(L);
fdac = 120e9; % GHz
fsim = fdac * 512 ; %oversampled simulation frequency
dt = 1/fsim; % sample time
nt = 4; % sampled signal length
omega = 2*pi*[(0:nt/2-1),(-nt/2:-1)]/(dt*nt) ; %angular frequency vector for optical signal of length nt
omega = 2*pi*fsim ; %angular frequency vector for optical signal of length nt
for rlz = 1:20000
db0 = [];
db1 = [];
delta_beta = [];
db0 = (rand(num_wave_plates,1)*2*pi - pi); % normal distr. between -pi <-> +pi
db1 = sqrt(3*pi/8)*(dgd/fsim)/ num_wave_plates .* omega; % linear increasing delta beta 1
%loop over waveplates
for n_wp = 1:length(db0)
delta_beta(n_wp,:) = (db1+db0(n_wp))./corr_len;
end
dphi = delta_beta*wp_len; %phase shift due to propagation constant = beta * L
if rlz == 1
figure()
plot(cumsum( delta_beta ));
end
dphi_end(rlz) = sum( dphi );
end
figure;
histogram((dphi_end),100);
deltaT = mean(abs(dphi_end));
%H = exp(-1j*dphi); % Filter to apply phase shift in freq. domain
end
PMD = 0.1 * ( 1e-12 / sqrt(1e3) ); %PMD s/sqrt(m)
L = 10e3; %m
sigma_dgd = PMD * sqrt(L); % in sec.
disp(['variance of DGD in pico seconds: ',num2str(sigma_dgd*1e12)]);
disp(['variance of DGD in pico seconds: ',num2str(sigma_dgd*1e12)]);
dgd = [];
for i = 1:4
dgd(i,:) = DGD(sigma_dgd,10000);
end
rms = PMD * sqrt(40e3) * 1e12 ; % in sec.
t = [0:0.01:100];
z2 = sqrt(2/pi).*(t.^2)/(rms.^3).*exp(-t.^2/(2*rms^2)) ;
figure;
histogram(dgd*1e12,1000,"Normalization","pdf");
hold on
plot(t(1:200),z2(1:200),'r');
hold off
function z = DGD(sigma_dgd,N)
%PMD = sqrt( (PMD*1e12)^2/3 );
% PMD = PMD * 1/( 1e-12 / sqrt(1e3) );
sigma_dgd = sigma_dgd*1e12;
x = [0:N];
almaga = randn(3,N)*sigma_dgd;
y = almaga.^2;
z0 = sqrt(sum(y,1));
disp(mean(z0));
disp(std(z0));
% disp(var(z0));
z1 = [];
for l = 0:0.1:N
z1 = [z1 sum(z0(1,:)>l & z0(1,:)<l+0.1) / N ];
end
z = z0 *1e-12;
rms = sigma_dgd;
t = [0:0.01:N];
z2 = sqrt(2/pi).*(t.^2)/(rms.^3).*exp(-t.^2/(2*rms^2)) ;
figure;
histogram(z0,1000,"Normalization","pdf");
hold on
plot(t(1:200),z2(1:200),'r');
return
end

59
Theory/Optical/PMD_maxwellshape/pmd_maxwell.m Normal file
View File

@@ -0,0 +1,59 @@
PMDcoeff = 0.1 * ( 1e-12 / sqrt(1e3) ); %PMD s/sqrt(m)
L_40 = 40e3; %m
mean_dgd = PMDcoeff * sqrt(L_40); % in sec. E(tau) == PMD value
disp(['defined mean of DGD in pico seconds: ',num2str(mean_dgd*1e12)]);
dgd_40km = [];
for i = 1:4000
dgd_40km(i,:) = getDGDrealization(mean_dgd,1);
end
disp(['simulated mean of DGD in pico seconds: ',num2str(mean(dgd_40km)*1e12)]);
%%% 2 test for segmented link %%%%
L_10 = 10e3; %m
mean_dgd = PMDcoeff * sqrt(L_10); % in sec. E(tau) == PMD value
dgd_10km = [];
for j = 1:4
for i = 1:1000
dgd_10km(i,j) = getDGDrealization(mean_dgd,1);
end
end
a = sum(dgd_10km,2) ;
%%% calc maxwell curve
mean_dgd = PMDcoeff * sqrt(L_40) * 1e12 ; % in sec.
q = sqrt(pi/8) * mean_dgd;
t = [0:0.01:100];
maxwell = sqrt(2/pi).* (t.^2)/(q.^3) .* exp(-t.^2/(2*q^2)) ;
%%% plot everything
figure;
histogram(dgd_40km*1e12,1000,"Normalization","pdf",'EdgeColor','none','FaceColor','b','FaceAlpha',0.6);
hold on
histogram(a*1e12,1000,"Normalization","pdf",'EdgeColor','none','FaceColor','g','FaceAlpha',0.6);
plot(t(1:200),maxwell(1:200),'r','LineWidth',3);
xline(mean_dgd,'LineWidth',3,'Color','magenta')
hold off
function z = getDGDrealization(mean_dgd,N)
mean_dgd = mean_dgd; %to pico seconds
q = sqrt(pi/8) * mean_dgd;
threenormaldists = randn(3,N)*q; %3x normal dist around [-1,1] with a deviation of mean_dgd
y = threenormaldists.^2;
z = sqrt(sum(y,1));
end

69
Theory/Optical/calcFWM/FWM_products.m Normal file
View File

@@ -0,0 +1,69 @@
w0 = [1290:2:1290+15*2]';
w0 = [1290:2:1290+15*2]';
w0 = [ 1302 1304 1306 1308]';
%w0 = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16]';
m = 0.5*(numel(w0)^3 - numel(w0)^2);
a = [1,1,1]';
w = w0;
for o = 2:3
p = nchoosek(w,3);
q = [];
parfor i = 1:size(p,1)
q_ = perms(p(i,:));
q = [q;q_];
end
p = q;
%p = unique(q,"rows");
w_ = p(:,1) + p(:,2) - p(:,3);
a_ = ones(size(w_)).* 1/o;
w = [w ; w_];
a = [a ; a_];
% w = unique(w);
% a = unique(a);
m(end+1) = 0.5*(numel(w)^3 - numel(w)^2);
end
figure(11)
hold on
lambda = min(w):max(w);
gen = sum(lambda == w,1);
gen(gen==0) = NaN;
stem(lambda,gen,"filled",'LineWidth',1,'Marker','o','MarkerSize',2,'LineStyle',':')
initial = sum(lambda == w0,1);
initial(initial==0) = NaN;
stem(lambda,initial,"filled",'LineWidth',1.5,'MarkerSize',5,'Marker','^');
xlabel('Wavelength');
ylabel('number of FWM products');
grid minor
legend('Generated Products', 'Initial Channel Position')
AxesMain = gca;
fig = gcf;
fontsize(AxesMain,8,"points")
fig.Units = "centimeters";
fig.Position = [2 2 8.5 7];

43
Theory/Optical/calcFWM/JLT_statistics_laser_and_zdw.m Normal file
View File

@@ -0,0 +1,43 @@
%% Laser Offset Statistics
figure
for i = 1
res = 1.7e6;
n_chann = 16;
df_T_exact = (-n_chann/2+0.5:n_chann/2).* 200e9;
for key = 1:100
laser_frequency_imperfection(key,:) = res .* round(randn(1,n_chann)*i*100);
df_T(key,:) = df_T_exact + laser_frequency_imperfection(key,:);
end
hold on
histogram(laser_frequency_imperfection.*1e-6,100,"Normalization","probability","EdgeColor","none","FaceAlpha",0.3,'DisplayName',['Std. Dev.: ',num2str(mean(std(laser_frequency_imperfection))*1e-6),' MHz']);
xlabel('Laser Frequency Offset in MHz');
ylabel('Probability');
title(['Laser deviations from exact grid.'])
end
%% ZDW Statistics
for k = 1:1000
% Set parameters
meanUniformMin = 1309;
meanUniformMax = 1315;
meanValue = 1310;
sigma = 2;
% Seed the random number generator (assuming Mersenne Twister)
rng(k);
% Generate normally distributed random numbers
randomNumbers(k,:) = normrnd(meanValue, sigma, [n_chann, 1]).';
end
figure;
histogram(randomNumbers,100,"Normalization","probability","EdgeColor","none");
xlabel('ZDW in nm');
ylabel('Probability');

65
Theory/Optical/calcFWM/analytical_calculation_paper.m Normal file
View File

@@ -0,0 +1,65 @@
%%FWM analysis from "Analytical Calculation of the Number of
%%Four-Wave-Mixing Products in Optical Multichannel Communication Systems"
N_ = [4,8,16];
df_hz = 400e9;
center_nm = 1310;
figure()
for i = 1:length(N_)
vec = (2*N_(i)-1:-1:2-N_(i)) -(N_(i)/2+0.5);
channelplan_hz = nm2hz(center_nm) + (vec * df_hz) ;
channelplan_nm = hz2nm(channelplan_hz);
[Mndg,Mdg] = getProducts(N_(i));
total(i) = sum(Mndg) + sum(Mdg);
subplot(1,length(N_),i)
xline(calcWavelengthPlan(N_(i), df_hz, center_nm));
hold on
stem(channelplan_nm,(Mndg+Mdg),'filled','LineWidth',1,'Marker','o','MarkerSize',2)
stem(channelplan_nm,(Mdg),'filled','LineWidth',1,'Marker','none');
ylim([0,100])
xlim([1260, 1365]);
grid off
xlabel('O-band wavelength region in nm');
ylabel('Number of FWM products');
title([num2str(N_(i)),' ch.'])
end
function [Mndg,Mdg] = getProducts(N)
s = abs(2-N-1) ;
for n = 2-N:2*N-1
if n<-N
Mdg(n+s) = 0;
elseif (-N <= n)&&(n <= 0)
Mdg(n+s) = N - ceil((N-n)/2);
elseif (1 <= n)&&(n <= N)
Mdg(n+s) = N - 1 - floor(n/2) - ceil((N-n)/2);
elseif (N < n)&&(n <= 2*N)
Mdg(n+s) = N - floor(n/2);
elseif n > 2*N
Mdg(n+s) = 0;
end
if n<-N
Mndg(n+s) = 0;
elseif (-N <= n)&&(n < 1)
Mndg(n+s) = ceil((N^2 + n^2 - 2*N - 2*n + 2*N*n)/4);
elseif (1 <= n)&&(n <= N)
Mndg(n+s) = ceil(((N^2 - 6*N - 2*n^2 + 2*n + 4)/4) + floor((N*n)/2));
elseif (N < n)&&(n <= 2*N)
Mndg(n+s) = floor(N^2 + n^2 /4 - N*n);
elseif n > 2*N
Mndg(n+s) = 0;
end
end
end

23
Theory/Optical/calcFWM/calcFwmEfficiency.m Normal file
View File

@@ -0,0 +1,23 @@
function [eta, deltaBeta] = calcFwmEfficiency(f_i, f_j, f_k, f_0, Ds, alphaDbPerKm, Lkm)
% FWM efficiency including phase mismatch and attenuation.
% alphaDbPerKm is the power attenuation in dB/km, Lkm is the fiber length in km.
deltaBeta = calcPhaseMatching(f_i, f_j, f_k, f_0, Ds);
alphaNpPerM = alphaDbPerKm .* log(10) ./ 10 ./ 1e3;
Lm = Lkm .* 1e3;
denominator = alphaNpPerM.^2 + deltaBeta.^2;
term1 = alphaNpPerM.^2 ./ denominator;
loss_term = 1 - exp(-alphaNpPerM .* Lm);
if abs(alphaNpPerM) < eps
term2 = 4 .* sin(deltaBeta .* Lm ./ 2).^2 ./ max((alphaNpPerM .* Lm).^2, eps);
else
term2 = 1 + 4 .* exp(-alphaNpPerM .* Lm) .* sin(deltaBeta .* Lm ./ 2).^2 ./ (loss_term.^2);
end
eta = term1 .* term2;
end

48
Theory/Optical/calcFWM/calcFwmPower.m Normal file
View File

@@ -0,0 +1,48 @@
function [P_fwm, eta, deltaBeta, Leff] = calcFwmPower( ...
f_i, f_j, f_k, f_0, Ds, alphaDbPerKm, Lkm, ...
P_i, P_j, P_k, gammaWInvKmInv, degeneracyFactor)
% Calculate FWM power for a fiber with attenuation and phase mismatch.
%
% Inputs:
% f_i, f_j, f_k, f_0 : frequencies in Hz
% Ds : dispersion slope in ps / (nm^2 km)
% alphaDbPerKm : attenuation in dB/km
% Lkm : fiber length in km
% P_i, P_j, P_k : launch powers in W
% gammaWInvKmInv : nonlinear coefficient in 1/(W km)
% degeneracyFactor : typically 3 for degenerate FWM, 6 for non-degenerate
if nargin < 8 || isempty(P_i)
P_i = 1;
end
if nargin < 9 || isempty(P_j)
P_j = P_i;
end
if nargin < 10 || isempty(P_k)
P_k = 1;
end
if nargin < 11 || isempty(gammaWInvKmInv)
gammaWInvKmInv = 1;
end
if nargin < 12 || isempty(degeneracyFactor)
degeneracyFactor = 1;
end
[eta, deltaBeta] = calcFwmEfficiency(f_i, f_j, f_k, f_0, Ds, alphaDbPerKm, Lkm);
alphaNpPerM = alphaDbPerKm .* log(10) ./ 10 ./ 1e3;
Lm = Lkm .* 1e3;
gammaWInvMInv = gammaWInvKmInv ./ 1e3;
if abs(alphaNpPerM) < eps
Leff = Lm;
else
Leff = (1 - exp(-alphaNpPerM .* Lm)) ./ alphaNpPerM;
end
P_fwm = degeneracyFactor .* eta .* ...
(gammaWInvMInv .* Leff).^2 .* ...
P_i .* P_j .* P_k .* ...
exp(-alphaNpPerM .* Lm);
end

20
Theory/Optical/calcFWM/calcPhaseMatching.m Normal file
View File

@@ -0,0 +1,20 @@
function deltaBeta = calcPhaseMatching(f_i, f_j, f_k, f_0, Ds)
% Approximate phase mismatch for degenerate FWM close to the ZDW.
% Inputs are frequencies in Hz.
% f_i, f_j : pump frequencies (equal in the degenerate case)
% f_k : signal frequency
% f_0 : zero-dispersion frequency
% Ds : dispersion slope in ps / (nm^2 km)
c = physconst('LightSpeed');
pump_frequency = 0.5 .* (f_i + f_j);
lambda_zdw_m = c ./ f_0;
dispersion_slope_si = Ds .* 1e3;
deltaBeta = -(2 .* pi .* lambda_zdw_m.^4 ./ c.^2) .* ...
dispersion_slope_si .* ...
(pump_frequency - f_0) .* ...
(pump_frequency - f_k).^2;
end

104
Theory/Optical/calcFWM/scriptFWM.m Normal file
View File

@@ -0,0 +1,104 @@
clear;
clc;
% Sweep the degenerate pump frequency around its nominal wavelength.
pump_detuning_hz = (-800:0.01:800) .* 1e9;
% Degenerate FWM setup: two pump photons at f_p and one signal at f_s
% generate an idler at f_i = 2*f_p - f_s.
pump_wavelength_nm = 1310;
signal_wavelength_nm = 1308;
zdw_wavelength_nm = 1310;
f_pump_nominal = wavelength2frequency(pump_wavelength_nm, 'nm');
f_signal_scalar = wavelength2frequency(signal_wavelength_nm, 'nm');
f_zdw_scalar = wavelength2frequency(zdw_wavelength_nm, 'nm');
f_pump = f_pump_nominal + pump_detuning_hz;
f_signal = f_signal_scalar .* ones(size(f_pump));
f_zdw = f_zdw_scalar .* ones(size(f_pump));
f_idler = 2 .* f_pump - f_signal;
% Fiber parameters
dispersion_slope_ps_nm2_km = 0.07;
attenuation_db_per_km = 0.21;
fiber_length_km = 10;
% Launch powers and nonlinear coefficient
pump_power_dbm = 10;
signal_power_dbm = 10;
pump_power_w = dbm2watt(pump_power_dbm);
signal_power_w = dbm2watt(signal_power_dbm);
gamma_w_inv_km_inv = 1.3;
degeneracy_factor = 3;
[P_fwm, eta, delta_beta, L_eff_m] = calcFwmPower( ...
f_pump, f_pump, f_signal, f_zdw, ...
dispersion_slope_ps_nm2_km, attenuation_db_per_km, fiber_length_km, ...
pump_power_w, pump_power_w, signal_power_w, ...
gamma_w_inv_km_inv, degeneracy_factor);
f_pump_thz = f_pump .* 1e-12;
f_zdw_thz = f_zdw_scalar .* 1e-12;
f_signal_thz = f_signal_scalar .* 1e-12;
idler_power_dbm = 10 .* log10(max(P_fwm, realmin) ./ 1e-3);
figure;
tiledlayout(2,1);
ax1 = nexttile;
plot(ax1, f_pump_thz, eta, 'LineWidth', 2);
hold(ax1, 'on');
xline(ax1, f_zdw_thz, '--r', 'ZDW', 'LineWidth', 1.2, ...
'LabelOrientation', 'horizontal', 'LabelVerticalAlignment', 'bottom');
xline(ax1, f_signal_thz, '--k', 'Signal', 'LineWidth', 1.2, ...
'LabelOrientation', 'horizontal', 'LabelVerticalAlignment', 'middle');
ylabel(ax1, 'FWM efficiency');
grid(ax1, 'on');
title(ax1, 'FWM Efficiency and Idler Power versus Pump Frequency');
ax2 = nexttile;
plot(ax2, f_pump_thz, idler_power_dbm, 'LineWidth', 2);
hold(ax2, 'on');
xline(ax2, f_zdw_thz, '--r', 'ZDW', 'LineWidth', 1.2, ...
'LabelOrientation', 'horizontal', 'LabelVerticalAlignment', 'bottom');
xline(ax2, f_signal_thz, '--k', 'Signal', 'LineWidth', 1.2, ...
'LabelOrientation', 'horizontal', 'LabelVerticalAlignment', 'middle');
xlabel(ax2, 'Pump frequency (THz)');
ylabel(ax2, 'FWM idler power (dBm)');
grid(ax2, 'on');
fprintf('Pump wavelength : %.3f nm -> %.6f THz\n', ...
pump_wavelength_nm, f_pump_nominal .* 1e-12);
fprintf('Signal wavelength : %.3f nm -> %.6f THz\n', ...
signal_wavelength_nm, f_signal_scalar .* 1e-12);
fprintf('ZDW wavelength : %.3f nm -> %.6f THz\n', ...
zdw_wavelength_nm, f_zdw_scalar .* 1e-12);
fprintf('Pump launch power : %.2f dBm -> %.4g W\n', ...
pump_power_dbm, pump_power_w);
fprintf('Signal launch power : %.2f dBm -> %.4g W\n', ...
signal_power_dbm, signal_power_w);
fprintf('Peak FWM efficiency : %.4g\n', max(eta));
fprintf('Peak FWM idler power : %.4g W (%.2f dBm)\n', ...
max(P_fwm), 10 .* log10(max(P_fwm) ./ 1e-3));
fprintf('Effective fiber length : %.4f km\n', L_eff_m ./ 1e3);
fprintf('Idler wavelength range : %.3f nm to %.3f nm\n', ...
min(frequency2wavelength(f_idler, 'nm')), max(frequency2wavelength(f_idler, 'nm')));
fprintf('Max |delta beta| : %.4g 1/m\n', max(abs(delta_beta)));
function wavelength = frequency2wavelength(frequency, outputUnit)
c = physconst('LightSpeed');
wavelength = c ./ frequency;
switch lower(outputUnit)
case 'm'
case 'nm'
wavelength = wavelength .* 1e9;
otherwise
error('Unsupported output unit "%s". Use "m" or "nm".', outputUnit);
end
end
function power_w = dbm2watt(power_dbm)
power_w = 1e-3 .* 10.^(power_dbm ./ 10);
end

124
Theory/Optical/calcFWM/validate_fwm.m Normal file
View File

@@ -0,0 +1,124 @@
%% Validate the analytical FWM product count against a brute-force reference
% The paper counts channel combinations, not only unique output frequencies:
% non-degenerate: i < j, k ~= i, k ~= j, n = i + j - k
% degenerate: i == j, k ~= i, n = 2*i - k
clear;
clc;
N_values = [4, 8, 16];
plot_N = 8;
fprintf('Validating analytical FWM product count from Goebel and Hanik (2008)\n');
for idxN = 1:numel(N_values)
N = N_values(idxN);
fprintf('\nN = %d\n', N);
[Mndg_ana, Mdg_ana] = getProducts(N);
[Mndg_brute, Mdg_brute, n_values] = getProductsBruteForce(N);
diff_ndg = Mndg_ana - Mndg_brute;
diff_dg = Mdg_ana - Mdg_brute;
fprintf(' Analytical total : %d\n', sum(Mndg_ana) + sum(Mdg_ana));
fprintf(' Brute-force total : %d\n', sum(Mndg_brute) + sum(Mdg_brute));
if all(diff_ndg == 0) && all(diff_dg == 0)
fprintf(' Match : yes\n');
else
fprintf(' Match : no\n');
fprintf(' Non-degenerate diff: %s\n', mat2str(diff_ndg));
fprintf(' Degenerate diff : %s\n', mat2str(diff_dg));
end
if N == plot_N
plotComparison(n_values, Mndg_ana, Mdg_ana, Mndg_brute, Mdg_brute, N);
end
end
function [Mndg, Mdg, n_values] = getProductsBruteForce(N)
n_values = (2 - N):(2*N - 1);
Mndg = zeros(size(n_values));
Mdg = zeros(size(n_values));
for idx = 1:numel(n_values)
n = n_values(idx);
% Degenerate products: two identical pumps and one different channel.
for i = 1:N
k = 2*i - n;
if isValidChannel(k, N) && (k ~= i)
Mdg(idx) = Mdg(idx) + 1;
end
end
% Non-degenerate products: unordered pump pair plus one third channel.
for i = 1:N
for j = (i + 1):N
k = i + j - n;
if isValidChannel(k, N) && (k ~= i) && (k ~= j)
Mndg(idx) = Mndg(idx) + 1;
end
end
end
end
end
function tf = isValidChannel(channel_idx, N)
tf = (channel_idx >= 1) && (channel_idx <= N) && (channel_idx == round(channel_idx));
end
function plotComparison(n_values, Mndg_ana, Mdg_ana, Mndg_brute, Mdg_brute, N)
figure;
subplot(1,2,1);
stem(n_values, Mndg_ana + Mdg_ana, 'filled', 'LineWidth', 1, 'Marker', 'o', 'MarkerSize', 2);
hold on;
stem(n_values, Mdg_ana, 'filled', 'LineWidth', 1, 'Marker', 'none');
title(['Analytical (N=', num2str(N), ')']);
xlabel('Product index n');
ylabel('Number of FWM products');
legend('Total', 'Degenerate');
grid on;
subplot(1,2,2);
stem(n_values, Mndg_brute + Mdg_brute, 'filled', 'LineWidth', 1, 'Marker', 'o', 'MarkerSize', 2);
hold on;
stem(n_values, Mdg_brute, 'filled', 'LineWidth', 1, 'Marker', 'none');
title(['Brute force (N=', num2str(N), ')']);
xlabel('Product index n');
ylabel('Number of FWM products');
legend('Total', 'Degenerate');
grid on;
end
function [Mndg, Mdg] = getProducts(N)
s = abs(2 - N - 1);
for n = 2 - N:2*N - 1
if n < -N
Mdg(n + s) = 0;
elseif (-N <= n) && (n <= 0)
Mdg(n + s) = N - ceil((N - n)/2);
elseif (1 <= n) && (n <= N)
Mdg(n + s) = N - 1 - floor(n/2) - ceil((N - n)/2);
elseif (N < n) && (n <= 2*N)
Mdg(n + s) = N - floor(n/2);
elseif n > 2*N
Mdg(n + s) = 0;
end
if n < -N
Mndg(n + s) = 0;
elseif (-N <= n) && (n < 1)
Mndg(n + s) = ceil((N^2 + n^2 - 2*N - 2*n + 2*N*n)/4);
elseif (1 <= n) && (n <= N)
Mndg(n + s) = ceil(((N^2 - 6*N - 2*n^2 + 2*n + 4)/4) + floor((N*n)/2));
elseif (N < n) && (n <= 2*N)
Mndg(n + s) = floor(N^2 + n^2/4 - N*n);
elseif n > 2*N
Mndg(n + s) = 0;
end
end
end

18
Theory/Optical/calcFWM/wavelength2frequency.m Normal file
View File

@@ -0,0 +1,18 @@
function frequency = wavelength2frequency(wavelength, inputUnit)
% Convert wavelength to optical frequency.
% Supported units: m, nm.
c = physconst('LightSpeed');
switch lower(inputUnit)
case 'm'
wavelength_m = wavelength;
case 'nm'
wavelength_m = wavelength .* 1e-9;
otherwise
error('Unsupported input unit "%s". Use "m" or "nm".', inputUnit);
end
frequency = c ./ wavelength_m;
end