3. Ray Tracing in Detail

3.1 Introduction

3.2 Sequential Tracing and Light Path

3.2.1 Light Source Ordering Example
3.2.2 Component Ording Example

3.3 Nonsequential Tracing and Resonate

3.3.1 Three-mirror Resonator
3.3.2 A Subtle Substitute for Resonate

3.4 Rays and Surfaces

3.4.1 How Rayica Generates Rays

3.4.2 Component Surfaces

Faces and Edges
Rolling Your Own Edge
Thick Components and Internal Reflections

3.5 Special Considerations with PropagateSystem

3.5.1 QuickSurfaceSort
3.5.2 Confine
3.5.3 Unconfine
3.5.4 Modular Systems

3.6 TurboTrace

3.6.1 Basics of TurboTrace

How to Enter Symbolic Parameters in Rayica

3.6.2 Lens Optimization
3.6.3 ConstructMeritFunction

Merit based on summation of optical path differences
Calculation of the beam centroid coordinates
Intensity-weighted beam-coordinate calculation
Rayica's built-in temporary variables

3.6.4 Wavefront Analysis

Two-Dimensional Analysis with Interpolation
Three-Dimensional Analysis with Zernike Polynomials

3.1 Introduction

Certain details of Rayica's ray tracing methods commend themselves to study.  By knowing which methods to employ when, you can optimize both your time and Rayica's results.  By knowing how Rayica performs ray-tracing, you can avoid various pitfalls associated with specific types of optical systems.

A large number of optical systems do not demand such fine-tuning.  Rayica handles them well using default methods.  Other optical systems, however, require more careful thought and some knowledge of Rayica's internal workings.

Most optical systems are sequential in nature.  An entrance aperture accepts incoming light, which then propagates from one element to the next, finally emerging from an exit aperture or terminating at an image plane.  A distinguishing feature of such systems is that a given ray intercepts a given surface only once.  All that need concern you about these systems is the order in which components are listed.  This order should match the intended light path of the system.

Other systems do not behave in this straightforward manner.  A ray can intercept the same surface more than once.  For example, a laser cavity involves two opposing mirrors.  The same ray may bounce between the mirrors many times before emerging from the laser.  In an optical fiber, the same ray may bounce off the inner surface many times before finally emerging from the end.

These systems require special handling.  Rayica makes a fundamental distinction between sequential and nonsequential systems.  This chapter discusses how and when Rayica invokes nonsequential ray-tracing.  It also covers the speed optimizations possible in Rayica, including TurboTrace.  These optimizations depend heavily upon the nature of the ray trace.  Speed is critically important for lens optimization problems.

Another category of systems that merit study of Rayica's ray-tracing behavior are modular optical systems.  Such systems facilitate design reuse.  You can precompute the fixed parts of a system in order to accelerate design iterations on the variable part.

Sometimes a component presents many faces to the external environment, such that light can enter from any of several bounding surfaces; only one of these surface illuminations has any bearing on the analysis.  Light entering through other surfaces is irrelevant to the design.  This situation also benefits from a knowledge of Rayica's ray-tracing details.  You can restrict Rayica's analysis to the entrance surface of interest.

None of these scenarios presents any difficulty to Rayica, but all of them require a measure of forethought tied to an understanding of the Rayica features that handle them.

Go to list of topics

3.2 Sequential Tracing and Light Path

Sequential ray-tracing is the default mode in Rayica.  The defining characteristic of this mode is that Rayica determines the order of surface intersections in advance, by inspection of input, without respect to physical geometry.  The system listing given by the user input fixes the order of intersection tests.  A given element may only transmit rays to components following it in the list.  Rayica checks the possible intersections in order, and does not backtrack.  Sequential ray-tracing is the most efficient mode of operation for this very reason.  Speed is the primary reason for using it.

3.2.1 Light Source Ordering Example

Sequential tracing demands that your system enumeration matches the light path through the system.  You must order the elements of your system in correct sequence.  Otherwise you will not obtain a correct trace.

Here is a correct system with two sources and two lenses.  The system elements appear in proper order for sequential ray tracing.  The sources precede the lenses, and the lenses appear in light-path order. Each beam refracts through each lens.

In[2]:=

Sequential tracing works correctly when the system list is properly ordered.

A seemingly minor alteration to the previous input will cause the sequential trace to fail.  If redSource follows L1 in the system enumeration, then redSource rays will ignore L1.  Only components listed after redSource can affect its rays.  You are not changing physical geometry here, merely reordering input.  The system list contains exactly the same elements, all in exactly the same geometric positions, the order of the Mathematica input has changed.

In[5]:=

L1 now fails to refract the lower beam from redSource because the system list does not follow light-path order.

The lower beam from redSource acts as though L1 does not exist.  The beam refracts only through lens L2.  The lenses have not moved, and the beam's point of origin has not changed, but its computed light path has changed.  This trace is therefore incorrect.

Go to list of topics

3.2.2 Component Ordering Example

Components too must appear in proper order for sequential ray tracing.  The system below has one source and four components, including an invisible Boundary.

In[8]:=

The prism functions correctly when the system list order mimics the light path.

Pay particular attention to the prism.  The only system elements that can transmit rays to the prism are those preceding it in the system list.  Observe what happens when you list mirror before prism.  The same elements still precede prism, but with the addition of mirror to the set.

In[15]:=

Putting mirror ahead of prism causes the sequential trace to fail.  The prism does not refract the beam.

The prism no longer affects the rays.  The trace would be identical if the prism was removed altogether.  The same elements still precede prism in the system enumeration, but the trace has failed.  Why?

The trouble arises from sequential processing.  In this kind of tracing, Rayica considers the system enumeration as a checklist of candidates for intersection or ray generation.  There is no backtracking.  Once Rayica reaches a certain point in the checklist, it does not return to previous items.

In this case the source is the first element, so Rayica generates its rays.  Rayica is then finished with source and checks it off the list.  The candidates are now, in order, {lens, mirror, prism, boundary}.  Rayica checks for intersections with lens, finds them, and checks lens off the list.  Now the remaining candidates are, in order, {mirror, prism, boundary}.  Rayica checks for intersections with mirror, finds them, and checks mirror off the list.  Rayica fails here to consider prism as a candidate.  This failure is an artifact of sequential processing.  From a processing viewpoint, the mirror is first in line, and therefore has priority.

Now the rays have physically bypassed prism and bounced off mirror.  Rayica has only {prism, boundary} remaining as candidates to the list.  The rays from mirror do not intersect prism, so Rayica checks prism off the list.  This action leaves only {boundary} on the list.  The rays from mirror do indeed intersect boundary.  Rayica finds their intersection and checks boundary off the list.  Now the checklist is empty, so Rayica terminates the trace.  (You generally always want a Boundary object as the last element in your system list.)

There are two ways to repair this failed trace.  One is to correct the system order to follow the physical light path.  Yet in many systems, the light path is more complicated than this simple example suggests.  A more general solution is Rayica's Resonate function.  Resonate instructs Rayica to perform nonsequential ray tracing for the set of components that are supplied as arguments.  In this mode, all surfaces are candidates for intersection at all times.

We wrap Resonate around the system components without changing their order.  Resonate accepts a list of components as its argument.  The argument here is the list {lens, mirror, prism}.

In[16]:=

The Rayica function Resonate enforces nonsequential ray tracing.  The trace is now correct.

Components inside the list argument to Resonate may appear in any order.  The Resonate function effectively becomes just another Component object in your system list.  The trace is sequential before and after the Resonate component.  Inside the Resonate component, it is nonsequential.

Go to list of topics

3.3 Nonsequential Tracing and Resonate

Sequential ray-tracing cannot handle all possible systems.  An important category of systems not amenable to sequential processing is that of resonating cavities.  A general resonating cavity is any component or system that confines light rays temporarily or indefinitely, not transmitting them directly.  This definition encompasses laser cavities, optical fibers, integrating spheres, and so on.  Such components and systems cause light to wander within their boundaries prior to escape.

