gnlse.m

function [Z, AT, AW, W] = gnlse(T, A, w0, gamma, betas, ...
                             loss, fr, RT, flength, nsaves)
% Propagate an optical field using the generalised NLSE
% This code integrates Eqs. (3.13), (3.16) and (3.17).
% For usage see the exampe of test_Dudley.m (below)
% Written by J.C. Travers, M.H. Frosz and J.M. Dudley (2009)
% Please cite this chapter in any publication using this code.
% Updates to this code are available at www.scgbook.info
% Minor bug fix: 16/03/2020

n = length(T); dT = T(2)-T(1); % grid parameters
V = 2*pi*(-n/2:n/2-1)'/(n*dT); % frequency grid
alpha = log(10.^(loss/10));    % attenuation coefficient

B = 0;
for i = 1:length(betas)        % Taylor expansion of betas
  B = B + betas(i)/factorial(i+1).*V.^(i+1);
end
L = 1i*B - alpha/2;            % linear operator

if abs(w0) > eps               % if w0>0 then include shock
    gamma = gamma/w0;    
    W = V + w0;                % for shock W is true freq
else
    W = 1;                     % set W to 1 when no shock
end

RW = n*ifft(fftshift(RT.'));   % frequency domain Raman
L = fftshift(L); W = fftshift(W); % shift to fft space

% === define function to return the RHS of Eq. (3.13)
function R = rhs(z, AW)
  AT = fft(AW.*exp(L*z));         % time domain field
  IT = abs(AT).^2;                % time domain intensity
  if (length(RT) == 1) || (abs(fr) < eps) % no Raman case
    M = ifft(AT.*IT);             % response function
  else
    RS = dT*fr*fft(ifft(IT).*RW); % Raman convolution
    M = ifft(AT.*((1-fr).*IT + RS));% response function
  end
  R = 1i*gamma*W.*M.*exp(-L*z);   % full RHS of Eq. (3.13)
end

% === define function to print ODE integrator status
function status = report(z, y, flag) % 
  status = 0;
  if isempty(flag)
    fprintf('%05.1f %% complete\n', z/flength*100);
  end
end

% === setup and run the ODE integrator
Z = linspace(0, flength, nsaves);  % select output z points
% === set error control options
options = odeset('RelTol', 1e-5, 'AbsTol', 1e-12, ...
                 'NormControl', 'on', ...
                 'OutputFcn', @report);
[Z, AW] = ode45(@rhs, Z, ifft(A), options); % run integrator

% === process output of integrator
AT = zeros(size(AW(1,:)));
for i = 1:length(AW(:,1))
  AW(i,:) = AW(i,:).*exp(L.'*Z(i)); % change variables
  AT(i,:) = fft(AW(i,:));           % time domain output
  AW(i,:) = fftshift(AW(i,:)).*dT*n;  % scale
end

W = V + w0; % the absolute frequency grid
end