1. Introduction

Determination of irradiance from extended sources in 3D is usually performed by means of Monte-Carlo raytracing, which becomes very computationally expensive when high accuracy is desired. Moreover these results are subject to stochastic noise, making the task of performing numerical optimisation on the optical surface parameters especially challenging [1].

In this paper a deterministic, quasi-analytic method is proposed which can compute the irradiance at specified target points quickly and without noise for a specific family of illumination devices. These two qualities lend themselves well both to the evaluation of merit functions in an optimisation process, where statistical noise is the cause of a great deal of issues [2], and for other analysis, e.g. determining the sensitivity of an optic’s performance with respect to tolerances of one or more of the system’s degrees of freedom.

The core notion is to determine the visible angular region of the source on the optical surface as seen by each target point, and to apply Stokes’ theorem to simplify the irradiance integral.

2. Method

In this paper a class of optical system will be considered which consists of two homogeneous optical media of different refractive indices, the interface between them being a smooth $({\mathcal{C}}^{k\ge 3})$ optical surface, along with a homogeneous Lambertian source immersed in the first medium, emitting light towards and through said surface, onto a target plane immersed in the second medium. Losses due to absorption and Fresnel reflection are neglected, and it is assumed that total internal reflection neither truncates nor perforates any visible angular regions.

2.1. Irradiance contour integral

The irradiance E at some target point t in a medium of refractive index n_t from a planar, simply-connected Lambertian source in a medium of refractive index n_r with basic radiance L (henceforth simply radiance) is given by:

In full generality, this set of rays could be incident from a topologically disconnected angular region (multiple patches) in the target space, or even multiply-connected (with holes). To avoid significantly complicating further analysis, a simply-connected angular region condition is presupposed. The contour of this region defines the edge-rays illuminating the target point.

By applying Stokes’ theorem to the point-to-surface view factor (≡ Ω^⊥/_π) [3], Eq. (1) can be replaced with a contour integral along the edge of the flashed portion of the optical surface, the beam-print. Condensed into vector notation, the projected solid angle in Eq. (1) can be rewritten:

By Fermat’s principle, the location of each s_c must be a stationary point (in u, v) of the optical path length (OPL) between each point r_c on the boundary of the source and the target point t. In this study, only one stationary point will be considered, namely the minimum. Therefore:

To summarise, given that the stated assumptions hold, evaluating the irradiance at a given target point simply requires finding the location of the edge of the source’s visible angular region with Eq. (3), performing the appropriate contour integral across it with Eq. (2), and finally substituting the projected solid angle Ω^⊥ into the general expression for irradiance with Eq. (1). This procedure is related to one developed for rotationally symmetric systems by Koch [4].

Albeit outside the scope of this paper, non-uniform sources, Fresnel losses and material losses could be accounted for with additional factors computed using e.g. Gaussian grid integration.

2.2. Discretisation scheme

The source edge is discretised into N points {r_i}. Associated with the target point is an angular region contour, a set {s_i} comprising a discretised subset of {s_c} as in Eq. (2). Each s_i is determined using its corresponding r_i, t by means of minimisation as in Eq. (3). The result of the minimisation process is the set {u_i, v_i}, and evaluating the surface across these gives {s(u_i, v_i)} ≡ {s_i}, the set of solved surface points forming the discretised contour.

As this minimisation occurs over the bounded region u, v ∈ [0, 1], each u_i, v_i pair can be found by first performing a coarse global optimisation pass (e.g. using “DIRECT”), then using this as the starting point for local optimisation (e.g. using “MMA”) for a refined result [6].

Evaluation of the contour integral for projected solid angle at t using Eq. (2) need not be performed using explicit numerical integration, but can instead with the hemisphere method [7], which approximates the resulting contour as a polygon, and explicitly converts the contour integral into a sum over the set of solved surface points:

Finally, substitution back into Eq. (1) gives the signed irradiance at each target point, where the sign results from both the winding direction around the beam-print contour and the orientation of each R with respect to the target normal N_t. For design tasks, the sign conveys potentially useful information about caustics and crossings [5], however when only performing forward simulations, the sign can be discarded, giving simply the irradiance.

2.3. Algorithm summary

2.4. Stopping criteria

Since it is generally not possible to determine Eq. (5) analytically, numerical optimisation must be used. Also needed are suitable stopping criteria, terminating the process once any is reached. For an iterative function f_n(x), suitable conditions at iteration n would be:

Since the OPL minima are non-zero, relative tolerances are more useful than absolute, and also Δf ≈ f″(Δx)²/₂ by Taylor expansion about the minimum [6]. Exact values chosen for this study were ϵ_f = ϵ_sys and ${\mathit{ϵ}}_x=\sqrt{{\mathit{ϵ}}_{\text{sys}}}$, with machine precision ϵ_sys ≈ 2.22 × 10⁻¹⁶.

2.5. Surface edges

For target points sufficiently far offset from a common optical axis, the edge of a beam-print may appear to pass off the edge of the optical surface, effectively truncating that target point’s view of the source. In these cases, some OPL minima will occur at points along the surface edge. If this is identified, or otherwise implicitly handled by the optimisation algorithm, the minimisation reduces to a simpler univariate optimisation in the OPL constrained along that edge line. The beam-print method will therefore still produce a valid irradiance value.

