1.

Introduction

The transport of light in strongly scattering turbid media, such as biological tissue, is generally modeled as a diffusion process, which is described by the diffusion equation for the fluence $\phi$.^1,2 In condensed scattering media, $\phi$ is dependent on the refractive index $n$, such that $\phi /n^2$ is conserved, enabling modeling of the diffusion of light in media with a spatially varying refractive index.³ However, a finite discontinuity in the refractive index, such as that which occurs at the boundary between two scattering media of differing refractive indices, gives rise to Fresnel reflection, such that $\phi /n^2$ is discontinuous, with a discontinuity $\mathrm{\Delta }\phi$ proportional to the diffuse radiant flux $J$ at the boundary:^4–7

The discontinuity in $\phi /n^2$ depends on the ratio of refractive indices at the boundary and is small for modest index ratios. Monte Carlo simulation of diffuse light transport across boundaries between turbid media with different refractive indices has shown that the error introduced when this correction is ignored is generally less than 10% for weakly absorbing scattering media, e.g., biological tissue illuminated with infrared radiation.⁶ An analytical solution for time-dependent diffusion between adjacent half-spaces presented previously⁸ supports this conclusion. However, the error is cumulative in multilayer scattering media and increases significantly for strongly forward-biased scattering and non-negligible absorption. A simple analytic method of assessing the error incurred in these circumstances when the discontinuity in $\phi /n^2$ is not taken into account would therefore be useful and is discussed below. Errors in modeling the diffuse fluence $\phi$ in turbid media can lead to systematic errors in diffuse transmittance and reflectance. They also introduce errors in scattering and absorption coefficients inferred from reflectance measurements and thus to errors in quantitative image reconstruction via diffuse optical tomography. Errors in estimating internal diffuse fluence may also impact on photodynamic therapy.

2.

Theory

2.1.

Interface Error

Using the definitions previously published,⁹

Eq. (2a)

$$\phi =1/2\int I(\xi )\mathrm{d}\xi \xi \in |-1,1|,$$

Eq. (2b)

$$\mathrm{J}=1/2\int I(\xi )\mathrm{d}\xi \xi \in |-1,1|,$$

where $I(\xi )$ is the angular intensity distribution and $\xi$ is the direction cosine with respect to the positive $z$-axis, we find the mean cosine of the radiant intensity distribution $\langle \xi \rangle =\int I(\xi )\xi \mathrm{d}\xi /\int I(\xi )\mathrm{d}\xi =J/\phi$. Upon dividing Eq. (1) by ${\phi }_1$, we find the relative error

Eq. (3a)

$$\mathrm{\Delta }\phi /{\phi }_1={(n_2/n_1)}^2-{\phi }_2/{\phi }_1=C(n_2/n_1)J_1/{\phi }_1=C(n_2/n_1)\langle {\xi }_1\rangle .$$

Similarly,

Eq. (3b)

$$\mathrm{\Delta }\phi /{\phi }_2={(n_2/n_1)}^2-{\phi }_1/{\phi }_2=C(n_2/n_1)J_2/{\phi }_2=C(n_2/n_1)\langle {\xi }_2\rangle .$$

Thus the error in applying the approximate boundary condition is directly proportional to the mean cosine $\langle \xi \rangle$ of the angular intensity distribution at the boundary. This is illustrated in Fig. 1, where the ratio $J/\phi$ versus refractive index ratio $n$ for a given relative error $\mathrm{\Delta }\phi /\phi$ (left-hand axis) and the relative error $\mathrm{\Delta }\phi /\phi$ versus $n$ for a given ratio $J/\phi$ (right-hand axis) are plotted. Thus the larger the index ratio $n$, the smaller the ratio $J/\phi$ required for a given relative error $\mathrm{\Delta }\phi /\phi$ and the larger the relative error $\mathrm{\Delta }\phi /\phi$ for a given ratio $J/\phi$ i.e., for a given angular radiance distribution $I(\xi )$ with mean cosine $\langle \xi \rangle$.

Fig. 1

The ratio $J/\phi$ versus refractive index ratio $n$ for a given relative error in diffuse fluence $\mathrm{\Delta }\phi /\phi$(left-hand axis) and the relative error $\mathrm{\Delta }\phi /\phi$ versus $n$ for a given ratio $J/\phi$ (right-hand axis).

