The scalar version of the SSPROP solves the scalar nonlinear Schrödinger equation:


using the symmetrized split-step Fourier method. It allows for almost arbitrary specification of the dispersion and loss, and it includes a simple approximation for the delayed Raman response time and the optical self-steepening effect.


A summary of the syntax and usage can be obtained from Matlab by typing “help ssprop” or “help sspropc“.

The compiled mex file (sspropc) can be invoked from Matlab using one of the following forms:

u1 = sspropc(u0,dt,dz,nz,alpha,betap,gamma);
u1 = sspropc(u0,dt,dz,nz,alpha,betap,gamma,tr);
u1 = sspropc(u0,dt,dz,nz,alpha,betap,gamma,tr,to);
u1 = sspropc(u0,dt,dz,nz,alpha,betap,gamma,tr,to,maxiter);
u1 = sspropc(u0,dt,dz,nz,alpha,betap,gamma,tr,to,maxiter,tol);

The last three arguments assume a default value if they are left unspecified.

sspropc may also be invoked with a single input argument, to specify options specific to the FFTW routines (discussed below):

sspropc -option

The corresponding Matlab m-file can be invoked using a similar syntax by replacing sspropc with ssprop. Note however, that the m-file version does not include the delayed Raman response time or self-steepening terms, and the m-file version does not allow for control of the FFTW options.

u1 = ssprop(u0,dt,dz,nz,alpha,betap,gamma);
u1 = ssprop(u0,dt,dz,nz,alpha,betap,gamma,maxiter);
u1 = ssprop(u0,dt,dz,nz,alpha,betap,gamma,maxiter,tol);

Input Arguments

vector (N) Input optical field, specified as a length-N vector time sequence. u0 represents the complex, slowly-varying envelope of the optical field. u0 should be normalized so that |u0|^2 is the optical power.
dt scalar The time increment between adjacent points in the vector u0.
dz scalar The step-size to use for propagation
nz scalar (int) The number of steps to take. The total distance propagated is therefore L = nz*dz
alpha scalar or vector (N) The linear power attenuation coefficient. Specifically, if the medium is excited with CW light, the power should fall with distance according to P(z) = P(0)exp(-α*z). Note: this implies that the fields decay as exp(-α*z/2).The loss coefficient alpha may optionally be specified as a vector of the same length as u0, in which case it is treated as vector that describes a wavelength-dependent loss coefficient α(ω) in the frequency domain. (The function wspace.m in the tools subdirectory can be used to construct a vector with the corresponding frequencies.)
betap vector A real-valued vector specifying the dispersion properties of the medium. The dispersion can be specified to any polynomial order by using a betap vector with the appropriate length. In most cases the nonlinear Schrödinger equation is solved with a slowly-varing envelope approximation in a reference frame that is moving at the group velocity. In this case, the first two elements of the betap vector (β0 and β1) should be set to zero. The first non-zero term, betap(3), would then represent β2, the group-velocity dispersion.The propagation constant can also be specified directly by replacing the polynomial argument betap with a vector of the same length as u0. In this case, the argument betap is treated as a vector describing propagation constant β(ω) in the frequency domain. (The function wspace.m in the tools subdirectory can be used to construct a vector with the corresponding frequencies.)
gamma scalar A real number that describes the nonlinear coefficient of the fiber, which is related to the mode effective area and the nonlinear refractive index n2.
tr scalar The Raman response time. This parameter is not available in the Matlab m-file version. (default = 0)
t0 scalar The optical cycle time (= 1/f). This parameter is not available in the Matlab m-file version. (default = 0)
maxiter scalar (int) The maximum number of iterations to make per step. If the solution does not converge to the desired tolerance within this number of iterations, a warning message will be generated. Usually this means that the chosen stepsize was too small. (default = 4)
tol scalar Convergence tolerance: controls to what level the solution must converge when performing the symmetrized split-step iterations in each step. (default = 10–5.)

Output Arguments

vector (N) Output optical field, specified as a length-N vector time sequence. Like u0, the output describes the slowly-varying envelope.


Several internal options of the routine can be controlled by separately invoking sspropc with a single argument:

sspropc -savewisdom
sspropc -forgetwisdom
sspropc -loadwisdom

The first command will save the accumulated FFTW wisdom to a file that can be later used. The second command causes sspropc to forget all of the accumulated wisdom. The last command forces FFTW to load the wisdom file from the current directory. The wisdom file (if it exists) is automatically loaded the first time sspropc is executed. The name of the wisdom file is “fftw-wisdom.dat” for the double-precision version of the program and “fftwf-wisdom.dat” for the single-precision version. This can be changed by recompiling the code. Note that the wisdom files are platform- and machine-specific. You should not expect optimal performance if you use wisdom files that were generated on a different computer.

The following four commands can be used to designate the planner method used by the FFTW routines in subsequent calls to sspropc.

sspropc -estimate
sspropc -measure
sspropc -patient
sspropc -exhaustive

The default method is patient. These settings are reset when the function is cleared or when Matlab is restarted.

These options are only available in the compiled version of the routine.


Units and Dimensions: The dimensions of the input and output quantities are arbitrary, as long as they are self consistent. For examle, if |u0|2 has dimensions of Watts and dz has dimensions of meters, then the nonlinearity parameter gamma should be specified in W-1m-1. Similarly, if dt is given in picoseconds, and dz is given in meters, then the dispersion polynomial coefficients betap(n) should have dimensions of ps(n-1)/m. It is of course possible to solve the normalized dimensionless nonlinear Schrödinger equation by setting some of the input terms to 1 or –1 as appropriate.

Periodicity:  SSPROP uses the FFT (DFT) to calculate the spectrum, which implies that the input and output signals are periodic in time. The periodicity is determined by the time increment and the length of the input vector, T = dt*length(u0). Because of the periodic boundary conditions used by the DFT, care must be taken to ensure that if the optical field at the edges of the window is not negligible it must be continuous in both magnitude and phase.

Self Steepening and Raman Effect: If the two parameters tr and to are not zero, SSPROP will emply a simple approximation of the self-steepening effect and the Raman self-frequency shift. This algorithm is significantly more complicated and slower than the conventional split-step method. These two effects usually do not need to be included in the simulation unless the pulsewidth is below 100 fs. Unless you need to include these effects, it is best to set both of these parameters to zero (the default value), in which case SSPROP will use a simpler and faster method.

Iterations and Tolerance: The last two optional parameters, maxiter and tol are related to the symmetrized split-step iteration algorithm. The algorithm uses a trapazoid integration equation to approximate the effect of the nonlinearity over a distance dz, but this approximation requires knowledge of the field at the subsequent distance-step. This problem is solved by using an iterative approach. maxiter represents the maximum number of iterations performed per step, and tol is a positive dimensionless number that tells the algorithm what level of convergence is required before the iteration stops.