Background
The scalar version of the SSPROP solves the scalar nonlinear Schrödinger equation:
using the symmetrized splitstep 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 selfsteepening effect.
Syntax
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 mfile can be invoked using a similar syntax by replacing sspropc with ssprop. Note however, that the mfile version does not include the delayed Raman response time or selfsteepening terms, and the mfile 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
u0

vector (N)  Input optical field, specified as a lengthN vector time sequence. u0 represents the complex, slowlyvarying 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 stepsize 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 wavelengthdependent 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 realvalued 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 slowlyvaring 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 nonzero term, betap(3), would then represent β_{2}, the groupvelocity 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 n_{2}. 
tr  scalar  The Raman response time. This parameter is not available in the Matlab mfile version. (default = 0) 
t0  scalar  The optical cycle time (= 1/f). This parameter is not available in the Matlab mfile 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 splitstep iterations in each step. (default = 10^{–5}.) 
Output Arguments
u1

vector (N)  Output optical field, specified as a lengthN vector time sequence. Like u0, the output describes the slowlyvarying envelope. 
Options
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 “fftwwisdom.dat” for the doubleprecision version of the program and “fftwfwisdom.dat” for the singleprecision version. This can be changed by recompiling the code. Note that the wisdom files are platform and machinespecific. 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.
Notes
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^{1}m^{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^{(n1)}/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 selfsteepening effect and the Raman selffrequency shift. This algorithm is significantly more complicated and slower than the conventional splitstep 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 splitstep 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 distancestep. 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.