Examples (2D and 3D)

Vectorizing over bundle angles

The Quickstart tutorial showed tracing a single bundle through a system. Now suppose there is a need to trace multiple types of bundles through a system all at once. Consider a system of two parabolic surfaces which have bundles incident at different angles. One could simply do the procedure in the tutorial three different ways, but this would be somewhat verbose. Instead this can be done more efficiently by defining functions over the variables which are of interest and then calling them on vectors of the different variable values.

Below a setup is implemented very similarly to the tutorial, but a function bund_ang(angle) is defined to originate a bundle of fixed spatial origin but variable y-direction-tangent angle.

using GeometricalOptics, Plots

surfslist = [aspherical, aspherical, zplane]
coeffslist = [[5., -6., -1.], [1., 2., -1.], [8.]]      # two parabolic surfaces
nlist = [1., -1., 1.]                               # index sign change means reflection

optsta_instance = opticalstack(coeffslist, surfslist, nlist)

bund_ang(angle) = bundle([0.], (-0.5:0.5:0.5).+1, 0., angle, 0.)

Then vectors to vary for the system are defined. Three different angles are of interest, and for visiblity it is best to plot in different colors and only plot surfaces along with a single set of rays (as the system does not change otherwise). Then a Plots.jl plot pane is initialized so that it can be passed into all of the vectorized cases of analysis and they can be plotted together.

angles = [0., 0.028, 0.04]                    # a list of three angles
colors = [:red, :green, :blue]                # each of which is colored differently
issurfplotvec = [true, false, false]          # only plot surfaces along with first bundle
plot_together_in = Plots.plot()               # must put this into the system plotter to get all angles.

Finally a function stack_render is defined to accept all of the varying parameters, and then the vectors are all passed into the function to return a vector of plots. The vector of plots are linked to the same original Plots.plot instance, so calling the last element of the vector of plots will show all three instances together.

stack_render(angle, pc, isp) = optsta_instance(bund_ang(angle); rend = "YZ", color = pc, numsurfpoints = 100, plobj = plot_together_in, issurfplot = isp)

plot_together_out_vec = stack_render.(angles, colors, issurfplotvec)
plot_together_out_vec[3]                    # the last entry in the vector of plots has all bundles

And an analogous procedure generates all of the spot plots together. A function is defined to set a now square bundle for a variable y-angle, and then another function plots all of the spots together on the final plane.

splot_init = Plots.plot()

bund2_ang(angle) = bundle((-0.5:0.2:0.5), (-0.5:0.2:0.5).+1, 0., angle, 0.)

function optsta_render_splot(angle, pc, isp)
    traceout = optsta_instance(bund2_ang(angle))
    splot(traceout, pl_splot = splot_init, color = pc)
end

splot_vec = optsta_render_splot.(angles, colors, issurfplotvec)
splot_vec[3]                                # the last entry in the vector of plots has all bundles

Advantageous start conditions and fixing a broken raytrace

Broadly, there are few ways to break the raytrace analysis. One way is to attempt to violate physics or basic geometry. Another way is due to extremely steep surface or ray angles with respect to z. This should be tested. Still another ays is to not provide proper start conditions for numerical raytrace methods. There is a utility to address this latter case.

A few other things should be tested:

• Rays must be known to converge on the surface.
  • Don't choose a function that rolls off before the ray can hit it; that cannot be fixed.
• The surface function must be designed to return a number at the ray incident position.
  • Example: built-in zernike is defined to return NaN when outside the normalization radius, so it is designed to fail in certain cases.

Using a modification of the above example, it is possible to "break" the raytrace. Take the same system for analysis as above, but change the parameters of the first asphere to be positioned and focused differently. Note that a few rays traced through the first surface are traced erratically. This is due to improper convergence of numerical raytracing.

using GeometricalOptics, Plots

surfslist = [aspherical, aspherical, zplane]
coeffslist = [[4., -4., -1.], [1., 2., -1.], [8.]]      # two parabolic surfaces
nlist = [1., -1., 1.]                               # index sign change means reflection

optsta_instance=opticalstack(coeffslist, surfslist, nlist)

bund_ang(angle) = bundle([0.], (-0.5:0.5:0.5).+1, 0., angle, 0.)

angles = [0., 0.028, 0.04]                     # a list of three angles
colors = [:red, :green, :blue]                 # each of which is colored differently
issurfplotvec = [true, false, false]           # only plot surfaces along with first bundle
plot_together_in = Plots.plot()                # must put this into the system plotter to get all angles.

stack_render(angle,pc,isp) = optsta_instance(bund_ang(angle); rend = "YZ",color = pc,ydom = -2:0.1:2, plobj = plot_together_in, issurfplot = isp)

plot_together_out_vec = stack_render.(angles, colors, issurfplotvec)
plot_together_out_vec[3]                    # the last entry in the vector of plots has all bundles
title!("broke the trace :(")

The spots show the erratic rays as well:

splot_init = Plots.plot()

bund2_ang(angle) = bundle((-0.5:0.2:0.5), (-0.5:0.2:0.5).+1, 0., angle, 0.)

