Extended adding-doubling method for fluorescent applications

Sven Leyre; Guy Durinck; Bart Van Giel; Wouter Saeys; Johan Hofkens; Geert Deconinck; Peter Hanselaer

doi:10.1364/OE.20.017856

1. Introduction

The behavior of light propagating through a fluorescent material is the subject of investigation in a number of fields. In the bio-medical field, fluorescence methods are valuable tools in diagnostics [1–3]. In lighting applications, a phosphor mixture is used to create white light from blue LEDs or from uv emission in a gas discharge [4–6]. In photovoltaics, fluorescent coatings are used to obtain a better match between the solar spectrum and the spectral response of the solar cell [7–9].

The optimization of the fluorescent material for lighting or photovoltaic applications is usually done either empirically [4, 10, 11], using Monte Carlo simulations [12–14] or by applying a simplified one dimensional model [15, 16]. The empirical method requires a lot of effort and time and allows usually only one parameter at the time to be investigated. Monte Carlo simulations take up a lot of computation time, since many rays must be traced in order to reduce the stochastic noise within acceptable margins [17]. The 1-D models usually depart from the well-known Lambert-Beer law which is not well suited to take scattering into account.

In this paper, a fast and accurate method to estimate the wavelength dependent angular distribution of light scattered and emitted by a fluorescent material is described. The proposed method is an extension of the existing adding-doubling algorithm for non-fluorescent samples [18].

Any macroscopic material can be considered as the conjunction of a number of plane parallel layers of the same material. If the reflection and transmission characteristics of one layer are known, the reflection and transmission characteristics of a combination of such layers can be calculated with the adding-doubling method. This method has been applied already in 1862 by Stokes [19]. The transmission and reflection characteristics of a plane parallel “single-scatter” layer can be calculated by solving the radiative transfer equation.

The adding-doubling method is commonly used to determine the angular light distribution of materials with volume scattering [20–22]. Since fluorescence can be interpreted as inelastic volume scattering, it should be possible to extend this method to determine the characteristics of materials showing both elastic and inelastic scattering.

An extended Kubelka Munk theory for fluorescence has been used to determine the reflection and transmission characteristics for the individual layers. However, the Kubelka Munk theory has some restrictions: an ideal diffuse irradiation is required, the flux propagating in the material must be diffuse, Fresnel reflections at the surface are not taken into account and the correlation between the concentration of the diffusing or fluorescent particles and the absorption and scatter constants is not straightforward [23–25]. Rosema et al. explored the possibility of the adding-doubling method for describing fluorescent layers using a two flux approach [26]. This approach returns total reflection and transmission but no angular information about the scattered radiation.

In this paper we exploit the full potential of the adding-doubling algorithm by using an n-flux model and starting from the radiative transfer equation for light propagation in a fluorescent layer to determine the reflection and transmission characteristics of a “single-scatter” layer. The n-flux model subdivides the flux propagating through a material into n fluxes each propagating within a conical segment with polar angle θ. The optical properties defined in the radiative transfer equation, in particularly absorption and scattering coefficient, are linearly dependent on the concentration [27]. This is advantageous if concentration is an optimization parameter. Furthermore, the use of the n-flux model provides angular information on the scattered and converted light.

2. General approach

For the adding-doubling method the incident flux is divided into a number of channels. The channels are conical segments, each labeled with the polar angle θ. The conical segments are represented in Fig. 1 .

Fig. 1 Subdivision of the flux into n conical segments

Download Full Size | PDF

Channel 1 contains the light travelling in directions lying between the normal on the sample and polar angle θ₁, channel 2 contains the light travelling in directions lying between θ₁ and θ₂ and so forth. The flux in the positive or the negative Z direction can be represented by a vector, where each element represents the amount of light within the conical segment. The vector contains n elements, equal to the number of channels the flux is divided in.

The reflected flux in each conical segment originating from the incident flux in each other conical segment can be written in matrix notation as given by Eq. (1):

[\begin{matrix} I^{-} (θ_{1}) \\ ⋮ \\ I^{-} (θ_{n}) \end{matrix}] = [\begin{matrix} R (θ_{1}, θ_{1}) & \dots & R (θ_{1}, θ_{n}) \\ ⋮ & ⋮ \\ R (θ_{n}, θ_{1}) & \dots & R (θ_{n}, θ_{n}) \end{matrix}] x [\begin{matrix} I^{+} (θ_{1}) \\ ⋮ \\ I^{+} (θ_{n}) \end{matrix}]

Herein is I⁺(θ_j) the incident flux in channel j propagating in the positive Z direction, I^-(θ_i) the reflected flux in channel i, propagating in the negative Z direction and R(θ_i, θ_j) gives the reflection from channel j into channel i. In a similar way a transmission matrix can be composed. The angle θ can take values from 0 to 90°.

3. Adding-doubling

3.1 Adding-doubling for non-fluorescent slabs

The adding-doubling method for slabs where only scattering and absorption occurs is extensively described in literature [20–22]. If the reflection (R) and transmission (T) matrices of two given slabs at a particular wavelength are known, the new R and T matrices of a slab consisting of the two slabs joined together can be calculated with the adding method. The doubling method is similar as the adding method for the combination of two identical slabs. For one slab the following equations can be written for each wavelength with Eq. (2) and Eq. (3):

I_{1}^{+} = R_{10} I_{1}^{-} + T_{01} I_{0}^{+}

I_{0}^{-} = R_{01} I_{0}^{+} + T_{10} I_{1}^{-}

Here, I_x^± is the flux within the channels in vector notation, the plus and minus sign denote the propagation direction along the positive and negative Z-axis respectively. The subindex represents the location of an interface starting from 0. R_xy respectively T_xy is the reflection matrix and transmission matrix for a flux incident on interface x of the slab defined by interface x and y. This is schematically shown in Fig. 2 .

Fig. 2 Schematic representation of the vectors and matrices for a single non-fluorescent slab

Download Full Size | PDF

For a second slab (slab 1-2), analogous equations can be written, given by Eq. (4) and Eq. (5):

I_{2}^{+} = R_{21} I_{2}^{-} + T_{12} I_{1}^{+}

I_{1}^{-} = R_{12} I_{1}^{+} + T_{21} I_{2}^{-}

After manipulation of the equations of slabs 0-1 and 1-2, the reflection and transmission matrices of the joined slab are given in Eq. (6) and Eq. (7).

R_{20} = R_{21} + T_{12} {(E - R_{10} R_{12})}^{- 1} R_{10} T_{21}

T_{02} = T_{12} {(E - R_{10} R_{12})}^{- 1} T_{01}

With E the unity matrix. A more extensive description can be found in literature [20–22].

3.2 Adding-doubling for fluorescent slabs

