Differential optics for illumination design in the presence of caustics

Alternatively to using time-consuming Monte-Carlo simulations the irradiation at a target plane can also be calculated by differential optics methods. In the case of caustics, these methods yield to infinite irradiance and its results are not directly comparable to those of MonteCarlo simulations. In this paper, a differential based algorithm for an on axis point source and a rotationally symmetric optical system is presented, which yields to the same results as a Monte-Carlo simulation. However, the differential optical methods are about three orders of magnitude faster than the latter one, thus allowing fast trial and error design of such kind of illumination systems. An applet is presented that uses sliders to change the shape of the lens and other properties of the illumination system whereas the irradiance profile is nearly immediately perceived. For beginners in the field, this does not only accelerate the design process itself but also the learning process is improved considerably. Some extensions and special cases are shortly discussed. [DOI: http://dx.doi.org/10.2971/jeos.2012.12018]


INTRODUCTION
Design of illumination systems incorporates regularly Monte-Carlo simulations tracing of the order of a million of optical rays.This is a robust but time consuming method to evaluate the irradiance profile at a target plane.Even for simple systems consisting of a point source and a single lens simulation times of the order of one minute have to be accepted.Targeting for a specified irradiance distribution, an approach by trail and error is therefore impeding for beginners.Additionally, for these simple illumination systems of standard optical elements there is, in contradiction to imaging optical systems, a lack of textbooks which allow to get a feeling of the effect of such systems.In this latter field, due to the fact that only of the order of hundred rays are sufficient to evaluate a typical target function like the geometrical MTF or similar measures, trial-and-error methods (called here in reference to [2] 'design by sliders') can be efficient to get a feeling of the problem.
Alternative ways to the above mentioned Monte-Carlo method of calculating the irradiance at a target plane, which will be presented here, are based a) on tracing local wavefronts associated with each ray and b) calculating analytically the local ray density at the target plane.Both methods use differentials as its base and called here differential optics methods.The advantages of these methods are that the irradiance can be calculated at much higher spatial density, with much lower computing effort and without any statistical noise compared to Monte-Carlo methods.However, at caustics, where the wavefront is folded, both differential optics methods yield to infinite irradiance.This unrealistic result has led to the fact that this methods are considered inappropriate for most of the illumination designers who deal with systems where caustics are not avoidable from the scratch.Only very few specialists have reported of successful use of the first mentioned methods in real world problems (where caustics could be avoided) [5].
Infinite irradiance at caustics is indeed computationally correct in the case where a point source with infinite radiance (but finite intensity) is used and diffraction is neglected.In Monte Carlo methods, based on the same assumptions, this outcome cannot be observed because always a receiving element with finite area (called bin) is used and thus this local effect is averaged out.In order to make the outcome of differential optical methods comparable to Monte Carlo methods, the same averaging at the receiver has to be realized.In this paper an appropriate algorithm is proposed and its convergence is discussed.Finally, an application based on these methods is presented that features design by sliders for a illumination system consisting of single lens, where the shape factor of the lens (and other properties) can be changed by the user in 'realtime' while observing the irradiance profile output.It is further pointed out that the methods are also powerful for the analysis and design of advanced illumination systems.

