Core design

The BMO package is intended to provide optical simulation capabilites with as much "out-of-the-box" comfort as possible, meaning that users should not have to worry about e.g. providing the correct optical sequence and exact alignment of objects. The following five design principles are core assumptions of the underlying API:

1. Optical interactions are decoupled from the underlying geometry representation

2. Optical elements are closed volumes or must mimic as such (exceptions apply, e.g. coatings)

3. Elements should be easily moveable and have working interactions for most angles of incidence

4. Without additional knowledge, tracing is performed non-sequentially

5. With additional knowledge, tracing is performed sequentially

Intersect-Interact-Repeat-Loop

The first two principles will be elaborated upon in more detail in the Geometry representation section. For the latter two design decisions, the following high-level solver schematic can be used to abstract the steps that are performed when calling solve_system! with an input system and beam:

This scheme is loosely referred to as the Intersect-Interact-Repeat-Loop and consists of the following steps:

1. Calculate the closest Intersection between a ray/beam and the objects within the system

2. Calculate the optical Interaction that occurs at the surface or within the volume of the element

3. Attach the next part of the ray chain or beam tree

4. Use the new information to repeat 1.

Once this procedure has been completed, the alignment of the system or other time-dependent optical properties (e.g. the phase of a Gaussian beamlet) can be updated. A new call of solve_system! solves the beam again from its start. This is described in more detail in the section: Tracing systems.

The next sections will focus on the Intersection and Interaction steps.

Intersections

Calculating intersections between straight lines, i.e. rays, and surfaces is a central challenge for every geometrical optics simulator. This must be done with high numerical precision, since many optical effects are sensitive on the order of the wavelength of the light under consideration with respect to position and direction [26, p. 265 ff]. In order to define this mathematically or algorithmically, many different methods exist [27]. The first question is, how is the geometry of the problem defined. This topic is treated in the Geometry representation section. The second question concerns then the algorithm or equation that allows to calculate the point of intersection between a ray and the surface of the element. This function is called intersect3d and is, at its core, defined for each shape and ray:

BeamletOptics.intersect3d Method

intersect3d(shape::AbstractShape, ::AbstractRay)

Defines the intersection between an AbstractShape and an AbstractRay, must return an Intersection or nothing. The default behavior for concrete shapes and rays is to log a warning (No intersect3d method defined for: ...) and to indicate no intersection, that is nothing, which will inform the tracing algorithm to stop. Refer to the Intersection documentation for more information on the return type value.

source

Regardless of the underlying concrete implementation, each call of intersect3d must return nothing or the following type:

BeamletOptics.Intersection Type

Intersection{T}

Stores data calculated by the intersect3d method.

Fields:

• object: a Nullable reference to the AbstractObject that has been hit. The intersect3d methods of shapes leave it nothing, the solver sets it

• shape: a Nullable reference to the AbstractShape of the object that has been hit. Provide it: composite objects and the check whether the rays of a beamlet hit the same shape identify the hit part by this reference

• t: length of the ray parametrization in [m]

• n: normal vector at the point of intersection

source

Since an optical element can consist of multiple joint shapes, the return type must store which specific part of the object was hit: the shape of the Intersection is that part, and the solver compares and looks up the hit part by this reference, e.g. to decide whether the rays of a beamlet hit the same shape. The object of the Intersection is set by the solver, the intersect3d method of a shape leaves it nothing.

The following diagram shows how trace_all reaches the shapes of the objects and groups of a System, including the Bounding sphere pretest, whose table is described below. A StaticSystem has no groups while tracing, since it flattens them when it is built.

The shape backend is selected by dispatch on the type of the shape, see Geometry representation.

Bounding sphere pretest

Since the tracing is non-sequential, every ray is tested against every shape of the system, and most of these tests are misses. A shape can therefore provide a sphere that encloses it. A ray that does not pass through this sphere is a miss without a call of intersect3d. This is optional: a shape without a sphere is tested as described above.

BeamletOptics.bounding_sphere_of Method

bounding_sphere_of(shape::AbstractShape)
bounding_sphere_of(object::AbstractObject)
bounding_sphere_of(group::AbstractObjectGroup)

Computes a sphere that encloses the shape, object or group in its current pose and returns it as a SingleBoundingSphere or MultiBoundingSphere, or returns NoBoundingSphere() if there is none.