At excitation wavelengths, the regular adding-doubling can be applied to obtain the reflection and transmission characteristics at each individual excitation wavelength interval. At an emission wavelength, the adding-doubling method described above must be adjusted for the contribution of the fluorescence. Indeed, the flux in forward and backward direction at an emission wavelength is also dependent on the incident flux at the excitation wavelengths.

For fluorescent slabs, the relevant wavelength ranges of excitation and emission are subdivided into respectively N wavelength intervals Δλ^X_i and L wavelength intervals Δλ^M_j, with i going from 1 to N and j going from 1 to L. Each wavelength region Δλ^X_i can have a contribution to the flux at emission wavelength interval Δλ^M_j, as schematically shown in Fig. 3 .

Fig. 3 Schematic representation of the contribution of all excitation wavelength intervals to one emission wavelength interval

Download Full Size | PDF

The matrix equations Eqs. (2) and (3) for a particular emission wavelength interval Δλ^M_j must be adapted and are given by Eq. (8) and Eq. (9):

\begin{matrix} I_{1}^{+} = R_{10} I_{1}^{-} + T_{01} I_{0}^{+} \\ + \sum_{i = 1}^{N} R_{10}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{1}^{-} (Δ λ_{i}^{X}) + T_{01}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{0}^{+} (Δ λ_{i}^{X}) \end{matrix}

\begin{matrix} I_{0}^{-} = R_{01} I_{0}^{+} + T_{10} I_{1}^{-} \\ + \sum_{i = 1}^{N} R_{01}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{0}^{+} (Δ λ_{i}^{X}) + T_{10}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{1}^{-} (Δ λ_{i}^{X}) \end{matrix}

R_xy and T_xy are the reflection and transmission matrices at the selected emission wavelength, representing scattering and absorption. The dependence on wavelength has been omitted from the equations for readability reasons. The matrices R^c_xy(Δλ^X_i, Δλ^M_j) and T^c_xy(Δλ^X_i, Δλ^M_j) represent the conversion of light from an excitation wavelength interval Δλ^X_i towards a particular emission wavelength interval Δλ^M_j for each incident and scattered conical segment.

Equation (8) is graphically illustrated in Fig. 4 . If no fluorescence occurs, all elements of R^c_xy and T^c_xy are identical to zero and Eqs. (8) and (9) become equal to Eqs. (2) and (3).

Fig. 4 Schematic representation of the vector-matrix equation for a fluorescent slab for one excitation wavelength Δλ^X_i (solid lines) and one emission wavelength Δλ^M_j (dashed lines). The equations represented are Eq. (2) for the excitation wavelength interval and Eq. (8) for the emission wavelength interval.

Download Full Size | PDF

Similar equations can be written for a second layer 1-2. The reflection and transmission matrices for the joined layer 0-2 are given in Eqs. (10-13). Details about the calculations can be found in appendix A.

R_{20} = R_{21} + T_{12} {(E - R_{10} R_{12})}^{- 1} R_{10} T_{21}

T_{02} = T_{12} {(E - R_{10} R_{12})}^{- 1} T_{01}

\begin{matrix} R_{_{20}}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) = \\ T_{12} (^{E - R_{10} R_{12}) - 1} [\begin{array}{l} R_{10} (\begin{array}{l} R_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {\begin{cases} {[E - R_{10}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) R_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M})]}^{- 1} \\ R_{10}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) T_{21}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) \end{cases}} \\ + T_{21}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) \end{array}) \\ + R_{10}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) R_{10}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M})]}^{- 1} T_{21} (Δ λ_{i}^{X}) \end{array}] \\ + R_{21}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) + T_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{10} (Δ λ_{i}^{X}) R_{12} (Δ λ_{i}^{X})]}^{- 1} R_{10} (Δ λ_{i}^{X}) T_{21} (Δ λ_{i}^{X}) \end{matrix}

\begin{matrix} T_{_{02}}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) = \\ T_{12} (^{E - R_{10} R_{12}) - 1} {\begin{cases} R_{10} R_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{10} (Δ λ_{i}^{X}) R_{12} (Δ λ_{i}^{X})]}^{- 1} T_{01} (Δ λ_{i}^{X}) \\ + R_{10}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{12} (Δ λ_{i}^{X}) R_{10} (Δ λ_{i}^{X})]}^{- 1} R_{12} (Δ λ_{i}^{X}) T_{01} (Δ λ_{i}^{X}) \\ + T_{01}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) \end{cases}} \\ + T_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{10} (Δ λ_{i}^{X}) R_{12} (Δ λ_{i}^{X})]}^{- 1} T_{01} (Δ λ_{i}^{X}) \end{matrix}

4. Initial matrices

4.1 Initial matrices for non-fluorescent slabs

In the previous section, it has been described how the transmission and reflection matrices for a stack of layers are obtained when the initial R and T matrices are known for each layer. In this section we explain how the initial R and T matrices can be determined. For non-fluorescent applications this is already extensively described in literature [28]. A short summary will be given. The initial matrices are calculated for thin “single-scatter” layers by solving the radiative transfer equation.

It is assumed that the sample is a uniform infinite plane parallel slab of finite thickness. The incident beam and the scattering behavior of the material are assumed to be azimuthally symmetric. The calculations start from the time-independent azimuthally-averaged radiative transfer equation in the case of unpolarized radiation in plane-parallel slab geometry without sources given by Eq. (14) [28]:

ν \frac{ϑ I (ν)}{ϑ z} = - (μ_{a} + μ_{s}) I (ν) + μ_{s} \int_{- 1}^{1} P (ν, ν') I (ν') d ν'

With I the flux within a selected channel, ν the cosine between the average propagation direction in the selected channel and the normal on the slab (ν can take values from −1 to 1). P(ν,ν’) is the azimuthally averaged phase function (normalized to 1), µ_a and µ_s are respectively the absorption and scattering coefficient. The phase function models the spatial redistribution of scattered light.

A commonly used phase function is the Henyey-Greenstein phase function [29]. This phase function was originally developed to model diffuse radiation in the galaxy, but it is also commonly used for medical [30] and biological [31] applications. The Henyey-Greenstein phase function P_HG(Υ) is given in Eq. (15), where Υ is the cosine of the angle between incident and scattered direction.

P_{H G} (ϒ) = \frac{1}{2} \frac{1 - g^{2}}{{(1 + g^{2} - 2 g ϒ)}^{\frac{3}{2}}}

Herein g represents the anisotropy factor. Using the Henyey-Greenstein phase function as a fixed phase function for all materials, a non-fluorescent material can be described with three optical parameters: µ_a, µ_s, g.

Wiscombe et al. extensively described the solution of the radiative transfer equation given in Eq. (14) [28]. After discretization and integration over a slab 0-1 with thickness Δd, Eq. (14) can be rewritten as Eq. (16) and Eq. (17):

(E + α) I_{1}^{+} = (E - α) I_{0}^{+} + β (I_{0}^{-} + I_{1}^{-})

(E + α) I_{0}^{-} = (E - α) I_{1}^{-} + β (I_{0}^{+} + I_{1}^{+})