The discontinuity in diffuse fluence $\phi$ quantified by Eq. (1) also implies a discontinuity in the mean cosine $\langle \xi \rangle$ of the radiance distribution $I(\xi )$ at the boundary, viz., $\langle {\xi }_1\rangle \ne \langle {\xi }_2\rangle$. To find the magnitude of the error in $\phi$ for a specific case requires numerical evaluation of the boundary fluxes.⁶ However, an estimate (lower-bound) can be made in terms of the mean cosine of the asymptotic angular radiance ${\langle \xi \rangle }_{\mathrm{as}}$ given by the following equation published previously:¹⁰

Eq. (4)

$${\langle \xi \rangle }_{\mathrm{as}}=(1-a)/{\mu }_{\mathrm{as}}$$

for scattering albedo $a={\mu }_s/({\mu }_s+{\mu }_s)={\mu }_s/{\mu }_e$, where ${\mu }_s$ is the scattering coefficient, ${\mu }_a$ is the absorption coefficient, ${\mu }_e={\mu }_s+{\mu }_a$ is the extinction coefficient, and ${\mu }_{\mathrm{as}}$ is the asymptotic attenuation coefficient. In the $\delta \text{-}{P}_1$ approximation,² ${\mu }_{\mathrm{as}}$ is replaced by ${\mu }_{\mathrm{eff}}$, the effective attenuation coefficient. ${\mu }_{\mathrm{eff}}$ and the mean cosine ${\langle {\xi }^{\prime }\rangle }_{\mathrm{as}}$ are determined by the reduced scattering coefficient ${\mu }_s^{\prime }$ and the absorption coefficient ${\mu }_a$, as follows:

Eq. (5a)

$${\mu }_{\mathrm{eff}}={[{\mu }_a/D^{\prime }]}^{1/2}={[3{\mu }_a({\mu }_a+{\mu }_s^{\prime })]}^{1/2}={(3{\mu }_a{\mu }_{\mathrm{tr}})}^{1/2}$$

Eq. (5b)

$${\langle {\xi }^{\prime }\rangle }_{\mathrm{as}}={\mu }_{\mathrm{eff}}{D}^{\prime }={[{\mu }_a/3({\mu }_a+{\mu }_s^{\prime })]}^{1/2}={({\mu }_a/3{\mu }_{\mathrm{tr}})}^{1/2},$$