More specifically, from the viewpoint of Rayica, a resonating cavity is any component or system for which it is possible that the same ray may hit the same surface (a) more than once or (b) from more than one other surface.  Proper analysis of a resonating cavity requires nonsequential ray tracing.  You instruct Rayica to apply nonsequential methods by means of the Resonate function.  This function creates a new Component object out of the constituents of the cavity system.

3.3.1 Three-mirror Resonator

The following mirror cavity demonstrates the need for the Resonate function.  First we trace the system without calling Resonate.  The first attempt will be a sequential ray trace.  The source is a SingleRay.  Note that the mirrors have no thickness.  In Rayica terms, each of them is a single surface, not a double surface.

In[17]:=

Three mirrors in a ring arrangement constitute a resonating cavity.  The trace fails in this case because Resonate was not used.

Instead of circulating between the mirrors, the ray terminated at the Boundary.  Not instructed otherwise, Rayica assumed that the trace was sequential.  The splitter was the first element encountered.  One transmissive and one reflective ray were created at the splitter.  Rayica then searched the system list for candidate components which might intersect these secondary rays.  The only component appearing in order after splitter was Boundary.  Rayica therefore intersected the two secondary rays with Boundary and terminated the trace.  Physically, the reflected secondary ray should have struck M1.  However owing to its position in the system list, M1 was no longer a candidate for intersection.

A natural temptation would be to reorder the system list to place M1 in sequence after splitter.

In[24]:=

The sequential ray trace does a little better with a revised system order.  However this trace is still incorrect because there should be a second encounter with the beam splitter.

The reordering has helped somewhat.  The ray correctly bounces off all three components.  The problem is that the ray should encounter the beam splitter a second time after its reflection from M2.  However it does not because this is still a sequential trace, and Rayica considers the trace complete once the system list is exhausted.  After M2, all that is left in the system list is boundary.  Rayica therefore intercepts the ray from M2 at the boundary and terminates the trace.

What you need is a way to ask Rayica to trace nonsequentially.  That is, Rayica should not ignore a component just because it was already hit once or appears in a certain input order.  The Resonate function causes Rayica to treat its arguments as a single Component with the property of a resonating cavity.

In[26]:=

Now we have a true resonating cavity by means of Resonate.  This trace is correct.

Resonate properly handles this cavity.  Rays circulate within the cavity and occasionally make their way out through the beam splitter.

Any time you ask Rayica to trace a resonating system, there is a question of when to terminate the trace.  In principle it could go on forever.  Rayica creates a fresh Ray object at each new intersection.  Some of the Ray objects do leave the cavity through the splitter.  However, many more are simultaneously created as reflected rays.  These could potentially circulate endlessly.

There are two options used by Rayica to regulate such traces: GenerationLimit and ThresholdIntensity.

One way that Rayica terminates such a trace is by restricting the number of Ray creation events.  The maximum number is given by GenerationLimit, an option to AnalyzeSystem.  The default value of the limit is moderately large, but it can be set to any value.

In[28]:=

The effect of lowering the ray GenerationLimit for the resonating cavity.

A second way that Rayica terminates such a trace is by monitoring the Intensity values of the circulating rays.  The minimum Intensity value permitted for a Ray is given by ThresholdIntensity, also an option to AnalyzeSystem.  For lossy systems, such as this cavity example, this provides a more natural ending to the trace.

Go to list of topics

3.3.2 A Subtle Substitute for Resonate

Sometimes you need Resonate even for mundane geometries.  If the light path backtracks on itself, you generally need Resonate.  However a subtle trick can allow you to trace some of these systems sequentially.  This trick accomplishes the same results as Resonate but without the speed penalty.  It requires careful analysis of the light path.  The idea is to include the same component multiple times in your system list.

Here is a system of one source and three components traced sequentially (without Resonate).  Sequential tracing cannot correctly analyze this system as given.

Even a simple geometry may demand nonsequential ray tracing.  The reflected rays in this sequential trace do not refract in a backward direction.

The trace has failed.  The reflected rays should have refracted through the lens a second time after reversing direction at the mirror.  Instead, the rays ignored the lens after the reflection.  We now call Resonate to advise Rayica that the system includes a nonsequential section.  This will repair the trace, but at the cost of slower execution associated with nonsequential tracing.

In[8]:=

Nonsequential tracing with Resonate handles this case properly.  Nonsequential tracing, however, incurs a speed penalty.

Resonate is generally safe under all circumstances; it has correctly traced this system.  However, it has also increased the execution time.  With a little care, we can attain the same results using sequential ray-tracing.  The light path crosses the lens two times.  Therefore, sequential tracing will work if we enter the lens two times in the system list.  One lens entry comes before mirror, and the other comes after mirror.

For speed optimization, you can trace resonating systems sequentially if you pay careful attention to the light path and system list.

This trick is rather subtle and not recommended until you are well acquainted with Rayica.  It requires a careful, accurate assessment of the light path for your system, and fastidious attention to the order of your system list.  By contrast, Resonate accepts components in any order.

Go to list of topics

3.4 Rays and Surfaces

Rayica records information about ray traces by means of Ray objects.  Rayica creates these objects dynamically to monitor the progress of the trace.  When the trace is complete, Ray objects serve as data repositories for ReadRays, the main function for system examination.  Ray generation also has ramifications for modular systems, discussed in Section 3.6.

Surfaces constitute a major topic in Rayica.  Certain issues surrounding component surfaces may influence your ray-tracing results.  Moreover, you can customize Rayica's behavior with various surface-related function calls.

3.4.1 How Rayica Generates Rays

Any ray trace spawns a certain number of Ray objects.  Each Source generates its particular set of Ray objects.  These propagate to the surfaces of various Component objects in the system.  At every intersection of a Ray with a Component, Rayica determines what action the Component has on the light.  Common interactions include reflection, refraction, beam splitting, and absorption.

Rayica embeds intercept data within the Ray object.  This data includes angle of incidence and point of intercept.  However, recording the entire optical action of a Component necessitates more software activity than a mere alteration of one Ray.  A solitary Ray cannot embody all of the relevant information.  Among other issues, Ray can only record a single straight path, not a complete, piecewise-linear path.  Rayica therefore creates multiple Ray objects to capture the total optical effect.

When Ray intersects a PlanoConvexLens, for example, Rayica registers the intercept location inside Ray, then creates another Ray to follow a refracted (but straight) path through the lens.  This second Ray encounters the back surface of the lens, where it refracts in turn.  Rayica registers its intercept location and creates a third Ray that leaves the lens.  So three distinct Ray objects record the overall light path.  This case is illustrated below.

In[30]:=

Rayica creates Ray objects 2 and 3 to record the optical effect of lens L on Ray number 1.

Two Rays were created where one existed before.

When Ray intercepts a Mirror, the result is similar.  Rayica registers the intercept and creates a reflected Ray to follow the new light path.  There are two Ray objects in this case.  A half-silvered mirror generates both a reflected and a transmitted Ray.

Go to list of topics

3.4.2 Component Surfaces

A primary task of Component objects is to maintain information about surfaces.  In general, you may view Component as a set of surfaces.  Component can have anywhere from one to dozens of surfaces.  The surfaces may be physically internal or external to the component.  Each surface has a parametric surface description, a surface normal function, and various other properties.

Surfaces within Component are numbered.  The number of a surface is called its SurfaceNumber.  There is no special meaning to surface numbers beyond their mutual uniqueness within Component.  Most Rayica components have a canonical surface ordering which assumes some default placement of the component.  This numbering system is merely convention; it has no effect on ray tracing.  Rays may enter through any external surface, whatever its number.

The reason to be concerned about surface numbers is for extraction of post-trace data.  When Ray intersects a surface, Rayica records certain identifiers into Ray that pertain to the intersection.  These include ComponentNumber, for the particular component within the system, and SurfaceNumber, for the particular surface within the component.  You often use these numbers in conjunction with ReadRays to isolate specific intercepts for examination.

Go to list of topics

Faces and Edges