With E the unity matrix and α and β are defined by Eq. (18) and Eq. (19):

α = M^{- 1} (E - μ_{s} \frac{Δ d}{2} p^{+ +} c)

β = M^{- 1} μ_{s} \frac{Δ d}{2} p^{+ -} c

M, c and p are matrices defined by respectively Eq. (20), Eq. (21) and Eq. (22):

M (i, j) = ν_{i} δ_{i j}

c (i, j) = c_{i} δ_{i j}

p_{i j}^{+ \pm} = P (ν_{i}, \pm ν_{j})

With ν_i the value of ν in conical segment i. The weights in matrix c are determined by the used quadrature scheme. Prahl et al. described the use of a Radau and Gaussian quadrature scheme for the discretization [18]. Further manipulation of Eqs. (16) and (17) and identification with Eqs. (2) and (3) results in the initial matrices given by Eq. (23) and Eq. (24) [28]:

R_{10} = 2 Γ β {(E + α)}^{- 1}

T_{01} = 2 Γ - E

With Γ defined by Eq. (25):

Γ = {[E + α - β {(E + α)}^{- 1} β]}^{- 1}

Similar calculations can be done to obtain R₀₁ and T₁₀, but these matrices will be identical to the ones given by Eqs. (23) and (24).

4.2 Initial matrices for fluorescent slabs

At the excitation wavelength and other non-emission wavelengths, the initial R and T matrices for a fluorescent slab are calculated as described above for non-fluorescent applications. For the calculation of the initial R and T matrices at an emission wavelength interval, an additional contribution of the fluorescence has to be taken into account [32].

Equation (26) gives the radiative transfer equation for light within the emission wavelength interval Δλ^M_j with central wavelength λ^M_j.

\begin{array}{l} ν \frac{ϑ I (ν)}{ϑ z} = - (μ_{a} + μ_{s}) I (ν) + μ_{s} \int_{- 1}^{1} P (ν, ν') I (ν') d ν' \\ + w (Δ λ_{j}^{M}) \sum_{i = 1}^{N} {μ_{a} (Δ λ_{i}^{X}) Q E (Δ λ_{i}^{X}) \frac{λ_{j}^{M}}{λ_{i}^{X}} [\int_{- 1}^{1} \frac{1}{2} I (ν', Δ λ_{i}^{X}) d ν'] Δ λ_{i}^{X}} \end{array}

In Eq. (26) I(ν) is the channel flux at the selected emission wavelength interval, I(ν’, Δλ^X_i) represents the channel flux at the excitation wavelength interval Δλ^X_i. The third term on the right hand side of Eq. (26) represents the contribution of each excitation wavelength interval Δλ^X_i (with central wavelength λ^X_i) to the selected emission wavelength interval.

In Eq. (26) µ_a and µ_s are respectively the absorption and scattering coefficient for the selected emission wavelength interval, µ_a(Δλ^X_i) is the absorption coefficient at excitation wavelength interval Δλ^X_i. QE(Δλ^X_i) expresses the quantum efficiency the excitation wavelength interval. The wavelength ratio is used to convert the number of photons to powerflux.

Since, the emission of the fluorescence is isotropic, the value for the phase function for fluorescence has been taken as 1/2 [32]. Only a part of the fluorescent light generated by each excitation wavelength interval Δλ^X_i will be emitted within the selected wavelength interval Δλ^M_j (schematically shown in Fig. 5 ). A weight factor w(Δλ^M_j) has been introduced to describe the fraction of light generated by an excitation at Δλ^X_i and emitted within the selected emission wavelength region. The coefficient w(Δλ^M_j) is normalized so that the integration over the entire emission wavelength region is 1 with Eq. (27):

Fig. 5 Schematic representation of the contribution of one excitation wavelength interval to all emission wavelength intervals, the contribution to each emission wavelength region is denoted with the factor w(Δλ^M_j).

Download Full Size | PDF

\sum_{j = 1}^{L} w (Δ λ_{j}^{M}) = 1

The same steps to solve the radiative transfer equation at non-emission wavelengths have been used to solve Eq. (26). After discretization and integration over a slab 0-1 with thickness Δd, Eq. (26) can be rewritten as Eq. (28) and Eq. (29) (see appendix B):

\begin{matrix} (E + α) I_{1}^{+} = (E - α) I_{0}^{+} + β (I_{0}^{-} + I_{1}^{-}) \\ + \sum_{i = 1}^{N} γ_{i} [I_{0}^{+} (Δ λ_{i}^{X}) + I_{1}^{+} (Δ λ_{i}^{X}) + I_{0}^{-} (Δ λ_{i}^{X}) + I_{1}^{-} (Δ λ_{i}^{X})] \end{matrix}

\begin{matrix} (E + α) I_{0}^{-} = (E - α) I_{1}^{-} + β (I_{0}^{+} + I_{1}^{+}) \\ + \sum_{i = 1}^{N} γ_{i} [I_{0}^{+} (Δ λ_{i}^{X}) + I_{1}^{+} (Δ λ_{i}^{X}) + I_{0}^{-} (Δ λ_{i}^{X}) + I_{1}^{-} (Δ λ_{i}^{X})] \end{matrix}

The matrices α and β are still defined by Eq. (18) and Eq. (19). The matrix γ_i is defined by Eq. (30):

γ_{i} = w (Δ λ_{j}^{M}) μ_{a} (Δ λ_{i}^{X}) Q E (Δ λ_{i}^{X}) \frac{λ_{M}}{λ_{X i}} \frac{Δ d}{4} p_{x} c

The matrix p_x is defined by Eq. (31):

p_{x} (i, j) = 1

Further manipulation of Eqs. (28) and (29) and identification with Eqs. (8) and (9) results in the initial matrices given by Eq. (32), Eq. (33) and Eq. (34):

R_{10} = 2 Γ β {(E + α)}^{- 1}

T_{01} = 2 Γ - E

\begin{matrix} R_{10}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) = \\ T_{01}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) = \\ [Γ β {(E + α)}^{- 1} + E] γ_{i} [E + T (Δ λ_{i}^{X}) + R (Δ λ_{i}^{X})] \end{matrix}

T(Δλ^X_i) and R(Δλ^X_i) are the reflection and transmission matrices calculated at the excitation wavelength region Δλ^X_i, calculated with the method for non-fluorescent slabs. The matrix Γ is still defined by Eq. (25). A more extensive description of the derivation of the equations can be found in appendix B.

Similar steps can be followed to calculate R₀₁, T₁₀, R^c₀₁(Δλ^X_i, Δλ^M_j) and T^c₁₀(Δλ^X_i, Δλ^M_j), the results will be the same as Eqs. (32-34) since the slab will behave identical for light entering the slab in the opposite direction.

