BeamLab 2.0

The all-new BeamLab 2.0 is here! This is the first major version in the history of BeamLab since our initial release 1.0 with some exciting new features such as:

  • Non-uniform simulation grids and new 'periodic' boundary conditions for the Mode Solver Toolbox
  • Extended waveguide rotation and anisotropy simulation options
  • Automatic detection of waveguide symmetries
  • New parameters FieldType and IntensityType for choosing between simulations based on either the 'E' (electric) or 'H' (magnetic) field and also defining the field type as well as any combination of the three field components x,y, and z for intensity and power evaluations
  • Cross-platform installer for Windows, macOS, and Linux

Read on to learn what’s new or directly jump to the release notes.

Non-uniform simulation grids

Simulations based on the finite difference method require the calculation domain to be discretized and enclosed by a boundary. If one wishes to calculate for example the field distribution or propagation constant of a waveguide mode, the refractive index distribution of the waveguide’s cross-section needs to be defined on a finite grid and inserted into the wave equation together with boundary conditions that force any field components to zero at the location of the domain boundary. The accuracy of the simulation is thus determined by (a) how small the grid’s mesh size can be made in the area where the refractive index changes and most of the mode’s energy is concentrated, and (b) how well that area can be isolated from the domain boundaries to keep their influence on the result to a minimum. It is also clear that with a uniform grid spacing, a performance trade-off exists when trying to use as fine a grid as possible and at the same time push the boundaries as far out as possible.

To improve the accuracy of mode calculations without having to significantly sacrifice performance, we introduced in the Mode Solver Toolbox of this release the possibility to define a non-uniform grid which allows one to place a fine grid in the area of interest, i.e., where most energy of the mode is concentrated, and a coarse grid in the area that only serves as a buffer zone to keep the boundaries sufficiently far away.

For better illustration of the various possibilities on how to change the grid structure, the following example shows the calculation results of the fundamental mode of a rectangular waveguide with a size of 4 µm × 2 µm and an index difference of 0.05 between core and cladding based on three different grid types at a wavelength of 1.55 µm.

Uniform grid:

gridPoints = [80 48];
gridSize = [20 12];

Non-uniform grid with two different constant mesh sizes:

gridPoints = [20 12];
gridSize = [5 3];

options.ExpandedGridSize = [20 12];
options.ExpandedGridPoints = [40 24];
options.ExpandedGridOrder = 0;

Non-uniform grid with cubic increase in mesh size:

gridPoints = [20 12];
gridSize = [5 3];

options.ExpandedGridSize = [20 12];
options.ExpandedGridPoints = [40 24];
options.ExpandedGridOrder = 3;

The first figure represents the case where the whole calculation domain (20 µm × 12 µm in size with a total of 80 × 48 = 3840 grid points defined by the parameter gridSize and gridPoints) is filled up with the standard uniform grid. In the remaining two figures the uniform grid is confined to an area of only 5 µm × 3 µm but the size of the calculation domain is still the same although with a largely reduced number of grid points of only 40 × 24 = 960 (defined now by the new parameters ExpandedGridSizeand ExpandedGridPoints) that are distributed across both the uniform and non-uniform grid. The change in mesh size inside the non-uniform grid is defined by a polynomial function whose order can be set by the parameter ExpandedGridOrder. A value of 0 corresponds to a constant mesh size, while a value of, e.g., 3 corresponds to an increase in mesh size according to a cubic function. In the latter case, the transition from uniform to non-uniform grid is much smoother, which further increases the calculation accuracy in the mode’s vicinity, while the mesh size rapidly increases towards to the domain boundaries. However, since the mode carries nearly no energy near the boundaries, the effect of the coarse grid on the calculation accuracy there is negligible.

Periodic boundary conditions

For optical devices that consist of a large number of periodic elements such as photonic crystals (PhCs), it is sensible to take advantage of this periodicity and investigate only the modal properties of one characteristic element, i.e., the unit cell, of the periodic array. To do so, the boundary conditions play an important role as they have to make sure that not only the field from the boundary at the opposite side is properly reproduced but possibly can also have some arbitrary phase shift. Starting with this release, periodic boundary conditions can also be used in the Mode Solver Toolbox. The relevant parameters are BoundaryX and BoundaryY that need to be set to 'periodic', as well as BoundaryXPhase and BoundaryYPhase that can take on arbitrary values as phase shift between opposite boundaries.

The following example shows the calculation of the fundamental mode at the wavelength of 1.55 µm of a PhC that consists of air holes with a diameter of 8 µm arranged on a hexagonal lattice with a lattice constant of 10 µm in a medium with a refractive index of 1.5. The figure on the left shows the unit cell’s refractive index distribution and the one in the center displays the intensity distribution of the fundamental mode. For reference, we added to the right another intensity distribution that was obtained by expanding the area of the calculation domain by a factor of 9 (i.e., 3 times in x- and y-direction, respectively), again with the periodic boundary conditions applied only to the domain boundaries but resulting in exactly the same mode confirming that the much more performant calculation on the single unit cell accurately determines the fundamental mode.

Refractive index distribution (left) and fundamental eigenmode (center and right) of a 2D photonic crystal with a hexagonal lattice

Extended waveguide rotation and anisotropy simulation options

In this release we also improved the way how to use anisotropy in the presence of twist in various waveguides. The amount of (geometrical) twist along waveguides can be set by the two optional waveguide parameters Rotation and RotationEnd which allow one to twist a waveguide by an arbitrary angle in either a clockwise or counterclockwise fashion. Hitherto, if the waveguide exhibited both anisotropy and a twist, the orientation of this anisotropy was synchronized with that twist. Starting with this release we have added the possibility to rotate the orientation of the anisotropy independent of the twist using the new optical parameters AnisotropyRotation and AnisotropyRotationEnd. As their names imply, these parameters allow one to change only the orientation of the anisotropy along the waveguide independent of the waveguide’s (geometrical) orientation or twist. As before, the orientation of the anisotropy still remains synchronized with the twist, if these parameter are not explicitly defined.