In extreme cases the source may become entirely truncated. This manifests as all solved {u_i, v_i} lying clustered together on the same edge or corner of the surface, unless the topological assumptions made in Sec. 2.1 are violated, which is explicitly neglected. This is a result of the OPL minimum’s trough being sufficiently deep (or broad) that it dominates even within u, v ∈ [0, 1]. The resulting irradiance value is therefore, correctly, zero.

3. Simulation example

3.1. Specification of test geometry

A sample system was constructed in the form of Fig. 1 consisting of a single planar source, a single refractive optical surface, and a dummy target plane. The source is a disk of 1mm radius in the z = 0 plane, with Lambertian emittance and a radiant (monochromatic) flux of 1W into the entire hemisphere. The surface is a bicubic B-spline, with a profile similar to actual freeforms in multiple applications, such as street lighting, defined over a roughly rectangular region, with a width and breadth of 5mm and 10mm respectively, and a height of approximately z = 4mm. A rendering is shown in Fig. 2. The target is a 40-by-80mm rectangle in the z = 40mm plane.

 

Fig. 1 A planar Lambertian source is immersed in an initial optical medium, separated by a freeform surface from a second medium, wherein a target point lies whose irradiance value is to be determined. The visible angular region (or beam-print) of the source on the surface, as seen by this target point, is projected into the target surface’s local space.

Download Full Size | PDF

 

Fig. 2 On the left is a wireframe render of the test profile used for verification. Of particular significance is the saddle point around the centre. On the right is a top-down view of the surface, showing a subset of solved beam-prints between the test source and target, whose sizes and shapes are determined by their relative orientation to and size of the target surface.

Download Full Size | PDF

All elements are mutually aligned and centred on the z-axis about x, y = 0. The source is immersed in a dummy glass-like material (n_r = 1.5), and the target in air or vacuum (n_t = 1.0).

3.2. Comparison with Monte-Carlo raytracing

The irradiance distributions through the test surface across the target plane were determined using both LightTools 8.5 [8] Monte-Carlo raytracing and the quasi-analytic beam-print approach. For the former, the target plane was broken into 20-by-40 bins, with 2.5 × 10⁷ rays traced. For the latter, the source boundary was discretised into N points, and the target plane into a 20-by-40 grid, corresponding to the centres of the raytracing bins. Both sets of calculations were performed on a typical Windows server (Intel Xeon X5570 @ 2.93GHz, 15GB RAM).

The beam-print simulation was repeated for increasing values of N from 20 to 200. Figure 3 shows that the irradiance values even at only N = 20 are visually equivalent to LightTools. This agreement is quantified in Table 1, demonstrating convergence by N = 200 using:

 

Fig. 3 Comparison plot of the target plane irradiance through the test surface. On the left, unsmoothed results from Monte-Carlo raytracing, which are very noisy even at 2.5 × 10⁷ rays. On the right, unsmoothed results from the beam-print method at N = 20, which are simultaneously in agreement with and less noisy than those from Monte-Carlo raytracing.

Download Full Size | PDF

Table 1. The execution time increases linearly with N. Results at N = 200 are considered to have converged to the quasi-analytic result. Comparison to LightTools is underneath, along with LightTools’ own estimate for the average Monte-Carlo statistical noise.

View Table

3.3. Discussion

For the beam-print method there are two primary sources of error to consider. The first is the accuracy of the discretised contour integral. Only with sufficiently many points along the source edge will the sum in Eq. (4) be a good approximation of the true integral in Eq. (2). For sources with a convex contour, the tendency will be to underestimate. The second is the tolerances used for OPL minimisation, with details given in Sec. 2.4. The result is that, for applicable cases, convergence with the beam-print method can be achieved quickly, and with relatively few rays.

A noteworthy property is the independence of the irradiance calculations at each target point. Although this provides an obvious path for parallelisation, it should be noted that this on the other hand prevents the ability to dynamically adjust the target points without recomputation.

While the considered design case and assumptions may seem to limit the scope of this work, three arguments support its interest. Firstly, the model is very similar to many single-chip LEDs with dome-type encapsulation, and thus behave closely. Secondly, for more complex light sources (multi-chip, inhomogeneous, non-Lambertian etc.), a first-approximation model as proposed here is necessary for initial designs and analysis, from which further optimisation, more complex models, or refinement using ray-sets can be developed. Thirdly, the formulation in Eqs. (3–4) provides a crucial extension to a related method published by Koch [4]. By further application of Fermat’s principle, the derivative of the irradiance with respect to the parameters of the surface representation can be deduced, providing quasi-analytical gradient information of enormous interest to various design techniques, from optimisation [9] to Newton-Raphson solution [10]. Although outside the scope of this paper, this is the subject of another currently in preparation.

4. Conclusions

Presented in this paper has been a fast, deterministic method for the evaluation of irradiance from an extended Lambertian source through a smooth $({\mathcal{C}}^{k\ge 3})$ single-surface refractive optical system.

Correctness has been demonstrated by comparison with Monte-Carlo raytracing results for a sample system, converging to a % Error, defined in Eq. (7), of 1.06%, well within LightTools’ average error estimate of 1.63%, at the same time producing an almost noiseless distribution.

This discretised method can, with trivial modification, be applied to a single reflective instead of refractive surface, and is expected to be extensible to multiple sequential surfaces, by expanding the OPL function in Eq. (3), and by including a u, v pair for each surface. Additionally, this formulation provides a basis to readily evaluate quasi-analytically the derivative of the irradiance values with respect to the surface parameters, the subject of further work.