Visualizing Component objects as organized collections of surfaces, instead of contiguous solids, more directly captures Rayica's behavior.  Think of multiple surfaces floating in space and you have the idea.  Rayica never defines the volume occupied by a component, but only selected faces.  The selected faces are those important for optical design.

This convention can break down; it fails only under peculiar conditions that seldom occur in practice.  Component edges typically matter very little.  Lenses are not designed to pass light through their edges, but through their faces.  Rayical properties of lens edges are mostly irrelevant to optical design work.  For this reason, most built-in Rayica components lack edge definitions.  Thus a PlanoConvexLens has only two surfaces, one planar and one convex; Rayica omits the lens circumference.  As far as Rayica is concerned, a PlanoConvexLens is a pair of infinitely thin surfaces with open space between them.

You may construct a highly contrived optical system to demonstrate how this convention works, and where it breaks down.  Following is a PlanoConvexLens turned by 90 degrees.  This rotation orients the faces away from the optic axis.  Two light sources illuminate the lens edge-on, down the optic axis: one exactly, the other through the curved face.  The light shining exactly edge-on passes straight through, unaffected, between the two faces of the lens.

In[91]:=

The standard PlanoConvexLens has no optical edge.  The edge-on (green) light travels between the faces without refraction.  The slightly off-sides light does refract, because it interacts with one of the faces.

Rayica's graphical display indicates an edge for aesthetic reasons only.  From a ray-tracing standpoint, this PlanoConvexLens is a pair of refractive surfaces floating in space; it has no edge at all.  This fact is proven by the trace above.  Only those rays striking a lens face experience refraction.  The remaining rays pass between the two  faces without disturbance.

Whenever edges concern your system, you should carefully explore the implications of Rayica's default treatment.

Go to list of topics

Rolling Your Own Edge

In some rare cases, you may need to create modified components with proper optical edges.  Rayica allows such definitions.  They require some facility with surface data  and a bit of creativity.  This section demonstrates how to accomplish the edge-completion task for PlanoConvexLens above.  This topic is somewhat advanced and may be skipped without penalty.

We can mate the two faces of PlanoConvexLens by means of CylindricalLens, which is shaped and placed in a highly unusual configuration.  The CylindricalLens will constitute the edge of the PlanoConvexLens.  In other words, the two main faces of CylindricalLens will form a single contiguous edge for PlanoConvexLens.  The result will be a completely closed Component that properly represents the full physical behavior of a plano-convex lens.

To produce this modified component, you must obtain the parameters of CylindricalLens.  The cylindrical lens will be rotated 90 degrees relative to the plano-convex lens.  Each face of CylindricalLens will cover a 180-degree arc.  The two radii of curvature will be identical but of opposite polarities, indicating opposite directions of curvature.  The curvature radius is one-half the aperture diameter of PlanoConvexLens, or 50/2.

If you view PlanoConvexLens edge-on, the outline of its edge is rectangular.  So the face-on aperture of CylindricalLens must also be rectangular.

In[131]:=

This side view of PlanoConvexLens reveals the rectangular aperture shape required for CylindricalLens, which will form its edge.

A list of two numbers specifies a rectangular aperture.  The height of the rectangle is the diameter of the PlanoConvexLens (50 units).  The width of the aperture is more difficult to ascertain.  It is the distance between the rim of surface 1 and the flat of surface 2.  Finding this distance requires some digging into the data of PlanoConvexLens.  The thickness of PlanoConvexLens is measured from the peak of its convex face to its planar face.  The width of the CylindricalLens aperture must equal this total thickness, less the thickness of the convex surface only.  Use SurfaceFunction of the convex face (surface 1) to compute the aperture width.

In[72]:=

Out[73]=

Out[74]=

Resonate will now consolidate edgelens with the previous lens definition to create newlens as a single component.  edgelens is nudged into position by nested Move calls.

In[140]:=

A test of the new lens proves it out.  This optical system is identical to the previous case but has lens replaced by newlens.

In[142]:=

This customized PlanoConvexLens has an edge.  The edge-on (green) light now refracts.

You must ensure that both the edge and the original component are made of the same material.  Since this material was unspecified in the present case, it defaulted to BK7 glass for both the plano-convex and the cylindrical lens.

Go to list of topics

Thick Components and Internal Reflections

All surfaces in Rayica are infinitely thin.  They have no inherent thickness.  Components, in turn, are composed of surfaces.  A thickcomponent can differ from its thin equivalent only by the number of surfaces used in the model.  For instance, Mirror can be modeled as thin, with a single reflective surface, or as thick, with both a reflective and a refractive surface.

Other components maintain no such distinction between thick and thin alternatives.  Where available, thickness is usually an extra parameter.  The following table shows how this extra parameter works in several common cases.  The left column shows the thin versions, and the right column shows the thick equivalents.

Out[9]//TableForm=

The thickness parameter: single-surface components (left) and double-surface components (right).

Many components automatically incorporate multiple surfaces and are always considered thick.  This class includes standard lenses and prisms.  Two exceptions are Screen and ThinLens; these components have only one surface.

Realize that thick models can experience internal reflections.  These reflections will depend on the incidence angles of your light sources.  Moreover, there is a Boolean option for lenses called FresnelReflections that instructs Rayica to track energy losses and first-surface reflections.

Here is a simple trace showing the energy intensity of Ray propagating through a lens.  FresnelReflections is an option for the lens.

In[4]:=

With FresnelReflections→False, the ray propagates through the lens without losing energy.

The light energy is constant before, within, and after the lens.  Under this default behavior, the light loses no energy and experiences no reflection.  We modify the default behavior by setting FresnelReflections→True.  Rayica will now track first-surface reflections and energy losses.

In[5]:=

With FresnelReflections→True, Rayica monitors the energy loss and surface reflections of the propagating ray.

In this case, about 88% of the light energy still follows the original light path; the rest is scattered.  Note the internal reflections.

Sequential ray-tracing works as you might expect when there are internal reflections.  The reflected rays propagate through the system list in reverse order, starting at their originating component.  That is why the example above includes Boundary at both ends of the system list.  The first Boundary captures backward-propagating rays.  The last Boundary captures the forward-propagating rays.  It is generally good practice and always safe to put Boundary at both ends of your system list.

Go to list of topics

3.5 Special Considerations with PropagateSystem

This section covers a number of special considerations with the use PropagateSystem and AnalyzeSystem. These considerations were especially important with previous versions of Rayica before the advent of TurboTrace and its associated functions. Today, however, their understanding has become less critical since TurboTrace does not utilize them. As such, this section can be skipped without penalty.

Go to list of topics

3.5.1 QuickSurfaceSort

When propagating rays within a component, Rayica must decide in which order to check surfaces for intersection, and which of all possible intersections is correct.  This problem is similar to that of sorting components at the system level.  At the lower level of an isolated component, the internal surfaces must be sorted in some kind of priority order for intersection processing.

Rayica uses geometry to make intelligent guesses about which surfaces should be checked first.  The actual details of Rayica's workings in this regard are complex.  However, Rayica presents an option called QuickSurfaceSort that controls the procedure in a general sense.  This option applies to AnalyzeSystem and PropagateSystem.  Note, however, that the QuickSurfaceSort option does not apply to TurboTrace and TurboPlot applications, since they do not directly employ QuickSurfaceSort at all (other for gaining scout trace information from PropagateSystem). In general, TurboTrace and TurboPlot always uses the QuickSurfaceSort->False setting.

In previous versions of Rayica, PropagateSystem used QuickSurfaceSort→True as its default setting in order to maximize its ray-trace speed. Since the advent of TurboTrace, however, speed considerations for PropagateSystem have been supplanted by the concern for roboust, trouble-free use and QuickSurfaceSort→False is now the default setting. When you set QuickSurfaceSort→True, Rayica employs an efficient, but occasionally fallible, method for sorting surfaces within components.  The method compares center-points of surfaces to sort them from nearest to farthest.  Sometimes you must set QuickSurfaceSort→False to obtain a correct trace.  This setting causes Rayica to perform actual surface intersections to determine the sort order.  Naturally, this extra calculation incurs a speed penalty, but it is more robust and Rayica has available closed-form intersection formulas for certain built-in surfaces like spheres, planes, and conics.