The following three videos show several examples of geometrical twist and rotating anisotropy in conjunction with an elliptical waveguide. The elliptical waveguide has a core with a width of 8 µm along the major and 4 µm along the minor axis. The index difference between the core and cladding is 0.5. Furthermore, twist is imposed by assuming that the core is rotated by 90 degrees over a distance of 1 mm, and index anisotropy is imposed by assuming that the waveguide consists of an uniaxial medium exhibiting a 1% higher refractive index along the extraordinary axis which is parallel to the x-axis at the input plane. The take-away conclusion of this simulation is that the perturbation induced by the core’s elliptical shape alone is usually not sufficient to entail a stable rotation of a linear input polarization in conjunction with (geometrical) twist.

Video of the light propagation in elliptical waveguides with and without geometrical twist and index anisotropy
Left: with both twist and index anisotropy following the twist
Center: with twist but without index anisotropy
Right: without twist but with index anisotropy rotating by the same angle as the twist shown in the other two figures

Exciting new demos

In this release, we have also added a total of 12 new demos. One of these shows the focussing mechanism of a micro ball lens when a Gaussian beam with a beam waist of 40 µm propagates through the lens with a diameter of 90 µm.

Three dimensional refractive index distribution of the ball lens (left) and intensity distribution in the x-z plane of a Gaussian beam propagating through the ball lens (right)

Another one shows how orthogonally polarized beams that initially do not interfere with each other suddenly change into interfering beams when passing through a polarizer located at the distance z = 0.35 mm. The polarizer is emulated by a thin element of zero thickness that completely removes the field’s x-component of the incoming beams.

Video of the two orthogonally polarized Gaussian beams propagating through a polarizer at a horizontal crossing angle of 10°

Interested in more examples? Please also check out our examples page.

Release notes

New features:

  • Add cross-platform installers for Windows, macOS, and Linux.
  • Add new function activatelicense which copies BeamLab license files to the ‘licenses’ directory.
  • Add new parameter FieldType which allows one to calculate 'E' (default) or 'H' fields in bpmsolver and modesolver.
  • Add new parameter IntensityType which defines the field type and field component(s) for intensity and power evaluations. If IntensityType is not set, it will assume the value 'Exy' in case FieldType is set to 'E', and 'Hxy' in case FieldType is set to 'H'. The field’s z-component can now also be taken into account by setting the parameter IntensityType to e.g. 'Exyz' or 'Hxyz'.
  • Add possibility to use a nonuniform grid in modesolver. The new parameters ExpandedGridPoints and ExpandedGridSize specify the total grid points and total size of the calculation area including both the inner uniform and outer nonuniform grid. The new parameter ExpandedGridOrder specifies the order (degree) of the polynomial describing the increase in mesh size towards the outer border of the nonuniform grid. The new parameter ExpandedGridDisplay is a switch for displaying the calculation area encompassing only the uniform grid or total grid, i.e., both the uniform and nonuniform grid.
  • Add new 'periodic' boundary condition type to BoundaryX and BoundaryY for modesolver. The amount of phase shift between opposite boundaries in x- and y-direction can be set with the new parameters BoundaryXPhase and BoundaryYPhase, respectively.
  • Add automatic detection of index symmetries for bpmsolver and modesolver.
  • Add new parameter CalculationCenter which allows one to specify the coordinates of the center of the calculation window.
  • Add new unit 'nm' to the parameter LengthUnit.
  • Add new parameters Index3DSmoothingLevel and Index3DSmoothingWidth to indexplot and beamset for specifying level and width, respectively, of the refractive index gradation at abrupt refractive index changes in the original refractive index matrix generated for the 3D refractive index contour.
  • Add possibility to include a complex effective refractive index in the output data structure of dispersionsolver.
  • Add possibility to select modes via parameter ModeSelect in output of dispersionsolver.
  • Add new parameter ModeSelect to dispersionplot for selecting mode(s) when displaying dispersion curves.
  • Add possibility to directly view documentation of a specific parameter with beamlabdoc. For example, the documentation page of PowerTrace can be displayed by executing the command beamlabdoc PowerTrace.
  • Add new function thinpolarizer which emulates an ideal polarizer of zero thickness for converting any input polarization to a linear output polarization at a specified angle.
  • Add possibility to use a nonlinear refractive index in waveguide function plc.
  • Add new parameter Rotation to multicore which allows one to rotate the whole multicore waveguide structure by means of a single parameter entry.
  • Add new parameters AnisotropyRotation and AnisotropyRotationEnd to waveguide functions homogeneous, singlecore, multicore, plc, rib, customwaveguide2D, and customwaveguide3D.
  • Add new parameter LensWidth to thinlens for defining the extension of the lens in x- and y-directions.
  • Add check for the occurrence of phase undersampling in thinlens.
  • Add new parameter Quiver to waveguide functions homogeneous, singlecore, multicore, plc, rib, customwaveguide2D, and customwaveguide3D to allow switching on and off quivers in individual sections.
  • Add new parameter BendPlaneAngle to rib waveguide function.
  • Add new parameter RefractiveIndex to gaussinput which allows one to generate a Gaussian input beam in a medium with a refractive index other than 1.
  • Add new parameter ModeSlicesXY to modeinput which allows one to define the mode location for modeinput independent of a modesolver calculation.
  • Add new parameter Angle to custominput which adds a phase tilt to the input field to have the input beam propagate at a specified angle.
  • Add new parameter CustomDatabase to getmaterial and getindex which allows one to define custom BeamLab Material Databases as text files.
  • Add new parameter LengthUnit to getindex which allows one to specify the unit of wavelength for refractive index calculations.
  • Add new parameters SlopesX and SlopesY to splinepath which allow one to specify the endpoint slopes of the spline in x- and y-direction, respectively.

New demos:

  • Add demo rib_modes_nonuniform_grid which calculates the field distributions and effective indices of the eigenmodes of a rib waveguide using a nonuniform grid.
  • Add demo photonic_crystal_modes which calculates the effective indices and field distributions of the eigenmodes of a photonic crystal with a hexagonal lattice geometry using periodic boundary conditions.
  • Add demo photonic_crystal_fiber_endlessly_singlemode which calculates the modes and dispersion characteristics of an endlessly single-mode photonic crystal fiber.
  • Add demo polarization_converter which shows how orthogonal beams that do not interfere with each other suddenly change to interfering beams by means of a polarizer.
  • Add demo total_internal_reflection which simulates the propagation of two Gaussian beams in a thick slab waveguide demonstrating total internal reflection.
  • Add demo freespace_talbot which shows how to generate fractal patterns of subimages of a Gaussian input beam akin to the Talbot effect.
  • Add demo twisted_waveguide which simulates the behavior of the mode field and its polarization state when an anisotropic waveguide is twisted along the propagation direction.
  • Add demo slot_bend_modes which calculates the field distributions and effective indices of the first three eigenmodes of a bent slot waveguide.
  • Add demo slab_modes_2d which calculates the field distributions and effective indices of the eigenmodes of a two-dimensional step-index slab waveguide with high index contrast.
  • Add demo ball_lens which simulates the propagation of a Gaussian beam through a ball lens.
  • Add demo ronchi_grating which simulates the beam propagation through a Ronchi phase grating.
  • Add demo grin_selfimaging which shows self-imaging in a graded-index multimode fiber.