The initial matrices R₁₀ and T₀₁, describing the behavior of incident light at the emission wavelength, are identical to the matrices in defined Eqs. (23) and (24). Indeed, fluorescence is independent on the absorption and scattering characteristics at the emission wavelengths. The R^c and T^c matrices describing the contribution due to fluorescence are equal for transmission and reflection. This meets the expectations since the emission due to fluorescence is assumed to be isotropic in the single-scatter layer. If there is no fluorescence (QE = 0), the matrices with superindex c become zero and the equations become equal to the equations valid for the non-fluorescent case. The formulas for the non-fluorescent case can thus be considered as special cases of the general formulas taking fluorescence into account.

5. Validation

To validate the extended adding-doubling method for fluorescent applications, the predictions obtained with this method have been compared to simulation results obtained by Monte Carlo ray tracing with TracePro 7.1 (Lambda Research).

The model in the Monte Carlo simulation consists of a slab with thickness of 1 mm and refractive index of 1.5. In the adding-doubling method, the geometry of the material is assumed to be a plane-parallel slab of infinite size with finite thickness. In the Monte Carlo simulation, the size of the sample has been chosen large enough such that the rays escaping through the edges will be negligible. The incident light is chosen perpendicular onto the surface of the slab.

Only absorption by arbitrarily chosen fluorescent particles has been assumed. Gaussian absorption and emission spectra with respectively central wavelength 450 and 600 nm and σ = 12 and 20 nm have been used.

The scattering properties of the fluorescent particle are determined by the scattering coefficient and the asymmetry factor. The spectral characteristics of these parameters have been chosen arbitrarily. The Henyey-Greenstein phase function was used for both calculations and simulations. The spectrum of the incident light has been chosen to be constant between 380 and 780 nm. For the calculations, a wavelength step of 10 nm, corresponding to 41 wavelengths, has been selected .

The angular distribution of the scattered light calculated with the adding-doubling method was validated by calculating the scattered intensities I(θ_s) for three wavelengths: one excitation wavelength (450 nm), one emission wavelength (600 nm), and a wavelength not participating in the fluorescence (750 nm).

For the adding-doubling calculations, 32 internal channels were used, 16 of them give rise to light escaping from the sample. The light within the other 16 channels is internally reflected. In Fig. 6 , the scattered intensities for incident angle θ_i = 0° obtained with the Monte Carlo simulations and with the adding- doubling method are shown for the three selected wavelengths.

Fig. 6 Scattered intensities for θ_i = 0° obtained through Monte Carlo simulations (marks) and calculated with the adding-doubling method (lines).

Download Full Size | PDF

In Fig. 6 a good agreement between the results obtained from the Monte Carlo simulations and those calculated with the extended adding-doubling method can be observed. The average deviation between the two methods is 2.7%. However, the total time required for the Monte Carlo simulations, tracing 500 000 rays per wavelength, was over 3 hours. Using this number of rays per wavelength, the relative standard deviation of the scattered intensity distribution due to statistical noise was approximately 1%.

The adding-doubling calculations, using 32 channels, took approximately 25 seconds. The extended adding-doubling code was implemented in Matlab 7.10 and run on the same PC as the TracePro simulations.

In Fig. 7 the total hemispherical reflection and transmission is shown for each wavelength. The deviation between the Monte Carlo simulations and the adding-doubling calculations is found to be within 2%.

Fig. 7 Total reflection and transmission obtained with Monte Carlo simulations (marks) and adding-doubling calculations (lines)

Download Full Size | PDF

6. Conclusion

By extending the theory of the adding-doubling method, a fast and accurate method to obtain the spectral and angular transmission and reflection characteristics of a fluorescent slab has been developed. With this method, accurate reflection and transmittance profiles for infinite plane parallel slabs with finite thickness, optically flat surfaces and an azimuthally symmetric incident beam can be obtained with the same accuracy, but over 400 times faster than with Monte Carlo simulations.

This extended adding-doubling method can be used to get a quick yet accurate estimate of the transmission and reflection characteristics of a simplified sample geometry for fluorescent applications. Due to the fast algorithm, the extended adding-doubling for fluorescence is particularly suitable for optimization purposes.

Appendix A

The matrix equations for a particular emission wavelength interval for a fluorescent slab 0-1 are given by Eq. (35) and Eq. (36):

\begin{matrix} I_{1}^{+} = R_{10} I_{1}^{-} + T_{01} I_{0}^{+} + \\ \sum_{i = 1}^{N} R_{10}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{1}^{-} (Δ λ_{i}^{X}) + T_{01}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{0}^{+} (Δ λ_{i}^{X}) \end{matrix}

\begin{matrix} I_{0}^{-} = R_{01} I_{0}^{+} + T_{10} I_{1}^{-} \\ + \sum_{i = 1}^{N} R_{01}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{0}^{+} (Δ λ_{i}^{X}) + T_{10}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{1}^{-} (Δ λ_{i}^{X}) \end{matrix}

R_xy and T_xy are the reflection and transmission matrices at the selected emission wavelength, representing scattering and absorption. The dependence on wavelength has been omitted from the equations for readability reasons.

The matrices R^c_xy(Δλ^X_i, Δλ^M_j) and T^c_xy(Δλ^X_i, Δλ^M_j) represent the conversion of light from an excitation wavelength interval Δλ^X_i towards a particular emission wavelength interval Δλ^M_j for each incident and scattered conical segment. Similar as with the adding-doubling method at non-emission wavelengths, the same equations can be written for a slab 1-2, given by Eq. (37) and Eq. (38):

\begin{matrix} I_{2}^{+} = R_{21} I_{2}^{-} + T_{12} I_{1}^{+} \\ + \sum_{i = 1}^{N} R_{21}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{2}^{-} (Δ λ_{i}^{X}) + T_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{1}^{+} (Δ λ_{i}^{X}) \end{matrix}

\begin{matrix} I_{1}^{-} = R_{12} I_{1}^{+} + T_{21} I_{2}^{-} + \\ \sum_{i = 1}^{N} R_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{1}^{+} (Δ λ_{i}^{X}) + T_{21}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{2}^{-} (Δ λ_{i}^{X}) \end{matrix}

Inserting Eq. (38) in Eq. (36) results in Eq. (39):

\begin{matrix} I_{1}^{+} = {(E - R_{10} R_{12})}^{- 1} \\ {\begin{cases} R_{10} [T_{21} I_{2}^{-} + \sum_{i = 1}^{N} R_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{1}^{+} (Δ λ_{i}^{X}) + T_{21}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{2}^{-} (Δ λ_{i}^{X})] \\ + T_{01} I_{0}^{+} + \sum_{i = 1}^{N} R_{10}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{1}^{-} (Δ λ_{i}^{X}) + T_{01}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{0}^{+} (Δ λ_{i}^{X}) \end{cases}} \end{matrix}

I₁^-(Δλ^X_i) can be written in function of I₂^-(Δλ^X_i) and I₀⁺(Δλ^X_i) given by Eq. (40):

