6. Experimenting with LensLab

Introduction

Loading

6.1 Building a Black Hole

6.2 Using a Prism to Create Rainbows and More

6.3 Making an Angular Amplifier

6.4 Modeling a Cassegrain Telescope

Introduction to Chapter 6

In this chapter, we use LensLab for performing some advanced optical experiments. In Section 6.1, we construct a model of a Cassegrain telescope. Then, in Section 6.2, we use the Resonate genetic building block to construct a reflective optical tunnel. Later, in Section 6.3, we evaluate prism dispersion of white light. Finally in Section 6.4, we use refraction for building a novel angular amplifier.

Loading LensLab

Make sure that the LensLab package is located either in the home directory, or on a directory path recognized by Mathematica for packages. The LensLab package is named LensLab.m and located in the LensLab directory, and the LensLab package is loaded with the following expression.

In[1]:=

Needs["LensLab`LensLab`"]

This loading process should only take a few seconds. In addition to being loaded as a package, the LensLab.m file is formatted as a Mathematica notebook. The LensLab source code is made accessible so that you can develop new functions of your own by studying LensLab's built-in functions. This is particularly helpful when you wish to model new component ideas in LensLab. However, you should receive permission from Optica Software before distributing any user-created functions that are derived from the LensLab source code. The unauthorized distribution of LensLab-derived code may be a LensLab license agreement violation or a copyright infringement.

Go to list of topics

6.1 Building a "Black Hole"

This example demonstrates how you can use LensLab to blend traditional mathematical functions with optics to create new and interesting optical components. Here we define a new component that acts as a reflective tapered tube, similar in nature to an optical fiber. To give this new component a tube shape, we use the hyperbolic function.

In[2]:=

S[x_,y_] = 100/(x^2 + y^2);

A simple plot of 100/x^2 reveals its basic nature.

In[3]:=

Plot[100./x^2,{x,-5,5},PlotRange->{0,100}];

Note the singularity at x = 0. This singularity must be carefully treated when defining the new component in LensLab. Since our component will be three-dimensional, we next plot the three-dimensional hyperbolic function.

In[4]:=

Plot3D[S[x,y],{x,-5,5},{y,-5,5},
    PlotRange->{0,100}];

Having seen the basic nature of the hyperbolic function, we are nearly ready to fully define the component function. The pole in our function will be masked by putting a hole at its center. In addition, the option FunctionCenter->0 will be used to redefine the surface function at S[0,0] to be zero, instead of infinity. We give the new component function the name BlackHole for its geometric likeness to the gravitation field around a black hole in space. As discussed in Chapter 5, we define BlackHole by calling a series of genetic building blocks, each of which has a specific job in defining the behavior of the optical component. In this example we use Resonate to allow the same surface to be used for multiple interactions with the ray.

In addition to Resonate, we use CustomMirror for the basic mirror structure and Hole to put a hole in the component center. BlackHole takes arguments for the size of its mouth and the size of its hole at the center.

In[5]:=

Clear[BlackHole];

