ssfs.m

   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