\begin{matrix} I_{1}^{-} (Δ λ_{i}^{X}) = {[E - R_{12} (Δ λ_{i}^{X}) R_{10} (Δ λ_{i}^{X})]}^{- 1} \\ [T_{21} (Δ λ_{i}^{X}) I_{2}^{-} (Δ λ_{i}^{X}) + R_{12} (Δ λ_{i}^{X}) T_{01} (Δ λ_{i}^{X}) I_{0}^{+} (Δ λ_{i}^{X})] \end{matrix}

R_xy(Δλ^X_i) and T_xy(Δλ^X_i) are the reflection and transmission matrices at excitation wavelength interval Δλ^X_i. In a similar way, I₁⁺(Δλ^X_i) can be written in function of I₂^-(Δλ^X_i) and I₀⁺(Δλ^X_i), given by Eq. (41):

\begin{matrix} I_{1}^{+} (Δ λ_{i}^{X}) = {[E - R_{10} (Δ λ_{i}^{X}) R_{12} (Δ λ_{i}^{X})]}^{- 1} \\ [T_{01} (Δ λ_{i}^{X}) I_{0}^{+} (Δ λ_{i}^{X}) + R_{10} (Δ λ_{i}^{X}) T_{21} (Δ λ_{i}^{X}) I_{2}^{-} (Δ λ_{i}^{X})] \end{matrix}

Inserting Eqs. (39), (40) and (41) in Eq. (37) results in Eq. (42):

\begin{matrix} I_{2}^{+} = R_{21} I_{2}^{-} \\ + T_{12} {(E - R_{10} R_{12})}^{- 1} [\begin{array}{l} R_{10} (T_{21} I_{2}^{-} + \sum_{i = 1}^{N} {\begin{cases} R_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{10} (Δ λ_{i}^{X}) R_{12} (Δ λ_{i}^{X})]}^{- 1} \\ [T_{01} (Δ λ_{i}^{X}) I_{0}^{+} (Δ λ_{i}^{X}) + R_{10} (Δ λ_{i}^{X}) T_{21} (Δ λ_{i}^{X}) I_{2}^{-} (Δ λ_{i}^{X})] \\ + T_{21}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{2}^{-} (Δ λ_{i}^{X}) \end{cases}}) \\ + T_{01} I_{0}^{+} + \sum_{i = 1}^{N} {\begin{cases} R_{10}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{12} (Δ λ_{i}^{X}) R_{10} (Δ λ_{i}^{X})]}^{- 1} \\ [T_{21} (Δ λ_{i}^{X}) I_{2}^{-} (Δ λ_{i}^{X}) + R_{12} (Δ λ_{i}^{X}) T_{01} (Δ λ_{i}^{X}) I_{0}^{+} (Δ λ_{i}^{X})] \\ + T_{01}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{0}^{+} (Δ λ_{i}^{X}) \end{cases}} \end{array}] \\ + {\sum_{i = 1}^{N} \begin{array}{l} R_{21}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) I_{2}^{-} (Δ λ_{i}^{X}) + \\ T_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{10} (Δ λ_{i}^{X}) R_{12} (Δ λ_{i}^{X})]}^{- 1} [T_{01} (Δ λ_{i}^{X}) I_{0}^{+} (Δ λ_{i}^{X}) + R_{10} (Δ λ_{i}^{X}) T_{21} (Δ λ_{i}^{X}) I_{2}^{-} (Δ λ_{i}^{X})] \end{array}} \end{matrix}

After further manipulation we obtain Eq. (43):

\begin{matrix} I_{2}^{+} = [R_{21} + T_{12} {(E - R_{10} R_{12})}^{- 1} R_{10} T_{21}] I_{2}^{-} \\ + [T_{12} {(E - R_{10} R_{12})}^{- 1} T_{01}] I_{0}^{+} \\ + \sum_{i = 1}^{N} [\begin{array}{l} T_{12} {(E - R_{10} R_{12})}^{- 1} (\begin{array}{l} R_{10} {\begin{cases} R_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{10} (Δ λ_{i}^{X}) R_{12} (Δ λ_{i}^{X})]}^{- 1} \\ R_{10} (Δ λ_{i}^{X}) T_{21} (Δ λ_{i}^{X}) + T_{21}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) \end{cases}} \\ + R_{10}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{12} (Δ λ_{i}^{X}) R_{10} (Δ λ_{i}^{X})]}^{- 1} T_{21} (Δ λ_{i}^{X}) \end{array}) \\ + R_{21}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) + T_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{10} (Δ λ_{i}^{X}) R_{12} (Δ λ_{i}^{X})]}^{- 1} R_{10} (Δ λ_{i}^{X}) T_{21} (Δ λ_{i}^{X}) \end{array}] I_{2}^{-} (Δ λ_{i}^{X}) \\ + \\ \sum_{i = 1}^{N} (\begin{array}{l} T_{12} {(E - R_{10} R_{12})}^{- 1} {\begin{cases} R_{10} R_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{10} (Δ λ_{i}^{X}) R_{12} (Δ λ_{i}^{X})]}^{- 1} T_{01} (Δ λ_{i}^{X}) \\ + R_{10}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{12} (Δ λ_{i}^{X}) R_{10} (Δ λ_{i}^{X})]}^{- 1} R_{12} (Δ λ_{i}^{X}) T_{01} (Δ λ_{i}^{X}) \\ + T_{01}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) \end{cases}} \\ + T_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{10} (Δ λ_{i}^{X}) R_{12} (Δ λ_{i}^{X})]}^{- 1} T_{01} (Δ λ_{i}^{X}) \end{array}) I_{0}^{+} (Δ λ_{i}^{X}) \end{matrix}

The R and T matrices at the selected emission wavelength for a layer consisting of two joined layers are finally given by Eqs. (44)-(47):

R_{20} = R_{21} + T_{12} {(E - R_{10} R_{12})}^{- 1} R_{10} T_{21}

T_{02} = T_{12} {(E - R_{10} R_{12})}^{- 1} T_{01}

\begin{matrix} R_{20}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) = \\ T_{12} {(E - R_{10} R_{12})}^{- 1} (\begin{array}{l} R_{10} {\begin{cases} R_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{10} (Δ λ_{i}^{X}) R_{12} (Δ λ_{i}^{X})]}^{- 1} \\ R_{10} (Δ λ_{i}^{X}) T_{21} (Δ λ_{i}^{X}) + T_{21}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) \end{cases}} \\ + R_{10}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{12} (Δ λ_{i}^{X}) R_{10} (Δ λ_{i}^{X})]}^{- 1} T_{21} (Δ λ_{i}^{X}) \end{array}) \\ + R_{21}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) + T_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{10} (Δ λ_{i}^{X}) R_{12} (Δ λ_{i}^{X})]}^{- 1} R_{10} (Δ λ_{i}^{X}) T_{21} (Δ λ_{i}^{X}) \end{matrix}

