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.

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.