The next example demonstrates how QuickSurfaceSort affects the outcome of a trace.  It is borrowed from the section on resonating cavities.  The difference here is that we give the mirrors thickness.  This change means that Rayica must sort the surfaces within each mirror component for purposes of intersection processing.  If a ray can intersect two surfaces, Rayica must decide which of the two it will strike.

In[48]:=

With QuickSurfaceSort→True, the trace fails in this particular geometry.

The trace is wrong.  Evidently, Rayica checked for an intersection with the back surface of M2 before the front surface.  This outcome shows that Rayica's default surface processing failed for this particular geometry.

Rayica produced the wrong outcome in spite of Resonate.  The Resonate directive ensures that all surfaces will be checked; it does not affect the order in which they are checked.  To change this order, you must use QuickSurfaceSort.  When this option is False, Rayica uses a more exact, but slower, method of sorting surfaces inside components.

With QuickSurfaceSort→False, Rayica uses a more exact method of sorting surfaces within components.  This trace is correct.

It is always safe to set QuickSurfaceSort→False, but this setting, just like Resonate, incurs a speed penalty.

Go to list of topics

3.5.2 Confine

Components can receive light from any exposed surface.  Sometimes concern focuses exclusively upon light that enters through a particular surface.  This surface might be termed an entry port for the component.  This situation immediately suggests a ray-tracing optimization; Rayica should follow only those rays that enter through this surface, ignoring the rest.  The Confine function instructs Rayica to consider only specific surfaces as candidates for light entry into the component.  Please note, however, that the following discussions are not used in TurboTrace and TurboPlot, which do not directly employ confined port information at all.  As such, this section only applies to the use of AnalyzeSystem and PropagateSystem.

Here is a custom function that creates a multi-sided prism component.  The output of this function is a set of adjoining LensSurface components consolidated by Resonate into a single, unified component.  As created, such a prism has no preferred entry port.

This is the appearance of a 9-sided PolygonPrism.

A custom prism.

The following trace shows how the prism refracts light with no entry port defined.  All exposed surfaces that happen to be facing the light will refract it.

The prism refracts indiscriminately with no privileged entry ports defined.

Suppose that surface 3 is an entry port; light impinging on other surfaces is considered irrelevant to the design.  Then you must notify Rayica to simplify the trace.  The prism is identical, but you Confine the light to enter through surface 3.

Confine causes the prism to receive light only through designated surfaces.

Light that would normally refract through other surfaces skips past the prism altogether, terminating at the boundary.  Light allowed into the prism through the entry port (surface 3) reflects internally and eventually exits.  Once inside the prism, light can exit through any surface.

You may define more than one surface as an entry port.  Instead of a single number, use list notation, Confine[prism,{3,4,5}], for example.

Go to list of topics

3.5.3 Unconfine

Sometimes you must eliminate entry port restrictions.  The component will then receive light from all exposed surfaces.  The Rayica function Unconfine accomplishes this task; it effectively cancels the action of Confine.  For instance, you might Unconfine the confinedPrism that was defined above:

Unconfine reverses the action of Confine.  The prism no longer has preferred entry ports.

The unconfinedPrism behaves as the original prism.  Light enters through all exposed surfaces.

Go to list of topics

3.5.4 Modular Systems

Rayical design, like other forms of engineering, often benefits from a block diagram approach.  A block diagram consolidates various details into conceptual blocks that perform high-level functions.  The blocks have defined interfaces, but their internal workings need not concern the designer.  Modular blocks ease system design tasks.

In terms of ray-tracing, the most important aspect of such a block is that  it can be traced ahead of time.  Rayica features an option called ModularSystem for this purpose.  This boolean option causes Rayica to create extra rays at the system output, poised for tracing into the next module.  This mechanism essentially turns a modular optical system into a light source for downstream modules.  Rayica need not recompute the first module when you trace the overall system.  This computational savings translates into better performance.

Under normal conditions ModularSystem→False, and Rayica terminates ray creation at the final surfaces.  The last intersection surfaces kill the rays.  However, when ModularSystem→True, Rayica traces through the final surfaces and leaves the final rays ready for tracing into downstream modules.  These rays encode intensity, direction, and intersection points at the output surfaces of the modular system.

A very simple modular system comprises SingleRay and Screen.  ReadRays will extract geometry data for each ray created by this system.

In[27]:=

Out[27]//MatrixForm=

The system has only one ray. You know this fact because two data per ray were requested, start and end, and two data were obtained.  Now set ModularSystem→True.

In[26]:=

Out[25]//MatrixForm=

The system now has two rays.  This fact is not evident from the graphics.  However, ReadRays reports two rays: one from the origin to the screen, and another essentially sitting on the screen, having equal start and end points.  This extra ray is the output of this modular system.  Rayica has primed an exit ray with data concerning the final intersection of the modular system.

You may assign any type of component as the output of a modular system.  There can be any number of output components.  Some make more sense than others.  A typical alteration made when converting systems to modular format is the replacement of a Boundary with a ClearBoundary.  The difference between these components is that a ClearBoundary does not absorb the light but merely intersects it.  For purposes of modularity, Rayica can record the intersection data there without altering the light itself (as would a Boundary).

The modular system concept is best demonstrated by example.  Here is a practical modular system, a beam-expanding spatial filter.  Note the use of a PinHole and the final Screen.

In[33]:=

A beam-expanding spatial filter becomes a modular system building-block.

Now create a collimator by feeding the output of spatialFilter into a collimating lens.  In doing so, there is no need for light sources.  The light source was defined inside the spatialFilter (a different approach will be shown below).  Furthermore, this new system will itself be modular.  This method follows the block-design approach.  Each system has a specific purpose.

In[34]:=

The spatialFilter module placed within a collimator system, itself also modular.

You may wish to study internal reflections within the collimator.  Create a system comprising the collimator and another off-axis spatialFilter.  A Boundary delimits both ends of the system list in order to catch both forward- and backward-propagating rays.

In[36]:=

The collimator and spatialFilter modules functioning within a system under study.

The Mathematica variables collimator and spatialFilter constitute system modules that can serve in any capacity.  They both contain pre-traced light rays.  In the example above, the only rays Rayica had to trace were those emerging from the output of the off-axis spatialFilter.

When speed is of less concern than flexibility, you can define modular systems lacking such pre-traced rays.  The basic idea is to compose a Mathematica list of components (and possibly sources as well) that serves as a modular system.  The difference here is that the list is just a list, not an OpticalSystem.  The collimator defined above is an OpticalSystem.

In[39]:=

Out[39]=

Yet, you may freely redefine Mathematica variables, and we do so to demonstrate another method of modular design.  This method does not depend on the ModularSystem option in any way.  The ModularSystem option applies only to OpticalSystem objects.  The following represents a completely different modular design approach.

In[40]:=

The same system built using a different modular method.

The trace is identical but consumes more time, for Rayica must trace the light through each and every component.  There is no pre-trace data.  Some computational savings was achieved by removing the output components from the modules.  The spatialFilter previously had Screen, and the collimator had ClearBoundary.  These components were not required since ModularSystem was not in use.

The question of whether to use the ModularSystem approach or the type shown revolves around the subject of light inputs.  The ModularSystem option effectively turns upstream modules into light sources for downstream modules.  It only makes sense to define an effective light source of this nature if its own light source is fixed.  Therefore, you should use the ModularSystem method if the input light to the module is fixed.  When, instead, the input light sources will be moved about, no pre-trace is really possible.

Go to list of topics

3.6 TurboTrace

TurboTrace speeds Rayica's ray-tracing performance up to 100 times beyond standard routines (and >1000 times in special cases).  The suite of TurboTrace functions is, perhaps, the most important branch of Rayica.  As a whole, TurboTrace and its allied functions constitute an alternative ray-tracing system within Rayica.  Yet this system merely builds upon concepts already presented in this chapter.

3.6.1 Basics of TurboTrace

