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