Wavica Users Guide

1. Introduction

2. Symbolic Calculations

2.1 Overview

Introduction
How to Enter Symbolic Parameters
How SymbolicTrace works
Light sources and SymbolicTrace
Options of SymbolicTrace
Working with SymbolicTrace
Multiple Light Paths
Infinite Precision Calculations

2.2 Symbolic Optimization

Infinite-Conjugate Spherical Lens Design
Finite-Conjugate Spherical Lens Design
Doublet Lens Design

2.3 Geometric Intensity Analysis

Set-up System
Preliminary Characterization
Inverse Position
Geometric Intensity Analysis

3. ABCD Matrix Calculations

3.1 Overview
3.2 SphericalMirror Examples

3.3 PlanoConvexLens Examples

Simple lens
Off-axis lens
Cylindrical lens

3.4 SphericalLens Examples
3.5 ThinLens, and Custom ABCD Input
3.6 Grating Examples

4. Gaussian Beam Calculations

4.1 Overview

4.2 Reflected Lens Example

4.3 Reflected ThickLens Example

4.4 Off-Axis Beam Example

4.5 CylindricalLens Example

4.6 Prism-Lens Example

5. Imaging Calculations

5.1 Overview
5.2 FindLensParameters
5.3 PupilFunction
5.4 Zernike Polynomials and Seidel Aberrations
5.5 PointSpreadFunction
5.6 OpticalTransferFunction
5.7 GeometricPointSpreadFunction

6. Interference and Wavefront Calculations

6.1 Overview

6.2 One-Dimensional Wavefronts

Example 1
Example 2
Example 3
Example 4
Example 5

6.3 Two-Dimensional Wavefronts

Example 1
Example 2
Example 3

6.4 Wavefront Propagation with Gaussian Wavelets

Introduction to Wavelet Analysis
Define LineOfWavelets
Define FindField
Collimated Beam of Wavelets
One Wavelet with Lens
Two Wavelets with Lens
Three Wavelets with Lens
Wavelet analysis of an Imaging System

1. Introduction

Welcome to Wavica, the first optical design package to utilize analytic solutions for optical systems. With its revolutionary approach to optical design, Wavica taps into the symbolic powers of Mathematica. Once loaded, Wavica fully integrates its capablities into the Rayica design system.

Copyright © 2000-2004 by Optica Software, Urbana Il. All Rights Reserved.

Loading Wavica

The Wavica package (the Wavica directory) should be located in the same directory as the Rayica package (the Rayica and RayicaTools  directories). Wavica is loaded with the following command. Wavica automatically loads the Rayica package before loading itself.

In[1]:=

Needs["Wavica`Wavica`"]

Go to list of topics

2. Symbolic Calculations: SymbolicTrace

2.1 Overview

Introduction

The principle function responsible for symbolic calculations is SymbolicTrace. SymbolicTrace creates the symbolic solution to an optical system. In its output, SymbolicTrace returns a host of symbolic parameters in the form of rules.  Each parameter has rule name identified with it, and the most important ones are: SymbolicRayEnd, SymbolicRayTilt, SymbolicIntensity, SymbolicSurfaceCoordinates, and SymbolicOpticalLength. Nearly all important functions in the Wavica package depend on SymbolicTrace in some fashion. Functions like GaussianTrace and FindABCDMatrix depend entirely on first order SymbolicTrace results. Therefore, it is best to first study SymbolicTrace in order to acquire a good understanding of these two functions. FindLensParameters uses SymbolicTrace for pupil calculations and finding paraxial information such as the lens principle points, magnification, and the lens focal length. Other functions such as PointSpreadFunction (PSF) and OpticalTransferFunction (MTF) use SymbolicTrace indirectly through FindLensParameters. However, it is not necessary to study SymbolicTrace before learning about these latter functions.

In[2]:=

?SymbolicTrace

How to Enter Symbolic Parameters

The basic Rayica system permits symbolic entry of most system parameters. These parameters may then be used to construct symbolic variables within the system. 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[3]:=

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

Out[3]=

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[4]:=

symbolicSource = LineOfRays[40, SymbolicWaveLength->λ, SymbolicIntensity->i, NumberOfRays->4]

Out[4]=

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[5]:=

LineOfRays[{s,10}]

Out[5]=

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

In[6]:=

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

Out[6]=

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[7]:=

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

Out[7]=

Once a system has been defined with symbolic parameters, you can use Wavica's SymbolicTrace function to calculate the symbolic solution to the optical system.

In[8]:=

TableForm[SymbolicTrace[symbolicLensSystem,{{f>0},{y,0,2}}, FinalFormat->False]]

Out[8]//TableForm=

SymbolicSurfaceBoundary→{50,50}

SymbolicOffAxis→{0,0}

SymbolicTranslationVector→{x,0,0}

SymbolicSurfaceNormalFunction→{1,0,0}

At the same time, you can use AnalyzeSystem to conduct a numeric ray-trace of the system.

In[9]:=

AnalyzeSystem[symbolicLensSystem,PlotType->TopView]

Out[9]=

When you include a Gaussian-beam light source in a system, you can use FindABCDMatrix to calculate the symbolic ABCD matrix of the system.

In[10]:=

symbolicGaussianSystem =
    {symbolicGaussianSource,
    Move[symbolicLens, {d,50}],
    Move[Screen[50],{{x,200},{y,0}}]}

Out[10]=

In[11]:=

FindABCDMatrix[symbolicGaussianSystem,{f>0}, MatrixForm->True]

Out[11]//MatrixForm=

If you use the ShowGaussian->True option with AnalyzeSystem, you can see the Gaussian wave-front behavior superimposed on the normal numeric ray-trace.

In[12]:=

AnalyzeSystem[symbolicGaussianSystem,PlotType->TopView,ShowGaussian->True]

Out[12]=

Go to list of topics

How SymbolicTrace works

SymbolicTrace initially works by tracing a single chief ray through the optical system and constructing symbolic expressions that represent light propagation through the optical surfaces in the system. First, SymbolicTrace uses Rayica to make a numerical trace of the chief ray through the optical system. From this chief ray, the surface sequence is determined as well as positions and directions of the chief ray at each surface. In addition to the numerical results, Rayica also provides any user-specified symbolic parameter information about the system to SymbolicTrace. SymbolicTrace then uses the numerical values of the initial ray trace as boundary conditions to determine the symbolic solution required in symbolic calculations. In particular, many surface functions have more than one intercept solution. In the case of a spherical surface, for example, there are two intercept solutions for the same spherical function, but SymbolicTrace must choose to work with only one of the two possible intercept solutions. In such an event, SymbolicTrace uses the solution determined from the numerical trace (where only numbers are present and no variables are present) to dictate its choice of symbolic ray-surface intercept solution.

SymbolicTrace performs its calculations in two stages. In the first stage, SymbolicTrace determines the symbolic solution at each optical surface. In this stage, the surface functions are calculated in the local coordinate system of the particular surface. After each local surface result is calculated, it is then transformed in the global coordinate system before progressing to the second stage of operation. Finally in the second stage of operation, SymbolicTrace constructs a single, nested solution of the optical system by cascading together the symbolic results of each individual surface.

In addition to constructing a symbolic solution, SymbolicTrace can use symbolic Taylor's series expansions during the course of its calculations in order to reduce the size and complexity of the resulting solution. For example, SymbolicTrace can give the paraxial solution to system by taking the first order Taylors series expansion of the result. Such a solution is useful for determining the abcd matrix representation of the system and for calculating Gaussian beam propagation through the system. However, in other cases, SymbolicTrace can use a higher order Taylor series expansion in order to maintain a higher degree of accuracy in the final results. In some cases, SymbolicTrace may return a fully unexpanded symbolic solution that completely models the system with machine-precision, or even infinite precision. In cases where Taylor's series expansions are used, the numerical surface intercepts from the ray-trace solution provide the fixed points for the symbolic Taylor's series calculations.

Go to list of topics

Light Sources and SymbolicTrace

In order to function, SymbolicTrace requires the presence of at least one light source as well as one or more optical surfaces. However, SymbolicTrace does not use light sources in the same way as PropagateSystem. In particular, only the chief ray of the light source is utilized and any other rays from the light source are ignored by SymbolicTrace. To illustrate this point, consider the following three examples.

In[13]:=

raySystem = {
    Move[SingleRay[],{-50,{y,0}}],
    PlanoConvexLens[{f,100},50,10],
    Move[Screen[50],100]};
AnalyzeSystem[raySystem,PlotType->TopView];
SymbolicOpticalLength/.SymbolicTrace[raySystem,{{y,0,3},{f>0}}]

Out[15]=

In[16]:=

linearSystem = {
    Move[LineOfRays[20],{-50,{y,0}}],
    PlanoConvexLens[{f,100},50,10],
    Move[Screen[50],100]};
AnalyzeSystem[linearSystem,PlotType->TopView];
SymbolicOpticalLength/.SymbolicTrace[linearSystem,{{y,0,3},{f>0}}]

Out[18]=

In[19]:=

wedgeSystem = {
    Move[WedgeOfRays[20],{-50,{y,0}}],
    PlanoConvexLens[{f,100},50,10],
    Move[Screen[50],100]};
AnalyzeSystem[wedgeSystem,PlotType->TopView];
SymbolicOpticalLength/.SymbolicTrace[wedgeSystem,{{y,0,3},{f>0}}]

Out[21]=

In all three cases, even though the three different light sources have different characteristics, the resulting symbolic solution is exactly the same. This is because the chief ray from each source follows the same trajectory through the optics. Therefore, when SymbolicTrace is involved, the type of light source is not as important as the starting position and direction of its chief ray.

Go to list of topics

Options of SymbolicTrace

SymbolicTrace has many special option settings, many of which are seldom changed by the user. However, one of the best ways to understand SymbolicTrace is to study its options. These many options are present in order to provide complete user control over the symbolic solution creation process. See the Help Browser listing for SymbolicTrace to learn more about its options. In addition to SymbolicTrace, many of these options are shared by other functions in Wavica including: FindABCDMatrix, GaussianTrace, and GaussianPlot.

In[22]:=

Options[SymbolicTrace]

Out[22]=

The SeriesOrder option is one of the most important ways of controlling the form of the resulting symbolic solution. SymbolicTrace uses SeriesOrder->1 as its default option setting. SeriesOrder->1 returns the first order Taylor's series symbolic solution to the optical system. SeriesOrder->Infinity returns the complete, unapproximated symbolic solution to the system. Sometimes, however, SeriesOrder->Infinity may result in a solution that is too large for your computer system's memory, possibly consuming hundreds of megabytes of space. Therefore, SeriesOrder must be used with some care and respect.

In[23]:=

?SeriesOrder

In[24]:=

opticalsystem = {
Move[LineOfRays[20], {-50,{y,0}}],
PlanoConvexLens[100,50,10,DesignWaveLength->.532],
Move[Screen[50],{d,100}]};
AnalyzeSystem[opticalsystem,PlotType->TopView];
SymbolicTrace[opticalsystem] (*SeriesOrder->1*)
ByteCount[%]

Out[26]=

Out[27]=

In[28]:=

SymbolicTrace[opticalsystem,SeriesOrder->2]
ByteCount[%]

Out[28]=

Out[29]=

In[30]:=

SymbolicTrace[opticalsystem, SeriesOrder->Infinity, MakeFloatingPoint->True]

In[31]:=

ByteCount[%]

Out[31]=

Go to list of topics

Working with SymbolicTrace

In addition to using the SeriesOrder option, the power expansion of a particular user-specified variable may be directly given. In such a case, the resulting solution order corresponds directly with the specified order. This is because the expansion takes place after the surface solutions are combined. If desired, the Taylor series expansion process can be closely monitered with the RunningCommentary->All option setting.

In[32]:=

SymbolicTrace[opticalsystem, {{y,0,2}}, RunningCommentary->All]

Out[32]=

In[33]:=

ByteCount[%]

Out[33]=

SymbolicTrace can take, within its second argument, a list of assumptions that get passed internally to the Simplify operations. These assumptions are used to further simplify the resulting expression. Here is the SymbolicTrace operation without any assumptions given.

In[34]:=

opticalsystem = {
Move[LineOfRays[20], {-50,{y,0}}],
PlanoConvexLens[{f,100},50,10,DesignWaveLength->.532],
Move[Screen[50],{d,100}]};
SymbolicTrace[opticalsystem]

Out[35]=

Here is the same operation with some assumptions included about the user variables. The RunningCommentary->Simplify setting is used here to moniter the internal simplification steps.

In[36]:=

SymbolicTrace[opticalsystem, {f>0, d∈Reals}, RunningCommentary->Simplify]

Out[36]=

In some cases, the second argument of SymbolicTrace may include both assumptions as well as specifications for Taylor series expansion of symbolic system variables.

In[37]:=

SymbolicTrace[opticalsystem, {{y,0,2},{f,100,1},{f>0, d∈Reals}}]

Out[37]=

Go to list of topics

Multiple Light Paths

In addition to using the SymbolicTrace with linearly ordered systems, it is also possible to work with systems involving multiple light paths. Each light path produces a different solution in the output. Here is an example that uses two light sources. In general, every distinct light source produces a distinct light path. In this event, the two light sources produce two distinct symbolic solutions.

In[38]:=

twoSourceSystem = {
    Move[LineOfRays[20],{-50,{y,0}}],
    Move[LineOfRays[20],{-50,{y,0}},5],
    PlanoConvexLens[100,50,10],
    Move[Screen[50],{d,100}]};
AnalyzeSystem[twoSourceSystem,PlotType->TopView];
SymbolicTrace[twoSourceSystem]

Out[40]=

In some cases, the optical arrangement may produce multiple light paths from a single light source. Here is an example that uses a beam splitter to create two light paths.

In[41]:=

beamSplitterSystem = {
    Move[LineOfRays[20],{-50,{y,0}}],
    PlanoConvexLens[100,50,10],
    Move[BeamSplitter[{50,50},50],50,-45],
    Move[Screen[50],{d1,100}],
    Move[Screen[50],{50,{d2,50}},90]
    };
AnalyzeSystem[beamSplitterSystem,PlotType->TopView];
SymbolicTrace[beamSplitterSystem]

Out[43]=

Go to list of topics

Infinite Precision Calculations

SymbolicTrace is also capable of producing solutions with infinite precision. This happens automatically when the input parameters are all pure integers, integer fractions, or symbols. Here is an example of a reflective spherical mirror calculation.

In[105]:=