TurboTrace incorporates compiled functions plus certain structural and numerical optimizations to maximize tracing speed.  Compiled functions utilize a special bytecode system that bypasses the standard Mathematica interpreter.  These bytecodes embody arithmetic and branching operations, much like assembly-language opcodes.  (See The Mathematica Book for more information.)

The implementation details of TurboTrace are transparent and largely unimportant for proper usage.  The major concern with this function is to understand its basic principles of operation.  The interfaces to the turbo system are practically identical to those of the standard system.  In general, you should use TurboTrace whenever possible.  The following equivalency table shows the primary interfaces to the turbo system. Note: this table should be taken conceptually and not literally.

TurboTrace works with the same kind of system descriptions as the standard functions.  However, TurboTrace differs from them in one important respect, its assumption about surface encounters. In particular, TurboTrace offers two different ways of determining optical surface encounters during a ray trace. These two methods are selected with the TurboTrace option settings: SequentialTrace and ScoutTrace.  As its default setting, TurboTrace uses SequentialTrace -> False and ScoutTrace -> Automatic. With this setting, TurboTrace traces the rays through optical surfaces in similar fashion as PropagateSystem, described previously. SequentialTrace -> False provides the most robust tracing environment.  However, in order to achieve maximum the ray-trace speed, SequentialTrace -> True is required. With SequentialTrace -> True, TurboTrace sorts the optical surfaces in your system by means of a scout trace.  This preliminary trace involves just a handful of scout rays but fixes the order of surface encounters for all other rays.  Internally, TurboTrace calls PropagateSystem to compute the scout rays.  These yield an ordered list of surface encounters.  Armed with this list, TurboTrace constructs a compiled function for ray tracing called, RayTraceFunction. Once invoked, RayTraceFunction sends the mass of remaining rays through ordered list of surface encounters, under the assumption that all rays encounter the same surfaces, in the same order, as the scout rays.  We refer to these main rays as turbo rays.

With SequentialTrace -> True, the fact that TurboTrace uses PropagateSystem to determine surface encounters means that there is no mystery about the path.  Everything that applies to PropagateSystem, as described in previous sections, applies to the scout trace.  Fundamentally, a scout trace can fail.  If the scout trace is wrong, so is the entire TurboTrace. Sometimes TurboTrace fails even after a successful scout trace.  This situation arises when the turbo rays have a wide crosswise dispersion relative to the optics.  A paraxial scout ray within a wide beam may follow the system quite well, while some stray edge ray may not.  Nonetheless, TurboTrace wrongly assumes that it does.  Therefore you must exercise caution with TurboTrace and SequentialTrace -> True.  You can often remedy these kinds of problems by manual specification of the scout rays.  For example, an automatic scout ray traced into a Cassegrainian telescope would be paraxial and thus strike the back of the secondary mirror, yielding incorrect scout data.  In this case you would define scout rays that bypass the central obstruction of the secondary mirror.  Alternatively, you can simply use the SequentialTrace -> False setting.

In addition to the SequentialTrace option, the auxiliary ScoutTrace option determines whether a scout trace is used in TurboTrace. With SequentialTrace -> True, a scout trace always occurs. However, with SequentialTrace -> False, a scout trace does not occur with the default ScoutTrace -> Automatic setting, but does occur with ScoutTrace -> True. In this later setting, only those surfaces that intersect the scout ray get included into the constructed RayTraceFunction.

Before we can examine TurboTrace in greater detail, we first need to explore how symbolic parameters are specified in Rayica.

Go to list of topics

How to Enter Symbolic Parameters in Rayica

The basic Rayica system permits symbolic entry of most system parameters. These parameters may then be used to construct symbolic variables within the system that later become the input parameters to the compiled RayTraceFunction. In many cases, the symbolic parameter "piggy-backs" together with a numerical value in the form of an ordered list. In general, Rayica requires that all symbolic parameters are also accompanied by a numerical value. Here is an example of symbolic parameter entry within PlanoConvexLens.

In[27]:=

symbolicLens = PlanoConvexLens[{f,100},{a,50},{t,10}]

Out[27]=

In other cases, the symbolic parameter is given in the form of an option. This occurs in particular for numerical option parameters. However, in this case, the symbolic value does not piggy-back with the numeric option, but instead gets its own option name that is prefixed with the word "Symbolic" in front of the corresponding numerical option name.  Two such options are SymbolicWaveLength and SymbolicIntensity. An exception is NumberOfRays, which must be a numeric at all times.

In[28]:=

LineOfRays[10, SymbolicWaveLength->λ, SymbolicIntensity->i, NumberOfRays->4]

Out[28]=

Most types of light source functions cannot accept symbolic values for its built-in parameters. The LineOfRays source given above is one such example. Here, the linewidth parameter must be a numeric value. In error message is generated otherwise.

In[30]:=

LineOfRays[{s,10}]

Out[30]=

An exception is the GaussianBeam source, which can take direct symbolic parameters.

In[29]:=

symbolicSource = GaussianBeam[{w,10},{div,.001}, SymbolicWaveLength->λ, NumberOfRays->4]

Out[29]=

The Move function accepts symbolic parameters as long as a numeric value is also given with each one. Angles may also be specified in a similar fashion. Remember, that any angular information will use units of Degrees for Rayica's Move functions. This is in contrast with Mathematica's built-in trigonometric functions, which takes radians by default.

In[32]:=

symbolicLensSystem =
    {symbolicSource,
    Move[symbolicLens, 50],
    Move[Screen[50],{{x,200},{y,0}},{α,10}]}

Out[32]=

Once a system has been defined with symbolic parameters, you can use TurboTrace to make compiled ray-trace code that uses these symbolic values as input parameters.

In[33]:=

turbotrace = TurboTrace[symbolicLensSystem]

Out[33]=

Here TurboTrace has generated a fully-compiled, RayTraceFunction that is returned along with a list of rules that describe the RayTraceFunction construction. RayTraceFunction contains all necessary information to carry out a complete ray trace operation. It uses all of the symbolic variables defined previously in symbolicLensSystem as input parameters. The SymbolicValues rule holds default numeric settings for each symbolic parameter. This was extracted from the numeric values originally given for each parameter. The EmbedRays, EmbedThresholdIntensity, and EmbedGenerationLimit options tell us whether each of these parameter types have been embedded within the compiled ray-trace code or whether the parameter has been left as an input variable to RayTraceFunction.

Go to list of topics

3.6.2 Lens Optimization

Iterative design procedures involving repeated traces will benefit substantially from TurboTrace.  Lens optimization falls into this category.  This section presents an example of lens optimization.  The objective is to find optimal surface curvatures for a lens focusing planar waves onto a screen at a fixed distance.

The optimization parameters are the radii of curvature for the front and back lens surfaces. Begin with each set to 100 mm, opposite polarities indicating opposite directions of curvature.  A Baffle fixes the desired focal plane. Here we assign symbolic parameters, r1 and r2 to the radius-of-curvature parameters of the front and back lens surfaces.

In[95]:=

A system in need of optimization.  The light does not focus optimally at surface B owing to unoptimized radii of curvature.

At this point, the variable sys is neither an OpticalSystem nor a TurboSystem.  It is merely a list of components.  You now require a means to alter the lens between traces, monitor its performance, and find the optimal solution for parameters, r1, and r2.  With this specific need in mind, Rayica already has a built-in function, called OptimizeSystem, which does this very thing. Next we demonstrate OptimizeSystem. In order to achieve the fastest answer, we employ SequentialTrace->True.

In[97]:=

soln = OptimizeSystem[sys, SequentialTrace->True]

Out[97]=

On the surface OptimizeSystem appears to be very simple. Underneath, however, OptimizeSystem is managing a number of complex activities. First, OptimizeSystem calls the ConstructMeritFunction, which, in-turn, internally calls TurboTrace with the necessary settings to construct a compiled merit function, called a RayTraceFunction, that performs a self-contained ray-trace of the components given by sys. Then, in its default mode, OptimizeSystem passes this RayTraceFunction to either FindMinimum or NMinimize, in order to determine an optimal solution. The resulting answer for the symbolic parameters is given in the returned SymbolicValues rule. Here, the NumberOfCycles tells how many iterations were required to find the optimal solution and FinalMerit indicates the optimal merit result.  Realize that FinalMerit is technically a demerit result because larger values are worse.