Enhancements:

  • Improve figure handles structures returned by bpmsolver, bpmplot, modesolver, modeplot, inputplot, indexplot, and all thin element waveguide functions.
  • Improve the way how IndexScannerStep and Index3DStep are calculated if not defined by the user.
  • Improve performance of dispersionsolver.
  • Improve performance and memory usage of gaussinput.
  • Improve Anisotropy parameter in singlecore, multicore, and customwaveguide2d. The orientation of the anisotropy in the x-y plane of these three waveguide functions is now also adjusted according to the amount specified by the Rotation parameter.
  • Improve power parameter assignment in thin element and fast element functions.
  • Improve getmaterial and getindex. When the first parameter is specified as a material structure, it is now possible to override the Delta value.
  • Improve the display of refractive index contours so that the index contour lines lie exactly at the section boundary.
  • Improve display of 3D refractive index contour plots such that the index contour begins and ends at exactly the section begin and end.
  • Improve command-line interface output and formatting.
  • Improve display of warnings in the command-line interface.
  • Improve figure titles.
  • Improve demo scripts.
  • Improve documentation. Add description of output structures of bpmsolver, modesolver, and dispersionsolver to documentation.

Parameter changes:

  • Remove parameter PowerType. The field type and field component(s) for intensity and power evaluations are now defined via the new parameter IntensityType.
  • Rename parameter CoreRotation to Rotation in singlecore.
  • Rename parameter WaveguideRotation to Rotation in customwaveguide2d.
  • Remove parameters Rotation and RotationEnd from waveguide function homogeneous. One can now define the rotation of the anisotropy via the new parameters AnisotropyRotation and AnisotropyRotationEnd.
  • Remove parameter InnerCladding from singlecore. One can now create an arbitrary number of core zones by defining coreWidth and coreIndex as cell arrays. Shift and ShiftEnd can now be defined for each core zone individually.
  • Change unit of parameter Phase in input field functions from degree to radian.

Bug fixes:

  • Fix bug in bpmsolver when SlicesXY contains the value 0 and the first section is a thin element.
  • Fix bug in bpmsolver for figure handles of x-z and y-z slice plots. If, e.g., both SlicesXZ and SlicesXZGraphType contain multiple values, the figure handles of all figures are now returned in a 2D structure array.
  • Fix bug in bpmsolver where only the first figure handle of multiple stacked slices plots corresponding to multiple SlicesXYGraphType values was returned.
  • Fix bug in bpmsolver when using a SlicesYZGraphType of type 'Int3D'.
  • Fix bug in Shift parameter of thinaperture.
  • Fix bug in gaussinput for 2D simulations in the y-z plane.
  • Fix bug in gaussinput when using a nonzero value of Angle in combination with certain parameter configurations.
  • Fix bug in install function when beamfunctions is not on the MATLAB path.
  • Fix bug with respect to index contour line generation.
  • Fix bug in documentation where the Help browser header overlaps in-page anchors.
  • Various minor bug fixes.

BeamLab 1.5

The all-new BeamLab 1.5 is here! This release contains some exciting new features as well as performance and stability improvements. Read on to learn what’s new or directly jump to the release notes.

Optical imaging simulations in BeamLab

Part of the BeamLab software is based on the Beam Propagation Method (BPM) which is typically used to evaluate the evolution of optical fields in free space and waveguides with rather simple input field distributions similar to those of Hermite- or Laguerre-Gaussian beams. However, in combination with complex images (e.g., the face of a wolf as shown in the example below) as input distribution, it can also be used to analyze a wide range of imaging systems incorporating lenses, apertures, optical elements that correct for aberrations or diffract the beam, and many more. BeamLab calculates the exact field distribution at any arbitrary plane between the object and image plane, and thus provides you with some unique and in-depth insights on how the imaging system works and how to optimize it.

In this respect, BeamLab 1.5 contains a couple of new features which make it very easy for you to simulate such imaging systems. For creating an input field distribution from an image file, we have added the function imageinput. It supports scaling, shifting, rotating, flipping (left-right or up-down), and inverting the image via the optional parameters Width, Shift, Rotation, FlipLeftRight, FlipUpDown, and Inverse, respectively. For example, the image wolf.png can be imported and rotated by 180° via the following simple command:

inputField = @(b) imageinput(b,'wolf.png','Rotation',180);

By using images as input field distributions you will be able to simulate imaging systems such as the following 4f optical system which consists of two lenses with identical focal lengths f  located at z = f and z = 3f:

Video of the light propagation in a two-lens system (4f optical system)

More slices options

Next, let’s assume you wish to add an aperture in the Fourier plane (2f-plane) of the 4f optical system to simulate an optical low-pass filter, i.e., a filter that blocks the higher spatial frequencies. In BeamLab 1.5, it is now easier than ever to display the aperture’s transmittance and the x-y intensity distribution (x-y intensity slice) directly behind the aperture by using the new optional parameters TransmittanceSlicesXY and SlicesXYSectionEnd, respectively:

apertureWidth = [50 50];
options.TransmittanceSlicesXY = true;
options.SlicesXYSectionEnd = true;

waveguide{3} = @(b) thinaperture(b,apertureWidth,options);

Transmittance of the thinaperture element (left) and xy intensity slice directly behind the aperture (right)

Large-scale optics

Starting with this release, you have now also more control over the amount and type of parameters that are returned in the output structure at the end of each simulation. Until the previous release, the distributions of intensity, each field component, and refractive index were always returned, resulting in a potential memory bottleneck for very large matrices, i.e., large values of gridPoints. The newly added parameter OutputParameters now allows you to work with larger matrices by reducing the output data to exactly those components that you need most. For example, if you are only interested in the intensity and the x-component of the electrical field of your simulation, you can now simply add the following line of code to your simulation file:

options.OutputParameters = {'Int','Ex'};

In combination with the powerful fasthomogeneous function introduced in BeamLab 1.3, this feature comes in handy if you wish to simulate, besides optical waveguides and micro-optics, also beam propagation through rather large-scale optical elements with lateral dimensions, e.g., in the centimeter range, that may require very large values of gridPoints.

Performance mode

For your convenience, the demo m-files that come bundled with BeamLab, are not optimized for speed but rather for making it as easy as possible for you to learn how to use BeamLab. We have added a new parameter called performanceMode at the begin of selected simulation files, which is a switch for optimizing the simulation performance. Its default value is false, but when set to true, the refractive index scan and field monitor functions will be turned off all together independent of the settings in the parameters IndexScanner, Index3D, and Monitor, resulting in a significantly reduced computation time.

For general tips on how to reduce the overall computation time of your BeamLab simulation, please refer to the simulation performance documentation page by executing the command beamlabdoc simulation_performance. If you want to optimize the performance of your scripts, do not hesitate to contact us — we are happy to help you with the implementation of your optical problem.

Exciting new demos

In this release, we have also added four new imaging demos. Among those is the imaging of a pure phase object through a 4f imaging system using the principle of phase contrast microscopy. By modulating the DC (zero) spatial frequency component in the Fourier plane at z = 2f with a quarter-wave phase plate, the inherently invisible phase distribution in the object plane can be converted to a visible intensity distribution in the image plane. The following two figures show the intensity and phase distributions in the object and image planes at z = 0 mm and z = 8 mm, respectively.

Intensity (left) and phase distribution (right) in the object plane

Intensity (left) and phase distribution (right) in the image plane

Interested in more examples? Please also check out our examples page.

Release notes

New features:

  • Add new function imageinput, which creates an input field distribution from an image file. It supports scaling, shifting, rotating, flipping (left-right or up-down), and inverting via the optional parameters Width, Shift, Rotation, FlipLeftRight, FlipUpDown, and Inverse, respectively.
  • Add new parameter OutputParameters that allows one to define which parameters are saved in the output structure at the end of the simulation.
  • Add new parameter PerformanceMode that allows one to turn off refractive index scan and field monitor functions all together and thus optimize the performance of a beam propagation simulation. This new parameter has been added to all relevant demo scripts.
  • Add possibility to calculate a virtual bend in multicore, plc, customwaveguide2d, and customwaveguide3d waveguides.
  • Add new parameter SlicesXYSectionEnd to all waveguide functions to plot and store the field distribution appearing immediately behind a waveguide section.
  • Add new parameter TransmittanceSlicesXY to thin element functions to plot and store the transmittance of a thin element.
  • Add new parameter ApertureInverse to thinaperture that allows one to quickly generate and use also the inverse of the aperture function.
  • Add new parameter ModeIndex that allows one to define an effective index value in whose vicinity the modes should be calculated.
  • Add new parameter Polarization to modeinput to allow the usage of a Stokes (polarization) vector also in conjunction with (semi-vectorial) modes as input field.
  • Add SymmetryX and SymmetryY as optional parameters to modeinput.
  • Add possibility to set the optional parameters SlicesXYGraphType, SlicesXYScale, SlicesXYRange, and SlicesXYView individually for inputplot.
  • Add information on the effective mode area to the output data structure of modesolver.
  • Add new optional parameters MonitorScale and MonitorRange for adjusting scale and range of monitor plots independently from x-y slice plots.
  • Add the possibility to use different scale and range values for each subplot in x-y, x-z, and y-z slice plots and monitor plots.
  • Add function beamlabdemos, which directly opens the examples pages of the BeamLab documentation in the Help browser.

New demos:

  • Add demo imaging_single_lens, which simulates the imaging of an object through a single lens with magnification M.
  • Add demo imaging_single_lens_aberration, which simulates the imaging of an object through a single lens with aberrations defined via a Zernike phase screen.
  • Add demo imaging_4f_system, which simulates the imaging of an object through a 4f optical system with different filter types.
  • Add demo imaging_phase_contrast, which simulates the imaging of a pure phase object through a 4f optical system using the principle of phase contrast microscopy.

Enhancements:

  • Improve calculation speed for fundamental modes in modesolver.
  • Improve automatic assignment of ModeNumber and ModeSelect parameters.
  • Improve bpmsolver, indexplot, and stacked slices plots for propagation structures containing fasthomogeneous sections.
  • Improve beamoverlap to allow for input fields with different number of field components.
  • Improve parsing of input parameters in beamset, inputplot, and indexplot.
  • Improve 2D graphs of 2D simulations.
  • Improve install function.
  • Improve output structures.
  • Improve titles and axis labels.
  • Improve error and warning messages.
  • Improve command-line output of modesolver calculations.
  • Improve demo scripts.
  • Improve documentation.

Parameter changes:

  • Rename parameter PowerTraceType to PowerTrace.
  • Rename parameter ShiftType to ShiftTransition in singlecore, multicore, plc, and rib waveguide functions.
  • Rename parameter ShiftFunction to ShiftTransitionFunction in singlecore, multicore, plc, and rib waveguide functions.
  • Remove parameter Monitor from thin element functions. For plotting (and storing) the transmittance of a thin element, one can now use the new parameter TransmittanceSlicesXY instead.

Bug fixes:

  • Fix bug when using a function handle for indexFunction2D or indexFunction3D in customwaveguide2d or customwaveguide3d, respectively.
  • Fix bug with respect to Anisotropy parameter.
  • Fix bug with respect to FontSizeAxes and FontName not working correctly for 3D index contour plots.
  • Various minor bug fixes.

BeamLab 1.4

Say hi to BeamLab 1.4 — another major release with some amazing new features and enhancements! This time we enhanced BeamLab to work with any input polarization state and a wide variety of anisotropic materials. Read on to learn what’s new (jump to release notes).

Arbitrary input polarizations