opticalSystem = {
    Move[LineOfRays[20],{0,{y,0}}],
    Move[SphericalMirror[{r,-100},50,10],25],
    Move[Screen[50],{d,-25}]};
AnalyzeSystem[opticalSystem,PlotType->TopView];
SymbolicTrace[opticalSystem, {{y,0,5},{r<0}},SymbolicRefractiveModels->{Air->1}]

Out[107]=

When refraction is involved, the index model needs to be substituted for a pure fraction of integer numbers within SymbolicTrace in order to maintain infinite precision. This example uses 3/2 for the BK7 refractive index and 1 as the value of Air. For this, SymbolicTrace takes the SymbolicRefractiveModels as an option.

In[47]:=

?SymbolicRefractiveModels

In[99]:=

opticalSystem = {
    Move[LineOfRays[20],{-50,{y,0}}],
    SphericalLens[{r,50},Infinity,50,10],
    Move[Screen[50],{d,100}]};
AnalyzeSystem[opticalSystem,PlotType->TopView];
SymbolicTrace[opticalSystem, {{y,0,5},{r>0}},SymbolicRefractiveModels->{BK7->3/2,Air->1}]

Out[101]=

In this example, the original component medium, BK7, is replaced with 3/2, and surrounding intrinsic medium, Air, is replaced with 1. We can also specify the level of precision with the Precision option of SymbolicTrace.

Go to list of topics

2.2 Symbolic Optimization

In this section, we will consider how to use Wavica for symbolic optimization purposes. In particular, we will examine several different examples that illustrate some ways that symbolic optimization can be carried out with Wavica. Unlike numeric optimization procedures, symbolic raises the possibility of the finding the global minimum automatically.

Infinite-Conjugate Spherical Lens Design

In this example, we find the optimal parameters for spherical lens system that focusses a plane wave onto a point at a surface. We first use SymbolicTrace to calculate a symbolic solution of a symbolic ray transversing through the system for different lens curvatures (represented by symbolic parameters c1 and c2) and different transverse y positions of ray entry into the lens. First we define the optical system. In order to make the solution more clean, we assign integer values to the refractive index parameters of the system.

In[51]:=

Clear[y,c1,c2];

In[52]:=

opticalsystem = {
    Move[LineOfRays[10, IntrinsicMedium->1],{-20,{y,0}}],
    SphericalLens[{1/c1,35},{1/c2,-55},30,10,ComponentMedium->3/2],
    Move[Screen[30],50]};

Here is a numeric ray-trace of the system. Note that we have the trace is not optimized but instead uses the numeric assignments made to the optical system during its initial definition.

In[53]:=

plot = TurboPlot[opticalsystem,PlotType->TopView, ShowArrows->True, EmbedRays->False];

Next, we use SymbolicTrace function to calculate the third order symbolic solution to the optical system. Although SymbolicTrace calculates several types of symbolic results, we will only use the SymbolicSurfaceCoordinates for optimization purposes. Here, MakeFloatingPoint retains the pure integer values of the numeric system parameters.

In[54]:=

Py = (SymbolicSurfaceCoordinates/.
    SymbolicTrace[opticalsystem,{{y,0,3}, {c1>0,c2<0,y∈Reals}},
        MakeFloatingPoint->False])[[1]]

Out[54]=

Now we are ready to calculate a symbolic figure of merit for the system performance. In particular, we want to base the performance of the system on every y value of the input ray across a finite aperture runs between -5 and 5. This is accomplished by integrating the square of the focal intercept solution given by Py, above.

In[55]:=

meritFunction = Expand[Integrate[Expand[Py^2],{y,-5,5}]]

Out[55]=

Here is a plot of the resulting merit function for diffferent lens curvatures.

In[56]:=

Plot3D[meritFunction,{c1,0,.1},{c2,-.1,0}];

Using minimization techniques in Calculus, we can find the symbolic global optimization for this lens configuration:

In[57]:=

Dc1 = D[meritFunction,c1]

Out[57]=

In[58]:=

Dc2 = D[meritFunction,c2]

Out[58]=

Next we use NSolve to find the roots of the these two equations. These roots contain all possible solutions to third order system model.

In[59]:=

extrema = NSolve[{Dc1==0,Dc2==0},{c1,c2}]

Out[59]=

Note that many of the roots contain imaginary components. In addition, some of the solutions are local minima and do not represent the global solution to the system. However, by shifting through the different roots and comparing their results we can determine the global solution to the system.

In[60]:=