\begin{matrix} T_{02}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) = \\ T_{12} {(E - R_{10} R_{12})}^{- 1} {\begin{cases} R_{10} R_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{10} (Δ λ_{i}^{X}) R_{12} (Δ λ_{i}^{X})]}^{- 1} T_{01} (Δ λ_{i}^{X}) \\ + R_{10}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{12} (Δ λ_{i}^{X}) R_{10} (Δ λ_{i}^{X})]}^{- 1} R_{12} (Δ λ_{i}^{X}) T_{01} (Δ λ_{i}^{X}) \\ + T_{01}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) \end{cases}} \\ + T_{12}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) {[E - R_{10} (Δ λ_{i}^{X}) R_{12} (Δ λ_{i}^{X})]}^{- 1} T_{01} (Δ λ_{i}^{X}) \end{matrix}

Appendix B

The calculation of the initial matrices starts from the radiative transfer equations for light within a particular emission wavelength interval Δλ^M_j with central wavelength λ^M_j given by Eq. (48):

\begin{matrix} ν \frac{ϑ I (ν)}{ϑ z} = - (μ_{a} + μ_{s}) I (ν) + μ_{s} \int_{- 1}^{1} P (ν, ν') I (ν') d ν' \\ + w (Δ λ_{j}^{M}) \sum_{i = 1}^{N} {μ_{a} (Δ λ_{i}^{X}) Q E (Δ λ_{i}^{X}) \frac{λ_{j}^{M}}{λ_{i}^{X}} [\int_{- 1}^{1} \frac{1}{2} I (ν', Δ λ_{i}^{X}) d ν'] Δ λ_{i}^{X}} \end{matrix}

In Eq. (48) I(ν) is the channel flux at the selected emission wavelength interval, I(ν’, Δλ^X_i) represents the channel flux at the excitation wavelength interval Δλ^X_i (with central wavelength λ^X_i). In Eq. (48) µ_a, µ_s and P(ν,ν’) are respectively the absorption and scattering coefficient and phase function for the selected emission wavelength interval, µ_a(Δλ^X_i) is the absorption coefficient at excitation wavelength interval Δλ^X_i. QE(Δλ^X_i) expresses the quantum efficiency the excitation wavelength interval. The wavelength ratio is used to convert the number of photons to powerflux. A weight factor w(Δλ^M_j) has been to describe the fraction of light generated by an excitation at Δλ^X_i and emitted within the selected emission wavelength region.

The integrals in Eq. (48) are discretized using quadrature schemes for [0,1] with points 0 → µ_M and corresponding weights c_i, which results in Eq. (49):

\begin{matrix} \pm ν_{i} \frac{ϑ I (\pm ν_{i})}{ϑ z} = - (μ_{a} + μ_{s}) I (\pm ν_{i}) + μ_{s} \sum_{j = 1}^{M} p_{i j}^{+ -} c_{j} I (\mp ν_{j}) + p_{i j}^{+ +} c_{j} I (\pm ν_{j}) \\ + w (Δ λ_{j}^{M}) \sum_{i = 1}^{N} {\begin{matrix} μ_{a} (Δ λ_{i}^{X}) Q E (Δ λ_{i}^{X}) \frac{λ_{j}^{M}}{λ_{i}^{X}} \\ [\sum_{j = 1}^{M} \frac{1}{2} c_{j} I (\pm ν_{j}, Δ λ_{i}^{X}) + \frac{1}{2} c_{j} I (\mp ν_{j}, Δ λ_{i}^{X})] Δ λ_{i}^{X} \end{matrix}} \end{matrix}

In Eq. (49) the cosine ν can take values from 0 to 1, p is the phase function in matrix notation defined by Eq. (50):

p_{i j}^{+ \pm} = P (ν_{i}, \pm ν_{j})

Since Eq. (48) can be written for all propagation directions, it can be rewritten in vector-matrix notation given by Eq. (51):

\begin{matrix} \pm M \frac{ϑ I^{\pm}}{ϑ z} = - (μ_{a} + μ_{s}) I^{\pm} + μ_{s} (p^{+ +} c I^{\pm} + p^{+ -} c I^{\mp}) \\ + w (Δ λ_{j}^{M}) \sum_{i = 1}^{N} {\begin{cases} μ_{a} (Δ λ_{i}^{X}) Q E (Δ λ_{i}^{X}) \frac{λ_{j}^{M}}{λ_{i}^{X}} \frac{1}{2} \\ (p_{x} c I^{\pm} (Δ λ_{i}^{X}) + p_{x} c I^{\mp} (Δ λ_{i}^{X})) Δ λ_{i}^{X} \end{cases}} \end{matrix}

The matrices M, c and p_x are defined by respectively Eq. (52), Eq. (53) and Eq. (54):

M (i, j) = ν_{i} δ_{i j}

c (i, j) = c_{i} δ_{i j}

p_{x} (i, j) = 1

I⁺ and I^- represent a channel flux at the selected emission wavelength interval in vector notation with the angle between the propagation direction and the positive Z-axis respectively smaller and larger than 90°. I⁺(Δλ^X_i) and I^-(Δλ^X_i) represent the fluxes at excitation wavelength interval Δλ^X_i.

Integration Eq. (51) over a slab 0-1 with thickness Δd results in Eq. (55):

\begin{matrix} \pm M (I_{1}^{\pm} - I_{0}^{\pm}) = - (μ_{a} + μ_{s}) Δ d 〈 I_{\frac{1}{2}}^{\pm} 〉 + μ_{s} Δ d (p^{+ +} c 〈 I_{\frac{1}{2}}^{\pm} 〉 + p^{+ -} c 〈 I_{\frac{1}{2}}^{\mp} 〉) \\ + w (Δ λ_{j}^{M}) \sum_{i = 1}^{N} {\begin{cases} μ_{a} (Δ λ_{i}^{X}) Q E (Δ λ_{i}^{X}) \frac{λ_{j}^{M}}{λ_{i}^{X}} \frac{1}{2} Δ d \\ (p_{x} c 〈 I_{\frac{1}{2}}^{\pm} (Δ λ_{i}^{X}) 〉 + p_{x} c 〈 I_{\frac{1}{2}}^{\mp} (Δ λ_{i}^{X}) 〉) Δ λ_{i}^{X} \end{cases}} \end{matrix}

<I^±_1/2> is defined by Eq. (56):

〈 I_{\frac{1}{2}}^{\pm} 〉 = \frac{1}{Δ z} \int_{Z 0}^{Z 1} I^{\pm} d z

The “diamond initialization” (DI) is used as the thin layer approximation [28]. The DI assumes a linear decrease of the channel flux within the thin layer. Using the DI approximation, Eq. (55) can be rewritten as Eq. (57):