BlackHole[aperture_,holeaperture_,opts___] :=
Block[{options},
    options = Flatten[{opts,Options[CustomMirror]}];
    Hole[
        Resonate[
            CustomMirror[
                Function[100./(#1^2 + #2^2)],
                aperture,
                "BlackHole",
                SurfaceLabel -> OtherShape,
                FunctionCenter -> 0.,
                options],
            "BlackHole"
        ],holeaperture,
    options]];

Notice that we use the CustomMirror options to define BlackHole. This eliminates the need to create a special options listing for BlackHole. First let's just make a three-dimensional plot of BlackHole.

In[7]:=

DrawSystem[BlackHole[5,2.5,GraphicDesign->Wire],
    Boxed->False];

Next we view the mouth of BlackHole.

In[8]:=

ShowSystem[%,PlotType->FrontView];

Finally we test out the effects of a single ray sent into the enlarged mouth of BlackHole.

In[9]:=

sys =
DrawSystem[{
    Move[Ray[],15,30],
    BlackHole[5,2.5,GraphicDesign->Sketch],
    Boundary[{0,-100,-100},{100,100,100}]},
    PlotType->TopView];

For a better view, the image might need to be expanded using the mouse cursor to drag the corner of the graphics frame. We see that the ray bounced around inside the tube, and did not leave the other side. Instead, it seems to have turned around and bounced back out the entrance. We next see the process in better detail by creating an animation. The basic Mathematica code used below is typical for making LensLab animations. Be forewarned that limited computer memories may become taxed by the number of images required for animations.

To see how many frames we need, we check the number of Ray objects in the system.

In[10]:=

Length[sys[[1]]]

Out[10]=

We use RayChoice->{IntersectionNumber->i} to display a different ray segment in each frame.

In[11]:=

Do[ShowSystem[
        sys,
        RayChoice->{IntersectionNumber->i},
        PlotType->TopView,
        PlotRange->{{-1,65},{-3,14}}],{i,1,63}];

In the previous ray tracing, we sent the ray into the component defined by BlackHole with an initial angle of 30 degrees. Let us now use a wedge of rays having angles shallower then 30 degrees. For more robust ray-tracing of some custom surfaces, you can use the SurfaceRayIntersections -> Solve option setting. In this case, the Solve function is used to find the symbolic solution of the surface intersection point while the default setting of SurfaceRayIntersections -> Automatic uses the FindRoot function. FindRoot can sometimes fail because sometime finds only a local minimum and may miss the global minimum for the surface intersection point.

In[12]:=

sys =
DrawSystem[
    {Move[WedgeOfRays[30, NumberOfRays->6],{15,0,0}],
    BlackHole[5,2,GraphicDesign->Sketch,SurfaceRayIntersections->Solve],
    Boundary[{0,-100,-100},{120,100,100}]},
PlotType->TopView];

This time all of the rays made it through the tube without turning around. Next we create an animation for this system. First we must determine how many ray bounces occurred inside the tube. This time the number of frames is not simply given by the number of Ray objects divided by the number of rays. Instead, we use the following expression to check how many ray bounces occurred inside the tube.

In[13]:=

Length[Union[ReadRays[sys,IntersectionNumber]]]

Out[13]=

Again, we use RayChoice->{IntersectionNumber->i} to display different ray segments in each frame.

In[14]:=

Do[ShowSystem[
        sys,
        RayChoice->{IntersectionNumber->i},
        PlotType->TopView,
        PlotRange->{{0,120},{-17,17}}], {i,1,21}];

Note that the shallower angles exit more quickly than the steeper angles. Next we use ReadRays to measure the total optical path length of the different rays.

In[15]:=

ReadRays[sys, OpticalLength, ComponentNumber -> 2]

Out[15]=

We also measure the difference in optical path length using the built-in Max and Min functions of Mathematica.

In[16]:=

(Max[#]-Min[#])&[%]

Out[16]=

This example only touches on the possibilities of building and analyzing novel optical components with LensLab.

Go to list of topics

6.2 Using a Prism to Create Rainbows and More

LensLab can be used to show the dispersive effects of prisms. Here we will do a simple experiment by sending a ray of "white" light through a typical prism (BK-7 glass). Our ray is actually six rays having wavelength distributions equally spaced between 400 nm and 700 nm. Since the index of refraction is a nonlinear function of wavelength, we can see the relationship between spatial distribution and wavelength on a flat optical surface intersecting with the rays.

In[17]:=

DrawSystem[{
    RainbowOfRays[{.4,.7}, NumberOfRays->6],
    Move[Prism[{60,135,60},100],{100,-50},-15],
    Boundary[{0,-200,-200},{300,200,200}]}];

In[18]:=

ShowSystem[%,PlotType->TopView];

Next we look at where the rays intersect with the boundary.

In[19]:=

ShowSystem[%,PlotType->Surface,
    RayChoice->{ComponentNumber->2}];

Unfortunately, in this picture, the spot size of the rays is too small to see the relationship between wavelength and position. Let's put a screen in the system that has an aperture size slightly larger than the spot size. We can determine the appropriate screen dimensions and placement coordinates by using the cursor to take measurements on the above two pictures. This is done by first pointing to an image, pressing the mouse button, and at the same time pressing the Command key. We see that the intersection spot ranges from -139 to -155 in the horizontal. Using measurements from PlotType->TopView shown above, we place the screen at {296, -145}.

In[20]:=

DrawSystem[{
    RainbowOfRays[{.4,.7}, NumberOfRays->6],
    Move[Prism[{60,135,60},100],{100,-50},-15],
    Move[Screen[155.-139.],{296,-145}],
    Boundary[{0,-200,-200},{300,200,200}]},
    PlotType->Surface,RayChoice->{ComponentNumber->2}];

We now see clearly how the positions are a nonlinear function of wavelength (which is why diffraction gratings are preferred over prisms for making careful spectroscopic measurements).

We next use Mathematica to make a plot of the horizontal component of the intersection points, and then fit a curve to these data points. First, let's trace through the system using 16 rays instead of 6.

In[21]:=

colorsys = DrawSystem[{
    RainbowOfRays[{.4,.7},NumberOfRays->16],
    Move[Prism[{60,135,60},100],{100,-50},-15],
    Move[Screen[155.-139.],{296,-145}],
    Boundary[{0,-200,-200},{300,200,200}]},
    PlotType->Off];

We then use ReadRays to read out ray parameter values for wavelength and position at the ray/screen intersection. First, let's check out again how to use ReadRays. To do this, we use the built-in help function, ?ReadRays.

In[22]:=

?ReadRays

Next we must extract the horizontal component of each ray/screen intersection, creating a list of numbers for plotting. Before we do this, we need to review Options[Ray] to learn the parameter name for surface intersections.

In[23]:=

Options[Ray]

Out[23]=

SurfaceCoordinates seems to be the right parameter we want. To be sure, we again use the on-line help function ?SurfaceCoordinates.

In[24]:=

?SurfaceCoordinates

Finally, we use ReadRays with SurfaceCoordinates to extract the horizontal values. These numbers are measurements of distance along the screen object.

In[25]:=

horizontalValues =
    ReadRays[colorsys,SurfaceCoordinates, ComponentNumber->2]

Out[25]=

Since we are only concerned with the first surface coordinate, we use Map with First to extract the first values.

In[26]:=

horizontalValues = Map[First[#]&,horizontalValues]

Out[26]=

We make the numbers start from 0.0 by subtracting the minimum offset from these values.

In[27]:=

horizontalValues = horizontalValues-Min[horizontalValues]

Out[27]=

Now we use ListPlot to generate a plot of these values.

In[28]:=

ListPlot[horizontalValues];

Next, we create a second list, containing the wavelength of each intersecting ray. By multiplying by 1000, we get units of nanometer.

In[29]:=

wavelengths =
    ReadRays[colorsys,WaveLength,ComponentNumber->2]1000

Out[29]=

Finally, we combine the two lists to make ordered pairs of numbers.

In[30]:=

displacements = Transpose[{wavelengths,horizontalValues}]

Out[30]=

Now we make a plot containing intersection distance as a function of wavelength. We make a solid curved line using PlotJoined->True.

In[31]:=

measuredPlot =
ListPlot[displacements,PlotJoined->True];

We now take things a step further by using the Fit function to come up with a calibrated function for our physical system. Let's first check out how to use Fit.

In[32]:=

?Fit

Now we do a quadratic fit.

In[33]:=

Fit[displacements,{1,x,x^2},x]

Out[33]=

Now we plot the fitted function.

In[34]:=

fittedPlot = Plot[%,{x,400,700}];

Finally, let's compare this graph with the graphed values "measured" using LensLab.

In[35]:=

Show[measuredPlot,fittedPlot];

We see that our fitted function doesn't quite match the measured values. If we wanted to, we could go back and try a different fitting function, perhaps involving a sin term or a higher-order polynomial, until a better match was found.

Go to list of topics

6.3 Making an Angular Amplifier

Using LensLab, it is possible to quickly explore novel ideas. This example illustrates a possible method for making an angular amplifier that takes an input signal containing a small angular deviation, and produces an output signal having an amplified angular deviation. Using the properties of Snell's law, this amplifier works by sending a beam of light through a glass/air interface such that the input beam, coming from the high refractive index side, hits the interface very near the critical angle for the interface. The result is that the input angle gets amplified by the interface.

In[36]:=

DrawSystem[{
    Move[WedgeOfRays[2, OpticalMedium -> BK7, NumberOfRays->6],{0,-25},40],
    Move[LensSurface[{100,100}],{50,0,0}],
    Boundary[{-100,-100,-50},{150,100,50}]},
    PlotType->TopView];

The previous diagram demonstrates the basic system. We characterize the system more accurately by creating an input/output model of the system behavior. First of all, let's take more data points by increasing the number of rays used in our example. This time, we won't make a plot, but instead use PropagateSystem and store the results in system.

In[37]:=

system = PropagateSystem[{
    Move[WedgeOfRays[2, OpticalMedium -> BK7, NumberOfRays->16],{0,-25},40],
    Move[LensSurface[{100,100}],{50,0,0}],
    Boundary[{-100,-100,-50},{150,100,50}]}];

We must extract from system only the input rays and output rays of the system. This is done using RaySelect.

In[38]:=

inputrays = RaySelect[system,IntersectionNumber->1]

Out[38]=

In[39]:=

outputrays = RaySelect[system,IntersectionNumber->2]

Out[39]=

Now we will extract the angle information from the input rays and output rays. For this we use the ray parameter RayTilt.

We extract the ray object angle by using Map and ReplaceAll, indicated by /. in the following expression. In addition ArcTan is employed for converting the direction vector into an angle. We store this result in outputangles.

In[40]:=

outputangles = Map[(N[ArcTan[(RayTilt/.#)[[2]]/(RayTilt/.#)[[1]]]/Degree])&, outputrays]

Out[40]=

Note that we could get the same result by using ReadRays and taking more steps. Since we are interested in relative angular displacements, rather than absolute displacements, we now subtract off the angular minimum using Min.

In[41]:=

outputangles = outputangles - Min[outputangles]

Out[41]=

We now plot the angular displacements.

In[42]:=

ListPlot[outputangles];

Similarly, we extract the input angular variation, storing the results in inputangles.

In[43]:=

inputangles = Map[(N[ArcTan[(RayTilt/.#)[[2]]/(RayTilt/.#)[[1]]]/Degree])&, inputrays]

Out[43]=

Again, we subtract off the angular minimum.

In[44]:=

inputangles = inputangles - Min[inputangles]

Out[44]=

Next we combine the input and output variables into a list of ordered pairs.

In[45]:=

transfercurve = Transpose[{inputangles,outputangles}]

Out[45]=

We now plot the output angle as a function of the input angle.

In[46]:=

measuredPlot = ListPlot[transfercurve,PlotJoined->True];

Next we create a model of this system. Initially, we apply a linear fit to the data points.

In[47]:=

fittedfunction = Fit[transfercurve,{1,x},x]

Out[47]=

Here we see that our amplifier has a gain of nearly 6.

Now we plot the fitted function.

In[48]:=

fittedPlot = Plot[fittedfunction,{x,0,2}];

Finally, let's compare this graph with the "measured" results.

In[49]:=

Show[measuredPlot,fittedPlot];

Clearly, our system is not terribly linear. Next we calculate and plot the error between the measured results and the fitted results. First we create a list of numbers using our fitted function.

In[50]:=

fittedyvalues = Table[fittedfunction,{x,0,2,2./15}]

Out[50]=

Next, we use ListPlot to plot the difference between the fitted and actual values.

In[51]:=

ListPlot[Transpose[{inputangles,
    (fittedyvalues - outputangles)}],PlotJoined->True];

As a further refinement to this example, we make a new Mathematica function that takes any optical system and returns a fitted function to describe the system gain. We will call this function ModelGain. This function is constructed by piecing together all of the steps performed separately above.

In[52]:=

ModelGain[opticalsystem_, outputIntersectionNumber_,
    fittingfunction_, options___]:=
Block[{system,inputrays,outputrays,outputangles,
    inputangles,fittedfunction,    transfercurve,inputrange,
    fittedplot,measuredplot},
(* Propagate optical system, if system has not already
been propagated *)
system = PropagateSystem[opticalsystem];
inputrays = RaySelect[system,IntersectionNumber->1];
outputrays = RaySelect[system,
    IntersectionNumber->outputIntersectionNumber];
outputangles = Map[(N[ArcTan[(RayTilt/.#)[[2]]/
        (RayTilt/.#)[[1]]]/Degree])&,outputrays];
outputangles = outputangles - Min[outputangles];
inputangles = Map[(N[ArcTan[(RayTilt/.#)[[2]]/
            (RayTilt/.#)[[1]]]/Degree])&,inputrays];
inputangles = inputangles - Min[inputangles];
inputrange = Max[inputangles];
transfercurve = Transpose[{inputangles,outputangles}];
(* Now plot measured parameters *)
measuredplot = ListPlot[transfercurve,PlotJoined->True];
fittedfunction = Fit[transfercurve,fittingfunction,x];
(* Now plot fitted parameters *)
fittedPlot = Plot[fittedfunction,{x,0,inputrange}];
(* Show both plots together *)
Show[measuredplot,fittedPlot];
(* Plot deviation between measured and fitted values *)
ListPlot[Transpose[{inputangles,(Table[fittedfunction,
{x,0,inputrange,N[inputrange/(Length[inputangles]-1)]}]-
outputangles)}],PlotJoined->True];
(* Finally, output the fitted function *)
fittedfunction];

Using ModelGain, we change parameters and then quickly test the system. For now, let's continue to use the same system and find a quadratic fit model for the system.

In[53]:=

ModelGain[system,2,{x,x^2}]

Out[53]=

Notice that our quadratic amplifier model more closely follows the measured parameters. In this example, we have examined a system behaving as an angular amplifier. In addition, we have seen how new Mathematica functions are quickly created for performing intricate tasks. We could conceivably build a multistage amplifier for higher, more linear gain.

Go to list of topics

6.4 Modeling a Cassegrain Telescope

In this section, we model a Cassegrain telescope with LensLab, including the primary and secondary mirror objectives. For this example, we combine several component functions into a special function that we call CassegrainTelescope. Within CassegrainTelescope, we use CylinderGraphic to create the telescope housing shroud. As inputs to CassegrainTelescope, we use the f/number of the system, the diameter of the primary mirror, the diameter of the secondary mirror, and the desired position in space of the eyepiece focal point. We now define CassegrainTelescope.

In[54]:=

Clear[CassegrainTelescope];

inch = 25.4;

In[56]:=

CassegrainTelescope[fnumber_,bigdiameter_,
    smalldiameter_,eyepieceposition_,options___]:=
Block[{bigmirrorthickness,smallmirrorthickness,focallength1,
    focallength2,gamma,position,holediameter},
bigmirrorthickness = bigdiameter/10;
smallmirrorthickness = smalldiameter/10;
focallength1 = fnumber bigdiameter;
position = N[smalldiameter focallength1 / bigdiameter];
gamma = N[(focallength1 - position + eyepieceposition +
    bigmirrorthickness)/position];
focallength2 = -N[position gamma/(gamma - 1)];
holediameter = N[.75 smalldiameter];
{Move[ParabolicMirrorWithHole[focallength1, bigdiameter,holediameter,bigmirrorthickness, options],-focallength1,180],
Move[ParabolicMirror[focallength2,smalldiameter, smallmirrorthickness,options], -position],
Move[CylinderGraphic[bigdiameter,-N[.1 bigdiameter+.1 smalldiameter+focallength1-position],options], {-position + .05 bigdiameter,0,0}],
Boundary[{-N[1.25(focallength1+eyepieceposition)], -bigdiameter/2,-bigdiameter/2}, {1.5eyepieceposition,bigdiameter/2,bigdiameter/2}]}];

Next, we use DrawSystem to trace a circle of rays through CassegrainTelescope.

In[57]:=

DrawSystem[
Move[{Move[CircleOfRays[15 inch, NumberOfRays->6],15 inch,180],
    CassegrainTelescope[2.5,20 inch,4 inch,10 inch, OpenSide->{-.25 Pi,1.3 Pi}, GraphicDesign->Solid]},
    {0,0,90}],Boxed->False];

Finally, we use PlotType->Surface with RayChoice->{IntersectionNumber->2} to examine the ray intersection points on the secondary mirror surface.

In[58]:=

ShowSystem[%,PlotType->Surface, RayChoice->{IntersectionNumber->2}];

In Chapter 7, you can learn more about designing optical systems using LensLab and Mathematica.

Go to list of topics

Copyright Statement

LensLab is a trademark of Optica Software.
Mathematica ® is a registered trademark of Wolfram Research, Inc. All other product names mentioned are trademarks of their producers. Mathematica is not associated with Mathematica Policy Research, Inc. or MathTech, Inc.
Copyright ©1995-2005 by Optica Software, Urbana, Illinois, Champaign, Illinois.

All rights reserved. No part of this document may be reproduced, stored in a retrieval system, or transmitted in any form or by any means, electronic, mechanical, photocopying, recording or otherwise, without prior written permission of the author, Optica Software.

Optica Software is the holder of the copyright to the LensLab package software and documentation ("Product") described in this document, including without limitation such aspects of the Product as its code, structure, sequence, organization, "look and feel", programming language and compilation of command names. Use of the Product, unless pursuant to the terms of a license granted by Optica Software. or as otherwise authorized by law, is an infringement of the copyright.

Optica Software makes no representations, express or implied, with respect to this Product, including without limitations, any implied warranties of merchantability or fitness for a particular purpose, all of which are expressly disclaimed. Users should be aware that included in the terms and conditions under which Optica Software is willing to license the Product is a provision that the author, Optica Software and distribution licensees, distributors and dealers shall in no event be liable for any indirect, incidental or consequential damages, and that liability for direct damages shall be limited to the amount of the purchase price paid for the Product.

In addition to the foregoing, users should recognize that all complex software systems and their documentation contain errors and omissions. Optica Software shall not be responsible under any circumstances for providing information on or corrections to errors and omissions discovered at any time in this document or the package software it describes, whether or not they are aware of the errors or omissions. Optica Software does not recommend the use of the software described in this document for applications in which errors or omissions could threaten life, injury or significant loss.

Go to list of topics

Created by Mathematica  (November 4, 2005)