Beginning with this release, it is now possible to define arbitrary polarization states of input fields via a polarization vector that corresponds to the normalized Stokes vector [S1 S2 S3] with S0 = 1. When using multiple input fields, the Polarization vectors can even be defined independently from each other. For example, the following demo shows the behavior of two Gaussian beams when they cross each other in free space — in one case they have the same circular polarization defined by a Polarization vector of [0 0 1] and thus strongly interfere with each other, while in the other case they are orthogonally polarized with Polarization vectors of [0 0 1] and [0 0 -1] and do not interfere at all. With BeamLab’s simple user interface, all it takes to define these input fields is two lines of MATLAB® code:

Case 1 (same circular polarization):

inputField{1} = @(b) gaussinput(b,[30 30],'Shift',[ 30 0 0],'Angle',[-5 0],'Polarization',[0 0 1]);
inputField{2} = @(b) gaussinput(b,[30 30],'Shift',[-30 0 0],'Angle',[ 5 0],'Polarization',[0 0 1]);

Case 2 (orthogonal circular polarization):

inputField{1} = @(b) gaussinput(b,[30 30],'Shift',[ 30 0 0],'Angle',[-5 0],'Polarization',[0 0 1]);
inputField{2} = @(b) gaussinput(b,[30 30],'Shift',[-30 0 0],'Angle',[ 5 0],'Polarization',[0 0 -1]);

The following video shows the beam propagation of the two different sets of input fields:

Video of the propagation of two Gaussian beams with the same (top) and orthogonal (bottom) polarizations

Analysis of bent waveguides

We’ve also added the possibility to calculate a virtual bend in singlecore and rib waveguides. The following demo shows how the fundamental mode propagating in a silica fiber with a core index delta of 2% starts to radiate light off the core if configured as a fiber loop with a bending radius of 0.5 mm. The waveguide structure consists of 3 sections where the fiber loop (section 2) with a total length of 3.14 mm was inserted between two straight fiber sections each having a length of 1 mm.

Schematic of the fiber loop

Video of the intensity profile of the beam propagating in the fiber loop

Intensity distribution in the fiber loop (xz slice)

Propagation through anisotropic media

The possibility of defining arbitrary polarization states at the input makes it also interesting to study the change of a light wave’s polarization state when propagating through birefringent media. In the following demo, a circularly polarized Gaussian beam is first sent through a quarter-wave plate transforming the circular input polarization into a linear output polarization with the polarization plane tilted by 45°. After a short length of free-space propagation, a half-wave plate then rotates the polarization plane further from 45° to 90°, i.e., parallel to the y-axis. While the quarter-wave plate’s two principal axes are oriented in parallel to the x– and y-axis, the principle axes of the half-wave plate have been rotated by 67.5° in counter-clockwise direction to achieve the additional rotation of the polarization plane by 45°. Both wave plates have a basic refractive index of 1.5 and a thickness of 10 µm requiring a relative birefringence of 1% and 2% for the quarter-wave and half-wave plate, respectively. The wavelength used in this simulation is 0.6 µm. The whole propagation structure with five sections is defined by the following simple MATLAB® code:

indexFunction{1} = @(b) homogeneous(b,2,1);
indexFunction{2} = @(b) homogeneous(b,10,1.5,'Anisotropy',[1 1.01 1]);
indexFunction{3} = @(b) homogeneous(b,2,1);
indexFunction{4} = @(b) homogeneous(b,10,1.5,'Anisotropy',[1 1.02 1],'Rotation',67.5);
indexFunction{5} = @(b) homogeneous(b,2,1);

Video of the beam propagating through a quarter-wave and a half-wave plate

Interested in more examples? Please also check out our examples page.

Release notes

New features:

  • Add the possibility to define the input polarization via a polarization vector (i.e., normalized Stokes vector) that allows one to work with arbitrary input polarizations. In this respect, the Polarization parameter has been moved from beamset to the input field functions. Both semi-vectorial and full-vectorial BPM calculations can be performed with arbitrary input polarizations. For semi-vectorial mode calculations the desired polarization can be set via the new parameter ModePolarization.
  • Add information on the E-field’s polarization state at the x-y slices and output plane to the output data structure of a BPM calculation.
  • Add the possibility to calculate a virtual bend in singlecore and rib waveguides.
  • Add possibility to define the anisotropy of waveguides by means of a 3- or 5-element vector representing the relative components [xx,yy,zz] or [xx,yy,zz,xy,yx] of the index tensor.
  • Add new optional parameters AnisotropyCore, AnisotropySubstrate, and AnisotropyCladding to rib, replacing the single parameter Anisotropy. This allows one to define the refractive index anisotropy for core, substrate, and cladding individually.
  • Add two new options 'DvsLambda' and 'DvsV' to the parameter DispersionGraphType for plotting the dispersion parameter D as a function of wavelength and normalized frequency, respectively.
  • Add scale options 'lincustom' and 'logcustom' to parameters SlicesXYScale, SlicesXZScale, SlicesYZScale, and PowerTraceScale. The new scale options allow one to define a linear or logarithmic scale, respectively, normalized to a custom value. In case of SlicesXYScale, SlicesXZScale, and SlicesYZScale, this custom value is defined by the new optional parameter SlicesCustomNorm and in case of PowerTraceScale, the custom value is defined by the newly added parameter PowerCustomNorm.
  • Add possibility to address simultaneously an arbitrary ensemble of cores with ModeCore and thus calculate the (super)modes defined by these cores.
  • Add new optional parameter ModeSelect that allows one to select specific modes (either in form of a scalar or a vector) out of a larger mode ensemble.
  • Add the possibility to define modesolver relevant parameters such as ModeCore, ModeNumber, ModePolarization, ModeSections, ModeSelect, and ModeSortingOrder directly as input parameters of modeinput. In addition, modeinput now also allows one to define an input distribution that consists of a linear combination of modes with different mode amplitudes by assigning a vector containing the relevant mode numbers to ModeSelect and a vector containing the appropriate mode amplitudes to ModeAmplitude.
  • Add new optional parameter Phase to all input field functions for specifying the input field’s phase in degrees.
  • Add possibility to display phase and index distributions in x-z and y-z slices.
  • Add possibility to display multiple x-z and y-z slice plots of different types.
  • Add new optional parameter QuiverLength that allows one to adjust the length of the quiver arrows.
  • Add new optional parameter ShowSectionTitles for displaying section titles inside the figure during monitor. The position inside the figure can be specified by the new parameter SectionTitlePosition.
  • Add option to beamlabdoc for searching the documentation for pages with words matching a specified expression.
  • Add demo fiber_loop, which simulates the propagation of the fundamental mode of a fiber through a fiber loop.
  • Add demo fiber_bend_modes, which calculates the first three eigenmodes of a fiber with a bending radius of 10 cm.
  • Add demo fiber_bend_loss, which calculates the bending loss of the first three eigenmodes of a fiber with a bending radius of 10 cm as a function of the wavelength.
  • Add demo quarterwave_plate, which simulates the propagation of a Gaussian beam with circular polarization through a quarter-wave plate.
  • Add demo noninterfering_beams, which simulates the propagation of two orthogonally polarized Gaussian beams that cross each other in free space.
  • Add demo fiber_taper, which simulates the propagation of the fundamental mode of a fiber whose size is changed by means of a fiber taper.

