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 (Ck3) 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 nt from a planar, simply-connected Lambertian source in a medium of refractive index nr with basic radiance L (henceforth simply radiance) is given by:

E=Lnt2Ω=Lnt2Ω[R^Nt]dΩ,
where ∬Ω 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, R^ denotes the unit-vector directions of these rays, Nt the unit normal at the target point, and Ω 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:

Ω=12C[Nt×(ts)ts2dsc]12[NtC(ts)×dscts2],
where sc denotes the location of the boundary curve C of the beam-print on the optical surface s(u, v), with local curvilinear coordinates 0 ≤ u, v ≤ 1, such that sc = s(uc, vc).

By Fermat’s principle, the location of each sc must be a stationary point (in u, v) of the optical path length (OPL) between each point rc 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:

uc,vc=argminu,v[0,1]{OPL(rc,s(u,v),t)},where OPL(r,s,t)=nrsr+ntts.

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 {ri}. Associated with the target point is an angular region contour, a set {si} comprising a discretised subset of {sc} as in Eq. (2). Each si is determined using its corresponding ri, t by means of minimisation as in Eq. (3). The result of the minimisation process is the set {ui, vi}, and evaluating the surface across these gives {s(ui, vi)} ≡ {si}, the set of solved surface points forming the discretised contour.

As this minimisation occurs over the bounded region u, v ∈ [0, 1], each ui, vi 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:

Ω=12[Nti=1NΓi],whereΓiAngle(Ri,Ri1)(Ri×Ri1)Ri×Ri1,
where {Ri} ≡ {tsi} denotes the set of vectors from each solved contour point to the target point, Angle(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 Nt. 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

  1. Represent border of planar source as N ordered points {ri}.
  2. Consider a target point t on a target surface with local normal Nt.
  3. For each i ∈ [1, N], find:
    ui,vi=argminu,v[0,1]{OPL(ri,s(u,v),t)},
    with OPL(r, s, t) defined in Eq. (3), and evaluate the set {si} = {s(ui, vi)}.
  4. Evaluate the discretised contour integral Ω=12i=1N[NtΓi], as in Eq. (4).
  5. Combine with the target medium’s refractive index nt and the source’s radiance 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 fn(x), suitable conditions at iteration n would be:

|Δfn||fn|<ϵf,Δxnxn<ϵx,
where Δfnfnfn−1, Δxnxnxn−1, and ϵf,x denote the relative tolerances in x and f respectively. The smaller these tolerances, the greater the accuracy of the solved ui, vi (thus 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 Δff″(Δx)2/2 by Taylor expansion about the minimum [6]. Exact values chosen for this study were ϵf = ϵsys and ϵx=ϵ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 {ui, vi} 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.

 figure: Fig. 1

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 | PPT Slide | PDF

 figure: Fig. 2

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 | PPT Slide | 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 (nr = 1.5), and the target in air or vacuum (nt = 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 × 107 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:

%Erro=(100%×)|ENEN=200||EN=200|,
where EN is substituted for the LightTools result as required. Also shown are results using LightToolsaccelerated mode, which approximates spline surfaces with a tessellation of planar sub-surfaces, yielding equally accurate results in much less time.

 figure: Fig. 3

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 × 107 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 | PPT Slide | PDF

Tables Icon

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.

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. (34) 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 (Ck3) 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).

References

  • View by:
  • |
  • |
  • |

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

2014 (1)

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]

2013 (1)

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

2007 (1)

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

2002 (1)

1992 (1)

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

1963 (1)

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

Bäuerle, A.

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

Berens, M.

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]

Cassarly, W. J.

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

Hirst, A.

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

Hottel, H. C.

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

Koch, D. G.

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

Loosen, P.

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]

Müller, G.

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]

Muschaweck, J.

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

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

Ries, H.

Sarofim, A. F.

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

Sparrow, E. M.

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

Stollenwerk, J.

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]

Völl, A.

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]

Wester, R.

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]

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

Adv. Opt. Technol. (1)

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

J. Heat Transf. (1)

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

J. Opt. Soc. Am. A (1)

Opt. Express (1)

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]

Proc. SPIE (2)

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

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

Other (4)

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

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

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

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

Cited By

OSA participates in Crossref's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (3)

Fig. 1
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.
Fig. 2
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.
Fig. 3
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 × 107 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.

Tables (1)

Tables Icon

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.

Equations (7)

Equations on this page are rendered with MathJax. Learn more.

E = L n t 2 Ω = L n t 2 Ω [ R ^ N t ] d Ω ,
Ω = 1 2 C [ N t × ( t s ) t s 2 d s c ] 1 2 [ N t C ( t s ) × d s c t s 2 ] ,
u c , v c = arg min u , v [ 0 , 1 ] { OPL ( r c , s ( u , v ) , t ) } , where OPL ( r , s , t ) = n r s r + n t t s .
Ω = 1 2 [ N t i = 1 N Γ i ] , where Γ i Angle ( R i , R i 1 ) ( R i × R i 1 ) R i × R i 1 ,
u i , v i = arg min u , v [ 0 , 1 ] { OPL ( r i , s ( u , v ) , t ) } ,
| Δ f n | | f n | < ϵ f , Δ x n x n < ϵ x ,
% Erro = ( 100 % × ) | E N E N = 200 | | E N = 200 | ,

Metrics