finalextrema = Select[Union[Chop[Transpose[{meritFunction/.extrema,extrema}]]],
    (#[[1]]==Re[#[[1]]]&&#[[1]]<1&&#[[1]]>0)&]

Out[60]=

Finally, we can verify the resulting solution by numerically tracing the resulting system and calculating the resulting focal point and spot size with FindFocusFast.

In[61]:=

FindFocusFast[TurboPlot[
    Move[LineOfRays[10,NumberOfRays->201,IntrinsicMedium->1],{-20,0}],
    plot,SymbolicValues->finalextrema[[1,2]], PlotType->TopView]]

Out[61]=

It is interesting to note that the calculated focal point position ( 49.9903) is slightly different from the final surface position (50). This is partly because the symbolic model uses an infinite number of rays in its integration calculation while TurboPlot is using only a finite number of rays. In addition, however, the symbolic model is a third order solution and, as such, does not perfectly model the actual system. Therefore, we will next try to improve the accuracy of our symbolic model by using a fifth order symbolic solution.

In[62]:=

Py = (SymbolicSurfaceCoordinates/.
    SymbolicTrace[opticalsystem,{{y,0,5}, {c1>0,c2<0,y∈Reals}}, MakeFloatingPoint->False])[[1]]

Out[62]=

Once again, we integrate across the system aperture to calculate a symbolic merit function for the system's performance.

In[63]:=

meritFunction = Expand[Integrate[Expand[Py^2],{y,-5,5}]]

Out[63]=

This time, since we already have the third-order global minimum for lens system, we will use its solution as a starting point for a search with FindMinimum on the fifth-order model.

In[64]:=

localminimum = FindMinimum[meritFunction,{c1,Evaluate[c1/.finalextrema[[1,2]]]},{c2,Evaluate[c2/.finalextrema[[1,2]]]}, AccuracyGoal->12]

Out[64]=

Since the third-order solution was already determined, the fifth-order result with FindMinimum is calculated very quickly. And we can, once again, verify the resulting solution by numerically tracing the resulting system and calculating the resulting focal point and spot size with FindFocusFast.

In[65]:=

FindFocusFast[TurboPlot[
    Move[LineOfRays[10,NumberOfRays->201,IntrinsicMedium->1],{-20,0}],
    plot,SymbolicValues->localminimum[[2]], PlotType->TopView]]

Out[65]=

When we trace the system with 101 rays, this time, the calculated focal point is much better than the third order model. However the result still does not completely match the focal surface position of the system model. However, when we increase the number of traced rays to 1001, the resulting focal position matches the final surface position precisely.

In[66]:=

FindFocusFast[TurboPlot[
    Move[LineOfRays[10,NumberOfRays->1001,IntrinsicMedium->1],{-20,0}],
    plot,SymbolicValues->localminimum[[2]], PlotType->TopView]]

Out[66]=

This final calculation indicates that the fifth-order symbolic solution is a close fit to the actual system.

Go to list of topics

Finite-Conjugate Spherical Lens Design

In this section, we will consider the application of symbolic optimization to a finite-conjugate lens system. Finite-conjugate lens systems are considerably more difficult to accurately represent symbolically than a similar finite-conjugate system. In particular, the use of symbolic angular parameters requires a greater number of expansion terms in order order to obtain the same modelling accuracy. This can significantly increase the size of the symbolic solution as well as the amount of time required for the symbolic calculation. Let us now consider such an optical system that images a point source of light onto a surface. This is defined for the opticalsystem shown below.

In[67]:=

opticalsystem = {
    Move[WedgeOfRays[10,NumberOfRays->101, IntrinsicMedium->1],-90,{θ/Degree,0}],
    SphericalLens[{1/c1,35},{1/c2,-35},30,10, ComponentMedium->3/2],
    Move[Screen[30],60]};

In[68]:=

plot = TurboPlot[opticalsystem,PlotType->TopView, EmbedRays->False];

Next we use SymbolicTrace to determine the fifth order solution to the symbolic optical path length of the system as function of lens curvatures, c1 and c2, and input ray angle, θ.

In[69]:=

OL = SymbolicOpticalLength/.SymbolicTrace[opticalsystem, {{θ,0,5}, {c1>0,c2<0,θ∈Reals}}, ReportedParameters->SymbolicOpticalLength, MakeFloatingPoint->False]

Out[69]=

Here is a plot of the symbolic optical path length as a function of ray angle. This uses the default initial settings of lens curvature.

In[70]:=

oplot = Plot[OL/.{c1->1/35,c2->-1/35},{θ,-3. Degree, 3. Degree}];

We can now compare this symbolic result with a numeric trace of the same system to see how close the fifth-order symbolic model matches the actual system. To this end, we need to extract the numeric transfer function of the optical system via a numeric ray-trace of the default lens configuration. We use ReadTurboRays to extract the numeric data from the trace.

In[71]:=

data = ReadTurboRays[plot, {RayTilt,OpticalLength,IntersectionNumber,RaySourceNumber}];

In[72]:=

inputangles = Map[{ArcTan[#[[1,1]],#[[1,2]]],#[[4]]}&,Select[data,(#[[3]]==1)&]]

Out[72]=

Next, we use ListPlot to plot the optical transfer function of the system.

In[73]:=

dataplot = ListPlot[Map[{#[[1,1]],#[[2,2]]}&,Transpose[{inputangles,Select[data,(#[[3]]==3)&]}]]];

Finally, we can compare a plot of numeric ray-trace points with the result given by the symbolic equation.

In[74]:=

Show[dataplot,oplot];

From this, we see that the symbolic result follows the numeric behavior quite well, although the symbolic result begins to deviate from the numeric result after about .05 radians. For this reason, we will limit the absolute angular values of our symbolic model to less than .05 radians. Based on this knowledge, we can now construct a merit function that integrates between the -.05 and .05 radians of angular space.

In[75]:=

meritFunction=Expand[Integrate[Expand[D[OL,θ]^2],{θ,-5/100,5/100}]]

Out[75]=

In[76]:=

Dc1 = D[meritFunction,c1]

Out[76]=

In[77]:=

Dc2 = D[meritFunction,c2]

Out[77]=

In[78]:=

Show[{Plot3D[Dc1,{c1,0,.05},{c2,-.05,0},ColorFunction->(Red&), PlotPoints->50],
    Plot3D[Dc2,{c1,0,.05},{c2,-.05,0},ColorFunction->(Green&), PlotPoints->50],
    Plot3D[0,{c1,0,.05},{c2,-.05,0},ColorFunction->(Blue&), PlotPoints->50]},
    ViewPoint->{-1.309, -1.983, 2.409}];

In[79]:=

extrema = NSolve[{Dc1==0,Dc2==0},{c1,c2}]

Out[79]=

In[80]:=

Length[extrema]

Out[80]=

In[81]:=

D2c1 = D[meritFunction,{c1,2}];
D2c2 = D[meritFunction,{c2,2}];
Dc1c2 = D[meritFunction,c1,c2];

In[84]:=

finalextrema =
    Union[Select[extrema,
        ((    c1==Re[c1]&&c2==Re[c2]&&
            Abs[c1]<.1&&Abs[c2]<.1&&
            (D2c1*D2c2-Dc1c2^2)>0&&D2c1>0
            )/.#)&]]

Out[84]=

In[85]:=

Map[{#,FindFocusFast[TurboPlot[
    plot,SymbolicValues->#, PlotType->TopView]]}&,
    finalextrema]

Out[85]=

In[86]:=

basicsystem = {
    Move[WedgeOfRays[10,NumberOfRays->101, IntrinsicMedium->1],-90],
    SphericalLens[{1/c1,35},{1/c2,-35},30,10, ComponentMedium->3/2],
    Move[Screen[30],60]};

In[87]:=

sol = OptimizeSystem[basicsystem,SymbolicValues->localminimum[[2]],SequentialTrace->True]

Out[87]=

In[88]:=

FindFocusFast[TurboPlot[plot, sol, PlotType->TopView]]

Out[88]=

Go to list of topics

Doublet Lens Design

Now we define a 3 surfaced lens system, using r1 = 100 mm, r2 = -300 mm, and r3 = -100 as initial curvature starting points for the lens. Although important parameters to optimize, in this example we will fix the spacings between the different lens surfaces to be s1 = 15 mm and s2 = 5 mm.

In[89]:=

In[90]:=

opticalsystem= {Move[LineOfRays[45,NumberOfRays->6],{0,{y,0}}],optics};

Here is a plot of this initial system:

In[91]:=

In[92]:=

ModelRefractiveIndex[BK7][WaveLength->.532]

Out[92]=

In[93]:=

ModelRefractiveIndex[SF11][WaveLength->.532]

Out[93]=

In[94]:=

?Rationalize

In[95]:=

Rationalize[ModelRefractiveIndex[SF11][WaveLength->.532],.02]

Out[95]=

In[96]:=

Rationalize[ModelRefractiveIndex[BK7][WaveLength->.532],.02]

Out[96]=

In[97]:=

Py = (SymbolicRayEnd/.
    SymbolicTrace[opticalsystem,{{y,0,3},{c1>0,c2<0,c3<0,y∈Reals}}, NormalizeSurfaceNormal->False, ReportedParameters->SymbolicRayEnd, RunningCommentary->All, Simplify->True, FinalFormat->False,
    SymbolicRefractiveModels->{BK7->3/2, SF11->9/5}])[[2]];

In[98]:=

ByteCount[Py]

Out[98]=

In[99]:=

Py//InputForm

Out[99]//InputForm=

y*(1. - 0.8333333333333331*c2 + c1*(-6.385445898048419 + 4.164421237857663*c2) + c2*(-23.99353490705878 - 53.30100786862726*c3) + 63.961209442352725*c3 +
  c1*(-39.96767453529393 - 408.42084246788704*c3 + c2*(119.90302360588177 + 266.3614190007958*c3)) +
  (-79.9515118029409*c2*(-0.0017089009379546 + c3)*c3*(0.48764285560680765 + c3) + 66.62625983578405*c2^2*(0.007971515150067884 + c3)*(0.20283168300952217 + c3)*(0.7735647111781162 + c3) +
    c3^2*(-0.39975755901470456 + 31.980604721176363*c3) - 18.5072943988289*c2^3*(0.5169292315766226 + c3)*(1.2975315262148452 + 0.7658726324299371*c3 + c3^2) +
    c1*(-612.6312637018306*(-0.0018540657824575426 + c3)*c3*(0.09770927331916474 + c3) + 1420.5942346709107*c2*(-0.0008775492072184491 + c3)*(0.042660433177280084 + c3)*(0.5220313761158 + c3) -
      1091.3419250727047*c2^2*(-0.028920975829296026 + c3)*(0.24057370655264448 + c3)*(0.8093056267501366 + c3) + 277.4598114591622*c2^3*(0.5169292315766226 + c3)*
       (1.297531526214845 + 0.7658726324299366*c3 + c3^2)) + c1^2*(3911.923789821072*(0.009138592301390268 + c3)*(0.02786745725476436 + c3)*(0.1672043655172598 + c3) +
      5915.952801241696*c2^2*(-0.05172453852830349 + c3)*(0.2630887177228902 + c3)*(0.840467368893461 + c3) - 8362.445782588522*c2*(0.5558346769858834 + c3)*
       (0.0021322577247055046 + 0.0794664514099874*c3 + c3^2) - 1386.5514377910222*c2^3*(0.5169292315766221 + c3)*(1.2975315262148461 + 0.7658726324299368*c3 + c3^2)) +
    c1^3*(-10624.499029050985*c2^2*(-0.0305711690535631 + c3)*(0.25638873105418436 + c3)*(0.852412801250604 + c3) - 8326.459239063663*(0.19890623747890349 + c3)*
       (0.01399394675838518 + 0.11160464690486825*c3 + c3^2) + 16290.898511211512*c2*(0.550657007706981 + c3)*(0.010673489313858552 + 0.14883676342155*c3 + c3^2) +
      2309.6737019676048*c2^3*(0.5169292315766224 + c3)*(1.2975315262148457 + 0.7658726324299371*c3 + c3^2)))*y^2)

In[100]:=

Timing[meritFunction=ExpandAll[Integrate[Py^2,{y,-5,5}]];][[1]]

Out[100]=

In[101]:=

ByteCount[meritFunction]

Out[101]=

In[102]:=

meritFunction//InputForm

Out[102]//InputForm=

83.33333333333333 - 7725.520072223724*c1 + 179467.42328278074*c1^2 - 77244.79245903643*c1^3 + 2.686384088965151*^6*c1^4 - 172351.7135478436*c1^5 + 1.1990021059754131*^7*c1^6 -
4137.811373398686*c2 + 212408.9693749286*c1*c2 - 990384.8241276548*c1^2*c2 + 2.877853121571819*^6*c1^3*c2 - 1.832924021410707*^7*c1^4*c2 + 1.0966652503587645*^7*c1^5*c2 -
9.906804433304237*^7*c1^6*c2 + 51572.782218815875*c2^2 - 505938.4792684064*c1*c2^2 + 1.0086273901604646*^6*c1^2*c2^2 - 1.0268587890209488*^6*c1^3*c2^2 + 1.6684772775022842*^7*c1^4*c2^2 +
2.817089520420558*^7*c1^5*c2^2 + 1.3119207727802566*^8*c1^6*c2^2 - 36205.8495172385*c2^3 + 1.548086693855232*^6*c1*c2^3 - 1.7922109105472267*^7*c1^2*c2^3 + 9.820077896997187*^7*c1^3*c2^3 -
3.008567653501071*^8*c1^4*c2^3 + 6.531375734981651*^8*c1^5*c2^3 - 1.2994516382756152*^9*c1^6*c2^3 + 770621.7432028732*c2^4 - 1.5362783655369086*^7*c1*c2^4 + 1.2129844142151959*^8*c1^2*c2^4 -
5.371592311617185*^8*c1^3*c2^4 + 1.8072237355127683*^9*c1^4*c2^4 - 4.875158776606653*^9*c1^5*c2^4 + 6.73440111955497*^9*c1^6*c2^4 - 46180.930506681754*c2^5 - 2.7131400739199156*^6*c1*c2^5 +
8.509113756155948*^7*c1^2*c2^5 - 8.508519403212641*^8*c1^3*c2^5 + 3.823942437314061*^9*c1^4*c2^5 - 7.626641205205066*^9*c1^5*c2^5 + 4.909316428643413*^9*c1^6*c2^5 + 3.4395836307242787*^6*c2^6 -
1.0313190087246548*^8*c1*c2^6 + 1.28845403488178*^9*c1^2*c2^6 - 8.585064554984132*^9*c1^3*c2^6 + 3.217664264503948*^10*c1^4*c2^6 - 6.431860511710032*^10*c1^5*c2^6 +
5.356995302780664*^10*c1^6*c2^6 + 10660.201573725453*c3 - 561926.2879904569*c1*c3 + 3.232042333179159*^6*c1^2*c3 - 7.54908270463919*^6*c1^3*c3 + 5.8659256157446004*^7*c1^4*c3 -
2.8316283716550138*^7*c1^5*c3 + 3.118053702644916*^8*c1^6*c3 - 273376.3555487509*c2*c3 + 3.5270072652928173*^6*c1*c2*c3 - 1.8163712023192946*^7*c1^2*c2*c3 + 1.0387401599868505*^8*c1^3*c2*c3 -
4.267297179618991*^8*c1^4*c2*c3 + 6.529090244318188*^8*c1^5*c2*c3 - 2.8495182332229667*^9*c1^6*c2*c3 + 257170.6274145298*c2^2*c3 - 4.921858767786396*^6*c1*c2^2*c3 +
4.669001802120763*^7*c1^2*c2^2*c3 - 2.901936267208843*^8*c1^3*c2^2*c3 + 1.022685816112428*^9*c1^4*c2^2*c3 - 2.3512932946217184*^9*c1^5*c2^2*c3 + 7.525267614801429*^9*c1^6*c2^2*c3 -
2.755184192395447*^6*c2^3*c3 + 6.106903637952671*^7*c1*c2^3*c3 - 5.1286435120359313*^8*c1^2*c2^3*c3 + 2.1863328565487666*^9*c1^3*c2^3*c3 - 6.154170340105774*^9*c1^4*c2^3*c3 +
1.5944062341812832*^10*c1^5*c2^3*c3 - 2.8492292033634136*^10*c1^6*c2^3*c3 + 3.603258448096*^6*c2^4*c3 - 8.615513250476956*^7*c1*c2^4*c3 + 9.426746241244583*^8*c1^2*c2^4*c3 -
6.633128882906318*^9*c1^3*c2^4*c3 + 3.2001849285345905*^10*c1^4*c2^4*c3 - 9.219747786670241*^10*c1^5*c2^4*c3 + 1.1486681236666006*^11*c1^6*c2^4*c3 - 6.197219187493397*^6*c2^5*c3 +
1.8369785950465462*^8*c1*c2^5*c3 - 2.2685185092057443*^9*c1^2*c2^5*c3 + 1.5206329953466343*^10*c1^3*c2^5*c3 - 5.9338406067303085*^10*c1^4*c2^5*c3 + 1.2931062340317773*^11*c1^5*c2^5*c3 -
1.2328440240524448*^11*c1^6*c2^5*c3 + 1.7368207534800593*^7*c2^6*c3 - 5.2076543271439034*^8*c1*c2^6*c3 + 6.50605988381383*^9*c1^2*c2^6*c3 - 4.335035832788493*^10*c1^3*c2^6*c3 +
1.624762376004455*^11*c1^4*c2^6*c3 - 3.2477735736504156*^11*c1^5*c2^6*c3 + 2.7050194491725345*^11*c1^6*c2^6*c3 + 339920.29887983925*c3^2 - 4.4365865591270365*^6*c1*c3^2 +
2.6616648285585493*^7*c1^2*c3^2 - 1.7285563214606816*^8*c1^3*c3^2 + 6.866960902507284*^8*c1^4*c3^2 - 1.1848255367861927*^9*c1^5*c3^2 + 4.70224481602699*^9*c1^6*c3^2 - 629861.8127289944*c2*c3^2 +
2.1357319956643984*^7*c1*c2*c3^2 - 2.8862047090173084*^8*c1^2*c2*c3^2 + 1.8474059925133746*^9*c1^3*c2*c3^2 - 6.177751732729872*^9*c1^4*c2*c3^2 + 1.7047245294845005*^10*c1^5*c2*c3^2 -
4.314392694858968*^10*c1^6*c2*c3^2 + 4.55635698327138*^6*c2^2*c3^2 - 1.1860743021115477*^8*c1*c2^2*c3^2 + 1.1756655133946195*^9*c1^2*c2^2*c3^2 - 5.947476723052926*^9*c1^3*c2^2*c3^2 +
2.0446201743547565*^10*c1^4*c2^2*c3^2 - 6.528121134630501*^10*c1^5*c2^2*c3^2 + 1.2959353641169475*^11*c1^6*c2^2*c3^2 - 1.0493986835673371*^7*c2^3*c3^2 + 2.6740229931503546*^8*c1*c2^3*c3^2 -
2.9733098654365163*^9*c1^2*c2^3*c3^2 + 2.021279641666957*^10*c1^3*c2^3*c3^2 - 9.458579618924527*^10*c1^4*c2^3*c3^2 + 2.8014186238640594*^11*c1^5*c2^3*c3^2 - 3.795151684822233*^11*c1^6*c2^3*c3^2 +
3.0018147354434848*^7*c2^4*c3^2 - 9.954338473361474*^8*c1*c2^4*c3^2 + 1.3964445468930927*^10*c1^2*c2^4*c3^2 - 1.064279825807562*^11*c1^3*c2^4*c3^2 + 4.6536821418298663*^11*c1^4*c2^4*c3^2 -
1.1059116472581458*^12*c1^5*c2^4*c3^2 + 1.113866914082069*^12*c1^6*c2^4*c3^2 - 5.178558258595619*^7*c2^5*c3^2 + 1.638247097835697*^9*c1*c2^5*c3^2 - 2.153543748403822*^10*c1^2*c2^5*c3^2 +
1.5103063055335773*^11*c1^3*c2^5*c3^2 - 5.974594786251178*^11*c1^4*c2^5*c3^2 + 1.2664737386600254*^12*c1^5*c2^5*c3^2 - 1.1254345721407913*^12*c1^6*c2^5*c3^2 + 3.508191276491249*^7*c2^6*c3^2 -
1.0518902105966748*^9*c1*c2^6*c3^2 + 1.3141541798710089*^10*c1^2*c2^6*c3^2 - 8.756306522358801*^10*c1^3*c2^6*c3^2 + 3.2818453962212317*^11*c1^4*c2^6*c3^2 - 6.560153600346769*^11*c1^5*c2^6*c3^2 +
5.4638485953785016*^11*c1^6*c2^6*c3^2 + 16029.069407182913*c3^3 - 1.4221516248009166*^7*c1*c3^3 + 2.679699296267854*^8*c1^2*c3^3 - 1.797613696507981*^9*c1^3*c3^3 +
5.943738639562176*^9*c1^4*c3^3 - 1.776566664414159*^10*c1^5*c3^3 + 4.339850126561682*^10*c1^6*c3^3 - 8.345161805281889*^6*c2*c3^3 + 2.351038179268103*^8*c1*c2*c3^3 -
2.556763097958812*^9*c1^2*c2*c3^3 + 1.484041581685373*^10*c1^3*c2*c3^3 - 6.1730200989166855*^10*c1^4*c2*c3^3 + 2.1198282836744748*^11*c1^5*c2*c3^3 - 3.799324303120444*^11*c1^6*c2*c3^3 +
2.060074313558452*^7*c2^2*c3^3 - 5.710662387683784*^8*c1*c2^2*c3^3 + 6.933132162419197*^9*c1^2*c2^2*c3^3 - 5.164305811038682*^10*c1^3*c2^2*c3^3 + 2.6436742399217947*^11*c1^4*c2^2*c3^3 -
8.450677606308281*^11*c1^5*c2^2*c3^3 + 1.202486643400443*^12*c1^6*c2^2*c3^3 - 5.30139018161296*^7*c2^3*c3^3 + 1.8173228134829347*^9*c1*c2^3*c3^3 - 2.6658534210166084*^10*c1^2*c2^3*c3^3 +
2.1501438516796686*^11*c1^3*c2^3*c3^3 - 1.0065036020618214*^12*c1^4*c2^3*c3^3 + 2.5859417316610435*^12*c1^5*c2^3*c3^3 - 2.833116754128015*^12*c1^6*c2^3*c3^3 + 1.3528095960521603*^8*c2^4*c3^3 -
4.561656737135868*^9*c1*c2^4*c3^3 + 6.394470239371801*^10*c1^2*c2^4*c3^3 - 4.7804553226309796*^11*c1^3*c2^4*c3^3 + 2.0145013786973215*^12*c1^4*c2^4*c3^3 - 4.545281498309997*^12*c1^5*c2^4*c3^3 +
4.2955378831305464*^12*c1^6*c2^4*c3^3 - 1.4038328886608866*^8*c2^5*c3^3 + 4.453359477695695*^9*c1*c2^5*c3^3 - 5.868697984160185*^10*c1^2*c2^5*c3^3 + 4.116783021460412*^11*c1^3*c2^5*c3^3 -
1.6227309804615405*^12*c1^4*c2^5*c3^3 + 3.4103653698807036*^12*c1^5*c2^5*c3^3 - 2.9872371843791484*^12*c1^6*c2^5*c3^3 + 4.347360399890937*^7*c2^6*c3^3 - 1.3035052784107587*^9*c1*c2^6*c3^3 +
1.6285035195219995*^10*c1^2*c2^6*c3^3 - 1.0850839428197389*^11*c1^3*c2^6*c3^3 + 4.066871954703166*^11*c1^4*c2^6*c3^3 - 8.129360611110518*^11*c1^5*c2^6*c3^3 +
6.770816395822441*^11*c1^6*c2^6*c3^3 + 5.11736249224076*^6*c3^4 - 1.2940899742885114*^8*c1*c3^4 + 1.3467206530712311*^9*c1^2*c3^4 - 8.473966072647891*^9*c1^3*c3^4 +
4.212661300460848*^10*c1^4*c3^4 - 1.5420246310297424*^11*c1^5*c3^4 + 2.612293075165982*^11*c1^6*c3^4 - 1.6257511840092458*^7*c2*c3^4 + 5.412254600477614*^8*c1*c2*c3^4 -
8.577515557890011*^9*c1^2*c2*c3^4 + 8.238513448531758*^10*c1^3*c2*c3^4 - 4.827974475520635*^11*c1^4*c2*c3^4 + 1.5578972891303337*^12*c1^5*c2*c3^4 - 2.095388335895248*^12*c1^6*c2*c3^4 +
6.925691332739684*^7*c2^2*c3^4 - 2.6883150522065763*^9*c1*c2^2*c3^4 + 4.439884591471104*^10*c1^2*c2^2*c3^4 - 3.9690004109212604*^11*c1^3*c2^2*c3^4 + 2.01355666309574*^12*c1^4*c2^2*c3^4 -
5.4692924431941*^12*c1^5*c2^2*c3^4 + 6.190452662809575*^12*c1^6*c2^2*c3^4 - 2.0887671880898958*^8*c2^3*c3^4 + 7.539209481791954*^9*c1*c2^3*c3^4 - 1.1322701335249127*^11*c1^2*c2^3*c3^4 +
9.062118007933304*^11*c1^3*c2^3*c3^4 - 4.0785453195595693*^12*c1^4*c2^3*c3^4 + 9.790989007273338*^12*c1^5*c2^3*c3^4 - 9.797230617841871*^12*c1^6*c2^3*c3^4 + 2.841011380977121*^8*c2^4*c3^4 -
9.53526046182426*^9*c1*c2^4*c3^4 + 1.3284284920609125*^11*c1^2*c2^4*c3^4 - 9.842365899214786*^11*c1^3*c2^4*c3^4 + 4.093656706140883*^12*c1^4*c2^4*c3^4 - 9.06944990355551*^12*c1^5*c2^4*c3^4 +
8.367233517685572*^12*c1^6*c2^4*c3^4 - 1.7179713794516996*^8*c2^5*c3^4 + 5.4316208750132885*^9*c1*c2^5*c3^4 - 7.136283501189009*^10*c1^2*c2^5*c3^4 + 4.990275429780605*^11*c1^3*c2^5*c3^4 -
1.9599143374877312*^12*c1^4*c2^5*c3^4 + 4.1008762635135396*^12*c1^5*c2^5*c3^4 - 3.572678824220341*^12*c1^6*c2^5*c3^4 + 3.8475751613931306*^7*c2^6*c3^4 - 1.1536505075778613*^9*c1*c2^6*c3^4 +
1.4412860024467558*^10*c1^2*c2^6*c3^4 - 9.603395261465988*^10*c1^3*c2^6*c3^4 + 3.5993324864151556*^11*c1^4*c2^6*c3^4 - 7.194785591298389*^11*c1^5*c2^6*c3^4 +
5.992423583647064*^11*c1^6*c2^6*c3^4 - 570736.0928192678*c3^5 - 7.29072212817693*^7*c1*c3^5 + 2.6767927963567743*^9*c1^2*c3^5 - 3.564657151210681*^10*c1^3*c3^5 +
2.3205158920013574*^11*c1^4*c3^5 - 7.484708095843163*^11*c1^5*c3^5 + 9.610578503463467*^11*c1^6*c3^5 - 5.404116909913511*^7*c2*c3^5 + 2.3903352382702436*^9*c1*c2*c3^5 -
4.2701901027519356*^10*c1^2*c2*c3^5 + 3.974114470350878*^11*c1^3*c2*c3^5 - 2.0420962630364282*^12*c1^4*c2*c3^5 + 5.511090147081978*^12*c1^5*c2*c3^5 - 6.116187006466256*^12*c1^6*c2*c3^5 +
2.3111670554659694*^8*c2^2*c3^5 - 8.862401862039667*^9*c1*c2^2*c3^5 + 1.4020484628226752*^11*c1^2*c2^2*c3^5 - 1.1723047132665042*^12*c1^3*c2^2*c3^5 + 5.467531887689922*^12*c1^4*c2^2*c3^5 -
1.3493307841533984*^13*c1^5*c2^2*c3^5 + 1.3772109559723006*^13*c1^6*c2^2*c3^5 - 3.8321278669981676*^8*c2^3*c3^5 + 1.3612498971470133*^10*c1*c2^3*c3^5 - 2.0022543260992398*^11*c1^2*c2^3*c3^5 +
1.5613596632552058*^12*c1^3*c2^3*c3^5 - 6.809568025665291*^12*c1^4*c2^3*c3^5 + 1.575260232515778*^13*c1^5*c2^3*c3^5 - 1.5104233003475846*^13*c1^6*c2^3*c3^5 + 3.119125293745963*^8*c2^4*c3^5 -
1.0428474965370953*^10*c1*c2^4*c3^5 + 1.4461452643130386*^11*c1^2*c2^4*c3^5 - 1.064943657528863*^12*c1^3*c2^4*c3^5 + 4.393389680376872*^12*c1^4*c2^4*c3^5 - 9.629761363676475*^12*c1^5*c2^4*c3^5 +
8.763293975809722*^12*c1^6*c2^4*c3^5 - 1.2480281811308348*^8*c2^5*c3^5 + 3.948303808072467*^9*c1*c2^5*c3^5 - 5.190377345869161*^10*c1^2*c2^5*c3^5 + 3.630065773632401*^11*c1^3*c2^5*c3^5 -
1.424885981885035*^12*c1^4*c2^5*c3^5 + 2.976856685347068*^12*c1^5*c2^5*c3^5 - 2.586503005541738*^12*c1^6*c2^5*c3^5 + 1.9615411836751346*^7*c2^6*c3^5 - 5.881452310245018*^8*c1*c2^6*c3^5 +
7.347853473069412*^9*c1^2*c2^6*c3^5 - 4.895929128946633*^10*c1^3*c2^6*c3^5 + 1.8349840119268198*^11*c1^4*c2^6*c3^5 - 3.667990267390709*^11*c1^5*c2^6*c3^5 + 3.055011313955417*^11*c1^6*c2^6*c3^5 +
2.2829443712770708*^7*c3^6 - 8.74657066262634*^8*c1*c3^6 + 1.3962688439914501*^10*c1^2*c3^6 - 1.1887732216584016*^11*c1^3*c3^6 + 5.693135318961332*^11*c1^4*c3^6 -
1.4541283027798484*^12*c1^5*c3^6 + 1.5475429343702817*^12*c1^6*c3^6 - 1.1414721856385353*^8*c2*c3^6 + 4.214832964236605*^9*c1*c2*c3^6 - 6.475449711264695*^10*c1^2*c2*c3^6 +
5.2977937052167883*^11*c1^3*c2*c3^6 - 2.434021621874772*^12*c1^4*c2*c3^6 + 5.9534963121059*^12*c1^5*c2*c3^6 - 6.05560278666632*^12*c1^6*c2*c3^6 + 2.3780670534136152*^8*c2^2*c3^6 -
8.450792910750084*^9*c1*c2^2*c3^6 + 1.2482397116653235*^11*c1^2*c2^2*c3^6 - 9.808128148134647*^11*c1^3*c2^2*c3^6 + 4.323512368750169*^12*c1^4*c2^2*c3^6 - 1.0136291335540053*^13*c1^5*c2^2*c3^6 +
9.87326541304291*^12*c1^6*c2^2*c3^6 - 2.642296726015127*^8*c2^3*c3^6 + 9.022982014082117*^9*c1*c2^3*c3^6 - 1.2800108643851454*^11*c1^2*c2^3*c3^6 + 9.655407994743712*^11*c1^3*c2^3*c3^6 -
4.084517716372149*^12*c1^4*c2^3*c3^6 + 9.187648004947873*^12*c1^5*c2^3*c3^6 - 8.585448185254704*^12*c1^6*c2^3*c3^6 + 1.651435453759454*^8*c2^4*c3^6 - 5.410121329581757*^9*c1*c2^4*c3^6 +
7.363626541833797*^10*c1^2*c2^4*c3^6 - 5.330504428334047*^11*c1^3*c2^4*c3^6 + 2.1647477676009907*^12*c1^4*c2^4*c3^6 - 4.676639493504751*^12*c1^5*c2^4*c3^6 + 4.199404003657192*^12*c1^6*c2^4*c3^6 -
5.504784845864844*^7*c2^5*c3^6 + 1.726959633454063*^9*c1*c2^5*c3^6 - 2.253002416429701*^10*c1^2*c2^5*c3^6 + 1.564802015954863*^11*c1^3*c2^5*c3^6 - 6.103253617905876*^11*c1^4*c2^5*c3^6 +
1.2676488844486382*^12*c1^5*c2^5*c3^6 - 1.0954966966062238*^12*c1^6*c2^5*c3^6 + 7.645534508145615*^6*c2^6*c3^6 - 2.2924242921956587*^8*c1*c2^6*c3^6 + 2.8639861225801277*^9*c1^2*c2^6*c3^6 -
1.90829514140837*^10*c1^3*c2^6*c3^6 + 7.152250333483447*^10*c1^4*c2^6*c3^6 - 1.4296791929871854*^11*c1^5*c2^6*c3^6 + 1.190757278919808*^11*c1^6*c2^6*c3^6

In[103]:=

{Max[#], Min[#], Min[Abs[#]], Apply[Plus,#], Apply[Plus,Abs[#]], Apply[Plus,Abs[#]]/Length[#]}&[
    (Apply[List,meritFunction]/.{c1->.1,c2->-.1,c3->-.1})]

Out[103]=

In[104]:=

meritFunctionTruncated =
    Apply[Plus,Transpose[
    Select[
        Transpose[{(Apply[List,meritFunction]/.{c1->.1,c2->-.1,c3->-.1})//N,
            Apply[List,meritFunction]}]
    ,
        (Abs[#[[1]]]>300)&
    ]][[2]]];
ByteCount[meritFunctionTruncated]

Out[105]=

In[106]:=

{Max[#], Min[#], Min[Abs[#]], Apply[Plus,#], Apply[Plus,Abs[#]], Apply[Plus,Abs[#]]/Length[#]}&[
    (Apply[List,meritFunctionTruncated]/.{c1->.1,c2->-.1,c3->-.1})]

Out[106]=

In[107]:=

Dc1 = D[meritFunctionTruncated,c1];
ByteCount[Dc1]

Out[108]=

In[109]:=

Dc2 = D[meritFunctionTruncated,c2];
ByteCount[Dc2]

Out[110]=

In[111]:=

Dc3 = D[meritFunctionTruncated,c3];
ByteCount[Dc3]

Out[112]=

In[113]:=

extrema = NSolve[{Dc1==0,Dc2==0,Dc3==0},{c1,c2,c3}]

Out[113]=

In[114]:=

Length[extrema]

Out[114]=

In[115]:=

finalextrema =
    Union[Select[extrema,
        ((    c1==Re[c1]&&c2==Re[c2]&&c3==Re[c3]&&
            Abs[Re[c1]]<.1&&Abs[Re[c2]]<.1&&Abs[Re[c3]]<.1
        )/.#)&]]

Out[115]=

In[116]:=

finalextrema =
    Map[Apply[Rule,Transpose[{{c1,c2,c3},#}],2]&,Union[Re[{c1,c2,c3}/.Select[extrema,
        ((    Abs[Re[c1]]<.1&&Abs[Re[c2]]<.1&&Abs[Re[c3]]<.1
        )/.#)&]]]]

Out[116]=

In[117]:=

Map[{meritFunction/.#,FindFocusFast[AnalyzeSystem[{
    LineOfRays[10,NumberOfRays->6],
    ComponentRendering[{Move[
        SphericalLensSurface[(1/c1)/.#,50,ComponentMedium->3/2],100],
        Move[SphericalLensSurface[(1/c2)/.#,50,ComponentMedium->{3/2,9/5}],115],
        Move[SphericalLensSurface[(1/c3)/.#,50,ComponentMedium->9/5],120]}],
        Move[Screen[50],200]}, PlotType->TopView, ShowArrows->False]]}&,
        finalextrema]

Out[117]=

In[118]:=

Timing[localminimum = FindMinimum[meritFunction,{c1,.01},{c2,-.003},{c3,-.01},MaxIterations->200]]

Out[118]=

In[119]:=

FindFocusFast[AnalyzeSystem[{
    LineOfRays[10,NumberOfRays->21],
    ComponentRendering[{Move[
        SphericalLensSurface[1/c1/.localminimum[[2]],50,ComponentMedium->BK7],100],
        Move[SphericalLensSurface[1/c2/.localminimum[[2]],50,ComponentMedium->{BK7,SF11}],115],
        Move[SphericalLensSurface[1/c3/.localminimum[[2]],50,ComponentMedium->SF11],120]}],
    Move[Screen[50],200]}, PlotType->TopView, ShowArrows->False]]

Out[119]=

In[120]:=

FindFocusFast[{
    LineOfRays[10,NumberOfRays->101],
    ComponentRendering[{Move[
        SphericalLensSurface[1/c1/.localminimum[[2]],50,ComponentMedium->BK7],100],
        Move[SphericalLensSurface[1/c2/.localminimum[[2]],50,ComponentMedium->{BK7,SF11}],115],
        Move[SphericalLensSurface[1/c3/.localminimum[[2]],50,ComponentMedium->SF11],120]}],
    Move[Screen[50],200]}]

Out[120]=

Go to list of topics

2.3 Geometric Intensity Calculations

Set-up System

In[121]:=

localminimum = {c1→0.041742,c2→-0.00807904}

Out[121]=

In[122]:=

optics = {
    SphericalLens[1/c1/.localminimum,1/c2/.localminimum,30,20,ComponentMedium->3/2],
    Move[Screen[50],{d+50,50}]};

In[123]:=

opticalsystem = {Move[LineOfRays[10,NumberOfRays->21],{-50,{y,0}}],optics};

In[124]:=

tracedsystem = AnalyzeSystem[opticalsystem,PlotType->TopView];

In[125]:=

data = ReadRays[tracedsystem, {SurfaceCoordinates[[1]], IntersectionNumber}];

In[126]:=

indata = Select[data,(#[[2]]==1)&]

Out[126]=

In[127]:=

dataplot = ListPlot[Map[{#[[1,1]],#[[2,1]]}&,Transpose[{indata,Select[data,(#[[2]]==3)&]}]]];

Go to list of topics

Preliminary Characterization

In[128]:=

FindFocus[opticalsystem]

Out[128]=

In[129]:=

PSF[opticalsystem,NumberOfPoints->512,PaddingFactor->10];

In[130]:=

GPSF[opticalsystem,NumberOfRays->512];

Go to list of topics

Inverse Position

In[131]:=

symbolicy = (SymbolicSurfaceCoordinates/.SymbolicTrace[opticalsystem,
    {{y,0,5},{y∈Reals,d∈Reals}}, ReportedParameters->SymbolicSurfaceCoordinates])[[1]]

Out[131]=

In[132]:=

oplot = Plot[symbolicy/.{d->0},{y,-5, 5}];

In[133]:=

Show[dataplot,oplot];

In[134]:=

Dy = D[symbolicy,y]

Out[134]=

In[135]:=

Cderivative = Compile[{y,d},Evaluate[Dy]];

In[136]:=

ytransfer = Compile[{y,d},Evaluate[symbolicy]];

In[137]:=

ParametricPlot[Evaluate[{symbolicy,1/Abs[Dy]}/.d->0],{y,-5,5}];

In[138]:=

inversey = Solve[ys==symbolicy,y]

Out[138]=

In[139]:=

Plot[Evaluate[y/.inversey/.d->0],{ys,-.02,.02}];

Go to list of topics

Geometric Intensity Analysis

In[140]:=

(*Fast Version including aperture effects: presently=10*)

In[141]:=

Clear[GeometricIntensity];
GeometricIntensity =
    Function[{ys,d},
        Apply[Plus,
        Map[If[Abs[#]<=5.,1/(Abs[dummy2]/.y->#+.0001),0]&,
            Select[
                Chop[dummy1]
            ,
                (Head[#]=!=Symbol&&Im[#]==0)&
            ]
        ]]
    ]/.{dummy1->(y/.inversey),dummy2->Dy};

In[143]:=

GeometricIntensity//InputForm

Out[143]//InputForm=

Function[{ys, d},
Plus @@ (If[Abs[#1] <= 5., 1/(Abs[0.0088364155326911 - 0.02376778501234205*d - 0.0015757513747870293*y^2 - 0.00006062604240982259*d*y^2 - 3.6305315986759584*^-6*y^4 -
         1.3656859916723253*^-7*d*y^4] /. y -> #1 + 0.0001), 0] & ) /@
   Select[Chop[{Root[ys - 0.0088364155326911*#1 + 0.02376778501234205*d*#1 + 0.0005252504582623431*#1^3 + 0.000020208680803274196*d*#1^3 + 7.261063197351916*^-7*#1^5 +
         2.7313719833446508*^-8*d*#1^5 & , 1], Root[ys - 0.0088364155326911*#1 + 0.02376778501234205*d*#1 + 0.0005252504582623431*#1^3 + 0.000020208680803274196*d*#1^3 +
         7.261063197351916*^-7*#1^5 + 2.7313719833446508*^-8*d*#1^5 & , 2], Root[ys - 0.0088364155326911*#1 + 0.02376778501234205*d*#1 + 0.0005252504582623431*#1^3 +
         0.000020208680803274196*d*#1^3 + 7.261063197351916*^-7*#1^5 + 2.7313719833446508*^-8*d*#1^5 & , 3],
      Root[ys - 0.0088364155326911*#1 + 0.02376778501234205*d*#1 + 0.0005252504582623431*#1^3 + 0.000020208680803274196*d*#1^3 + 7.261063197351916*^-7*#1^5 + 2.7313719833446508*^-8*d*#1^5 & , 4],
      Root[ys - 0.0088364155326911*#1 + 0.02376778501234205*d*#1 + 0.0005252504582623431*#1^3 + 0.000020208680803274196*d*#1^3 + 7.261063197351916*^-7*#1^5 + 2.7313719833446508*^-8*d*#1^5 & ,
       5]}], Head[#1] =!= Symbol && Im[#1] == 0 & ]]

In[144]:=

Options[NIntegrate]

Out[144]=

In[145]:=

Off[NIntegrate::slwcon];
NIntegrate[GeometricIntensity[ys,0],{ys,-.04,.04},MaxRecursion->12]

Out[146]=

In[147]:=

On[NIntegrate::slwcon];

In[148]:=

Plot[GeometricIntensity[ys,0],{ys,-.03,.03}];

In[149]:=

paraxialfocus = d/.Solve[Evaluate[(Dy/.y->0)==0],d][[1]]

Out[149]=

In[150]:=

Plot[GeometricIntensity[ys,paraxialfocus],{ys,-.03,.03}];

In[151]:=

Plot[GeometricIntensity[0,d],{d,-.5,.5}];

In[152]:=

Plot[GeometricIntensity[0,d],{d,-.5,.5},PlotRange->{0,1000}];

In[113]:=

Note that the following plots are very time consuming to construct. Some systems may require 10 - 15 minutes for each plot. The results are well worth the wait however.

In[153]:=

Plot3D[GeometricIntensity[ys,d], {ys,-.05,.05}, {d,-.5, .5},PlotPoints->150,PlotRange->{0,500}];

In[154]:=

ContourPlot[GeometricIntensity[ys,d], {ys,-.05,.05}, {d,-.5, .5},PlotPoints->200,PlotRange->{0,500}, Contours->20, ColorFunction->Function[Hue[.65-#*(.65),1,#*.9+.1]]];

In[155]:=

DensityPlot[GeometricIntensity[ys,d], {ys,-.05,.05}, {d,-.5, .65}, PlotPoints->200, PlotRange->{0,260}, ColorFunction->Function[Hue[.65-#*(.65),1,#*.9+.1]], Mesh->False];

Go to list of topics

3. ABCD Matrix Calculations: FindABCDMatrix

3.1 Overview

In[156]:=

?FindABCDMatrix

In[157]:=

Options[FindABCDMatrix]

In[158]:=

opticalsystemreflected = {
GaussianBeam[40.,.01,SymbolicIntensity->100,SymbolicWaveLength->λ,NumberOfRays->6],
Move[PlanoConvexLens[{f,100},{a,50},{t,10},DesignWaveLength->.532],{d1,50}],
Move[Mirror[{b,20}],175,-45],
Move[Screen[50],{175,50},90]};

In[159]:=

AnalyzeSystem[opticalsystemreflected,PlotType->TopView];

In[160]:=

FindABCDMatrix[opticalsystemreflected,{f>0},MatrixForm->True]

In[161]:=

FindABCDMatrix[opticalsystemreflected,{f>0},MatrixForm->True,DecomposeABCD->True]

In[162]:=

FindABCDMatrix[opticalsystemreflected,MatrixForm->True,DecomposeABCD->True,NumericalResults->True]

Go to list of topics

3.2 SphericalMirror Examples

In[12]:=

opticalsystem = {
    Move[LineOfRays[40,IntrinsicMedium->Vacuum],{d1,-70}],
    SphericalMirror[{R,-100},50],
    Move[Screen[25],{d2,-50}]};

In[13]:=

AnalyzeSystem[opticalsystem,PlotType->TopView];

In[17]:=

result = FindABCDMatrix[opticalsystem,{R<0,d1<0,d2<0},ABCDConstruction->Horizontal]

In[18]:=

Simplify[result/.{d1->0,d2->0},{R<0}]

Go to list of topics

3.3 PlanoConvexLens Examples

In[19]:=

opticalsystem = {
Move[LineOfRays[40,SymbolicIntensity->100,SymbolicWaveLength->λ], {-d1,-50}],
PlanoConvexLens[{f,100},{a,50},{t,10},DesignWaveLength->.532],
Move[Screen[50],{d2+t,100}]};

In[20]:=

AnalyzeSystem[opticalsystem,PlotType->TopView];

In[21]:=

FindABCDMatrix[opticalsystem,{f>0},MatrixForm->True, DecomposeABCD->True, ABCDConstruction->Horizontal]

In[23]:=

FindABCDMatrix[opticalsystem,{f>0,a>0,d1>0,t>0,d2>0}, MatrixForm->True, DecomposeABCD->True,
    ComplexABCD->True, ABCDConstruction->Horizontal]

Go to list of topics

In[24]:=

offaxisopticalsystem = {
Move[LineOfRays[40,SymbolicIntensity->100,SymbolicWaveLength->λ], {-d1,-50}],
PlanoConvexLens[{f,100},{a,50},{t,20},DesignWaveLength->.532,OffAxis->{10,0}],
Move[Screen[50],{d2+t,100}]};

In[25]:=

AnalyzeSystem[offaxisopticalsystem,PlotType->TopView];

In[26]:=

FindABCDMatrix[offaxisopticalsystem,{f>0,a>0,d1>0,t>0,d2>0}, MatrixForm->True, DecomposeABCD->True, ComplexABCD->True]

In[27]:=

result0 = FindABCDMatrix[offaxisopticalsystem,{f>0,a>0,d1>0,t>0,d2>0}, MatrixForm->False, DecomposeABCD->True, ComplexABCD->True, ABCDConstruction->Horizontal]

In[28]:=

FindABCDMatrix[offaxisopticalsystem,{f>0,a>0,d1>0,t>0,d2>0}, MatrixForm->True, ABCDConstruction->Horizontal]

Go to list of topics

In[29]:=

gaussiansystem3D = {
GaussianBeam[{1,1},{.1,.1},SymbolicIntensity->100,SymbolicWaveLength->λ],
Move[PlanoConvexCylindricalLens[{f,100},{{a,50},{b,40}},{t,10},DesignWaveLength->.532],{d1,50}],
Move[Screen[{{h,30},{v,50}}],{d2+t,100}, TwistAngle->30]};

In[30]:=

gaussiansystem = {
GaussianBeam[1,.1,SymbolicIntensity->100,SymbolicWaveLength->λ],
Move[PlanoConvexCylindricalLens[{f,100},{{a,50},{b,40}},{t,10},DesignWaveLength->.532],{d1,50}],
Move[Screen[{{h,30},{v,50}}],{d2+t,100}, TwistAngle->30]};

In[31]:=

FindABCDMatrix[gaussiansystem3D,{f>0},MatrixForm->True]

In[32]:=

FindABCDMatrix[gaussiansystem,{f>0},MatrixForm->True]

In[47]:=

cylindricalsystem = {
Move[LineOfRays[40,SymbolicIntensity->100,SymbolicWaveLength->λ, IntrinsicMedium->Vacuum], {-d1,-50}],
PlanoConvexCylindricalLens[{f,100},{{a,50},{b,40}},{t,10},DesignWaveLength->1/2, ComponentMedium->3/2],
Move[Screen[{{h,30},{v,50}}],{d2+t,100}, TwistAngle->30]};

In[34]:=

AnalyzeSystem[cylindricalsystem];

In[48]:=

FindABCDMatrix[cylindricalsystem, {f>0}, MatrixForm->True]

In[49]:=

FindABCDMatrix[cylindricalsystem, {f>0}, MatrixForm->True, ABCDConstruction->Horizontal]

In[50]:=

FindABCDMatrix[cylindricalsystem, {f>0}, MatrixForm->True, ABCDConstruction->Vertical]

In[38]:=

FindABCDMatrix[cylindricalsystem, {f>0}, MatrixForm->True, DecomposeABCD->True]

In[39]:=

FindABCDMatrix[cylindricalsystem, {f>0}, MatrixForm->True, DecomposeABCD->True, ComplexABCD->True]

In[40]:=

FindABCDMatrix[cylindricalsystem, MatrixForm->True, NumericalResults->True]

In[41]:=

FindABCDMatrix[cylindricalsystem, MatrixForm->True, NumericalResults->True, ComplexABCD->True]

In[42]:=

FindABCDMatrix[cylindricalsystem, MatrixForm->True, NumericalResults->True, ComplexABCD->True, DecomposeABCD->True]

Go to list of topics

3.4 SphericalLens Examples

In[35]:=

opticalsystem = {Move[LineOfRays[40, IntrinsicMedium->Vacuum],{{-d1,-50},{y,0}},{s/Degree,0}],
SphericalLens[{R1,100},{R2,-100},50,{t,10},DesignWaveLength->1/2],
Move[Screen[50],{d2,100}]};

In[93]:=

AnalyzeSystem[opticalsystem,PlotType->TopView];

In[36]:=

FindABCDMatrix[opticalsystem,y,s,{R1>0,R2<0,d1>0,t>0,d2>0,n2>0,n1>0},
            SymbolicRefractiveModels->{Vacuum->n1,BK7->n2}, MatrixForm->False]

Out[36]=

In[37]:=

result = FindABCDMatrix[opticalsystem,y,s,{R1>0,R2<0,d1>0,t>0,d2>0,n2>0,n1>0},
            SymbolicRefractiveModels->{Vacuum->n1,BK7->n2}, MatrixForm->True]

Out[37]//MatrixForm=

In[38]:=

result/.{d1->0,t->0,d2->0}

Out[38]//MatrixForm=

In[39]:=

result/.{d1->0,d2->0}

Out[39]//MatrixForm=

In[40]:=

result/.{d1->d,t->0,d2->0}

Out[40]//MatrixForm=

Go to list of topics

3.5 ThinLens, and Custom ABCD Input

In[99]:=

thinsystem = {Move[LineOfRays[40],{-d1,-50}],
ThinLens[{f,100},50],
Move[Screen[50],{d2,100}]};

In[100]:=

AnalyzeSystem[thinsystem,PlotType->TopView];

In[101]:=

result = FindABCDMatrix[thinsystem,{f>0,d1>0,d2>0}, ABCDConstruction->Horizontal]

In[102]:=

result/.{d1->0,d2->0}

In[103]:=

result/.{d1->d,d2->0}

In[104]:=

thicksystem = {Move[LineOfRays[40],{-d1,-50}],
ThickLens[{f,100},50,{t,10}],
Move[Screen[50],{d2,100}]};

In[105]:=

AnalyzeSystem[thicksystem,PlotType->TopView];

In[106]:=

TurboPlot[thicksystem,PlotType->TopView];

In[107]:=

result = FindABCDMatrix[thicksystem,{f>0,d1>0,t>0,d2>0}, ABCDConstruction->Horizontal]

In[108]:=

result/.{d1->0,t->0,d2->0}

In[109]:=

result/.{d1->0,d2->0}

In[110]:=

result/.{d1->d,t->0,d2->0}

You can enter your own ABCD matrix expressions using Rayica's built-in ABCDOptic function.

In[111]:=

abcdsystem = {Move[LineOfRays[40],{-d1,-50}],
ABCDOptic[{{1,0},{{-1/f,-1/100},1}},{a,50}],
Move[Screen[50],{d2,100}]};

In[112]:=

AnalyzeSystem[abcdsystem,PlotType->TopView];

In[113]:=

result = FindABCDMatrix[abcdsystem,{f>0}, ABCDConstruction->Horizontal]

In[114]:=

result/.{d1->0,d2->0}

In[115]:=

FindABCDMatrix[abcdsystem,{f>0}, ABCDConstruction->Horizontal, ComplexABCD->True, DecomposeABCD->True]

ABCDOptic now also works with astigmatic systems.

In[116]:=

abcdsystem = {Move[GridOfRays[40],{-d1,-50}],
ABCDOptic[{{{1,0},{{-1/f1,-1/100},1}},{{1,0},{{-1/f2,-1/200},1}}},{a,50}],
Move[Screen[50],{d2,100}]};

In[117]:=

AnalyzeSystem[abcdsystem];

In[118]:=

result = FindABCDMatrix[abcdsystem,{f>0}]

In[119]:=

result/.{d1->0,d2->0}

Go to list of topics

3.6 Grating Examples

In[120]:=

?Grating

In[13]:=

gratingsystem = {Move[SingleRay[SymbolicWaveLength->λ, IntrinsicMedium->Vacuum],{-d1,-50}],
Grating[{g,200},50,GratingMedium->IntrinsicMedium,GratingThickness->0,ComponentMedium->IntrinsicMedium],
Move[Screen[50],{d2,100}]};

In[14]:=

sol = AnalyzeSystem[gratingsystem,PlotType->TopView];

In[124]:=

ReadRays[sol,RayEnd]

Out[124]=

In[11]:=

result = FindABCDMatrix[gratingsystem, {λ>0,g>0,d1>0,t>0,d2>0},MatrixForm->True, SymbolicRefractiveModels->{Vacuum->1}, ABCDConstruction->Horizontal, FilterTrace->False]

Out[11]=

In[10]:=

result/.{d1->0,d2->0}

Out[10]//MatrixForm=

In[12]:=

result = FindABCDMatrix[gratingsystem,{g>0,λ>0,d1>0,t>0,d2>0},MatrixForm->True, SymbolicRefractiveModels->{Vacuum->1}, ABCDConstruction->Horizontal, FilterTrace->True]

Out[12]//MatrixForm=

In[13]:=

FindABCDMatrix[gratingsystem, NumericalResults->True,MatrixForm->False, FilterTrace->False]

Out[13]=

Go to list of topics

4. Gaussian Beam Calculations: GaussianTrace

4.1 Overview

In[228]:=

?GaussianTrace

In[229]:=

Options[GaussianTrace]

In[230]:=

?GaussianPlot

In[231]:=

Options[GaussianPlot]

In[232]:=

?GaussianBeam

In[233]:=

Options[GaussianBeam]

In[17]:=

gratingsystem = {Move[GaussianBeam[5.,.01,SymbolicWaveLength->λ,NumberOfRays->3],{-d1,-50}],
Grating[{g,400},50,GratingMedium->IntrinsicMedium,GratingThickness->0,ComponentMedium->IntrinsicMedium],
Move[Screen[50],{d2,100}]};

In[18]:=

sol = AnalyzeSystem[gratingsystem,PlotType->TopView, ShowGaussian->True];

Go to list of topics

4.2 Reflected Lens Example

Set-up

In[24]:=

opticalsystemreflected = {
GaussianBeam[20.,.01,SymbolicWaveLength->λ,NumberOfRays->6],
Move[PlanoConvexLens[{f,100},{a,50},{t,10},DesignWaveLength->.532],{d1,50}],
Move[Mirror[{b,20}],175,-45],
Move[Screen[50],{175,50},90]};

In[21]:=

AnalyzeSystem[opticalsystemreflected, PlotType->TopView];

In[25]:=

FindABCDMatrix[opticalsystemreflected,MatrixForm->True,DecomposeABCD->True,NumericalResults->True]

In[27]:=

FindABCDMatrix[opticalsystemreflected,{f>0},MatrixForm->True,DecomposeABCD->True,NumericalResults->False,SymbolicRefractiveModels->{Air->1,BK7->3/2}]

Go to list of topics

GaussianTrace

In[238]:=

GaussianTrace[opticalsystemreflected]

Go to list of topics

GaussianPlot

In[29]:=

GaussianPlot[opticalsystemreflected];

In[30]:=

GaussianPlot[opticalsystemreflected, RenderedParameters->BeamCurvature];

In[31]:=

GaussianPlot[opticalsystemreflected, RenderedParameters->ComplexBeamParameter];

In[32]:=

finalresultreflected =
    GaussianPlot[opticalsystemreflected, Plot3D->True, PlotType->TopView]

Go to list of topics

ShowGaussian

In[33]:=

ShowSystem[opticalsystemreflected,PlotType->TopView,ShowGaussian->True];

In[34]:=

AnalyzeSystem[opticalsystemreflected,PlotType->TopView,ShowGaussian->True];

Go to list of topics

4.3 Reflected ThickLens Example

Set-up

In[35]:=

thinopticalsystemreflected2 = {
GaussianBeam[{20.,15.},{.01,.05},SymbolicWaveLength->λ,NumberOfRays->6],
Move[ThinLens[{{f1,100},{f2,150}},50],{d1,50}],
Move[Mirror[{b,20}],175,-45],
Move[Screen[50],{175,50},90]};

In[36]:=

AnalyzeSystem[thinopticalsystemreflected2];

In[37]:=

FindABCDMatrix[thinopticalsystemreflected2,MatrixForm->True,DecomposeABCD->True,NumericalResults->True,SymbolicRefractiveModels->{Air->1}]

In[38]:=

FindABCDMatrix[thinopticalsystemreflected2,MatrixForm->True,DecomposeABCD->True,NumericalResults->False,SymbolicRefractiveModels->{Air->1}]

Go to list of topics

GaussianTrace

In[39]:=

GaussianTrace[thinopticalsystemreflected2]

Go to list of topics

GaussianPlot

In[40]:=

GaussianPlot[thinopticalsystemreflected2]

Go to list of topics

ShowGaussian

In[41]:=

ShowSystem[thinopticalsystemreflected2,ShowGaussian->True,CreateStereoView->Crossed];

Go to list of topics

4.4 Off-Axis Beam Example

Set-up

In[61]:=

opticalsystemoffaxis = {
Move[GaussianBeam[{5.,5.},{.01,.01},SymbolicWaveLength->λ],{0,10}],
Move[PlanoConvexLens[{f,100},{a,50},{t,10},DesignWaveLength->.532],{d1,50}],
Move[Mirror[20],175,-45],
Move[Screen[50],{175,50},90]};

In[62]:=

AnalyzeSystem[opticalsystemoffaxis];

In[63]:=

FindABCDMatrix[opticalsystemoffaxis,MatrixForm->True,DecomposeABCD->True,NumericalResults->True]

In[64]:=

FindABCDMatrix[opticalsystemoffaxis,NumericalResults->True,MatrixForm->True]

Go to list of topics

GaussianTrace

In[65]:=

GaussianTrace[opticalsystemoffaxis]

Go to list of topics

GaussianPlot

In[66]:=

GaussianPlot[opticalsystemoffaxis];

In[68]:=

GaussianPlot[opticalsystemoffaxis,RenderedParameters->ComplexBeamParameter];

Go to list of topics

ShowGaussian

In[69]:=

ShowSystem[opticalsystemoffaxis, PlotType->TopView, ShowGaussian->True];

In[70]:=

ShowSystem[opticalsystemoffaxis, PlotType->SideView, ShowGaussian->True];

In[71]:=

ShowSystem[opticalsystemoffaxis, ShowGaussian->True];

Go to list of topics

4.5 CylindricalLens Example

Set-up

In[72]:=

cylindricalsystem = {
GaussianBeam[{5,10},{.01,.005},SymbolicWaveLength->λ],
Move[PlanoConvexCylindricalLens[{f,100},{{a,50},{b,40}},{t,10},DesignWaveLength->.532],{d1,50}],
Move[Screen[{{h,30},{v,50}}],{d2+t,100}, TwistAngle->30]};

In[73]:=

AnalyzeSystem[cylindricalsystem];

In[74]:=

FindABCDMatrix[cylindricalsystem,{f>0},MatrixForm->True,DecomposeABCD->True]

In[75]:=

FindABCDMatrix[cylindricalsystem,NumericalResults->True,MatrixForm->False]

Go to list of topics

GaussianTrace

In[76]:=

GaussianTrace[cylindricalsystem]

Go to list of topics

GaussianPlot

In[77]:=

GaussianPlot[cylindricalsystem];

In[78]:=

GaussianPlot[cylindricalsystem,RenderedParameters->ComplexBeamParameter];

Go to list of topics

ShowGaussian

In[79]:=

ShowSystem[cylindricalsystem, PlotType->TopView, ShowGaussian->True];

In[80]:=

ShowSystem[cylindricalsystem, PlotType->SideView, ShowGaussian->True];

In[81]:=

ShowSystem[cylindricalsystem, ShowGaussian->True];

Go to list of topics

4.6 Prism-Lens Example

Set-up

In[82]:=

prism = Prism[{60,50,60},50];
prismgaussiansystem = {
Move[GaussianBeam[{5,5},{.01,.01}],0,20],
Move[prism,{30,-20}],
Move[BiConvexLens[45,30,10],60],
Move[Screen[50],120]};

In[84]:=

AnalyzeSystem[prismgaussiansystem];

In[85]:=

FindABCDMatrix[prismgaussiansystem,MatrixForm->True,DecomposeABCD->True,NumbericalResults->True]

Go to list of topics

GaussianTrace

In[86]:=

GaussianTrace[prismgaussiansystem]

Go to list of topics

GaussianPlot

In[87]:=

GaussianPlot[prismgaussiansystem];

In[88]:=

GaussianPlot[prismgaussiansystem,RenderedParameters->ComplexBeamParameter];

Go to list of topics

ShowGaussian

In[89]:=

ShowSystem[prismgaussiansystem, PlotType->TopView, ShowGaussian->True];

In[90]:=

ShowSystem[prismgaussiansystem, PlotType->SideView, ShowGaussian->True];

In[91]:=

ShowSystem[prismgaussiansystem, ShowGaussian->True];

Go to list of topics

5. Imaging Calculations

5.1 Overview

In[106]:=

Wavica has several built-in functions userful for characterizing imaging systems. These functions include: FindLensParameters, PupilFunction, PointSpreadFunction, GeometricPointSpreadFunction, and OpticalTransferFunction. In addition to these names, many of these functions can also take an abbreviated alias that include: PF, PSF, GPSF, OTF and MTF. The following examples shall serve to introduce the basic use of these highlighted functions. Next, evaluate the examples to see the resulting answers.

In[122]:=

optics = {PlanoConvexLens[100,50,10],
    Move[Screen[1],{d+103,103}]};

In[123]:=

opticalsystem = {Move[LineOfRays[20],{-50,{y,0}}],optics};

In[124]:=

AnalyzeSystem[opticalsystem,PlotType->TopView];

In[125]:=

FindFocus[{Move[LineOfRays[20,NumberOfRays->21],-50],optics}]

Out[125]=

In[126]:=

lensparameters = FindLensParameters[opticalsystem,FindImagePoint->True]

Out[126]=

In[127]:=

pupilfunction = PupilFunction[lensparameters]

Out[127]=

In[128]:=

?PSF

In[129]:=

Options[PointSpreadFunction]

Out[129]=

In[130]:=

psf = PointSpreadFunction[pupilfunction,NumberOfPoints->512,RenderedParameters->PointSpreadFunction,PaddingFactor->8]

Out[130]=

In[131]:=

?MTF

In[132]:=

ModulationTransferFunction[psf]

Out[132]=

In[133]:=

MTF[opticalsystem,NumberOfRays->256,FindImagePoint->True,GeometricPointSpreadFunction->True]

Out[133]=

Go to list of topics

Define Systems

This section serves to define the various optical systems used in later in imaging calculations.

In[134]:=

sys3D =
    {PointOfRays[{10,10}, NumberOfRays -> 5],
    Move[PlanoConvexLens[{f1,100}, 50, 9, CurvatureDirection -> Back],90.5],
    Move[ApertureStop[50,30],100],
    Move[PlanoConvexLens[{f2,100}, 50, 9], 100.5],
    Boundary[208, GraphicDesign -> Off]};

In[135]:=

AnalyzeSystem[sys3D];

In[136]:=

sys =
    {WedgeOfRays[10, NumberOfRays -> 5],
    Move[PlanoConvexLens[{f1,100}, 50, 9, CurvatureDirection -> Back],90.5],
    Move[ApertureStop[50,30],100],
    Move[PlanoConvexLens[{f2,100}, 50, 9], 100.5],
    Boundary[208, GraphicDesign -> Off]};

In[137]:=

AnalyzeSystem[sys,PlotType->TopView];

In[138]:=

offaxis3D =
    {MoveDirected[PointOfRays[{10,10}],{-90.5,30},{0,0},SideOfObject->After],
    BiConvexLens[50, 50, 20],
    Boundary[108, GraphicDesign -> Off]};

In[139]:=

AnalyzeSystem[offaxis3D];

In[140]:=

colsys3D =
    {Move[PointOfRays[{10,10}],-100],
    PlanoConvexLens[100, 50, 10],
    Boundary[108, GraphicDesign -> Off]};

In[141]:=

AnalyzeSystem[colsys3D];

In[142]:=

colsys32 =
    {Move[PointOfRays[{10,10}, NumberOfRays->32],-100],
    PlanoConvexLens[100, 50, 10],
    Boundary[108, GraphicDesign -> Off]};

Go to list of topics

5.2 FindLensParameters

In[143]:=

?FindLensParameters

In[144]:=

Options[FindLensParameters]

Out[144]=

In[145]:=

FindLensParameters[sys]

Out[145]=

FindLensParameters can also be used with astigmatic imaging systems. For example, you can check the focal length along each astigmatic axis.

In[146]:=

FindLensParameters[offaxis3D,ABCDConstruction->Full3D,FindPupils->False]

Out[146]=

Here we can see that it has two focal lengths for each direction.

In[147]:=

FocalLength/.%

Out[147]=

In[148]:=

FindLensParameters[sys,FindImagePoint->True]

Out[148]=

In[149]:=

ImagePoint/.%

Out[149]=

In[150]:=

FindLensParameters[sys,FindImagePoint->True,FocalFraction->0]

Out[150]=

In[151]:=

ImagePoint/.%

Out[151]=

In[152]:=

FindLensParameters[sys,FindImagePoint->True,FocalFraction->1]

Out[152]=

In[153]:=

ImagePoint/.%

Out[153]=

In[154]:=

FindLensParameters[sys,FindImagePoint->True,FocalFraction->.5]

Out[154]=

In[155]:=

ImagePoint/.%

Out[155]=

Go to list of topics

5.3 PupilFunction

In[156]:=

?PupilFunction

In[157]:=

Options[PupilFunction]

Out[157]=

In[158]:=

PupilFunction[lensparameters]

Out[158]=

In[159]:=

PupilFunction[sys3D]

Out[159]=

In[161]:=

PupilFunction[sys3D,Plot2D->False]

Out[161]=

In[18]:=

PupilFunction[offaxis3D,Plot2D->False]

Out[18]=

In[19]:=

pupfun = PupilFunction[sys3D, FindImagePoint->True, FocalFraction->0]

Out[19]=

In[20]:=

PupilFunction[sys,FindPupils->False]

Out[20]=

In[26]:=

PSF[sys, NumberOfPoints->2048, NumberOfRays->256, FocalFraction->0]

Out[26]=

Go to list of topics

5.4 Zernike Polynomials and Seidel Aberrations

In[27]:=

?ZernikeFit

In[28]:=

?SeidelAberrations

In[73]:=

ZernikeFit may be called directly. SeidelAberrations is one of the results returned by ZernikeFit.

In[29]:=

ZernikeFit[ReadTurboRays[TurboTrace[colsys32, ReportedSurfaces->Last],{SurfaceCoordinates,OpticalLength}]]

Out[29]=

In[73]:=

When Compiled->False is given, the output is left uncompiled.

In[30]:=

ZernikeFit[ReadTurboRays[TurboTrace[colsys32, ReportedSurfaces->Last],{SurfaceCoordinates,OpticalLength}],Compiled->False]

Out[30]=

In[73]:=

ZernikeFit is also used internally by PupilFunction for three-dimesional sources when the option ZernikeFit->True is set.

In[31]:=

ZernikeFit/.Options[PupilFunction]

Out[31]=

In[73]:=

ZernikeFit is also used internally by PupilFunction for three-dimensional light sources. SeidelAberrations is also returned by PupilFunction.

In[32]:=

{OpticalPathDifference,ResidualFitError,SeidelAberrations}/.
    PupilFunction[offaxis3D,FindImagePoint->True,Plot2D->False,ZernikeFit->True]

Out[32]=

In[33]:=

{OpticalPathDifference,ResidualFitError,SeidelAberrations}/.
    PupilFunction[colsys3D,FindImagePoint->True,Plot2D->False,ZernikeFit->True]

Out[33]=

Go to list of topics

5.5 PointSpreadFunction

In[34]:=

?PointSpreadFunction

In[35]:=

Options[PointSpreadFunction]

Out[35]=

In[36]:=

PSF[pupfun,NumberOfPoints->128]

Out[36]=

In[38]:=

Options[Graphics3D]

Out[38]=

In[45]:=

Options[Plot3D]

Out[45]=

In[46]:=

PSF[pupfun, NumberOfPoints->128, Plot2D->False, PlotRange->All, PlotPoints->200, Mesh->False]

Out[46]=

In[58]:=

PSF[offaxis3D, NumberOfPoints->1024, NumberOfRays->32, Plot2D->False, FocalFraction->0, PlotPoints->100, ColorFunction->Automatic, Mesh->False]

Out[58]=

In[59]:=

PSF[sys3D, FindImagePoint->True, FocalFraction->0, NumberOfPoints->128]

Out[59]=

In[66]:=

PSF[sys, NumberOfPoints->2048, NumberOfRays->128]

Out[66]=

Go to list of topics

5.6 OpticalTransferFunction

In[63]:=

?OpticalTransferFunction

In[64]:=

?MTF

In[68]:=

MTF[sys, NumberOfPoints->2048, NumberOfRays->128]

Out[68]=

In[72]:=

MTF[sys, NumberOfPoints->2048, NumberOfRays->128, FindImagePoint->True]

Out[72]=

In[75]:=

OpticalTransferFunction[sys, NumberOfPoints->1024, NumberOfRays->2048, FindImagePoint->True, IntensityTransform->False]

Out[75]=

In[76]:=

ModulationTransferFunction[sys3D, FindImagePoint->True, FocalFraction->0, NumberOfPoints->128]//Timing

Out[76]=

In[77]:=

MTF[sys3D, FindImagePoint->True, FocalFraction->0, NumberOfPoints->128, Plot2D->False]

Out[77]=

In[81]:=

OpticalTransferFunction[offaxis3D, FindImagePoint->True, NumberOfPoints->512]

Out[81]=

Out[349]=

In[350]:=

(*Note NumberOfPoints->256 can run out of Kernal memory sometimes*)

In[78]:=

MTF[sys3D, FindImagePoint->True, FocalFraction->1, NumberOfPoints->256]//Timing

Out[78]=

In[80]:=

MTF[sys3D, FindImagePoint->True, FocalFraction->1, NumberOfPoints->256, Plot2D->False, PlotPoints->100, ColorFunction->Automatic, Mesh->False];

Go to list of topics

5.7 GeometricPointSpreadFunction

In[12]:=

?GeometricPointSpreadFunction

In[13]:=

Options[GeometricPointSpreadFunction]

Out[13]=

In[14]:=

GPSF[sys3D, FindImagePoint->True, FocalFraction->0]

Out[184]=

In[15]:=

GPSF[sys3D, FindImagePoint->True, FocalFraction->0, Plot2D->False]

Out[15]=

In[16]:=

GPSF[sys3D, FindImagePoint->True, FocalFraction->0, Plot2D->False]

Out[16]=

In[17]:=

GeometricPointSpreadFunction[offaxis3D, FindImagePoint->True, FocalFraction->0, Plot2D->False]

Out[17]=

In[187]:=

MTF[sys3D, FindImagePoint->True, FocalFraction->0, NumberOfPoints->128, GeometricPointSpreadFunction->True]

Out[187]=

In[185]:=

GPSF[sys, NumberOfPoints->512]

Out[185]=

In[186]:=

GPSF[sys, NumberOfPoints->512, FindImagePoint->True]

Out[186]=

Go to list of topics

6. Interference and Wavefront Calculations

6.1 Overview

In[3]:=

Wavica contains two tools used for interference and wavefront calculations. These functions are FindWaveFronts and FindInterference. Finally, at the end of this chapter, we introduce a method to calculate propagated optical wavefronts that uses Gaussian wavelets.

In[21]:=

?FindWaveFronts

In[22]:=

Options[FindWaveFronts]

Out[22]=

In[23]:=

?FindInterference

In[24]:=

Options[FindInterference]

Out[24]=

Go to list of topics

6.2 One-Dimensional Wavefronts

Example 1

In[25]:=

sys = {
MoveDirected[WedgeOfRays[10],{0,-50},{100,0},SideOfObject->After],
MoveDirected[WedgeOfRays[10],{0,50},{100,0},SideOfObject->After],
Move[Screen[100],100]};

In[26]:=

DrawSystem[sys,PlotType->TopView];

In[27]:=

wavefront = FindWaveFronts[sys]

Out[27]=

In[28]:=

wavefront = FindWaveFronts[sys,PlotDomain->{-.1,.1}]

Out[28]=

In[29]:=

wavefront = FindWaveFronts[sys,PlotDomain->{-.1,.1}]

Out[29]=

In[30]:=

FindInterference[wavefront,PlotDomain->{-.005,.005}];

In[31]:=

FindInterference[sys,PlotDomain->{-.005,.005}];

Go to list of topics

Example 2

In[32]:=

sys = {
MoveDirected[WedgeOfRays[2],{0,-2},{100,0},SideOfObject->After],
MoveDirected[LineOfRays[5],{0,2},{100,0},SideOfObject->After],
Move[Screen[100],100]};

In[33]:=

DrawSystem[sys,PlotType->TopView];

In[34]:=

wavefront = FindWaveFronts[sys]

Out[34]=

In[35]:=

FindInterference[wavefront,PlotDomain->{-.1,.1}];

Go to list of topics

Example 3

This example shows how to create each wave-front separately and interfere them together afterwards.

In[36]:=

Out[36]=

In[37]:=

In[38]:=

In[61]:=

In[65]:=

In[45]:=

In[49]:=

planesource = {Move[WedgeOfRays[10],{-640+115,185}],
            Move[L,{-640+115+100,185}],
            MoveReflected[M,{-640+115+100,185},{-640+115+100+100,185},{-425,0}],
            MoveReflected[M,{-640+115+100,185},{-640+115+100+100+75,185},{-425,0}]};

In[50]:=

planesource = {Move[WedgeOfRays[10],{-640+115,185}],
            Move[L,{-640+115+100,185}],
            MoveReflected[M,{-640+115+100,185},{-640+115+100+100+75,185},{-425,0}]};

In[66]:=

sys1 = Flatten[{planesource,
            Move[HOE,-425,MoveRelative->False]}];

In[69]:=

sys2 = {    reference,
            Move[Piston,20],
            Move[{Move[Baffle[{25,120},Labels→""],{0,.5 65+12.5}],
            Move[Baffle[{25,120},Labels→""],{0,-.5 65-12.5}]},5],
            Move[H,-4.5],
            Move[L1,-30,180],
            Move[L2,-80],
            Move[CL,-375],
            Move[HOE,-425,MoveRelative->False]};

In[70]:=

In[71]:=

wave2 = FindWaveFronts[sys2];

In[72]:=

wave1 = FindWaveFronts[sys1];

In[73]:=

FindInterference[Join[wave1,wave2]];

In[74]:=

FindInterference[Join[wave1,wave2],PlotDomain->{0,.01}];

Go to list of topics

Example 4

This example shows how to create and interfere two wave-fronts simultaneously.

In[3]:=

Out[3]=

In[4]:=

In[5]:=

In[10]:=

In[11]:=

In[12]:=

In[16]:=

planesource = {Move[WedgeOfRays[10],{-640+115,185}],
            Move[L,{-640+115+100,185}],
            MoveReflected[M,{-640+115+100,185},{-640+115+100+100+75,185},{-475,0}]};

In[17]:=

combinedsys =
        {    reference,
            Move[Piston,20],
            Move[{Move[Baffle[{25,120},Labels→""],{0,.5 65+12.5}],
            Move[Baffle[{25,120},Labels→""],{0,-.5 65-12.5}]},5],
            Move[H,-4.5],
            Move[L1,{-30,0},180],
            Move[L2,{-80,0},0],
            Move[CL,-375],
            planesource,
            Move[HOE,-475,MoveRelative->False]};

In[18]:=

combinedsys =
        {    reference,
            planesource,
            Move[Piston,20],
            Move[{Move[Baffle[{25,120},Labels→""],{0,.5 65+12.5}],
            Move[Baffle[{25,120},Labels→""],{0,-.5 65-12.5}]},5],
            Move[H,-4.5],
            Move[L1,{-30,0},180],
            Move[L2,{-80,0},0],
            Move[CL,-375],
            Move[HOE,-475,MoveRelative->False]};

In[19]:=

In[20]:=

FindWaveFronts[combinedsys]

Out[20]=

In[21]:=

FindInterference[combinedsys,PlotDomain->{0,.01}];

Go to list of topics

Example 5

In[275]:=

?GaussianBeam

In[541]:=

system =
    {
        Move[GaussianBeam[1.5,.001],-50],
        PlanoConcaveLens[-30,20,2,"L1",CurvatureDirection->Back],
        Move[PlanoConvexLens[160,30,6.5,"L2"],130],
        Move[BeamSplitter[{70,30},50,10],175,-45],
        Move[Mirror[50,10],250],
        Move[Mirror[50,10],{175,50},90.01],
        Move[BeamSplitter[{70,30},50,10,""],175,-45,GraphicDesign->Off],
        Move[Screen[50],{175,-50},90]
    };

In[542]:=

AnalyzeSystem[system,PlotType->TopView];

In[543]:=

wavefronts = FindWaveFronts[system]

In[544]:=

FindInterference[wavefronts];

In[545]:=

FindInterference[system]

Go to list of topics

6.3 Two-Dimensional Wavefronts

Example 1

In[17]:=

sys = {
MoveDirected[PointOfRays[{2,2}],{0,-2},{100,0},SideOfObject->After],
MoveDirected[GridOfRays[{5,5}],{0,2},{100,0},SideOfObject->After],
Move[Screen[100],100]};

In[18]:=

DrawSystem[sys,PlotType->TopView];

In[19]:=

wavefront = FindWaveFronts[sys,PlotDomain->{{-.5,-.1},{-.2,.2}}]

Out[19]=

In[34]:=

FindInterference[wavefront,PlotDomain->{{-.6,.1},{-.3,.3}}, PlotPoints->500]

Out[34]=

Go to list of topics

Example 2

In[21]:=

L1=PlanoConcaveLens[-30,20,2,"L1",CurvatureDirection->Back];
L2=PlanoConvexLens[160,30,6.5,"L2"];

In[23]:=

sys1 = {
Move[GridOfRays[{3,3}],-50],
L1,
Move[L2,130],
MoveReflected[Mirror[75],{130,0},{200,0},{70,75}],
MoveDirected[Screen[50],{70,75},{200,0}]};

In[24]:=

DrawSystem[sys1,PlotType->TopView,FormatType->OutputForm];

In[25]:=

sys2 = {
Move[GridOfRays[{3,3}],-50],
L1,
Move[L2,130],
MoveReflected[LensSurface[75],{130,0},{200,0},{70,75}],
Move[MoveReflected[Mirror[75],{130,0},{200,0},{70,75}],10],
MoveReflected[LensSurface[75],{130,0},{200,0},{70,75}],
MoveDirected[Screen[50],{70,75},{200,0}]};

In[26]:=

DrawSystem[sys2,PlotType->TopView,FormatType->OutputForm];

In[27]:=

wavefront1 = FindWaveFronts[sys1,ZernikeFit->True]//Timing

Out[27]=

In[28]:=

wavefront2 = FindWaveFronts[sys2]//Timing

Out[28]=

In[29]:=

FindInterference[Join[wavefront1[[2]],wavefront2[[2]]]]//Timing

Out[29]=

Go to list of topics

Example 3

In[30]:=

system =
    {
        Move[GaussianBeam[1.5,.01,FullForm->True],-50],
        PlanoConcaveLens[-30,20,2,"L1",CurvatureDirection->Back],
        Move[PlanoConvexLens[160,30,6.5,"L2"],130],
        Move[BeamSplitter[{70,30},50,10],175,-45],
        Move[Mirror[50,10],250],
        Move[Mirror[50,10],{175,50},90.01],
        Move[BeamSplitter[{70,30},50,10,""],175,-45,GraphicDesign->Off],
        Move[Screen[50],{175,-50},90]
    };

In[31]:=

AnalyzeSystem[system];

In[32]:=

interference = FindInterference[system]

Out[32]=

In[33]:=

Plot3D[Evaluate[(InterferenceFunction/.interference)][x,y],{x,-13.5,6.7},{y,-10.,10.},PlotPoints->100,Mesh->False];

Go to list of topics

6.4 Wavefront Propagation with Gaussian Wavelets

Introduction to Wavelet Analysis

Wavelet analysis is a very new concept that is still in its infancy. This section introduces the idea and uses it in several examples. At the moment, there are two special functions used here for wavelet analysis: a new light source, called LineOfWavelets that constructs a set of GaussianBeam wavelet functions, and a special function, called FindField that calculates the interference between the different Gaussian wavelets. Both LineOfWavelets and FindField will eventually become built-in functions of Wavica, along with several other functions that includes a point source function called WedgeOfWavelets, but for the moment, they are defined here externally until their final formats are fully worked out. The permanent versions of these two functions may be altered and their present formats may not get supported in the future. Presently, these functions must be taken as "experimental". Currently, all calculations are only carried out in the two-dimensions of the optical and transverse axes. In the near future, there will be additional support for full three-dimensional calculations with the GridOfWavelets and PointOfWavelets source functions to be developed as well.

Go to list of topics

Define LineOfWavelets (Evaluate This Subsection)

In this section, we define LineOfWavelets. You must evaluate the following two subsections before running any wavelet analysis examples.

In[517]:=

Clear[LineOfWavelets];

In[518]:=

Options[LineOfWavelets]:={
    PaddingFactor->1,
    NumberOfRays->11,
    WaveLength->.532,
    IntrinsicMedium->Vacuum};

In[519]:=

LineOfWavelets[linewidthin_,opts___]:=
Block[{divergence, y, options, beamspotsize, paddingfactor, numberofrays, wavelength, beamoffset, intrinsicmedium, refractiveindex, linewidth},
    options = Flatten[{opts,Options[LineOfWavelets]}];
    {paddingfactor, numberofrays, wavelength, intrinsicmedium} =
        {PaddingFactor, NumberOfRays, WaveLength, IntrinsicMedium}/.options;
    If[paddingfactor<1&&numberofrays>1,
        Print["Warning, LineOfWavelets will not be uniform."];
        Print["Set PaddingFactor greater or equal to 1 for uniform illumination."]
    ];
    If[numberofrays>1,
        linewidth = linewidthin;
        beamoffset = linewidth/(numberofrays-1)
    ,
        beamoffset = linewidthin;
        linewidth = 0
    ];
    beamspotsize = beamoffset*paddingfactor;
    refractiveindex = ModelRefractiveIndex[intrinsicmedium][options];
    divergence = 2*(wavelength*10^-3)/(Pi*refractiveindex*beamspotsize)+10^-10;
    Table[
        Move[GaussianBeam[beamspotsize,divergence,options],{0,y}],
        {y,-linewidth/2,linewidth/2,beamoffset}]
];

Go to list of topics

Define FindField (Evaluate This Subsection)

In this section, we define FindField. FindField returns a symbolic function that represents either the field intensity (for  IntensityTransform->True) or the complex optical field (IntensityTransform->False). You must evaluate this subsection (as well as the LineOfWavelets subsection) before running any wavelet analysis examples.

In[520]:=

Clear[FindField];

In[521]:=

Options[FindField]:=
    {IntensityTransform->True,
    CoherentSource->True,
    NormalizeOutput->True,
    Compiled->True};

In[522]:=

FindField[opticalsystem_,opts___] :=
Block[{rp, sp, xc, yc, xp, yp, tx, ty, sc, zo, wo, W, R, field, k, wl, ol, minOL, gaussianresults, ds, scale, intensity, results, sourcewaist, width, centers, widths, options, intensitytransform, coherentsource, normalizeoutput, compiled},
    options = Flatten[{opts,Options[FindField]}];
        {intensitytransform, coherentsource, normalizeoutput, compiled} =
            {IntensityTransform, CoherentSource, NormalizeOutput, Compiled}/.options;
    gaussianresults = GaussianTrace[opticalsystem, ABCDConstruction->Horizontal,
        NumericalResults->True, ReportedSurfaces->{-1}];
    If[Head[gaussianresults[[1]]]=!=List,gaussianresults={gaussianresults}];
    Clear[xp,yp];
    centers = widths = {};
    minOL = Sort[OpticalLength/.gaussianresults][[1]];
    scale = 1/Length[gaussianresults];
    field = Map[
        (
        {zo,wo,sc,wl,ol,intensity} =
            {BeamScaleLength,BeamWaist,WaistDistance,WaveLength*10^-3,OpticalLength,Intensity/100}/.#;
        sc = -sc;
        {tx,ty} = (ABCDRotationMatrix/.#)[[1]][[{1,2}]];
        {xc,yc} = (ABCDCenterPoint/.#)[[{1,2}]];
        centers = {centers,yc};
        k = 2 Pi/wl;
        ds = Dot[{xp,yp-yc},{tx,ty}];
        rpSq = (xp)^2+(yp-yc)^2-ds^2;
        ol = ol - minOL + ds;
        sp = sc + ds;
        W = wo*Sqrt[1+(sp/zo)^2];
        width = W/.{xp->0,yp->0};
        widths = {widths,width};
        R = sp*(1+(zo/sp)^2);
        energy = NIntegrate[((Exp[-rpSq/W^2])^2/W)/.xp->0,{yp,yc-3*width,yc+3*width},MaxRecursion->12];
        Sqrt[intensity/W]*
            Exp[-rpSq/W^2]*Exp[-I(k*ol-ArcTan[zo,sp])]*Exp[-I*k*rpSq/(2*R)]/Sqrt[energy]
        )&
    ,
        gaussianresults
    ];
    If[intensitytransform === True,
        If[coherentsource === True,
            field = Abs[Apply[Plus,field]]^2
        ,
            field = Apply[Plus,Map[(Abs[#]^2)&,field]]
        ];
        centers = Sort[Flatten[centers]][[{1,-1}]];
        width = Sort[Flatten[widths]][[-1]];
        If[normalizeoutput=!=False,
            field = field/NIntegrate[field/.xp->0,{yp,-3*width+centers[[1]],3*width+centers[[2]]}]
        ];
    ,
        If[normalizeoutput=!=False,
            field = scale*Apply[Plus,field]
        ,
            field = Apply[Plus,field]
        ]
    ];
    field = field/.{xp->#1,yp->#2};
    If[compiled===True,
        Apply[Compile,{{#1,#2}, field}]
    ,
        Apply[Function,{{#1,#2}, field}]
    ]
]

Go to list of topics

Collimated Beam of Wavelets

In[46]:=

waveletSystem = {
    LineOfWavelets[10,NumberOfRays->9,PaddingFactor->1],
    Move[Screen[50],50]}

Out[46]=

In[47]:=

ShowSystem[waveletSystem,ShowGaussian->True,PlotType->TopView];

In[48]:=

intensity = FindField[waveletSystem]

Out[48]=

In[49]:=

Plot[intensity[0,y],{y,-15,15},PlotRange->All];

In[50]:=

Plot3D[intensity[x,y],{x,0,50},{y,-15,15},PlotRange->All,PlotPoints->50];

Go to list of topics

One Wavelet with Lens

In[51]:=

waveletSystem = {
    Move[LineOfWavelets[10,NumberOfRays->1],-50],
    PlanoConvexLens[100,50,10],
    Move[Screen[50],103]};

In[52]:=

ShowSystem[waveletSystem,ShowGaussian->True,PlotType->TopView];

In[53]:=

intensity = FindField[waveletSystem];

In[54]:=

Plot[intensity[0,y],{y,-.05,.05},PlotRange->All];

In[55]:=

Plot3D[intensity[x,y],{x,0,.55},{y,-.04,.04},PlotRange->All,PlotPoints->50];

Go to list of topics

Two Wavelets with Lens

In[523]:=

waveletSystem = {
    Move[LineOfWavelets[10,NumberOfRays->2,PaddingFactor->.1],-50],
    PlanoConvexLens[100,50,10],
    Move[Screen[50],103]};

In[524]:=

ShowSystem[waveletSystem,ShowGaussian->True,PlotType->TopView];

In[525]:=

Timing[intensity = FindField[waveletSystem]][[1]]/60/Second*"Minute"

In[526]:=

Plot[intensity[0,y],{y,-.07,.07},PlotRange->All];

In[527]:=

Plot[intensity[-10,y],{y,-1,1},PlotRange->All];

In[61]:=

Plot3D[intensity[x,y],{x,-1,1},{y,-.05,.05},PlotRange->All,PlotPoints->100];

Go to list of topics

Three Wavelets with Lens

In[62]:=

waveletSystem = {
    Move[LineOfWavelets[10,NumberOfRays->3,PaddingFactor->.1],-50],
    PlanoConvexLens[100,50,10],
    Move[Screen[50],103]};

In[63]:=

ShowSystem[waveletSystem,ShowGaussian->True,PlotType->TopView];

In[64]:=

Timing[intensity = FindField[waveletSystem]][[1]]/60/Second*"Minute"

Out[64]=

In[65]:=

Plot[intensity[0,y],{y,-.07,.07},PlotRange->All];

In[66]:=

Plot[intensity[-10,y],{y,-1,1},PlotRange->All];

In[67]:=

Plot3D[intensity[x,y],{x,-2,2},{y,-.1,.1},PlotRange->All,PlotPoints->100];

Go to list of topics

Wavelet analysis of an Imaging System

In[243]:=

In this section, we calculate the wavefront propagation for a plane wave traveling through a lens. Here is an overall schematic of the system.

In[68]:=

ShowSystem[{
    Move[LineOfWavelets[10,NumberOfRays->1],-50],
    PlanoConvexLens[100,50,10],
    Move[Screen[50],103]},PlotType->TopView,ShowGaussian->True];

In[243]:=

Next, we define the system using 21 Gaussian wavelets to approximate a collimated beam of light. Because each Gaussian wavelet follows the symbolically calculated solution for the system, relatively few numbers of wavelets are required to accurately model the system's behavior.

In[69]:=

waveletSystem = {
    Move[LineOfWavelets[10,NumberOfRays->21],-50],
    PlanoConvexLens[100,50,10],
    Move[Screen[50],103]};

In[243]:=

Finally, we use FindField to construct a symbolic solution of the field intensity. This calculation will take about 10 minutes on many computers.

In[70]:=

Timing[intensity = FindField[waveletSystem]][[1]]/60/Second*"Minutes"
ByteCount[intensity]

Out[70]=

Out[71]=

In[243]:=

Here is plot of the intensity near the plane of focus. The vertical axis is intensity. The horizontal axis is measured in millimeters.

In[72]:=

Plot[intensity[0,y],{y,-.022,.022},PlotRange->All,PlotPoints->100];

In[243]:=

The vertical axis is normalized to have a unity area. We can check this result with NIntegrate. In fact, you could also use the Integrate to do this symbolically, but the calculation time would be prohibitive.

In[73]:=

NIntegrate[intensity[0,y],{y,-.1,.1}]

Out[73]=

In[243]:=

We now compare this plot with the point spread function that is calculated numerically for the same surface location. Although not identical, the two results are consistent with each other.

In[74]:=

psf = PointSpreadFunction/.PSF[{
    Move[LineOfRays[10],-50],
    PlanoConvexLens[100,50,10],
    Move[Screen[50],103]},NumberOfPoints->256,PaddingFactor->10];

In[243]:=

NIntegrate is used again to check its area.

In[75]:=

NIntegrate[psf[y],{y,-.1,.1}]

Out[75]=

In[243]:=

Here is an symbolic intensity plot at a different plane near the focus.

In[76]:=

Plot[intensity[-1,y],{y,-.1,.1},PlotRange->All,PlotPoints->100];

In[243]:=

Again, the numeric point spread function is consistent in form. Note that the numeric PSF shows eveidence of high frequency ringing that is absent from the symbolic solution. This is due the hard edge of the pupil function in the numeric PSF.

In[77]:=

psf = PointSpreadFunction/.PSF[{
    Move[LineOfRays[10],-50],
    PlanoConvexLens[100,50,10],
    Move[Screen[50],102]},NumberOfPoints->512,PaddingFactor->10];

In[243]:=

Again we try a symbolic intensity plot but, this time we move further away from the plane of focus.

In[78]:=

Plot[intensity[-10,y],{y,-1,1},PlotRange->All,PlotPoints->100];

In[243]:=

The numeric PSF is again consistent, but with a great deal of high frequency ringing present.

In[79]:=

psf = PointSpreadFunction/.PSF[{
    Move[LineOfRays[10],-50],
    PlanoConvexLens[100,50,10],
    Move[Screen[50],93]},NumberOfPoints->4096,PaddingFactor->10];

In[243]:=

We can find the plane of peak intensity by making a symbolic plot along the optical axis of the system.

In[80]:=

Plot[intensity[x,0],{x,-10,10},PlotRange->All];

In[243]:=

We can now compare plots of the symbolic and numeric intensity profiles at the plane of peak intensity. Again, the overall shapes are in agreement with each other.

In[81]:=

Plot[intensity[.133,y],{y,-.022,.022},PlotRange->All,PlotPoints->100];

In[82]:=

psf = PointSpreadFunction/.PSF[{
    Move[LineOfRays[10],-50],
    PlanoConvexLens[100,50,10],
    Move[Screen[50],103.133]},NumberOfPoints->1024,PaddingFactor->10];

In[243]:=

Finally, we make a three-dimensional plot of symbolic intensity between the optical and transverse axes. The peak indicates the paraxial focus of the system.

In[83]:=

Plot3D[intensity[x,y],{x,-2,2},{y,-.1,.1},PlotRange->All,PlotPoints->100];

In[243]:=

Here is a similar intensity plot using the ContourPlot function. The horizontal dimension is measured along the optical axis of the system. The plot has been a thresholded in order to observe the low-intensity details.

In[84]:=

ContourPlot[intensity[x,y], {x,-1, 1}, {y,-.05,.05}, PlotPoints->100, PlotRange->{0,50}, Contours->30, ContourLines->False, ColorFunction->Function[Hue[.65-#*(.65),1,#*.9+.1]]];

In[243]:=

Here is the same plot with all intensity values included.

In[85]:=

ContourPlot[intensity[x,y], {x,-1, 1}, {y,-.05,.05}, PlotPoints->200, PlotRange->All, Contours->30, ContourLines->False, ColorFunction->Function[Hue[.65-#*(.65),1,#*.9+.1]]];

Go to list of topics

Created by Mathematica  (November 10, 2005)