Enhancements:

  • Eliminate the need for the user-defined parameters FieldSymmetryX and FieldSymmetryY. Field symmetries are now set automatically according to the input field and index symmetries SymmetryX and SymmetryY.
  • Extend full-vectorial bpmsolver to work also with anisotropic homogeneous and singlecore waveguides that are rotated by some arbitrary angle.
  • Improve power area definition. The SmoothingLevel parameter is now also used for smoothing the rim of the power area.
  • Improve display of quiver arrows in x-y slice plots.
  • Improve performance of license manager.
  • Improve demo scripts.
  • Improve documentation.

Bug fixes:

  • Fix bug in dispersionsolver when using cores with index trenches.
  • Fix bug in indexplot when using fasthomogeneous sections.
  • Fix bug when modeling a tapered rib waveguide.
  • Fix license manager bug when using macOS or Linux.
  • Various minor bug fixes.

BeamLab 1.3

Today, we are pleased to announce the release of our next major version BeamLab 1.3 including numerous new and (hopefully also for your application) most attractive features. Throughout the development of this release, we placed emphasis mainly on one thing: performance.

Faster than ever

BeamLab is now faster than ever since we’ve significantly improved the performance of our bpmsolver for full-vectorial calculations and for calculations using a WideAngleOrder greater than 0, a new parameter which has also been introduced in this release. The following diagram shows a comparison of the execution times of a full-vectorial BPM simulation with our demo gradedindex_waveguide for three different use cases:

  1. No waveguide symmetry, i.e., SymmetryX = false and SymmetryY = false
  2. SymmetryX = true
  3. SymmetryX = true and SymmetryY = true

Thanks to the optimized algorithms, BeamLab 1.3 outperforms BeamLab 1.2 by up to a factor of 4 with respect.

Benchmark results using an Intel® Core™ i7-7500U CPU @ 2.70GHz, 16 GB RAM, and MATLAB® version 9.3 (R2017b)

Next to these changes to bpmsolver, we’ve also significantly improved the performance of the function gaussinput and the performance of our license manager.

Not yet fast enough?

In BeamLab 1.3, we’ve also added a new waveguide function called fasthomogeneous, which opens up completely new possibilities in terms of decreasing the execution time of your simulations. Unlike the conventional function homogeneous, the new function fasthomogeneous allows you to skip entire homogeneous propagation sections by calculating in almost a twinkle of an eye the resulting field at any desired z-position thanks to an innovative combination of the beam propagation method and the angular spectrum method. For illustrative purposes, the following two figures show how fasthomogeneous can be used to replace a part of a homogeneous propagation section in order to significantly decrease the overall execution time. Since fasthomogeneous does not calculate any intermediate steps between z = 0.3 mm and 0.7 mm, this part is depicted as a white area in xz and yz slices.

Conventional free-space beam propagation using the function homogeneous

Free-space beam propagation with a fasthomogeneous section

Like any other waveguide function, fasthomogeneous can also be used together with a complex refractive index for modeling optical attenuation or gain. Furthermore, by setting the optional parameter PhaseConjugation to false or true, it is possible to quickly obtain the outcome of a forward or backward (i.e., phase-conjugate) propagation. The following figure shows the case where in the first half of the simulation a Gaussian beam is amplified while being propagated forward, whereas in the second half the amplified beam is propagated backwards to its original input state. Note that backward propagation of an amplified beam involves attenuation which is automatically taken care of.

Forward and backward propagation as well as amplification and attenuation capabilities of fasthomogeneous

Exciting new demos

We’ve added two new demos, one of which shows the beam propagation in one core of a photonic lantern. The device consists of 6 cores surrounded by an inner and outer cladding tapered down to a multimode fiber where the inner and outer cladding of the photonic lantern take over the role of the core and cladding of the multimode fiber.

3D refractive index contour (left) and intensity profiles in the photonic lantern

Video of refractive index (left) and intensity of the beam propagating in the photonic lantern (right)

Interested in more examples? Please also check out our examples page.

Release Notes

New features:

  • Add waveguide function fasthomogeneous. This function calculates propagation through a homogeneous section in very short time. Furthermore, by using the optional parameter PhaseConjugation, one can also quickly obtain the outcome of a backward (phase-conjugate) propagation.
  • Add new possibilities for defining input fields. Input fields of bpmsolver are now defined via function handles or cell arrays of function handles, which allows one to efficiently define now also multiple and superimposed input fields. For backward compatibility, it is still possible to use the old definition via field structures generated when directly calling an input field function.
  • Add function inputplot, which displays the input field defined in a beamProblem. The input field to be displayed can also be specified as an input field function handle or cell array of input field function handles and passed to inputplot as second, optional parameter. Consequently, the Monitor parameters have been removed from all input field functions.
  • Add possibility to define complex refractive indices with negative or positive imaginary parts to all waveguide functions for modeling optical attenuation or gain. Until now, only negative imaginary parts were allowed.
  • Add parameter WideAngleOrder. This parameter allows one to specify the order of a Padé approximant up to order 3. A value of 0 corresponds to the Padé order (0,1) and represents paraxial approximation. Furthermore, the values 1, 2, and 3 correspond to the Padé orders (1,1), (2,2), and (3,3), respectively. The default value of WideAngleOrder is 0 except for homogeneous sections where BPM automatically adapts the algorithm for a WideAngleOrder of 1.
  • Add parameter ProfileExponent to the waveguide functions plc and rib. This parameter allows one to specify the exponent of the core’s permittivity distribution to simulate various types of graded-index profiles.
  • Add possibility to specify VideoResolution in dots per inch (DPI). VideoResolution can now alternatively be defined as a string containing '-r' and an integer value indicating the resolution in DPI. For example, '-r300' sets the output resolution to 300 DPI.
  • Add demo photonic_lantern, which simulates the beam propagation in a photonic lantern with vanishing cores.
  • Add demo three_beam_interference. This demo shows how to define multiple input fields (three Gaussian beams with different Shift and Angle values) that propagate simultaneously in free space.