\begin{matrix} \pm M (I_{1}^{\pm} - I_{0}^{\pm}) = - (μ_{a} + μ_{s}) Δ d \frac{1}{2} (I_{0}^{\pm} + I_{1}^{\pm}) \\ + μ_{s} Δ d \frac{1}{2} [p^{+ +} c (I_{0}^{\pm} + I_{1}^{\pm}) + p^{+ -} c (I_{0}^{\mp} + I_{1}^{\mp})] \\ + w (Δ λ_{j}^{M}) \sum_{i = 1}^{N} (\begin{array}{l} μ_{a} (Δ λ_{i}^{X}) Q E (Δ λ_{i}^{X}) \frac{λ_{j}^{M}}{λ_{i}^{X}} \frac{1}{4} Δ d \\ {\begin{matrix} p_{x} c [I_{0}^{\pm} (Δ λ_{i}^{X}) + I_{1}^{\pm} (Δ λ_{i}^{X})] \\ + p_{x} c [I_{0}^{\mp} (Δ λ_{i}^{X}) + I_{1}^{\mp} (Δ λ_{i}^{X})] \end{matrix}} Δ λ_{i}^{X} \end{array}) \end{matrix}

Equation (57) can be rewritten as Eq. (58) and Eq. (59):

\begin{matrix} (E + α) I_{1}^{+} = (E - α) I_{0}^{+} + β (I_{0}^{-} + I_{1}^{-}) \\ + \sum_{i = 1}^{N} γ_{i} [I_{0}^{+} (Δ λ_{i}^{X}) + I_{1}^{+} (Δ λ_{i}^{X}) + I_{0}^{-} (Δ λ_{i}^{X}) + I_{1}^{-} (Δ λ_{i}^{X})] \end{matrix}

\begin{matrix} (E + α) I_{0}^{-} = (E - α) I_{1}^{-} + β (I_{0}^{+} + I_{1}^{+}) \\ + \sum_{i = 1}^{N} γ_{i} [I_{0}^{+} (Δ λ_{i}^{X}) + I_{1}^{+} (Δ λ_{i}^{X}) + I_{0}^{-} (Δ λ_{i}^{X}) + I_{1}^{-} (Δ λ_{i}^{X})] \end{matrix}

E represents the unity matrix and α, β and γ_i are defined by respectively Eq. (60), Eq. (61) and Eq. (62):

α = M^{- 1} (E - μ_{s} \frac{Δ d}{2} p^{+ +} c)

β = M^{- 1} μ_{s} \frac{Δ d}{2} p^{+ -} c

γ_{i} = w (Δ λ_{j}^{M}) μ_{a} (Δ λ_{i}^{X}) Q E (Δ λ_{i}^{X}) \frac{λ_{j}^{M}}{λ_{i}^{X}} \frac{Δ d}{4} p_{x} c

Solving Eq. (59) to I₀⁺ and inserting in Eq. (58) results in Eq. (63):

\begin{matrix} (E + α) I_{1}^{+} = (E - α) I_{0}^{+} \\ + β {\begin{cases} {(E + α)}^{- 1} (E - α) I_{1}^{-} + {(E + α)}^{- 1} β I_{0}^{+} + I_{1}^{+} \\ + {(E + α)}^{- 1} \sum_{i = 1}^{N} γ_{i} [\begin{array}{l} I_{0}^{+} (Δ λ_{i}^{X}) + I_{1}^{+} (Δ λ_{i}^{X}) \\ + I_{0}^{-} (Δ λ_{i}^{X}) + I_{1}^{-} (Δ λ_{i}^{X}) \end{array}] + I_{1}^{-} \end{cases}} \\ + \sum_{i = 1}^{N} γ_{i} [I_{0}^{+} (Δ λ_{i}^{X}) + I_{1}^{+} (Δ λ_{i}^{X}) + I_{0}^{-} (Δ λ_{i}^{X}) + I_{1}^{-} (Δ λ_{i}^{X})] \end{matrix}

To further solve this equation we use the relationships at excitation wavelength given by Eqs. (2) and (3). For the initial homogeneous layer, the reflection and transmission characteristics are independent on the irradiation direction, thus R₁₀=R₀₁ and T₁₀=T₀₁. Using this identities in Eq. (63), Eq. (64) is obtained:

\begin{matrix} I_{1}^{+} = (2 Γ - E) I_{0}^{+} + 2 Γ β {(E + α)}^{- 1} I_{1}^{-} \\ + \sum_{i = 1}^{N} {\begin{cases} [Γ β {(E + α)}^{- 1} + E] γ_{i} \\ [E + T (Δ λ_{i}^{X}) + R (Δ λ_{i}^{X})] [I_{0}^{+} (Δ λ_{i}^{X}) + I_{1}^{+} (Δ λ_{i}^{X})] \end{cases}} \end{matrix}

With Γ defined by Eq. (65):

Γ = {[E + α - β {(E + α)}^{- 1} β]}^{- 1}

The R and T matrices for a fluorescent thin “single-scatter” layer are thus given by Eqs. (66)-(68):

R_{10} = 2 Γ β {(E + α)}^{- 1}

T_{01} = 2 Γ - E

\begin{matrix} R_{10}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) = \\ T_{01}^{c} (Δ λ_{i}^{X}, Δ λ_{j}^{M}) = \\ [Γ β {(E + α)}^{- 1} + E] γ_{i} [E + T (Δ λ_{i}^{X}) + R (Δ λ_{i}^{X})] \end{matrix}

Acknowledgment

The authors would like to thank 'Strategic Initiative Materials' in Flanders (SIM) and the Institute for Innovation through Science and Technology in Flanders (IWT) for their financial support through the Sol-Cap project within the SIBO program. Wouter Saeys is funded as a postdoctoral fellow of the Research Foundation – Flanders (FWO).

References and links

1. J. Wu, M. S. Feld, and R. P. Rava, “Analytical model for extracting intrinsic fluorescence in turbid media,” Appl. Opt. 32(19), 3585–3595 (1993). [CrossRef] [PubMed]

2. A. Liebert, H. Wabnitz, N. Zołek, and R. Macdonald, “Monte Carlo algorithm for efficient simulation of time-resolved fluorescence in layered turbid media,” Opt. Express 16(17), 13188–13202 (2008). [CrossRef] [PubMed]

3. A. J. Welch, C. Gardner, R. Richards-Kortum, E. Chan, G. Criswell, J. Pfefer, and S. Warren, “Propagation of fluorescent light,” Lasers Surg. Med. 21(2), 166–178 (1997). [CrossRef] [PubMed]

4. I. Seo, J. Jung, B. J. Oh, and K. Whang, “Improvement of luminance and luminous efficacy of mercury-free, flat fluorescent lamp by optimizing phosphor profile,” IEEE Trans. Plasma Sci. 38(5), 1097–1100 (2010). [CrossRef]

