Commit | Line | Data |
---|---|---|
5117fc58 AIL |
1 | function [x, y] = ssfs(xin, yin, D, lambda, z, dz, Tsamp, gamma, alpha) |
2 | %% Split-step Fourier solver, simulates chromatic dispersion, | |
3 | %% attenuation, Kerr effect, and power splitting / amplification. | |
4 | %% Params: | |
5 | %% - xin, yin: input waveform (x and y polarizations) | |
6 | %% - D: dispersion coefficient (ps / (nm km)) | |
7 | %% - lambda: wavelength (nm) | |
8 | %% - z: length of fibre (km) | |
9 | %% - dz: step size (km) | |
10 | %% - Tsamp: sampling time (s) | |
11 | %% - gamma: Non-linear coefficient (W^-1 / km) | |
12 | %% - alpha: attenuation (dB / km) | |
13 | %% Output: | |
14 | %% - x, y: output waveform (both polarizations) | |
15 | ||
16 | %% Convert everything to SI base units | |
17 | c = 299792458; % m/s | |
18 | D = D * 1e-6; % s/m^2 | |
19 | lambda = lambda * 1e-9; % m | |
20 | z = z * 1e3; % m | |
21 | gamma = gamma * 1e-3; % watt^-1 / m | |
22 | dz = dz * 1e3; % m | |
23 | alpha = alpha / 10 * log(10); % Np/km | |
24 | alpha = alpha * 1e-3; % Np/m | |
25 | ||
26 | stepnum = z / dz; | |
27 | ||
28 | %% Frequency response of CD | |
29 | n = length(xin); | |
30 | fs = 1 / Tsamp; | |
31 | omega = (2*pi * fs / n * [(0 : floor((n-1)/2)), (-ceil((n-1)/2) : -1)]).'; | |
32 | dispDFT = exp(-1j * omega.^2 * D * lambda^2 * dz / (4 * pi * c)); | |
33 | ||
34 | %% Convenient variables to reduce typing | |
35 | hhz = 1j * (8/9) * gamma * dz; % Kerr phase shift per power | |
36 | attn = -alpha * dz / 2; % attenuation | |
37 | ||
38 | %% Initial Kerr half step | |
39 | P = abs(xin) .^ 2 + abs(yin) .^ 2; | |
40 | x = xin .* exp(P .* hhz / 2 + attn / 2); | |
41 | y = yin .* exp(P .* hhz / 2 + attn / 2); | |
42 | ||
43 | for i = 1 : stepnum | |
44 | %% CD in frequency domain | |
45 | xDFT = fft(x); | |
46 | yDFT = fft(y); | |
47 | x = ifft(xDFT .* dispDFT); | |
48 | y = ifft(yDFT .* dispDFT); | |
49 | ||
50 | %% Kerr effect in time domain | |
51 | P = abs(x) .^ 2 + abs(y) .^ 2; | |
52 | x = x .* exp(P .* hhz + attn); | |
53 | y = y .* exp(P .* hhz + attn); | |
54 | ||
55 | %% Power split after 40km | |
56 | if i * dz == 40e3 | |
57 | splitnum = 1024; % energy factor, amplitude factor is sqrt of this | |
58 | x = x ./ sqrt(splitnum); | |
59 | y = y ./ sqrt(splitnum); | |
60 | %% Splitter loss - 1:4 coupler * 5 levels, 0.3 dB per level | |
61 | %% so total loss is 1.5 dB | |
62 | x = x ./ sqrt(10 ^ 0.15); | |
63 | y = y ./ sqrt(10 ^ 0.15); | |
64 | end | |
65 | end | |
66 | %% Final Kerr effect has overshot by half step, so cancel this | |
67 | P = abs(x) .^ 2 + abs(y) .^ 2; | |
68 | x = x .* exp(-P .* hhz / 2 - attn / 2); | |
69 | y = y .* exp(-P .* hhz / 2 - attn / 2); | |
70 | end |