Enhancements:

  • Improve performance of bpmsolver for full-vectorial calculations.
  • Improve performance of gaussinput.
  • Improve performance of license manager.
  • Improve the definition of graded-index core profiles in conjunction with core shapes other than circular.
  • Improve multicore waveguide code with respect to 2D BPM simulations.
  • Improve figure titles.
  • Improve demo scripts.
  • Improve documentation.

Bug fixes:

  • Fix license manager bug when using certain MATLAB® versions.
  • Fix bug when displaying the phase distribution during a BPM calculation.
  • Various minor bug fixes.

BeamLab 1.2.1

We’ve just released BeamLab 1.2.1! In this minor release, we’ve replaced the optional parameter VideoFileFormat with the new, more flexible parameter VideoProfile, which enables you to use all video profiles supported by the MATLAB® function VideoWriter:

  • 'Archival': Motion JPEG 2000 with lossless compression
  • 'Motion JPEG AVI' (default): AVI with Motion JPEG encoding
  • 'Motion JPEG 2000': Motion JPEG 2000
  • 'MPEG-4': MPEG-4 with H.264 encoding
  • 'Uncompressed AVI': Uncompressed AVI with RGB24 video
  • 'Indexed AVI': Uncompressed AVI with indexed video
  • 'Grayscale AVI': Uncompressed AVI with grayscale video

Furthermore, this new release contains also a couple of small bug fixes such as wrong titles in figures generated by modesolver and erroneous phase shifts of the field at the interface between two different media.

Release Notes

New features:

  • Add support for all video profiles supported by the MATLAB® function VideoWriter. For this purpose, the parameter VideoFileFormat has been replaced by the new parameter VideoProfile, which can be set to the values 'Archival', 'Motion JPEG AVI' (default), 'Motion JPEG 2000', 'MPEG-4', 'Uncompressed AVI', 'Indexed AVI', and 'Grayscale AVI'.

Enhancements:

  • Improve titles and axis labels.
  • Improve demo scripts.
  • Improve documentation.

Bug fixes:

  • Fix a bug that introduced an erroneous phase shift of the field at the boundary between two different media.
  • Fix a bug that was causing wrong titles in figures generated by modesolver.
  • Various minor bug fixes.

BeamLab 1.2

Say hi to BeamLab 1.2 — our biggest release yet! Read on to learn more about the most exciting new features and performance improvements (jump to release notes).

Enhanced waveguide design flexibility

In BeamLab 1.2, we’ve added a new class of functions called path functions. These functions allow the definition of complex 3D trajectories with just a few lines of MATLAB® code, making the definition of complex waveguide structures both easy and flexible. Instead of having to split up a complex waveguide structure into many different simpler waveguide sections, the entire structure can now be directly defined in a single section with a much shorter and simpler code as exemplified below for the case of our directional coupler demo:

options.ShiftType = 'custom';

points = [20 0 0;
    20 0 100;
    -2 0 1600;
    -2 0 5600;
    20 0 7100];
options.ShiftFunction{1} = @(z) linsinpath(z,points);

points = [-20 0 0];
options.ShiftFunction{2} = @(z) linsinpath(z,points);

Here, the function linsinpath is one of our new path functions which creates a 3D trajectory consisting of sinusoidal transitions, linear transitions, or combinations thereof, connecting an arbitrary number of points defined as a N x 3 matrix with each row corresponding to the x-, y-, and z-coordinates of a transition point in space.

One of our new demos showcasing this new feature is an integrated RGB beam combiner [1] where all three waveguide arms are defined by these path functions:

Intensity profiles in the RGB beam combiner (left: red input channel; center: green input channel; right: blue input channel)

The following video shows the beam propagation through the RGB beam combiner for each of the three color channels:

Video of red, green, and blue beams propagating in the RGB beam combiner

Simplified coding

Applications that make use of new BeamLab 1.2 functions such as path functions and the input field function modeinput can reduce the code length by more than 70%. For example, the waveguide definition of our directional coupler demo required 105 lines of code in BeamLab 1.1, but can be reduced to less than 30 lines of code in BeamLab 1.2, as exemplified in the code snippet shown above.

Exciting new demos

We’ve added a new Mach-Zehnder interferometer (MZI) demo with two input and two output arms showcasing xz slice plots of the real part of the propagating complex electric field (in contrast to the absolute value that is typically shown in BPM based plots). This is just one example of the many new plotting possibilities enabled by the newly added parameter SlicesXZGraphType in BeamLab 1.2:

Field profiles in the MZI (left: complete propagation structure; right: detail plot)

We also added a new demo that demonstrates deflection of a Gaussian beam in free space by a series of ideal prisms using the new waveguide function thinprism:

Video of beam deflection via multiple prisms

Interested in more examples? Please also check out our examples page.

Improved parameter orders

In BeamLab 1.2, we’ve also changed the parameter order of the waveguide functions multicore, plc, and rib to improve consistency throughout all waveguide functions. Along with this change, we’ve also replaced the required parameter shift by an optional parameter Shift in the plc waveguide function. The new parameter order allows you now to change your propagation structure more easily from one type of waveguide to another, e.g., from plc to rib.

BeamLab 1.1 and earlier:

multicore(beamProblem,len,coreWidth,coreNumber,ringRadius,coreIndex,claddingIndex);
plc(beamProblem,len,coreWidth,coreNumber,shift,coreIndex,claddingIndex);
rib(beamProblem,len,coreWidth,substrateWidth,substrateCenter,coreIndex,substrateIndex,claddingIndex);

Starting in BeamLab 1.2:

multicore(beamProblem,len,coreNumber,coreWidth,coreIndex,claddingIndex,ringRadius);
plc(beamProblem,len,coreNumber,coreWidth,coreIndex,claddingIndex);
rib(beamProblem,len,coreNumber,coreWidth,coreIndex,claddingIndex,substrateWidth,substrateCenter,substrateIndex);