5. J. P. You, N. T. Tran, Y. Lin, Y. He, and F. G. Shi, “Phosphor-concentration-dependent characteristics of white LEDs in different current regulation modes,” J. Electron. Mater. 38(6), 761–766 (2009). [CrossRef]

6. J. H. Park and J. H. Ko, “Optimization of the emitting structure of flat fluorescent lamps for LCD backlight applications,” J. Opt. Soc. Korea 11(3), 118–123 (2007). [CrossRef]

7. P. Chung, H. Chung, and P. H. Holloway, “Phosphor coatings to enhance Si photovoltaic cell performance,” J. Vac. Sci. Technol. A 25(1), 61–66 (2007). [CrossRef]

8. W. G. J. H. M. van Sark, “Enhancement of solar cell performance by employing planar spectral converters,” Appl. Phys. Lett. 87(15), 151117 (2005). [CrossRef]

9. E. Klampaftis and B. S. Richards, “Improvement in multi-crystalline silicon solar cell efficiency via addition of luminescent material to EVA encapsulation layer,” Prog. Photovolt. Res. Appl. 19(3), 345–351 (2011). [CrossRef]

10. D. Bera, S. Maslov, L. Qian, J. S. Yoo, and P. H. Holloway, “Optimization of the yellow phosphor concentration and layer thickness for down-conversion of blue to white light,” J. Disp. Technol. 6(12), 645–651 (2010). [CrossRef]

11. R. G. Young and E. G. F. Arnott, “The effect of phosphor coating weight on the lumen output of luorescent lamps,” J. Electrochem. Soc. 112(10), 982–984 (1965). [CrossRef]

12. W.-T. Chien, C.-C. Sun, and I. Moreno, “Precise optical model of multi-chip white LEDs,” Opt. Express 15(12), 7572–7577 (2007). [CrossRef] [PubMed]

13. Y. Shuai, N. T. Tran, and F. G. Shi, “Nonmonotonic phosphor size dependence of luminous efficacy for typical white LED emitters,” IEEE Photon. Technol. Lett. 23(9), 552–554 (2011). [CrossRef]

14. C. C. Chang, R.-L. Chern, C. C. Chang, C.-C. Chu, J. Y. Chi, J.-C. Su, I.-M. Chan, and J.-F. T. Wang, “Monte Carlo simulation of optical properties of phosphor-screened ultraviolet light in a white light-emitting device,” Jpn. J. Appl. Phys. 44(8), 6056–6061 (2005). [CrossRef]

15. D.-Y. Kang, E. Wu, and D.-M. Wang, “Modeling white light-emitting diodes with phosphor layers,” Appl. Phys. Lett. 89(23), 231102 (2006). [CrossRef]

16. J. Y. Chi, J.-S. Chen, C.-Y. Liu, C.-W. Chu, and K.-H. Chiang, “Phosphor converted LEDs with omni-directional-reflector coating,” Opt. Express 17(26), 23530–23535 (2009). [CrossRef] [PubMed]

17. J. Chen and X. Intes, “Comparison of Monte Carlo methods for fluorescence molecular tomography-computational efficiency,” Med. Phys. 38(10), 5788–5798 (2011). [CrossRef] [PubMed]

18. S. A. Prahl, M. J. C. van Gemert, and A. J. Welch, “Determining the optical properties of turbid mediaby using the adding-doubling method,” Appl. Opt. 32(4), 559–568 (1993). [CrossRef] [PubMed]

19. G. G. Stokes, “On the intensity of the light reflected from or transmitted through a pile of plates,” Proc. R. Soc. Lond. 11(0), 545–556 (1860). [CrossRef]

20. Z. Zhang, P. Yang, G. Kattawar, H.-L. Huang, T. Greenwald, J. Li, B. A. Baum, D. K. Zhou, and Y. Hu, “A fast infrared radiative transfer model based on the adding–doubling method for hyperspectral remote-sensing applications,” J. Quant. Spectrosc. Radiat. Transf. 105(2), 243–263 (2007). [CrossRef]

21. G. N. Plass, G. W. Kattawar, and F. E. Catchings, “Matrix operator theory of radiative transfer. 1: Rayleigh scattering,” Appl. Opt. 12(2), 314–329 (1973). [CrossRef] [PubMed]

22. P. J. Flatau and G. L. Stephens, “On the fundamental solution of the radiative transfer equation,” J. Geophys. Res. 93(D9), 11037–11050 (1988). [CrossRef]

23. H. G. Völz, Industrial color testing, Fundamentals and Techniques (Wiley-VCH, 2001)

24. S. N. Thennadil, “Relationship between the Kubelka-Munk scattering and radiative transfer coefficients,” J. Opt. Soc. Am. A 25(7), 1480–1485 (2008). [CrossRef] [PubMed]

25. W.-F. Cheong, S. A. Prahl, and A. J. Welch, “A review on the optical properties of biological tissues,” IEEE J. Quantum Electron. 26(12), 2166–2185 (1990). [CrossRef]

26. A. Rosema, W. Verhoef, J. Schroote, and J. F. H. Snel, “Simulating fluorescence light-canopy interaction in support of laser-induced fluorescence measurements,” Remote Sens. Environ. 37(2), 117–130 (1991). [CrossRef]

27. C. F. Bohren and D. Huffman, Absorption and Scattering of Light by Small Particles (John Wiley, 1983).

28. W. J. Wiscombe, “On initialization, error and flux conservation in the doubling method,” J. Quant. Spectrosc. Radiat. Transf. 16(8), 637–658 (1976). [CrossRef]

29. L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J. 93, 70–83 (1941). [CrossRef]

30. J. F. Beek, P. Blokland, P. Posthumus, M. Aalders, J. W. Pickering, H. J. C. M. Sterenborg, and M. J. C. Gemert, “In vitro double-integrating-sphere optical properties of tissues between 630 and 1064 nm,” Phys. Med. Biol. 42(11), 2255–2261 (1997). [CrossRef] [PubMed]

31. W. Saeys, M. A. Velazco-Roa, S. N. Thennadil, H. Ramon, and B. M. Nicolaï, “Optical properties of apple skin and flesh in the wavelength range from 350 to 2200 nm,” Appl. Opt. 47(7), 908–919 (2008). [CrossRef] [PubMed]

32. D. Yudovsky and L. Pilon, “Modeling the local excitation fluence rate and fluorescence emission in absorbing and strongly scattering multilayered media,” Appl. Opt. 49(31), 6072–6084 (2010). [CrossRef]

Extended adding-doubling method for fluorescent applications

Abstract

1. Introduction

2. General approach

3. Adding-doubling

3.1 Adding-doubling for non-fluorescent slabs

3.2 Adding-doubling for fluorescent slabs

4. Initial matrices

4.1 Initial matrices for non-fluorescent slabs

4.2 Initial matrices for fluorescent slabs

5. Validation

6. Conclusion

Appendix A

Appendix B

Acknowledgment

References and links

Cited By

Figures (7)

Equations (68)

Optics Express