## Abstract

Presented is a quasi-analytic method for irradiance evaluation through a single refractive surface from a single Lambertian source. The method is compared to Monte-Carlo raytracing for a sample system, producing in much less time an irradiance distribution equal to within the latter’s statistical noise. In addition to its interest to optical analysis, the method is also useful for optical design problems by allowing for fast, noise-free evaluation of merit functions, along with their derivatives.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 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*with basic radiance

_{r}*L*(henceforth simply radiance) is given by:

_{Ω}dΩ is understood to denote the solid angle integral over every ray direction impinging on the target point which, having passed through the optical surface, originated from the source, $\widehat{\mathbf{R}}$ denotes the unit-vector directions of these rays,

**N**

*the unit normal at the target point, and Ω*

_{t}^{⊥}the projected solid angle subtended by these rays into the target point’s local plane.

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:

**s**

*denotes the location of the boundary curve*

_{c}*C*of the beam-print on the optical surface

**s**(

*u*,

*v*), with local curvilinear coordinates 0 ≤

*u*,

*v*≤ 1, such that

**s**

*=*

_{c}**s**(

*u*,

_{c}*v*).

_{c}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**

*on the boundary of the source and the target point*

_{c}**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**

*} comprising a discretised subset of {*

_{i}**s**

*} as in Eq. (2). Each*

_{c}**s**

*is determined using its corresponding*

_{i}**r**

*,*

_{i}**t**by means of minimisation as in Eq. (3). The result of the minimisation process is the set {

*u*,

_{i}*v*}, and evaluating the surface across these gives {

_{i}**s**(

*u*,

_{i}*v*)} ≡ {

_{i}**s**

*}, the set of solved surface points forming the discretised contour.*

_{i}As this minimisation occurs over the bounded region *u*, *v* ∈ [0, 1], each *u _{i}*,

*v*pair can be found by first performing a

_{i}*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:

**R**

*} ≡ {*

_{i}**t**−

**s**

*} denotes the set of vectors from each solved contour point to the target point, Angle(*

_{i}**a**,

**b**) the angle between vectors

**a**and

**b**, and ⊕ the

*circular next*operator, offset 1.

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

- Represent border of planar source as
*N*ordered points {**r**}._{i} - Consider a target point
**t**on a target surface with local normal**N**._{t} - For each
*i*∈ [1,*N*], find:$${u}_{i},{v}_{i}=\underset{u,v\in [0,1]}{\mathrm{arg}\mathrm{min}}\{\text{OPL}({\mathbf{r}}_{i},\mathbf{s}(u,v),\mathbf{t})\},$$with OPL(**r**,**s**,**t**) defined in Eq. (3), and evaluate the set {**s**} = {_{i}**s**(*u*,_{i}*v*)}._{i} - Evaluate the discretised contour integral ${\mathrm{\Omega}}^{\perp}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{\displaystyle {\sum}_{i=1}^{N}[{\mathbf{N}}_{t}\cdot {\mathrm{\Gamma}}_{i}]}$, as in Eq. (4).
- Combine with the target medium’s refractive index
*n*and the source’s radiance_{t}*L*to yield the signed irradiance*E*at the target point, as in Eq. (1).

#### 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:

*f*≡

_{n}*f*−

_{n}*f*

_{n}_{−1}, Δ

**x**

*≡*

_{n}**x**

*−*

_{n}**x**

_{n}_{−1}, and

*ϵ*denote the

_{f,x}*relative*tolerances in

**x**and

*f*respectively. The smaller these tolerances, the greater the accuracy of the solved

*u*,

_{i}*v*(thus

_{i}*E*), but the longer the runtime, and the greater the risk of floating-point errors preventing convergence.

Since the OPL minima are non-zero, *relative* tolerances are more useful than *absolute*, and also Δ*f* ≈ *f*″(Δ*x*)^{2}/_{2} by Taylor expansion about the minimum [6]. Exact values chosen for this study were *ϵ _{f}* =

*ϵ*

_{sys}and ${\mathit{\u03f5}}_{x}=\sqrt{{\mathit{\u03f5}}_{\text{sys}}}$, with machine precision

*ϵ*

_{sys}≈ 2.22 × 10

^{−16}.

#### 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*} 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

_{i}*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.

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*= 1.0).

_{t}#### 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^{7} 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:

*E*is substituted for the

_{N}*LightTools*result as required. Also shown are results using

*LightTools*’

*accelerated mode*, which approximates spline surfaces with a tessellation of planar sub-surfaces, yielding equally accurate results in much less time.

#### 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.

## Funding

FP7 People: Marie-Curie Actions, PITN-GA-2013-608082, “ADOPSYS.”

## Acknowledgments

The research leading to these results has received funding from the European Union’s Seventh Framework Programme under the Marie Skłodowska-Curie grant agreement stated above. Responsibility for the content of this publication lies solely with the authors.

The authors extend their thanks to Tobias Schmidt for helpful discussions on the subject, and to Synopsys for the academic license of *LightTools* at UPM.

## References and links

**1. **R. Wester and A. Bäuerle, “Light Shaping for Illumination,” Adv. Opt. Technol. **2**, 301–311 (2013).

**2. **W. J. Cassarly, “Illumination merit functions,” Proc. SPIE **6670**, 66700K (2007). [CrossRef]

**3. **E. M. Sparrow, “A New and Simpler Formulation for Radiative Angle Factors,” J. Heat Transf. **85**(2), 81–87 (1963). [CrossRef]

**4. **D. G. Koch, “Simplified irradiance/illuminance calculations in optical systems,” Proc. SPIE **1780**, 17800F (1992).

**5. **H. Ries and J. Muschaweck, “Tailored freeform optical surfaces,” J. Opt. Soc. Am. A **19**(3), 590–595 (2002). [CrossRef]

**6. **Steven G. Johnson, “The NLopt nonlinear-optimization package,” http://ab-initio.mit.edu/nlopt.

**7. **H. C. Hottel and A. F. Sarofim, *Radiative Transfer*, (1969, McGraw-Hill).

**8. ** Synopsys Inc., “LightTools® 8.5,” https://optics.synopsys.com/lighttools.

**9. **R. Wester, G. Müller, A. Völl, M. Berens, J. Stollenwerk, and P. Loosen, “Designing optical free-form surfaces for extended sources,” Opt. Express **22**, 552–560 (2014). [CrossRef]

**10. **A. Hirst and J. Muschaweck, “Irradiance tailoring for extended sources in 3D by implicit integral equation solution,” OSA Tech. Dig. (online), FT4B.2 (2015).