For a shape, the method is an optional part of the AbstractShape interface. With the default NoBoundingSphere() every ray is tested via intersect3d. With a sphere, a solve tests a ray against the sphere first and skips the shape if it misses, which pays off for shapes with a costly intersect3d. The sphere is computed once per solve_system! and not stored in the shape, hence the method may take some time, e.g. look at all vertices of a mesh. The table holds spheres of floating-point shapes only: for other number types, e.g. dual numbers under ForwardDiff, the shape is tested without a sphere.

The sphere must enclose every point of the shape, for all its parameters: a hit outside of the sphere is lost without a warning. It should also be tight, since a ray that hits the sphere is tested as before. A safety margin is not needed, the solver adds its own tolerance. Draw the result with render!(ax, sphere) to look at it.

# a cylinder of the radius `r` and the height `h`, with the origin at the center of its base
function BeamletOptics.bounding_sphere_of(c::MyCylinder)
    return BeamletOptics.SingleBoundingSphere(c, Point3(0, c.h / 2, 0), sqrt(c.r^2 + (c.h / 2)^2))
end

For an object, the result follows from its shapes, see AbstractShapeTrait: the sphere of the shape of a SingleShape object, and for a MultiShape object or an object group the MultiBoundingSphere that encloses the spheres of all its parts, or NoBoundingSphere() if one part has none.

source

All Signed Distance Functions (SDFs) and Meshes have a bounding sphere. An object of several shapes and an object group have the sphere around the spheres of their parts, such that a ray that misses it is tested against none of them. The table holds spheres of floating-point shapes only: for other number types, e.g. dual numbers under ForwardDiff, the pretest is skipped and the shapes are tested directly.

The spheres are not stored in the shapes or objects. solve_system! computes them once per call, for the poses that the objects have at that time, via bounding_sphere_of(x) into a table, which the intersection code reads per ray via bounding_sphere_of(table, x). The solve of a beam group computes the table once for all its beams, and the solves of its beams keep it instead of computing it again. Outside of a solve there is no table, i.e. a direct call of intersect3d(object, ray) tests the object without a sphere.

The table lives in the task-local storage of the solving task, not in a ScopedValue, because the intersection code reads it for every object and ray and this read is several times cheaper. A task that is spawned within a solve does not see the table of its parent. It passes the table on by calling BeamletOptics.with_bounding_spheres(f, table), as the solve of a beam group does. The system must not change while a solve runs, since the spheres are those of the poses at the start. The types in the sphere code are asserted, and the methods for the parts of an object are kept separate from those for a shape, because inference would otherwise widen their argument types and the shapes would be tested by code compiled for abstract types.

BeamletOptics.bounding_sphere_of Method

bounding_sphere_of(table::BoundingSphereTable, x)
bounding_sphere_of(::Nothing, x)

Returns the sphere that the table holds for x, or NoBoundingSphere() if there is none. Nothing is computed.

source

BeamletOptics.SingleBoundingSphere Type

SingleBoundingSphere{T} <: AbstractBoundingSphere

A sphere that encloses one shape, see bounding_sphere_of.

Fields

• pos: the center in world coordinates [m]

• radius: the radius [m]

source

BeamletOptics.MultiBoundingSphere Type

MultiBoundingSphere{T} <: AbstractBoundingSphere

A sphere that encloses several bounding spheres, i.e. all parts of a MultiShape object or all objects of an object group, see bounding_sphere_of. A ray that misses it hits none of the parts, hence a solve tests it before the parts. It is tested against a ray like a SingleBoundingSphere and differs only in what it states: the sphere of several things.

Fields

• pos: the center in world coordinates [m]

• radius: the radius [m]

source

BeamletOptics.NoBoundingSphere Type

NoBoundingSphere <: AbstractBoundingSphere

States that a shape or object has no bounding sphere. It is then tested against every ray via intersect3d, see bounding_sphere_of.

It has no center and no radius: position and radius throw an ArgumentError. Code that needs them handles NoBoundingSphere by a method of its own, or tests for it first.

source

The spheres of a system can be shown via render!(ax, BeamletOptics.bounding_sphere_of(x)) for one sphere or the show_bounding_sphere keyword of render!.

Interactions

Optical interactions are performed after the point of intersection has been determined. The interact3d interface allows users to implement algorithms that calculate or try to mimic optical effects. The fidelity of the algorithm is effectively only limited by the amount of information that can be passed into the interact3d interface. The method is defined as follows:

BeamletOptics.interact3d Method

interact3d(::AbstractSystem, object::AbstractObject, ::AbstractBeam, ::AbstractRay)

Defines the optical interaction between an incoming/outgoing beam/ray of light and an optical element, must return an AbstractInteraction or nothing. The default behavior is that no interaction occurs, i.e. return of nothing, which should stop the system tracing procedure. Refer to the AbstractInteraction typedocs for more information on the return type value.

A MultiShape object receives this call for whichever of its parts was hit. Its method finds that part, e.g. via shape(intersection(ray)), and hands the beam to it, see CubeBeamsplitter. Parts that do not interact as one belong in an ObjectGroup or another AbstractObjectGroup instead, where each part interacts on its own.

source

As with the intersect3d method, a predefined return type must be provided in order to make the solve_system! interface work. The AbstractInteraction is used in order to create "building blocks" from which the output beam is constructed.

BeamletOptics.AbstractInteraction Type

AbstractInteraction

Describes how an AbstractBeam and an AbstractObject interact with each other. This data type stores information from the interact3d function and provides it to the solver. The solver can use this data to extend the AbstractBeam.

Implementation reqs.

Subtypes of AbstractInteraction must implement the following:

Fields

• hint: a nullable Hint for the solver. The field is required, since the solver reads it for every interaction, but it may be nothing

Beam data

It is required that concrete implementations of this type provide some form of data on how to extend the beam. For instance, refer to theBeamInteraction and GaussianBeamletInteraction.

source

An object whose interact3d method stores data in the object, e.g. a detector its hits, also implements initialize!, which discards this data before a new solve.

The interact3d return type holds one new beam segment per interaction, which the solver appends to the current beam. Further beams, e.g. the reflected part at a beamsplitter, are attached to the beam as children by the component, and the solver queues them after the current beam, see Tracing systems. The developer needs to take into account that after e.g. a lens surface air-to-glass interaction, the solver "forgets" that the next logical step is to immediately test against the lens again, since the most likely step will be the refraction at the glass-to-air surface. In order to alleviate this issue, the Hint type can be used.

Re-emitting components

Some components do not continue a beam where it ended, but hand it to another solver and emit the result elsewhere, e.g. a fiber that takes the light at its input facet and releases it at its output facet. Their interact3d method samples the incoming beamlet segment as a field, lets the other solver compute the field at the exit and starts new beamlets there. The exchange format is the PlaneField of OpticsBase.jl (see OpticsBase interoperability):

function BeamletOptics.interact3d(::AbstractSystem, fiber::MyFiber, agb::AstigmaticGaussianBeamlet, id::Int)
    f = PlaneField(agb, id; size, spacing, origin, axes)      # the segment that hit the fiber, on its input plane
    g = other_solver(f)                                       # PlaneField at the exit, with the full phase
    BeamletOptics.relaunch!(agb, [GaussianModeDecomposition(g)])
    return nothing
end

The new beamlets are attached with BeamletOptics.relaunch! instead of children!. Children added by children! continue the parent: its length, optical path and reference plane are counted on. A re-emitted beamlet already carries the full phase in its amplitude and starts at another place, so it must count its optical path from its own start.

BeamletOptics.relaunch! Function

relaunch!(beam::B, new_beams::AbstractVector{B}) where {B <: AbstractBeam}

Attaches new_beams as children of beam that start their own optical reference: they are traced and removed together with beam, but their length, optical path and phase are counted from their own start instead of being continued from beam. The new beams carry the full phase in their amplitude.

Use it in the interact3d of a component that re-emits light elsewhere, e.g. the output end of a fiber or a solver coupled through a field: the light leaves at another place than where beam ended, so the geometric continuation of children! (parent length, reference plane of a beamlet) does not apply.

The new_beams replace any children that beam already has. A beam that the component keeps and attaches again is solved from its first ray, like every beam.

source

PlaneField(beamlet, id; ...) is for components without a detector, e.g. in the interact3d of a component that hands the incoming beamlet to another solver: the segment that just hit the component is id. It builds a single hit of the segment and passes it to PlaneField(hits; ...), which also takes the hits of several segments.

OpticsBase.PlaneField Method