CALCULATING THE IRRADIANCE FOR SPHERICAL SYSTEMS BY THE DIFFERENTIAL Q-U-METHOD
Spherical surfaces are due to manufacturing and metrology reasons the standard geometry of optical elements.Most ofthe-shelf lenses are of this type.The focusing properties in terms of spot-diagrams of such lenses for point source objects are teached in nearly all optic courses and a variety of online accessible applets are available to get insight by trial and error.
In simple cases, design by sliders is faster than a systematic approach by aberration theory.Also for combinations of such lenses, e.g.triplets, such an approach has been proposed [1], [2].The author of this paper suggests that a comparable introductionary subject in illumination design is the analysis of the irradiance profile of a point source collimated by a spherical lens (and later on by a combination of these).However, to the knowledge of the author, this subject is not treated by textbooks and similar applets like in imaging applications are not available.The author himself has teached students and other newcomers in illumination design by using educational licenses of complex commercial software based on Monte Carlo methods from the scratch.This kind of teaching, supported by standard illumination ray tracing, is threatened by a number of facts, one of them is that students get disappointed by the high latency between input (e.g.change of the lens parameters) and output, i.e. the irradiance profile.In this chapter a method is presented that is suitable for rapid calculation of the irradiance profile of spherical systems illuminated by a point source, thus allowing trial-and-error methods and design by sliders for this subject.
The so called Q-U-method is an analytical method for tracing meridional rays.This known method is extended here by calculating the ray density based on differentials and from that in a straight forward way the irradiance.We repeat here shortly the Q-U-method and add the differential with respect to the starting ray angle ϕ: The starting equations for a meridional ray, emerging on axis with an angle ϕ 1 at distance s from the first surface vertex are: From these quantities, the angle of incidence is calculated by using the curvature c i of the i th surface: Refraction is described by the well know equation From this quantity, the ray angle after the surface is also straight forward to calculate: The parameter Q is now updated as follows: The derivative dQ i dϕ is again derived by standard math.Because the result is a bit lengthy, it is omitted here.The transfer from surface i to the next surface i + 1 at vertex distance d i+1 reads as: For the final height of the ray at the last surface k one gets Again, the derivative dy k dϕ and dz k dϕ are easy to derive but are a bit lengthy an omitted here.The height of the ray at some distance d behind the last surface is finally derived as follows: where all the quantities (except d ) are functions of ϕ It is now straight forward, however a bit lengthy, how to coagulate the ray density dh dϕ successively.
In the case of a single lens (k = 2) this differential can even be calculated by symbolic math software in close form.The author has done this with the Symbolic Math Toolbox of MAT-LAB®.From the symbolic expression, an executable function can be automatically generated, which has in this case about 150 lines of code.From this differential and the given intensity I (ϕ) of the point source the irradiance at the impact point h (ϕ) of the ray (with emerging ray angle ϕ) is finally calcu- lated by Together with Eq. ( 8) this results in a parametric description of point pairs (h (ϕ) , E (ϕ)).
In the limit ϕ → 0 (on axis) the above equation is not directly applicable.In this limit the ray height yields with ∆s = a + f = s − s H + f (the axial offset of the point source with respect to the focal point) and the locations of the principal planes s H , s H with respect to the surface vertices of the system: From this result one gets the irradiance in the limit ϕ → 0 (please note that this is not necessarily identical with the on axis irradiance due to other contributions, see below): Consider now the following simple example: a lambertian point source is located at the focal point of a plano-convex lens of focal length of 100 mm with thickness of 10 mm, a diameter of 50 mm and a refraction index of 1,5.The irradiance profile at a distance of 2000 mm from the lens is seeked (or just consider a scaled version of this example).An application of such a scaled version would be the starting design (based on a point source) of a spot light with a (obviously truncated) lambertian LED and a single spherical lens.A simple ray plot of the situation is shown in Figure 1.
The corresponding result for the irradiance based on 500 rays is displayed in Figure 2. From this figure the typical problems associated with this method described up to here are obvious: • The irradiance at the caustic approaches infinity (only due to the finite ray resolution, one gets here finite values).The peak is however quite narrow.
• At the caustic the wavefront is 'folded' inwards, as a consequence two irradiance values at the same height below the caustic are observable.
• The folded wavefront extends to negative height values which yield to a third irradiance value at positive height from the rays emitting downwards from the source (not shown in Figure 1).
• The density of irradiance values changes considerably with the ray height at the detector plane (better observable with less rays).Result of a Monte-Carlo simulation for comparison in blue.

