Home

Awesome

Readme and User Manual for BPM-Matlab

What is BPM-Matlab?

BPM-Matlab is a MATLAB-based numerical simulation tool for solving the paraxial Helmholtz equation using the Douglas-Gunn Alternating Direction Implicit (DG-ADI) method to efficiently model the electric field propagation using a Finite-Difference Beam Propagation Method (FD-BPM) in a wide variety of optical fiber geometries with arbitrary refractive index profiles.

Included is also a solver for the Fast Fourier Transform Beam Propagation Method, FFT_BPM.

You can find the latest version of BPM-Matlab at https://github.com/ankrh/BPM-Matlab.

License

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/.

How do I get set up?

Requirements

BPM-Matlab function files

All the functions needed for running BPM-Matlab are located in subfolders of the folder that the examples are in. If you set the folder containing the examples to your MATLAB working direction or you add it to the MATLAB path, all the functions will be automatically located. If you choose to keep your model files (see below) in a different folder to the example files, you will need to manually add the folder with the example files to your MATLAB path.

How do I use BPM-Matlab?

Compilation

The folders include all the executables necessary, so you don't need to compile anything. If, however, you want to change the routine in the FDBPMpropagator.c source code located in the folder "src"), you will need to recompile the mex-file. Check out the header of the source-file on how to do so. The source code is written in such a way that it can be compiled as either C or C++ code using either GCC, MSVC, clang or as CUDA source code using NVCC.

Model files

In BPM-Matlab, you set up your model in a single m-file. You can find a few examples to get you started. Once you're familiar with those, you can start to write your own model files based on those example files.

Segments

If the model contains multiple segments (for example, a stretch of straight fiber followed by a stretch of bent fiber), BPM-Matlab will add a new segment with each call to the FD_BPM() or FFT_BPM() solver functions, using as the input E-field the simulation output of the previous segment (see below). As such, you can stack segments one after the other by only re-defining the parameters that change. Check Example2.m for an example.

List and explanation of parameters

The below parameters apply to the FD_BPM solver. The FFT_BPM solver uses a subset of the below parameters.

Visualization parameters

These parameters affect how the visualizations will look. They will not affect the actual calculation of the E-field result.

If the user is on one of these MATLAB releases, BPM-Matlab will automatically downsample the plots to <= 500x500 (the simulation grid is unaffected, only the plots have lowered resolution). Ideally, the user should simply update MATLAB to a version in which Mathworks have fixed the bug (R2020a or newer).

Solver parameters

These parameters will affect how the solver will numerically describe the problem. You have to ensure yourself that the pixel size and z step size are small enough for the simulation to converge. This is typically done by doing a manual parameter scan of the resolution parameters.

Geometry parameters

These parameters describe the spatial layout of the structure you want to model.

Optical and material properties

Refractive index

The refractive index may be defined as either a 2D or 3D refractive index.

In the case of a 2D refractive index, it is assumed that it is unchanged throughout the segment, with the exception of any tapering and twisting you may have defined. For a 3D refractive index, the third dimension of the array corresponds to the different z coordinates along the segment.

Also be aware that after a segment has been run, the final refractive index profile is stored in the model (P.n), ready to be used automatically as the input to the next segment.

You are free to redefine the simulation grid from one segment to the next (for example, if one segment requires finer resolution than the others). In that case, the refractive index will automatically be interpolated from the old grid to the new grid using the interpn function without the need for user interaction.

There are two ways to define the refractive index:

Method 1, using initializeRIfromFunction()

This method is probably the most common and consists of defining a MATLAB function at the end of the model file that provides the 2D or 3D (possibly complex) refractive index profile as an output.

You run the function initializeRIfromFunction() as shown in example 1 (for 2D) or example 14 (for 3D) to make BPM-Matlab apply your function to the simulation grid in preparation for running the solver.

If you want to work with a 2D refractive index profile, the function you define must take 4 arguments as inputs: X, Y, n_background and nParameters. nParameters is a cell array that is often empty but is available for the user to use in whatever way they want. For example, a user might choose to pass in a core width, core refractive index etc. as different elements in the cell array nParameters.

If you want to work with a 3D refractive index, the function must take 5 arguments as inputs: X, Y, Z, n_background and nParameters.

The syntax of the initializeRIfromFunction() function is: P = initializeRIfromFunction(P,hFunc,nParameters,Nz_n); Here, hFunc is a MATLAB function handle, typically generated with the @ operator (see example 1). The parameter Nz_n is only required when running initializeRIfromFunction() to initialize a 3D refractive index. It describes how many z slices the function should be evaluated at (see example 14). Nz_n should be high enough to resolve all the structure in the z direction, but the higher it is the more memory the solver will need and the slower it will run.

Method 2, defining P.n manually

This method can be used if your refractive index is non-trivial and you want to load it from a data file, for example. In this method you define the refractive index object P.n manually. See example 12. P.n is a "BPMmatlab.refractiveIndexProfile" object with five properties:

The n property is a 2D or 3D array of values corresponding to the (possibly complex) refractive index. For 3D arrays, the first slice along the third dimension corresponds to the refractive index at z = 0, and the last slice corresponds to that at z = P.Lz.