OpticsBase.PlaneField(beam::BeamletOptics.AstigmaticGaussianBeamlet, id::Integer;
    size::NTuple{2,Int}, spacing::NTuple{2,Real}, origin, axes, progress = false)
OpticsBase.PlaneField(beam::BeamletOptics.GaussianBeamlet, id::Integer; kwargs...)

Samples the field of segment id of a traced beamlet beam on a plane, as PlaneField(detector; ...) does for the beamlets that hit a detector.

size and spacing are as for PlaneField(detector; ...); origin (in [m]) and axes (columns u, v, n) place the plane and have no default. The field is that of the segment's straight continuation through the plane (no further interaction), including the optical path from the source (parents included), and follows the conventions of PlaneField(detector; ...): no projection factor, H from the beamlet direction, a scalar GaussianBeamlet polarized along u, R = Inf, n the refractive index of the segment.

Both beamlet models can be sampled on any of their segments, also within a medium.

Warning

The scalar GaussianBeamlet takes no Fresnel losses at surfaces: within a medium its field carries the full power of the beamlet, while that of an AstigmaticGaussianBeamlet is reduced by the transmission of each surface it passed.

source

OpticsBase.PlaneField Method

OpticsBase.PlaneField(hits::AbstractVector{<:Union{GaussianBeamletHit, AstigmaticGaussianBeamletHit}};
    size::NTuple{2,Int}, spacing::NTuple{2,Real}, origin, axes, progress = false)

Samples the coherent field of a vector of beamlet hits on a plane. Build the hits with BeamletOptics.GaussianBeamletHit(beam, id) or BeamletOptics.AstigmaticGaussianBeamletHit(beam, id) (public, not exported).

Arguments

• hits: all GaussianBeamletHits or all AstigmaticGaussianBeamletHits, with a shared vacuum wavelength and refractive index.

• size (keyword, required): number of samples (nx, ny) along u, v.

• spacing (keyword, required): sample spacings (Δu, Δv) in [m].

• origin (keyword, required): center of the plane in global coordinates in [m].

• axes (keyword, required): 3×3 matrix with columns u, v, n.

• progress (keyword): show a progress bar. Default: false.

The field follows the conventions of PlaneField(detector; ...).

Throws an ArgumentError if hits is empty, if the hits are not beamlet hits, or if they do not share wavelength and refractive index.

source

Hints

As mentioned in the previous section, the BeamletOptics.Hint interface allows developers to manipulate the non-sequential solver algorithm into testing against a specific component and shape during the next cycle of the Intersect-Interact-Repeat-Loop. This interface has very high priority during intersection testing.

BeamletOptics.Hint Type

Hint

A Hint can be passed as part of an AbstractInteraction and will inform the tracing algorithm about which AbstractObject in the AbstractSystem will be hit next.

Info

The Hint does not need to result in a guaranteed Intersection. However, if the hinted shape is intersected, it will be immediately assumed as the correct global intersection.

Fields

• object: the object that might or will be intersected next. The solver calls interact3d for this object, hence for a part of a MultiShape object it is the composite itself, not the part

• shape: the underlying shape that will be intersected next, i.e. shape(object), relevant for multi-shape objects

Hint(object) uses shape(object), which is a tuple for a MultiShape object and fails. Pass the hit part explicitly there, Hint(object, part), with the composite as object and the shape that will be hit as part.

source

The main reason for this is the intersection ambiguity encountered at interfaces between air-tight component interfaces, e.g. Plate beamsplitters or cemented Doublet lenses. This is caused by the fact that for a ray with a starting point at this interface, technically both shapes are being "touched" at the same time. Additional program logic considering the ray direction of propagation can not always resolve this ambiguity. Therefore the task of providing additional information to the solver via the Hint interface is placed as a burden on the developer.

The diagram below shows this for the CubeBeamsplitter, whose interact3d method finds the hit part by identity, hands the beam to it and hints the coating; the DoubletLens hints the other lens in the same way.

The solver picks the interact3d method by the object of the hint, not by its shape. For a BeamletOptics.MultiShape component the hint therefore names the composite as object and the part that is hit next as shape, as Hint(cbs, shape(cbs.coating)) does above; Hint(object) alone fails for such an object, since shape(object) is a tuple.

Tip

The Hint interface is needed most for multishape objects with joint surfaces, but any component can use it: a single lens hints itself when a ray enters it or is totally internally reflected, since the next intersection lies on the same lens.