function optsta_render_splot(angle, pc, isp)
    traceout = optsta_instance(bund2_ang(angle))
    splot(traceout, pl_splot = splot_init, color = pc)
end

splot_vec = optsta_render_splot.(angles, colors, issurfplotvec)
splot_vec[3]                                # the last entry in the vector of plots has all bundles
title!("broke the trace :(")

This case can be fixed by using the optional argument init_prop_guesses to an instance of opticalstack. The optional argument can be passed in as a vector of same length as the number of surfaces, or as nothing. When nothing, all rays are first guessed to propagate zero distance to the next surface. When passed in as a vector, each value of the vector is the initial guess of signed propagation distance to the corresponding surface. Because the error arises in finding intersection of the rays with the second surface, the initial guess for the second surface is altered, and in particular made negative (due to reflection) and given magnitude close to the distance expected for the ray. Then all of the rays in a bundle will be traced separately, only with newly specified starting guesses for distances traversed between surfaces.

The issue with the raytrace is rectified by adding the vector of values to provide a better initial guess for the negative propagation to the second surface:

using GeometricalOptics, Plots

surfslist = [aspherical, aspherical, zplane]
coeffslist = [[4., -4., -1.], [1., 2., -1.], [8.]]      # two parabolic surfaces
nlist = [1., -1., 1.]                               # index sign change means reflection

optsta_instance = opticalstack(coeffslist, surfslist, nlist)

bund_ang(angle) = bundle([0.], (-0.5:0.5:0.5).+1, 0., angle, 0.)

angles = [0., 0.028, 0.04]                      # a list of three angles
colors = [:red, :green, :blue]                 # each of which is colored differently
issurfplotvec = [true, false, false]            # only plot surfaces along with first bundle
plot_together_in = Plots.plot()               # must put this into the system plotter to get all angles.

init_prop_guesses_choice = [0., -2., 0.]

stack_render(angle, pc, isp) = optsta_instance(bund_ang(angle); rend = "YZ", color = pc, ydom = -2:0.1:2, plobj = plot_together_in, issurfplot = isp, init_prop_guesses = init_prop_guesses_choice)

plot_together_out_vec = stack_render.(angles, colors, issurfplotvec)
plot_together_out_vec[3]                    # the last entry in the vector of plots has all bundles
title!("fixed the trace :)")

The spots behave now too:

splot_init = Plots.plot()

bund2_ang(angle) = bundle((-0.5:0.2:0.5), (-0.5:0.2:0.5).+1, 0., angle, 0.)

function optsta_render_splot(angle, pc, isp)
    traceout = optsta_instance(bund2_ang(angle); init_prop_guesses = init_prop_guesses_choice)
    splot(traceout, pl_splot = splot_init, color = pc)
end

splot_vec = optsta_render_splot.(angles, colors, issurfplotvec)
splot_vec[3]                                # the last entry in the vector of plots has all bundles
title!("fixed the trace :)")

Vectorizing over index, 3D analysis

The same sort of thing that was done for different bundle angles can also be done for different system configurations. Let's examine a few variations on this, plotting in 3D instead of 2D and using a custom bundle array defined to be a circular bundle footprint. First define the parameters of the system that are fixed. The bundle is fixed with a circular footprint defined by bundle_as_array. The array arguments are defined such that one array subscript corresponds to a radial distance from the central axis, and another subscript corresponds to an azimuthal angle around the axis.

using GeometricalOptics, Plots

surfslist = [aspherical, aspherical, zplane]
coeffslist = [[1., 6., -5.], [2., -6., -5.], [8.]]
# index list is no longer up here, as the material is varied in modeling.

rdom = 0.01:0.2: 1.
thdom = 0.:pi/10:2pi
r = rdom' .* ones(length(thdom))
th = ones(length(rdom))' .* thdom
bund_instance = bundle_as_array(r .* cos.(th), r .* sin.(th), zeros(size(r)), zeros(size(r)), 0.)

Then define the variable system parameters. The selected index values for the second material are in the vector nsel, and all three of these will be tested; results will be plotted in different colors on the same plot with surfaces only plotted alongside the first case.

nsel = [1.46, 1.5, 1.56]                # a list of three indices,
colors = [:red, :green, :blue]         # each of which is colored differently
issurfplotvec = [true, false, false]    # only plot surfaces with one of the bundles.
plot_lens_init = Plots.plot()       # must put this into the system plotter to get all cases.

Now make a function to vary the system over the desired values, and then call over the vectors defined above:

function optsta_render(n_u, pc, isp)
    nlist_chr = [1., n_u, 1.]
    optsta_instance = opticalstack(coeffslist, surfslist, nlist_chr)
    optsta_instance(bund_instance; rend = "3Dcirc", color = pc, numsurfpoints = 100, plobj = plot_lens_init, issurfplot = isp)
end

plot_lens_vec = optsta_render.(nsel, colors, issurfplotvec)
plot_lens_vec[3]                # only need to plot the last of the three elements in the vector; it caught all cases.