where ${\mu }_s^{\prime }={\mu }_s(1-g)$, $g$ is the scattering asymmetry, the transport coefficient ${\mu }_{\mathrm{tr}}={\mu }_a+{\mu }_s^{\prime }$ and $D=1/3{\mu }_{\mathrm{tr}}$ is the diffusion coefficient.¹ Thus ${\mu }_{\mathrm{eff}}\Rightarrow 0$, ${\langle {\xi }^{\prime }\rangle }_{\mathrm{as}}$ (isotropic radiance) when ${\mu }_a\Rightarrow 0$ (zero absorption). More precise evaluation of the effective attenuation coefficient ${\mu }_{\mathrm{eff}}$ (and of ${\langle {\xi }^{\prime }\rangle }_{\mathrm{as}}$, which is required for forward-biased scattering in absorbing media, involves higher moments of the phase function.¹¹ In the $\delta \text{-}{P}_3$ approximation, the effective attenuation coefficient is as follows:

Eq. (6a)

$${\mu }_{\mathrm{eff}}^{(3)}={[(\beta -{({\beta }^2-36\gamma )}^{1/2})/18]}^{1/2},$$

where $\beta =27{\mu }_a{{\mu }_t^{\prime }}^{(1)}+28{\mu }_a{{\mu }_t^{\prime }}^{(3)}+35{{\mu }_t^{\prime }}^{(2)}{{\mu }_t^{\prime }}^{(3)}$, $\gamma =105{\mu }_a{{\mu }_t^{\prime }}^{(1)}{{\mu }_t^{\prime }}^{(2)}{{\mu }_t^{\prime }}^{(3)}$, ${{\mu }_t^{\prime }}^{(m)}={\mu }_a+{\mu }_s(1-g^m)$, $m=1$, 2, 3^12,13 and the mean cosine of the asymptotic radiance

Eq. (6b)

$${\langle {\xi }^{\prime }\rangle }_{\mathrm{eff}}^{(3)}=(1-a^{\prime })/{\mu }_{\mathrm{eff}}^{(3)}\mathrm{.}$$

The dependence of the relative error $\mathrm{\Delta }\phi /\phi$ on scattering asymmetry $g$ is shown in Fig. 2 for scattering albedoes in the range $a\in |0,0.99|$ for accurate values of ${\mu }_{\mathrm{eff}}$.¹¹ It can be seen that $\mathrm{\Delta }\phi /\phi$ is only weakly dependent on scattering asymmetry for $g<0$ (backward-biased scattering), even for strong absorption ($a=0.2$, i.e., ${\mu }_a=4{\mu }_s$), while increasing rapidly for forward-biased scattering ($g>0$), approaching 10% for $g\ge 0.99$ when $a=0.9$. In the $\delta \text{-}P1$ approximation, the scattering asymmetry is reduced: $g^{\prime }\in |0,0.5|$ for $g\in |0,1|$, but so is the scattering albedo: $a^{\prime }={\mu }_s^{\prime }/({\mu }_s^{\prime }+{\mu }_a)$ via the reduced scattering coefficient ${\mu }_s^{\prime }={\mu }_s(1-g)$, potentially offsetting a reduced error in $\phi$. The error in diffuse fluence increases for interfaces with higher index ratios, $\mathrm{\Delta }\phi /\phi$ exceeding 20% for $g=0.95$ when $n=1.25$.

Fig. 2

Relative error in diffuse fluence $\mathrm{\Delta }\phi /\phi$ versus scattering asymmetry $g$ for a range of scattering albedoes $a$, with the error increasing rapidly for $g>0.9$. The limiting case for zero scattering ($a=0$) is indicated by the horizontal dashed line. The vertical dotted line marks the maximum asymmetry in the $\delta \text{-}{P}_1$ approximation ($g^{\prime }=0.5$). Index ratio $n=1.06$.

3.

Results

3.1.

Interface

The boundary condition Eq. (9a) has been evaluated analytically by Shendeleva for time-dependent diffusion in adjoining media with isotropic scattering, and it has been validated by Monte Carlo simulation.³ Validation of the analytic method presented herein is provided by comparison with the results of Ripoll and Nieto-Vesperinas, who evaluated the error using numerical methods.⁶ Figure 3 shows the relative errors in the diffuse fluences at the common boundary between two adjoining media versus the refractive index $n_2$ of a scattering medium adjoining an aqueous scattering medium ($n_1=1.333$), as calculated from the analytic formulae given above [Eqs. (12a), (12b), (12c), (12c′), (12d)]. The results show precise agreement with the numerical data (plotted points) published previously,⁶ confirming the validity of the simpler analytic method, which can therefore replace the previous numerical methods for rapid evaluation of the error in similar cases.

Fig. 3

Comparison of analytic results for diffuse fluence error $\mathrm{\Delta }\phi /\phi$ (curves) with numerical results (data points from Fig. 8 in Ref. 6) at the interface between an aqueous scattering medium ($n_1=1.333$) and a scattering medium with refractive index $n_2$ varied in the range of $1\le n_2\le 2$. Scattering parameters: ${\mu }_{s1}^{\prime }=15{\mathrm{cm}}^{-1}$, ${\mu }_{a1}=0.035{\mathrm{cm}}^{-1}$, ${\mu }_{s2}^{\prime }=10{\mathrm{cm}}^{-1}$, ${\mu }_{a2}=0.24{\mathrm{cm}}^{-1}$, $g=0.8$.

Figure 4 shows the dependence of the fractional errors $\mathrm{\Delta }\phi /\phi$ in the diffuse fluences versus scattering albedo $a$ for two adjoining media with disparate scattering parameters ($g=0$, 0.95), calculated in the $P_1$ and $P_3$ approximations to the diffusion parameters for Henyey-Greenstein scattering,¹⁶ with refractive index ratio $n=1.41/1.34=1.06$ (tissue/aq). Initially, the error increases rapidly with absorption ($a<1$) for $g=0$, reaching a broad maximum $\mathrm{\Delta }\phi /\phi \sim 6\%$ when $a\sim 0.6$; in contrast, the error for $g=0.95$ in the adjoining medium increases quasilinearly to $\sim 5\%$ when $a=0$. The results show that the $P_1$ approximation seriously underestimates the error in diffuse fluence (by $\sim 40\%$), whereas $P_3$ is $\le 10\%$ low compared with the accurate value and is preferred for analytic evaluation of $\mathrm{\Delta }\phi /\phi$. Overall, the error increases sharply when there is non-negligible absorption in a turbid medium with strongly forward-biased scattering.

Fig. 4

Relative error in diffuse fluence $\mathrm{\Delta }\phi /\phi$ versus scattering albedo $a$ on either side of the interface between two homogeneous media with disparate scattering parameters: $g=0$ (upper curves), $g=0.95$ (lower curves) for index ratio $n=1.06$. The filled squares (▪) are data points calculated with accurate values of the diffusion parameters $D$, $\kappa$.^10,11

3.2.

Multilayers

For diffusion of light in multiple layers of finite thickness, the diffuse fluence ${\phi }_k(z)$ in the $k$th layer may be expressed as previously described:¹⁵

Eq. (13)

$${\phi }_k(z)=a_k\exp ({\mu }_{\mathrm{eff}\text{-}k}z)+b_k\exp (-{\mu }_{\mathrm{eff}\text{-}k}z)$$

with a similar expression for ${\phi }_{k+1}(z)$ in the ($k+1$)th layer. The boundary conditions described previously,^4,7

Eq. (14a)

$${(n_{k+1}/n_k)}^2{\phi }_k(z_k)-{\phi }_{k+1}(z_k)=C(n_{k+1}/n_k)J_k(z_k)$$

Eq. (14b)

$$-D_k\nabla {\phi }_k(z_k)=-D_{k+1}\nabla {\phi }_{k+1}(z_k),$$

yield the simple recurrence relations (for $n=n_{k+1}/n_k>1$),

Eq. (15a)

$$a_{k+1}=1/2\{[n^2+1+C(n){D\mu }_{\mathrm{eff}}]a_k+[n^2-1-C(n){D\mu }_{\mathrm{eff}}]b_k\exp (-2{\mu }_{\mathrm{eff}}h)\}$$

Eq. (15b)

$$b_{k+1}=1/2\{[n^2-1+C(n){D\mu }_{\mathrm{eff}}]a_k\exp (2{\mu }_{\mathrm{eff}}h)+[n^2+1-C(n){D\mu }_{\mathrm{eff}}]b_k\},$$

when $D_{k+1}{\mu }_{\mathrm{eff}\text{-}k+1}=D_k{\mu }_{\mathrm{eff}\text{-}k}=D{\mu }_{\mathrm{eff}}$, and ${\mu }_{\mathrm{eff}\text{-}k+1}{h}_{k+1}={\mu }_{\mathrm{eff}\text{-}k}{h}_k={\mu }_{\mathrm{eff}}h$, where $h_k$, $h_{k+1}$ are the widths of the $k$th and ($k+1$)th layers, enabling the coefficients $a_{k+1}$, $b_{k+1}$ to be related $a_k$, $b_k$. The results for the approximate boundary condition [Eq. (10)] are obtained upon setting $C(n)=0$ in Eqs. (15a), (15b). Successive application of these relations yields the coefficients $a_k$, $b_k$ for all the layers involved, with appropriate boundary conditions chosen for the first and last.¹⁵ A parallel set of coefficients, $a_k^{\prime }$, $b_k^{\prime }$, for $C(n)=0$, enables direct comparison of the accurate and approximate fluxes ${\phi }_k$, ${\phi }_k^{\prime }$ in each layer and thus evaluation of the cumulative error for the multilayer system. This is illustrated in Fig. 5, with $\mathrm{\Delta }\phi /\phi =1.4\%$ at a single interface (for $a=0.995$, $g=0.95$, $n=1.1$) and the cumulative error increasing with the total number of layers, exceeding 30% for five layers when ${\mu }_s^{\prime }/{\mu }_a=1$. For multilayer media with higher index ratios, or for a larger number of layers, the cumulative error can easily exceed 100%.

Fig. 5

Relative error in diffuse fluence $\mathrm{\Delta }\phi /\phi$ versus ${\mu }_s^{\prime }/{\mu }_a$ in layer 1 (upper points) and layer 5 (lower points) of a five-layer medium on a half-space ($g=0.95$, index ratios $n=1.1$): open squares (□) represent $P_1$ (diffusion) values, filled squares (▪) represent accurate values of diffusion parameters $D$, $\kappa$.

3.3.

Perturbing Layer

The dependence of the diffuse reflectance of a layered medium on changes in the optical properties of a subsurface layer is of special interest.¹⁷ To illustrate this, consider a simple three-layer model comprising two-plane parallel layers supported on a semi-infinite layer (half-space) and vary the properties of the middle layer. The arrangement is sketched in Fig. 6. The optical properties are given in Table 1. The problem is analyzed via the equations for the diffuse fluences in the three regions:

Eq. (16a)

$$\text{Half-space}{\phi }_3(z)=b_3\exp (-{\mu }_{\mathrm{eff}\text{-}3}z)$$

Eq. (16b)

$$\text{Mid-layer}{\phi }_2(z)=a_2\exp ({\mu }_{\mathrm{eff}\text{-}2}z)+b_2\exp (-{\mu }_{\mathrm{eff}\text{-}2}z)$$

Eq. (16c)

$$\text{Surface-layer}{\phi }_1(z)=a_1\exp ({\mu }_{\mathrm{eff}\text{-}1}z)+b_1\exp ({\mu }_{\mathrm{eff}\text{-}1}z)-S_0/D\exp (-{\mu }_{\mathrm{tr}}z),$$

with the diffuse fluence ${\phi }_3(z)$ decaying exponentially in the half-space, the flux ${\phi }_2(z)$ in the midlayer having both exponential terms, and the flux ${\phi }_1(z)$ in the surface layer including the exponential source $S_0\exp (-{\mu }_{\mathrm{tr}}z)$ for plane parallel illumination at the surface.² The boundary conditions defined in Eqs. (14a) and (14b) are applied successively at the interfaces to evaluate the coefficients $a_k$, $b_k$ for $k=1$, 2, 3. The flux ${\phi }_1(z)$ in the surface layer is extrapolated to zero a finite distance beyond the surface equal to the linear extrapolation distance $z_b$ determined by the refractive index $n_1$.⁴ The radiant flux at the surface $J_s=-D_{1}d{\phi }_1/dz{|}_s=-D_1{\phi }_1(0)/z_b$. A set of coefficients $a_k^{\prime }$, $b_k^{\prime }$ for $C(n_{k+1}/n_k)=0$ at each interface yields the uncorrected surface flux $J_s^{\prime }$, allowing the relative error $\mathrm{\Delta }{J}_s/J_s$ (and hence the relative error in diffuse reflectance) to be determined. A representative set of curves showing the dependence of $\mathrm{\Delta }{J}_s/J_s$ on the ratio ${\mu }_s^{\prime }/{\mu }_a$ for selected values of the refractive index ratio $n=n_{21}$ is plotted in Fig. 7, with the relative error in surface flux $\mathrm{\Delta }{J}_s/J_s$ increasing with $n$ and ${\mu }_a$, reaching 30% for $n=1.5$ and ${\mu }_s^{\prime }/{\mu }_a=1(a^{\prime }=0.5)$, i.e., for ${\mu }_s/{\mu }_a=20$ ($a=0.95$, $g=0.95$). In general, the error is larger (or smaller) for higher (or lower) values of scattering asymmetry $g$ and also increases with absorption.

Fig. 6

Schematic of three-layer system comprising two finite layers on a semi-infinite substrate, subject to plane illumination normally incident on the first layer.

Fig. 7

Relative error in diffuse radiant flux $\mathrm{\Delta }{J}_s/J_s$ at the surface of a three-layer system (two finite layers on a semi-infinite substrate) versus the ratio ${\mu }_s^{\prime }/{\mu }_a$ in the middle layer (layer 2) with refractive index in the range $n=1.1$ to 1.5 (Table 1).

Table 1

Optical properties of the three-layer system.

Layer	${\mu }_s^{\prime }$/${\mathrm{mm}}^{-1}$	${\mu }_a/{\mathrm{mm}}^{-1}$	$d/\mathrm{mm}$	$n$
1	1	0.001	3	1.0
2	1	0.001 to 1.0	3	1.1 to 1.5
3	1	0.001	Inf.	1.0

4.

Discussion

The analytic method of asymptotic planar fluxes enables straightforward evaluation of the error in diffuse fluence $\phi$ without recourse to a Monte Carlo simulation. The magnitude of the error is readily found from the optical properties of the adjoining scattering media, namely, the refractive index ratio, the reduced scattering coefficients, the absorption coefficients, and scattering asymmetries. It is simply expressed via the product $C(n)\langle \xi \rangle$, where $C(n)$ is a monotonically increasing function of the index ratio $n=n_2/n_1$ (Ref. 4) and $\langle \xi \rangle$ is the mean cosine of the boundary radiance, a key result of the above analysis. This is approximated by the asymptotic mean cosine ${\langle \xi \rangle }_{\mathrm{as}}$ of the radiance far from the boundary, expressed in terms of scattering albedo $a$ and diffuse attenuation coefficient ${\mu }_{\mathrm{eff}}$. An equivalent formula for the error is $C(n){D\mu }_{\mathrm{eff}}$, where $D$ is the diffusion coefficient. The dependence of flux error on scattering asymmetry $g$ is of some interest, viz., for forward-biased scattering in turbid media with near-negligible absorption (typical of biological tissue in the near-IR⁶). The error increases rapidly as $g\Rightarrow 1$, but is virtually independent of scattering asymmetry for negative values $g<0$ (Fig. 2). For forward-biased scattering media with non-negligible absorption, e.g., biological media in the visible spectrum, the error becomes progressively less dependent on scattering asymmetry as absorption increases, ultimately becoming independent of $g$ in the limit of zero scattering albedo $a\Rightarrow 0$.

Having found a simple means of estimating the magnitude of the flux error at a boundary, accurate formulae for the diffuse fluence $\phi$ in adjoining media incorporating the correction term $C(n){D\mu }_{\mathrm{eff}}$ are obtained from the boundary condition Eq. (1). It is noted that the correction applies to the diffuse fluence distribution throughout the turbid medium, not simply at the boundaries. In the case of multilayer media, repeated application of the boundary conditions yields recurrence relations for the coefficients $a_k$, $b_k$ of the formula for the diffuse fluence. A corresponding set of coefficients $a_k^{\prime }$, $b_k^{\prime }$ obtained from the approximate boundary condition [with $C(n)=0$] enables (1) the flux error in successive layers to be obtained and (2) its parametric dependence on the optical properties of the layers to be investigated. This is of importance for comparison with experimental determination of diffuse reflectance, as well as for the inverse problem of determining optical constants from reflectance measurements.

The principal aim of the present work is to provide a simple analytic method for estimating the error incurred in using the approximate form of the boundary condition [Eq. (10)]. This has been checked against the computational results reported previously⁶ (Sec. 3.1 and Fig. 3) and illustrated with several examples relevant to biomedical optics. The work presented herein concerns analysis of the error involved in applying the approximate diffusion boundary condition [Eq. (10)] rather than the error in using diffusion theory per se. Thus a numerical evaluation would merely quantify the “error within the error,” whereas the analytic method provides a simple means of estimating its magnitude. The analytic approach was never intended to replace accurate radiative transfer computations where these are merited, e.g., Phillips and Jacques,¹⁷ but rather as a simple check on the diffusion approximation, e.g., the widely used $\delta \text{-}{P}_1$ formulation.²

5.

Conclusion

A simple analytic method of estimating the error involved in applying a commonly used approximate boundary condition for diffuse radiation in two adjoining scattering media with differing refractive indices has been presented. The method is based on asymptotic planar fluxes and enables the relative error to be readily evaluated without recourse to Monte Carlo simulation. Three examples of its application have been considered: (1) evaluating the error in calculating the diffuse fluences at a boundary between two media with differing refractive indices and dissimilar scattering properties, (2) the dependence of the relative error in diffuse fluence $\phi$ in a multilayer medium with discontinuous refractive index on the ratio of the reduced scattering coefficient to the absorption coefficient ${\mu }_s^{\prime }/{\mu }_a$, and (3) the dependence of the relative error in radiant flux $J_s$ at the surface of a three-layer medium on the ratio ${\mu }_s^{\prime }/{\mu }_a$ in the middle layer. In addition to its dependence on refractive index ratio $n=n_2/n_1$ via the function $C(n)$, the fluence error increases with scattering asymmetry $g$ in forward-biased scattering media and is cumulative in multilayered media with refractive index increments between layers. The same methodology can be applied to multilayer problems with cylindrical symmetry, with the system being first converted to planar geometry via a Hankel transform, to allow 1-D analysis (as here), followed by reconversion of the solution to cylindrical symmetry via an inverse Hankel transform.¹⁵