Lx and Ly are the side lengths that correpond to the first and second dimensions of the n array.

xSymmetry and ySymmetry have the same meaning as the ones described above for the solver itself, but in this object they describe which symmetry assumptions n should be interpreted with. For example, if xSymmetry = ySymmetry = 'Symmetry', the P.n array will be interpreted as being just the part of the refractive index that is in the first quadrant.

Initial electric field

Defining the initial electric field is very similar to defining the refractive index, see above.

As with the refractive index, the solver will store the final electric field of a segment in P.E, ready to be automatically used as the input to the next segment.

You are free to redefine the simulation grid from one segment to the next (for example, if one segment requires finer resolution than the others). In that case, the E-field will automatically be interpolated from the old grid to the new grid using the interpn function without the need for user interaction.

There are three ways to define the initial E-field:

Method 1, using initializeEfromFunction()

You define a MATLAB function at the end of the model file that provides the 2D (possibly complex) electric field profile as an output. You run the function initializeEfromFunction() as shown in example 1 to make BPM-Matlab apply your function to the simulation grid in preparation for running the solver.

The function must always take 3 arguments: X, Y and Eparameters. It must return a 2D array of electric field values corresponding to the X and Y coordinate locations.

Similar to nParameters, Eparameters is a cell array that is at the user's disposal to optionally pass additional arguments to the function.

Method 2, defining P.E manually

As with the refractive index, if your initial electric field is non-trivial and doesn't have a simple analytical expression that you can put into a function, you can define P.E manually by setting the properties yourself.

P.E is a "BPMmatlab.electricFieldProfile" object with seven properties:

The field property is the complex E-field 2D array in the initial z slice of the segment. Lx and Ly are the side lengths that correpond to the first and second dimensions of the field array. xSymmetry and ySymmetry have the same meaning as the ones described above for the solver itself, but here describes which symmetry assumptions field should be interpreted with. For example, if xSymmetry = ySymmetry = 'Symmetry', the field array will be interpreted as being just the part of the electric field that is in the first quadrant.

label and neff are used only for electric fields found by the findModes() function, see below. They are not necessary for the user to specify. label is a char array that contains a text description of the field, such as 'LP21' for a mode. neff is the effective refractive index of the mode. It isn't used by the solver but provided for the user's convenience.

Method 3, using findModes()

This method is perhaps the most common. BPM-Matlab includes a mode solver that can calculate the supported modes of the waveguide. The mode finder starts its search at effective refractive indices equal to n_0, so if the mode finder fails to identify the lower-order modes, try increasing n_0. The modes are not used for the propagation simulation itself. The mode finder will not return unguided modes. The modes will be labeled with LPlm designations if the refractive index is assessed to be radially symmetric, otherwise the modes will simply be labeled sequentially (Mode 1, 2, 3 etc.). To use the mode solver, find the modes supported by the waveguide (after you set all the other parameters of P, especially P.n) by using P = findModes(P,nModes); where nModes is the number of modes to try to find.

After nModes, three different Name,Value pairs can be added to the list of arguments. They can be any of the following:

findModes will upon completion set P.modes, an array of "BPMmatlab.electricFieldProfile" objects with each element correspondng to one mode.

After calling findModes(), you can set P.E equal to one of the elements of P.modes. For example P.E = P.modes(1) will set the initial electric field to the first found field, which is often the fundamental mode, depending on the choice of n_0.

A function is provided to make it easy to inject a field that is a superposition of modes. After running findModes(), you can run P.E = modeSuperposition(P,modeIdxs,coefficients) with the arguments

Invoking the solver

After you have defined your model, invoke P = FD_BPM(P); to call the actual simulation tool. During the simulation, the display will be updated according to P.updates.

Optional FFT-BPM solver

In addition to the FD_BPM solver, BPM-Matlab comes with an FFT-BPM solver. FD-BPM should be used for propagation through media with a non-uniform refractive index such as guiding structures, while FFT-BPM can be used for propagation through media with uniform refractive index. Check Example3.m for a comparison.

Unlike FD_BPM, FFT_BPM can step arbitrarily large steps in z without accumulating numerical artifacts, assuming that the beam does not diverge to the edge of the simulation window. If that is the case, however, you will need to put the step size low enough that the absorber layer around the window gets enough steps to get rid of the escpaing light without itself introducing numerical artifacts in the form of "reflections" going back into the simulation window.

Note that many, but not all the parameters listed above for FD_BPM are supported in FFT_BPM.

Multiple segments

Each call to P = FD_BPM(P); or P = FFT_BPM(P); will add a new segment to the simulation that uses as an input E-field and refractive index profiles the outputs of the previous segment, such that you can stack segments one after the other by only re-defining the parameters that changed. Check Example2.m for an example.

Contribution guidelines

If you want to report a bug or are missing a feature, shoot us an email, see below.

Who do I talk to?

This repository is maintained by the Technical University of Denmark "Biophotonics Imaging Group" as well as the "Diode Lasers and LED Systems" group. The main responsibles are Anders Kragh Hansen: ankrh@fotonik.dtu.dk and Madhu Veettikazhy: madve@dtu.dk