Tracing logic

Tracing systems

In the initial state, it is assumed that the problem consists of objects <: AbstractObject (in a system) and a beam <: AbstractBeam with a defined starting position and direction. No additional information is provided, and the specific path of the beam is not known beforehand. Consequently, brute force tracing of the optical system is required, involving testing against each individual element to determine the trajectory of the beam.

A call of solve_system! traces the beam tree breadth-first: it traces one beam after another and queues the children that the components attach, up to the keyword depth_max of solve_system! branching levels.

BeamletOptics.trace_system! Function

trace_system!(system::AbstractSystem, beam::Beam{T}; r_max = get_default_r_max()) where {T <: Real}

Trace a Beam through an optical system. The number of rays of the beam can be capped by r_max.

Tracing logic

The intersection of the last ray of the beam with any objects contained within the system is tested. If an object is hit, the optical interaction is calculated. If no interaction occurs or no further objects are hit, the tracing procedure is stopped.

Arguments

• system: The optical system through which the Beam is traced.

• beam: The Beam object to be traced.

• r_max: Maximum number of rays of the beam, which ends the tracing.

source

trace_system!(system::AbstractSystem, gauss::GaussianBeamlet{T}; r_max = get_default_r_max(), check_invariant = true) where {T <: Real}

Trace a GaussianBeamlet through an optical system. The number of rays of each component beam can be capped by r_max.

Tracing logic

The chief, waist and divergence beams are traced step-by-step through the system. For each intersection after a tracing_step!, the intersections are compared. If all rays hit the same target, the optical interaction is analyzed, else the tracing stops. Tracing also stops with a warning if the Lagrange invariant of a new segment leaves get_lagrange_tolerance of its start value.

Arguments

• system: The optical system through which the GaussianBeamlet is traced.

• gauss: The GaussianBeamlet object to be traced.

• r_max: Maximum number of rays of each component beam, which ends the tracing.

• check_invariant: enables or disables the check of the Lagrange invariant.

source

trace_system!(system, agb::AstigmaticGaussianBeamlet; r_max = get_default_r_max(), check_invariant = true, threshold = get_invariant_threshold())

Trace an AstigmaticGaussianBeamlet through an optical system. The number of rays of each component beam can be capped by r_max. All 9 component beams (chief + 8 parabasal) are traced in lockstep: the chief ray is traced first for each intersection, followed by the auxiliary rays. All rays must hit the same shape; otherwise tracing stops.

source

For a Beam, trace_system! runs the loop below until the beam has r_max rays, its last ray hits nothing or interact3d returns nothing; the methods for beamlets run one tracing_step! per component ray and one interact3d for the whole beamlet per step, and also stop if the rays hit different shapes or the optical invariant check fails.

This non-sequential mode is comparatively safe in determining the "true" beam path, but will scale suboptimally in time-complexity with the amount of optical elements.

Object order

Unlike with classic, surface-based ray tracers, the order in which objects are listed in the System object vector/tuple is not considered for the purpose of tracing.

CPU and GPU support

Parallelizing the execution of a solve_system! call on the CPU is straight-forward, since the solution is independent for each beam or ray and the solver can run on multiple threads: the beams of a beam group are solved in parallel, see below. Objects that are mutated during a solve, e.g. detectors like the Detector, need special consideration when implementing them, since multiple threads might be able to access the underlying memory, leading to race conditions. Specifically, this means ensuring e.g. atomic write and read access. The Detector takes a lock for its hits and, within a beam group solve, collects them in a buffer per chunk of beams.

With respect to GPU acceleration, this is not the case. Currently, all available implementations of solve_system! are highly branching algorithms which can not be implemented in a parallelized way easily. This will most likely require a specific new subtype of the BeamletOptics.AbstractSystem with determinable sequential properties. This is not a development goal as of the writing of this section.

The beams of a group are solved in chunks, and all beams of a group share one table of bounding spheres. It belongs to a task, so each task of the threads sets it for its beams. Each chunk collects its detector hits and stores them at its end. A free task takes the next chunk, because the beams of a group differ in cost (e.g. the inner rings of a source pass a lens, the outer ones miss it), and those of similar cost lie next to each other. Four chunks per thread balance most of the load. Each further chunk costs a channel access and a flush, which shows if all beams cost the same.