The final lens design incorporates the results of our search.  These are the curvatures.

Here is a TurboPlot of the optimized lens with the SymbolicValues option.  Note TurboPlot responds much more quickly when it is called a second time with a previous TurboPlot result. In fact, it responds much more quickly than working with AnalyzeSystem.

In[98]:=

TurboPlot[turbosys, PlotType->TopView, soln];

Now the light focuses optimally at surface B.

Finally, a comment about optimization in Mathematica: FindMinimum and NMinimize are merely two of a growing collection of optimization routines availableConstructMeritFunction for Mathematica.  FindMinimum and NMinimize are Mathematica's two primary built-in function for optimization, but they is not necessarily the most efficient or most effective. There also exists a number of third-party optimization packages for Mathematica.  Some packages offer the promise of global optimization while others offer greater efficiency and speed.  In general, different types of problems will require different types of optimization solutions.  Unfortunately, it is not possible to make a single optimization algorithm that performs optimally for all problems.  The trick is to find the one that works best for your particular system.  Nevertheless, FindMinimum and NMinimize both provide reasonable starting points because of the general nature of their design.

Go to list of topics

3.6.3 ConstructMeritFunction

The following examples demonstrate the use of ConstructMeritFunction to perform application-specific, high-level, ray-trace calculations. Here is Rayica's built-in explanation for ConstructMeritFunction:

?ConstructMeritFunction

ConstructMeritFunction constructs a single compiled function that combines ray-tracing and post-processing of the ray-traced data. You can use ConstructMeritFunction to perform optimizations and perturbation analysis on various optical parameters. Internally, ConstructMeritFunction is used by OptimizeSystem for such purposes. We will now use ConstructMeritFunction to obtain two different types of measurements on an imaging system. We will first measure the optical path differences of a point-source image and test how this changes with different lens positions, which gives us an indication of the optical phase aberrations present in the system. Then, we will measure the projected focal point position on a screen for different lens positions. First, however, we need to define an optical system.

For this optical system, we will include three symbolic parameters that we can manipulate later on. Two of these parameters, y1 and y2, correspond with the lateral positions of the two lenses and the third parameter, x, corresponds with the longitudinal position of the second lens.

In[25]:=

lenssystem =
    {WedgeOfRays[20, NumberOfRays->11],
    Move[PlanoConvexLens[100,50,10, CurvatureDirection->Back],{90,{y1,0}}],
    Move[PlanoConvexLens[100,50,10],{{x,140},{y2,0}}],
    Move[Screen[75],250]};

In[26]:=

AnalyzeSystem[lenssystem,PlotType->TopView];

Go to list of topics

Merit based on summation of optical path differences

We are now ready to use ConstructMeritFunction to calculate optical path differences. For this we need to use the MeritType option with ConstructMeritFunction.

?MeritType

From this we can see that to calculate optical path differences, we must use the MeritType -> OpticalLength setting with ConstructMeritFunction. Next, we call ConstructMeritFunction:

result = ConstructMeritFunction[lenssystem, MeritType->OpticalLength]

ConstructMeritFunction works by internally calling TurboTrace with the ReportedParameters and ReportedFunction options. Here, we can see that OpticalLength parameter is used in ReportedParameters. The ReportedFunction setting contains a small subroutine that uses the OpticalLength result. It first calculates the minimum optical length and uses it to find the optical path difference. It then calculates the sum of the squares across all of the optical path differences. From this a single number is returned that provides a measure of optical phase aberration. The RayTraceFunction, passed back from ConstructMeritFunction, is custom-built to measure the combined optical path differences. As such, we will assign the RayTraceFunction to the meritfunction variable:

meritfunction = RayTraceFunction/.result;

In order to obtain a measurement, we must assign three numeric values to the meritfunction for the three symbolic parameters: {x, y1, y2}:

meritfunction[130,0,0]

Here we see that for x = 130, y1 = 0, and y2 = 0, a single number is returned that corresponds with the combined optical path differences of the traced rays in the system. Next, we use the Plot function to construct a plot of the optical path variance as a function of the longitudinal lens position, x.

Plot[meritfunction[x,0,0],{x,120,180}];

From this plot, we can see that path variance is a minimum near x = 152. This gives us the best focus for the system. For the remaining calculations, we will fix the value of x to 152 and examine the effects of y1 and y2 on the measurement. For this, we use the Plot and ContourPlot functions on y1 and {y1,y2}.

Plot[meritfunction[152,y1,0],{y1,-5,5}];

ContourPlot[meritfunction[152,y1,y2],{y1,-5,5},{y2,-5,5}];

From these measurements, we can see that the phase aberrations are minimized with y1 = y2 = 0.

Go to list of topics

Calculation of the beam centroid coordinates