References

[1] A. Nakao et al., “Integrated waveguide-type red-green-blue beam combiners for compact projection-type displays,” Optics Communications, vol. 330, pp. 45-48, 2014.

Release Notes

New features:

  • Add new class of path functions. These functions allow the definition of complex 3D trajectories with just a few lines of MATLAB® code, making the definition of complex waveguide structures both easy and flexible. They can be used anywhere in BeamLab, where custom functions define a lateral offset in x- and y-direction as a function of the longitudinal z-value. For example, path functions can be used as ShiftFunction in the waveguide functions singlecore, multicore, plc, and rib to define complex trajectories along which waveguide cores are shifted. Path functions can also be used as custom PowerCenterFunction and PowerAreaTransitionFunction in a large number of waveguide functions, offering a high flexibility in terms of defining power calculation trajectories and areas, respectively.
  • Add path function linsinpath, which interpolates the lateral x- and y-values between transition points via sine functions, linear functions, or combinations thereof. The optional parameter TransitionType specifies the absolute length of the sinusoidal transition after a transition point and before the next one. A value of 0 corresponds to a purely linear shift between two transition points and a value of Inf (or greater than the transition length between two transition points) corresponds to a purely sinusoidal transition, respectively.
  • Add path function splinepath, which uses the MATLAB® function spline for interpolating the lateral x- and y-values between transition points via a cubic spline interpolation.
  • Add path function pchippath, which uses the MATLAB® function pchip for interpolating the lateral x- and y-values between transition points via a piecewise cubic Hermite interpolating polynomial (PCHIP).
  • Add function modeinput. This function generates an input field for a BPM simulation based on the eigenmode number ModeNumber of the waveguide core specified in ModeCore of waveguide section specified in ModeSections instead of having to use a separate modesolver call in conjunction with custominput. The function modeinput uses the cached output from a previous call if the waveguide function and other relevant parameters have not changed. Since only the waveguide from the section specified in ModeSections is relevant for modeinput, other sections can be modified and/or additional sections can be added or sections can be removed without having to recalculate the input field.
  • Add function getcommonvars. This new function automatically returns all common variables, which have been defined in a waveguide function up to the point where getcommonvars is executed. This eliminates the need to manually specify the common variables which should not be cleared in a subsequent execution of sectionclear.
  • Add waveguide function thinprism. This function emulates an ideal prism of zero thickness that deflects an input beam by a specified angle.
  • Add waveguide function thinzernike. This function creates a circular phase screen defined by a Zernike polynomial.
  • Add demo rgb_beam_combiner, which simulates the beam propagation in an integrated silica RGB beam combiner using the new path function linsinpath for defining the waveguide structure.
  • Add demo mach_zehnder_interferometer, which simulates the beam propagation through a Mach-Zehnder interferometer (MZI) with two input and two output arms.
  • Add demo rib_splitter, which simulates the beam propagation in a polymer rib waveguide splitter using the new path function pchippath for defining the waveguide structure.
  • Add demo freespace_deflection, which simulates the deflection of a Gaussian beam by several ideal prisms using the new waveguide function thinprism.
  • Add demo directional_coupler_performance, which performs the same simulation much faster as the original demo directional_coupler by omitting the index scanner and live monitor functionalities.
  • Add possibility to define SuperGaussDegree in gaussinput for x- and y-direction independently.
  • Add parameters SlicesXZGraphType and SlicesYZGraphType. Similar to the parameter SlicesXYGraphType, these two new parameters specify the plot types of x-z and y-z slice plots. These parameters now allow also the generation of 2D or 3D slice plots of the real part, imaginary part, amplitude, or phase of the field distribution, besides the intensity distribution.
  • Add parameter ModeCore. This parameter specifies a single core in a multi-core type waveguide function, such as multicore, plc, or rib, to be used for a mode evaluation. In conjunction with the input function modeinput, it can be used to select a specific core as the location of field input for BPM simulations, and thus eliminates the need for defining a separate dummy waveguide section that was hitherto used only for an input field specific mode evaluation.
  • Add parameter NumericalAperture to gaussinput to take into account the numerical aperture of the optical system generating the Gaussian input beam. When set to Inf (which corresponds to the default value), also beam waists beyond the classical resolution limit can be generated, e.g., beam widths smaller than the wavelength. Note however that the use of such small beam waists cause the beam to strongly diverge shortly behind the location of the waist.
  • Add new materials to material database.

Enhancements:

  • Change parameter order in waveguide functions multicore, plc, and rib to improve consistency throughout all waveguide functions.
  • Replace required parameter shift by optional parameter Shift in plc waveguide function.
  • Improve dispersionsolver to work also with complex refractive indices.
  • Improve input field functions.
  • Improve waveguide functions with respect to 2D calculations.
  • Improve memory efficiency in modesolver.
  • Remove numerical noise in imaginary part of modes.
  • Remove restriction of the video export functionality requiring MATLAB version 8.2 (R2013b) or later.
  • Use IndexRange also for all monitor and x-y slice plots.
  • Add possibility to define MonitorGraphType, SlicesXYGraphType, SlicesXZGraphType, SlicesYZGraphType, IndexSlicesXYGraphType, and DispersionGraphType as string.
  • Add possibility to set SlicesXYGraphType, SlicesXZGraphType, SlicesYZGraphType in bpmplot independent of previous choices.
  • Add possibility to set SlicesXYGraphType in modeplot independent of previous choices.
  • Add possibility to define CoreShapeFactor, Shift, and ShiftEnd as cell arrays in rib.
  • Add possibility to change the faceted shading line width via optional parameter LineWidth.
  • Add general tips how to reduce the overall computation time of BeamLab simulations to documentation.
  • Add subdirectories to demos directory.
  • Improve titles and axis labels.
  • Improve warning and error messages.
  • Improve demo scripts.
  • Improve documentation.

Bug fixes:

  • Fix bug with respect to semi-vectorial calculations using y-polarization.
  • Fix units of field and intensity distributions for 2D calculations.
  • Fix bug in plc when using materials from the material database.
  • Fix bug in rib waveguide function with respect to 2D calculations.
  • Fix bug in stacked slices plots when plotting the field’s phase.
  • Various minor bug fixes.