YIELDING COMPARABLE RESULTS TO MONTE-CARLO SIMULATIONS
The big difference to the result of a Monte-Carlo-simulation in the example above may lead to the consequence, that the proposed approach is not suitable.To convince users, the method should yield comparable results as vastly used Monte-Carlos simulations.In these latter ones, due to the unavoidable statistical noise, detector bins have a relative large size in the range of a 10th to a 100th of the detector size.The bins are typically (but not necessarily) equidistant.The irradiance output of Monte-Carlos simulations is an average value over the bin area.
In order to yield to these requested comparable results, the non-equidistant irradiance of Eq. ( 8) has to be averaged over the same bin area used in a Monte-Carlo simulation.In the algorithm used in this paper, the first step is to interpolate the irradiance value at the lower and upper bound h bin,lb , h bin,ub of the bin.We than use trapezoidal integration to calculate the flux in the bin (to be more precise: the flux in the ring with the same cross section as the bin) and average it over the corresponding bin area h bin = h bin,ub + h bin,lb 2 .
Care has however to be taken at the caustic, which is easily detectable from the results of Q-U-differentials (a sign change in dh (ϕ) dϕ.Trapezoidal integration including the irradi- ance value closest to the caustic can yield in a few cases, depending on the choice of rays, to drastically wrong results.It was heuristically assumed here, that the product E (h ) • h behaves like E c • h c 1 − h h c near the caustic at location h c .Assume, that the caustic is detected between ray with index i c and i c + 1.The parameters of the function above can then be estimated by the two values closest to the caustic: The irradiance contribution at the bin including the caustic is therefore calculated by where the first integral is calculated by trapezoidal integration and the second one analytically with parameters from Eq. ( 13).The same procedure is used for the contribution of the folded part of the wavefront.With this approach a robust behavior of the averaging is achieved, see below the details on the sensitivity.
It was also investigated if a better choice of the power m of the heuristic model 1 − h h c m (in the procedure described above, m = −1/2 was used) could improve the result.The power m can be estimated, if the three data points closest to the caustic are fitted to the model.However, an improvement could not be observed.
As in most cases an incoherent source is assumed (the limit of a point source is just a geometrical approximation) the integration has to be done for all parts of overlapping wavefronts separately and the resulting irradiance has to be added up.Doing so the result of Figure 3 is achieved for the introduced example.
It is obvious that the method proposed in this contribution based on differential optics achieves after the modification ex- plained above comparable absolute results also in the presence of caustics, justifying its relevance.

DISCUSSION
Comparing the proposed method with the result of a Monte-Carlo simulation the following facts should be realized: • The proposed method does not show statistical noise like Monte-Carlo methods.
• The speed of the proposed method, including all described evaluation, for Figure 3 (500 differential rays) is about 20 msec, whereas a Monte-Carlo simulation (2 000 000 randomly emitting rays without any symmetry considerations) takes about 30 sec.Thus the method is about 3 magnitudes faster, allowing for comfortable trial and error or design by slider methods.
• The proposed method results in an only transversal averaged profile of the irradiance.In Monte-Carlo methods 2D-bins are used, thus also some averaging in sagittal direction is realized.
• In this implementation, the method is limited to spherical systems with an on-axis point source.Extensions are discussed below.
As explained above, the averaged irradiance at the bin which includes the caustic may be sensitive depending on the way the averaging is implemented at this point.Employing the method proposed in this paper (i.e.Eq. ( 14)), Figure 4 shows this averaged value of the example depending on the number of rays used.
It is seen that the result converges rapidly, for more than 500 rays the change is less than 1%.For comparison, a Monte-Carlo simulation with 2 000 000 rays will result in about 1800 rays in the bin, therefore the statistical uncertainty in a standard implementation would be about 2%.

ILLUMINATION DESIGN BY SLIDERS
Based on the algorithm presented above, a MATLAB® applet was developed mainly for education purpose that allows to simulate the irradiance profile behind a spherical lens illuminated by an on-axis lambertian source (it is also possible to change the source model to a gerneralized lambertian source).By default, the source is located at the focal point of the lens ('ds = 0' in Figure 5).The user may change the source on axis position ds, the shape of the lens via the Coddingtion shape factor X (X = (c 1 + c 2 )/(c 1 − c 2 )) and the detector location d by sliders, see Fig. 5 (the slider for the detector location is located below the layout window and is therefore not displayed in this figure).
As output of the applet, the layout (see Figure 1) and the averaged irradiance profile, similar to Figure 3 is displayed.The output change is immediately perceived after an input is given, what would not be possible with Monte-Carlo based methods.
The applet allows to pre-design appropriate illumination tasks by trial and error using sliders.The user quickly experience the following facts: • Locating the source in the focal point of the lens will result in a bright center spot together with a bright ring in some distance of the center due to the caustic, regardless of the lens shape.The profile strongly differs from a homogeneous spot.Some kind of unexpected focusing of the rays emitted at larger angle due to the caustic is present.Some connection to the spherical aberration of lenses may be associated by a user familiar to aberration theory.
• In order to achieve a more homogeneous profile, the source has to be considerably defocused towards the lens.The bright center spot will disappear but the irradiance of the bright ring will be amplified.Only for a very divergent output ray fan a near to homogeneous (but very large) profile can be achieved.
• Changing the shape of the lens from plano-convex (X = −1) to convex-plan (X = 1) aggravates the situation.A shape near to plano-convex (X = −0.8 for n = 1.5, known perhaps from spherical aberration theory of a thin lens) is the best compromise.
• A pre-design for the position of the source and the shape of the lens can be quickly found.This pre-design has to be verified or detailed by simulations based on a extended light source model, and probably Monte-Carlo methods.However, the overall design process is considerably speeded up.
For newcomers in the field of illumination design, getting the same insight and results with standard Monte-Carlo based illumination design software would be very laborious and time-consuming.

ALTERNATIVE DIFFERENTIAL METHOD AND EXTENSION FOR ASPHERICAL SYSTEMS
In several papers [3]- [5] it has been shown, how differential geometry can be exploited to calculate the irradiance of a local wavefront associated with a ray.Rays can therefore be extended using the second derivatives of the surfaces in order to trace the associated local wavefronts and its irradiance.This method yields in the case of a spherical system with an on axis point source exactly the same result as described in Sec. 2 and shown in Figure 2. In the author's implementation, the speed of the latter method is only about a factor of 2 slower than the analytical method, thus still allowing illumination design by sliders.As already mentioned, the method is directly applicable to complex systems with more than just one element including also aspherics.
A common goal in illumination design is to homogeneously illuminate a target.Experimenting with the above introduced applet, quickly the question arises, which lens shape achieves this goal for a lambertian point source.From Eq. ( 9) it is obvious that the system with h (ϕ) = f • sin (ϕ), i.e. a system which fulfills the sine-condition, is the solution [6].In the case of a single element with two surfaces, this condition leads to a coupled set of differential equations [7] and its solution to highly aspherical surfaces.The method mentioned in the paragraph before is able to confirm this solution.
Another application of the last mentioned method is tolerancing of illumination systems, were caustics could be avoided.This is for example the case when a small, inhomogeneously irradiating source is combined with a free form surface (continuous in the second derivative) in order to achieve a predescribed irradiance on a target surface, see e.g.[8] and Figure 6.Tolerancing of such systems requires a lot more than ten simulations, each with a different system set up.In the case of a standard Monte-Carlo based approach, this is very time consuming and may lead to some discomfort of the designer.The approach suggested here, based on differential geometry and local wavefronts, can accelerate this phase dramatically.

OUTLOOK
The application above directly leads to the questions, if there is an approach in the case of caustics and off-axis sources.In that case the method described in Sec. 3 has to be extended from one to two dimensions.The author is currently working on an implementation of this extension.First results show, that this is basically feasible yielding to non-noisy irradiance distributions.However, the acceleration of simulation time is no longer as pronounced as in the one dimensional case.

P
FIG. 1 Layout (with different scaling in z and y) of a simple system (plano-convex lens with n = 1.5, d = 10, f = 100, h = 25, point source at the focal point with N A = 0.24) with corresponding rays (51 rays, equidistant in positive emitting angle).

FIG. 2
FIG. 2 Ray height and irradiance at the detector plane (d = 2000) for a lambertian point source (I(0) =1W/sr) calculated with the Q-U-method and its derivative for 500 rays, emitting equidistantly in positive direction, as shown in Figure 1 (black).

P
FIG. 3 Center bin location h bin and averaged irradiance after adding all overlapping parts of the wavefront (d =2000) for a (truncated) lambertian point source (I(0)=1W/sr) calculated with the Q-U-method and its derivative for 500 equidistant emitting rays shown (green).Result of a comparable Monte-Carlo simulation (2 000 000 rays, 33 bins at positive height shown) for comparison in blue.

FIG. 4
FIG.4Averaged irradiance at the bin which includes the caustic (h bin = 8 with bin size h bin,ub − h bin,lb = 0.5) for the system of Figure1with respect to the number of equidistant emitted rays n ray .For n ray ≥ 500 the change is less than approx.1%.The absolute value of the irradiance at n ray ≥ 500 is within the range of the standard deviation of a Monte Carlo simulation with an equivalent of a total of 30 million rays.

P
FIG. 5 Part of the graphical user interface (only the input) for simulating the irradiance profile behind a spherical lens illuminated by a lambertian point source.The shape of the lens can be changed via the Coddingtion shape factor X .Not shown is the slider that allows to change the detector position d .

P
FIG.6Plane symmetric irradiance distribution produced by a nonsymmetric small point source together with a free form mirror (only the part right to the axis of symmetry is displayed).Left: Monte-Carlo simulation (10 6 rays).Right: Irradiance calculated by differential optics with local wavefronts (1000 extended rays), without binning.The latter shows no noise, has higher spatial resolution and its calculation is about 2 orders of magnitude faster.To the above right, a caustic is present, which is not perceivable in the Monte-Carlo simulation due to its binning.