ConstructMeritFunction works by constructing and passing ReportedParameters and ReportedFunction options to TurboTrace. As shown in the previous example, ConstructMeritFunction uses the MeritType option to determine the construction of ReportedParameters and ReportedFunction. For such built-in operations, this simplifies the overall task. In some situations, however, ConstructMeritFunction does not have a ready made MeritType setting to perform a needed operation. In these instances, it is necessary to directly specify the ReportedParameters and ReportedFunction options rather than using MeritType. One such calculation is the measurement of the focus centroid coordinates on a surface. In this instance, we will use ReportedParameters -> SurfaceCoordinates in order to work with the SurfaceCoordinates information to recover the ray-intercept coordinates. In order to obtain the beam centroid position, we will need to calculate the average of the ray-intercept coordinates. This is accomplished with Function[Apply[Plus,#]/Length[#]] and passed in the ReportedFunction option.

result = ConstructMeritFunction[lenssystem, ReportedParameters -> SurfaceCoordinates,
            ReportedFunction -> Function[Apply[Plus,#]/Length[#]]]

In such instances, although you can still use ConstructMeritFunction, it is not actually necessary and you can, in fact, call TurboTrace directly for yourself:

result = TurboTrace[lenssystem, ReportedParameters -> SurfaceCoordinates,
            ReportedFunction -> Function[Apply[Plus,#]/Length[#]]]

Here, you can see that the returned RayTraceFunction is exactly the same with either ConstructMeritFunction or a direct call to TurboTrace. The RayTraceFunction is now custom-built to measure the average focus spot position. As such, we will assign the RayTraceFunction to the focus variable:

focus = RayTraceFunction/.result;

As in the previous example, in order to obtain a measurement, we must assign three numeric values to the focus for the three symbolic parameters: {x, y1, y2}:

focus[152,0,10]

This time, the resulting calculation returns a list of two numbers that correspond with the focal coordinates on the final surface. Next, we can use Plot to make a plot of the first coordinate as a function of y2 by setting x, and y1 to fixed values. We use First to extract the first coordinate of the result.

Plot[First[focus[152,0,y2]], {y2,-10,10}];

From this plot, we notice that something odd occurs after y2 = 8. To investigate this glitch further, we can call TurboPlot with these conditions.

In[27]:=

plot = TurboPlot[lenssystem,SymbolicValues->{x->152,y1->0,y2->10},PlotType->TopView];

Here we see that the lens has moved so far that a ray has missed the lens entirely! This, of course, skews the calculated focal positions. To see the effect of lens position in more detail, we can use Table to construct an animation of the lens movement:

In[28]:=

Table[TurboPlot[plot,SymbolicValues->{x->152,y1->0,y2->y},PlotType->TopView], {y,0,10,1}];

Double-click above to obtain an animation.

Go to list of topics

Intensity-weighted beam-coordinate calculation

In the previous example, we did not use the intensity information from the ray trace in order to calculate the beam spot centroid position. Each ray intercept was given an equal weighting regardless of its intensity value. However, there is no difficulty in including the intensity information into ReportedParameters and ReportedFunction. We can learn more about this process if we pass only ReportedParameters -> {SurfaceCoordinates, Intensity} to TurboTrace without including a ReportedFunction setting. This allows us to see how the parameters are initially passed to ReportedFunction.

In[32]:=

(RayTraceFunction/.TurboTrace[lenssystem, ReportedParameters -> {SurfaceCoordinates, Intensity}])[152,0,0]

Out[32]=

Here we see that the SurfaceCoordinates and Intensity parameters are returned in a nested list of three ordered numbers. This helps us to see how to construct a new ReportedFunction setting that uses this input structure to calculate the intensity-weighted centroid position. We will assign this function to the reportedfunction variable.

In[39]:=

reportedfunction =
Function[
     {TV1,TV2} = Apply[Plus,Map[(#[[3]]*{#[[1]],#[[2]]})&,#]];
     TV3 = Apply[Plus, #][[3]];
     {TV1,TV2}/TV3
];

Here we have used some intermediate variables {TV1, TV2, TV3} in order to break the ReportedFunction calculation into several simplified steps. First, we weighted the SurfaceCoordinate values by the Intensity and assign its combined result to the {TV1, TV2} variables. Next, we calculate the total integrated Intensity and assigned this value to TV3. Finally, we divided the weighted SurfaceCoordinate values by the total Intensity. We are now ready to call TurboTrace and construct our RayTraceFunction.

In[40]:=

result = TurboTrace[lenssystem, ReportedParameters -> {SurfaceCoordinates, Intensity},
    ReportedFunction -> reportedfunction]

Out[40]=

In[41]:=

focus = RayTraceFunction/.result;

In[42]:=

focus[152,0,5]

Out[42]=

Go to list of topics

Rayica's built-in temporary variables

Rayica uses a number of built-in temporary variables for its compiled ray-trace calculations. These variables are specially designed to work with RayTraceFunction. Unlike most variables in Mathematica, which can accept different types of formats, these variables only designed to take single real-number assignments. If a vector assignment is required, then a list of variable names must be used instead. Such variables start with the TV prefix followed by a numeric value that runs from 1 to 10: TV, TV1, TV2, TV3, TV4, TV5, TV6, TV7, TV8, TV9, TV10.

Go to list of topics

3.6.4 Wavefront Analysis

Design procedures involving large numbers of rays during a single trace will also benefit substantially from TurboTrace.  Wavefront analysis falls into this category.  The objective here is to find the optical path function at a surface in the system.  This in turn leads to determination of the phase front and, finally, complex field at the measurement surface.

This section presents two examples that use wavefront analysis.  In the first example, a two-dimensional analysis on a beam-expander/ collimator system is performed.  This allows you to introduce the basic techniques of wavefront determination.  In the second example, a three-dimensional wavefront analysis using the same optical components is performed.  The optical system of both examples has been chosen such that the geometric optical path length dominates the wavefront behavior and that diffraction due to aperture edges can be safely neglected.  You will find other, more advanced, wavefront analysis examples in the Help Browser that include diffraction effects, arbitrary precision calculations, and even interference between overlapping wavefronts.

Go to list of topics

Two-Dimensional Analysis with Interpolation

This first example performs the wavefront analysis in a two-dimensional plane.  This allows you to introduce basic methods before entering into a full 3D wavefront analysis.

First set up the system definitions and use AnalyzeSystem with a small number of rays to get a display of the basic system.

In[22]:=

The basic beam expanding collimator system. The input source is LineOfRays. The wavefront impinging on the measurement surface is not quite planar.

Here, two sources that differ only in the number of rays invoked are defined: the basicSource is used for visualization purposes only, and the denseSource is used in the actual wavefront analysis.  The optics consists of a tiny negative lens in series with a large plano-convex lens, followed by a boundary.  The analysis surface will be the far-right boundary of the system.  Implicated in this example are the physical units of millimeters.  The optics to exhibit a certain amount of residual phase curvature at the measurement surface are intentionally chosen.  This allows for a better demonstration of the basic method of wavefront analysis.

Next, trace more rays through the system.  If PropagateSystem is used, it would take a very long time to finish the trace, so TurboTrace is used instead.  TurboTrace normally keeps the last surface information and discards the findings of earlier surfaces.

In[29]:=

Out[29]=

You can now view the result with TurboPlot, but here only the last ray segments appear, since the previous ones have been automatically deleted during the trace.  You can otherwise use the ReportedSurfaces option with TurboTrace to specify which ray-surface results are retained.  The default setting is ReportedSurfaces->Last.  ReportedSurfaces->All would retain all ray-surface segments.

In[30]:=

A call to TurboPlot shows you only the last ray segments. The preceding ray-surface information was automatically discarded by TurboTrace to save memory.

In order to calculate the wavefront information, you need the optical length of each ray at the last surface.  Instead of calling ReadRays, use ReadTurboRays because it is much faster for turbo-rays (ReadTurboRays works only with TurboTrace results, although ReadRays can work with either TurboTrace or PropagateSystem).  You will need the position of each intersection point on the surface as well as the optical length, so ask for {SurfaceCoordinates[[1]],OpticalLength}.  Because the trace is two-dimensional, only the horizontal surface dimension is of interest, so take the first part of SurfaceCoordinates.

In[31]:=

Because opticalLength contained only the last ray-intersections, it was not necessary for us to give additional selection properties to ReadTurboRays.  In such cases, the third parameter entry of selection properties for ReadTurboRays is simply omitted and opticallength automatically holds only the last surface intersection data.

Using the built-in ListPlot function of Mathematica, plot the optical path length data.

In[32]:=

Optical path length as a function of the surface position. Note the discrete plot points because of the discrete data.

You are now ready to convert the optical length into wavelength units. However, in order to keep accurate machine precision numbers, first subtract off a constant optical length from the data. The resulting optical path difference still portrays the wavefront.

In[33]:=

Next plot the wavefront.

In[35]:=

The wavefront has units of wavelengths. The data is still discrete.

There are several interesting aspects to this wavefront result.  The wavefront has a range of 8000 wavelengths.  This is because the optical path length has a path difference of only 0.4 millimeters.  If the entire optical path length had been converted into wavelengths, then, with the minimum path length greater than 80 millimeters, this absolute wavefront would have taken on values greater than one million.

In[36]:=

Out[36]=

Later in this example, the wavefront values in a complex exponential function to calculate the complex field will be used.  If the absolute phase information had been converted, such results would have easily exceeded the limits of standard machine precision.  Nevertheless, together with Mathematica, Rayica presents some methods that permit such arbitrary precision calculations.  This is discussed further in the Help Browser.  For the moment, however, play it safe and work instead with optical path differences, since only the relative phase is important for this analysis.

Until now, the data has always taken the form of discrete points.  But in order to be really useful in phase-front and complex field calculations, a continuous model for our data is needed.  Mathematica provides the built-in Interpolation function for this very purpose.  First, supply Interpolation with the list of wavefront points and Interpolation returns a special continuous function that is matched to the discrete data.

In[117]:=

?Interpolation

In[37]:=

For two-dimensional data sets, Interpolation permits the use of non-uniform point spacings.  This feature is very important to get an accurate model of the wavefront, since perfectly spaced rays at the measurement surface cannot be guaranteed.  Now a continuous plot of the wavefront with OPDfunction can be made.

In[45]:=

A continuous wavefront plot with OPDfunction.

It is a simple matter to obtain the complex amplitude field, once you have the interpolated wavefront function.  Here a constant field strength is assumed.  Otherwise the ray Intensity values would be needed to construct a field strength function.  See the Help Browser for further information on such intensity-dependent calculations.

In[46]:=

Here we evaluate a single field point.

In[47]:=

OpticalField[5]

Out[47]=

You can also plot the complex field over a particular range.  Plot the real and imaginary values with two different colors.  When trying to plot out the full measurement surface, the plotted curves do not make any sense.  The plotted functions are oscillating too rapidly.

In[48]:=

The complex optical field plotted as real and imaginary values.  We cannot follow the plot at this scale.

Note the gaps present in the plot's top and bottom edges.  These gaps are merely artifacts in Mathematica's plotting process.  Mathematica's Plot routine is unable to keep up with the rapid oscillations of the plot, although you can try to increase the resolution with PlotPoints for better plot accuracy.  The actual data would appear as a solid rectangle of color at this scale.

Fortunately, because of the continuous nature of our interpolated phase function, you are able to choose a very tiny plot scale.  This time you can see all of the details in the plot.  In fact, you are able to extract complex field information with far better resolution than the original ray trace spacing!

In[49]:=

Note that the optical field values still oscillate quite rapidly over a small scale.  This time you can follow the plot. This plot range is only 0.3 millimeters across!

The complex field function is very important for all sorts of interference and diffraction calculations.  See the Help Browser for further discussion and examples.

Go to list of topics

Three-Dimensional Analysis with Zernike Polynomials

Now you examine the same optics using a three dimensional source, GridOfRays.  Again you are interested in obtaining a continuous optical path  difference (OPD) function for your system.  However, you shall soon discover  that accurate OPD modeling is more challenging in three dimensions.  For  readers becoming familiar with Rayica for the first time, this topic is more advanced and may be skipped  without penalty.

In[88]:=

The same beam expander/collimator optics from the previous example.  This  time the input source is a GridOfRays.

To obtain the optical path points, nest ReadTurboRays together with a TurboTrace call.  However, instead of calling TurboTrace directly, call TurboPlot to get a plot of the ray intersection points on the last surface.  TurboPlot works with TurboTrace in an analogous way that AnalyzeSystem works with PropagateSystem: initially, it runs the trace and, finally, it produces a graphical  display of the trace.  However, TurboPlot is restricted to the display of two-dimensional surface-spot diagrams.  You can also run ShowSystem with TurboTrace, but TurboPlot works much faster than the TurboTrace-ShowSystem combination.

In[75]:=

This TurboPlot spot diagram of the last surface shows the ray intercept points.  Note that the grid of ray-surface coordinates are not perfectly straight.  Such coordinate distortions become an important concern when modeling the continuous OPD function.

Next, find the optical path difference by subtracting the minimum path length and dividing by the wavelength.

In[76]:=

ListPlot3D can give you an approximate plot of the three-dimensional phase front.  Inside ListPlot3D, first partition the wavefront data into a matrix (since waveFront is presently a one-dimensional list of points).  Unfortunately, ListPlot3D must always assume a uniform grid spacing.  Therefore, you must discard the surface coordinate information from your data before plotting the OPD: use Transpose[..][[3]] to accomplish this.

In[79]:=

A three-dimensional plot of the discrete OPD.  The plot only approximates our phase front, since it does not take into account the ray-grid distortion at the surface.

You are ready to construct a continuous model for the optical path difference (OPD) function.  In the two-dimensional case, you had used the built-in Interpolation function for this purpose.  Although Interpolation does work with data sets of any dimension, for three-dimensional and higher-order data, Interpolation does not permit the specification of nonlinear grid lines.  Unfortunately, this missing feature is very important for you to get an accurate model of your wavefront, since the rays were not perfectly spaced on the measurement surface and the resulting grid lines are clearly curved.  Therefore you need to consider an approach that is different from Interpolation.  Fortunately, Mathematica provides you with other methods.

One such method involves constructing a least-squares fit of the data points to a linear combination of functions.  Because the optics are radially symmetric, you can use a linear combination of Zernike polynomials to model the wavefront.  The Zernike polynomials form a basis set of functions.  Once the basis coefficients are found by the least squares method, you can use them to construct a very accurate model of your OPD.

First, however, you must fabricate the Zernike polynomial functions in Mathematica.  By consulting Section 9.2.1 of the Principles of Optics book by Born and Wolf (1980), you can construct the required functions.  It will require several steps to arrive at the desired results.  Since the Zernike polynomials are so widely used in wavefront applications, it is well worth examining here in some detail.

According to Born and Wolf: "The circle polynomials of ZERNIKE are polynomials in two real variables X, Y, which, when expressed in polar coordinates (X = ρ Sin[θ], Y = ρ Cos[θ]), are of the form

=

where l ≤> 0 and n ≥ 0 are integers, n ≥ |l|, and n - |l| is even."

However, instead of using complex polynomials, V, it is suggested by the book to use the real polynomials, U, of the forms

and

Here, are called the radial polynomials (apparently related to Jacobi's polynomials and terminating hypergeometric series) and are given by

where m≥0 and n≥0 are integers, mn and n+m is even.

Next, equation (4) to define a Mathematica function R for the radial polynomials :

In[10]:=

The book also gives a table that shows the explicit form of the radial polynomial for the first few values of n and m. Here produce a similar Table in Mathematica:

In[30]:=

Out[30]//TableForm=

	n=0	n=1	n=2	n=3	n=4	n=5	n=6
m=0	1
m=1		ρ
m=2
m=3
m=4
m=5
m=6

The radial polynomials for m ≤ 6, n ≤ 6.

Now that you have developed the radial polynomial function, construct the real polynomials of Zernike given by equations (2) and (3) above.

In[149]:=

ZernikePolynomial only produces a single pair of polynomials for each n and m value.  For your purposes, however, you need to reconstruct the entire series of polynomials that run up to and include the given values of n and m.  For this, define the ZernikeList function.

In[150]:=

ZernikeList produces a flattened Table of ZernikePolynomial terms.  Select then removes any zeros present.  Here is a display of some terms.

In[151]:=

Out[151]=

Out[152]=

28 polar-coordinate functions are generated by ZernikeList for all terms running up to n = m = 6.

Finally, you are ready to construct a continuous OPD model.  Use the built-in Fit function of Mathematica to perform a least squares fit of the Zernike polynomials to your ray-trace wavefront data.

In[153]:=

?Fit

Because the Zernike polynomials are designed for use over the unit circle, you must renormalize your wavefront coordinates to occupy a unit circle. For this, you need the maximum radius of your wavefront coordinates.

In[154]:=

Out[154]=

Now use Fit to build your model.  Since the original data is given in rectangular coordinates, x and y, you must convert the wavefront coordinates for ZernikeList into normalized polar coordinates by the equations: and θ = ArcTan[x,y].  Lastly, Chop is employed to delete any near-zero coefficients returned by Fit.

In[155]:=

Out[156]=

The fitted Zernike polynomials of the OPD using terms up to n = m = 6.  Only six out of the original 28 terms remain in the final equation.  Chop has removed the 22 other near-zero terms.

In the previous result, Chop has served a very useful purpose of removing the near-zero terms. This prevents you from having to wade through vast quantities of polynomials when only a few of the terms are actually important to the wavefront model. This also allows you to start out with as many terms as you wish, since the near-zero terms are automatically deleted. You can also specify a second argument to Chop that manually sets its cutoff threshhold (it is normally ).

The amount of RMS residual error that was not fitted by the order of Zernike coefficients is given.

In[157]:=

Out[157]=

The residual error is considerably smaller than a single wavelength unit (i.e., one wavelength unit = 1). This indicates an excellent fit to the data. If you wanted even greater accuracy, you could simply add more Zernike polynomial terms to the fitting operation.

With an accurately fitted equation, you can at last plot a truly accurate picture of the OPD function.

In[163]:=

A truly accurate three-dimensional plot of the continuous OPD.  This time, the plotted function fully takes into account the ray trace grid distortion.

With a continuous wavefront function in hand, it is now a simple matter to obtain the complex amplitude field at the surface. Again in this example, for simplicity, assume a constant field strength.

In[159]:=

Here is a contour plot of the optical field phase near the surface center.

In[161]:=

A contour plot of the optical field phase.  The actual field function is smooth but zagged plot lines arise from finite plot points.

In imaging systems analysis, the optical field at the exit pupil of an imaging system (known as the pupil function) is very important, since both the point spread function (PSF) and modulation transfer function (MTF) are often derived from the pupil function.  Such calculations can be based on the methods presented here.  See the Help Browser for further discussion and examples.

Go to list of topics

Created by Mathematica  (November 19, 2004)