OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 17, Iss. 26 — Dec. 21, 2009
  • pp: 23423–23442
« Show journal navigation

Implicitly causality enforced solution of multidimensional transient photon transport equation

Chintha C Handapangoda and Malin Premaratne  »View Author Affiliations


Optics Express, Vol. 17, Issue 26, pp. 23423-23442 (2009)
http://dx.doi.org/10.1364/OE.17.023423


View Full Text Article

Acrobat PDF (947 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

A novel method for solving the multidimensional transient photon transport equation for laser pulse propagation in biological tissue is presented. A Laguerre expansion is used to represent the time dependency of the incident short pulse. Owing to the intrinsic causal nature of Laguerre functions, our technique automatically always preserve the causality constrains of the transient signal. This expansion of the radiance using a Laguerre basis transforms the transient photon transport equation to the steady state version. The resulting equations are solved using the discrete ordinates method, using a finite volume approach. Therefore, our method enables one to handle general anisotropic, inhomogeneous media using a single formulation but with an added degree of flexibility owing to the ability to invoke higher-order approximations of discrete ordinate quadrature sets. Therefore, compared with existing strategies, this method offers the advantage of representing the intensity with a high accuracy thus minimizing numerical dispersion and false propagation errors. The application of the method to one, two and three dimensional geometries is provided.

© 2009 Optical Society of America

1. Introduction

Many models have been developed to solve the one-dimensional photon transport equation (PTE), which assumes horizontally uniform plane-parallel media. The modeling of the three-dimensional photon transport is usually considered difficult, even for the steady state case, because it involves solving integro-differential equations of four or five variables [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

]. However, many diverse methods for solving the multidimensional steady state photon transport problem, which discard the time-dependency complications of the transient problem considered here, have been found in literature [2

2. K. F. Evans, “The spherical harmonics discrete ordinate method for three-dimensional atmospheric radiative transfer,” J. Atmos. Sci. 55, 429–446 (1998). [CrossRef]

, 3

3. G. L. Stephens, “Radiative transfer through arbitrarily shaped optical media. part I: A general method of solution,” -J. Atmos. Sci. 45, 1818–1836 (1988). [CrossRef]

].

The use of short pulse lasers in engineering and biomedical applications, such as optical tomography and laser ablation has motivated research in transient photon transport through biological media [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

]. Transient signals have the possibility to probe specific resonances in biological media, which far exceed the capabilities rendered by steady state methods. Most of the existing transient models had been developed for the one-dimensional PTE; but there are a few models developed for the two-dimensional PTE [4

4. Z. M. Tan and P. F. Hsu, “Transient radiative transfer in three-dimensional homogeneous and non-homogeneous participating media,” J. Quant. Spectrosc. Radiat. Transfer 73, 181–194 (2002). [CrossRef]

]. Only recently have researchers started to work on models for the three-dimensional transient photon transport in scattering media such as biological tissue. Also, there has been renewed interest in having sensors implanted in tissue for monitoring and characterization [5

5. C. C. Handapangoda, M. Premaratne, D. M. Paganin, and P. R. D. S. Hendahewa, “Technique for handling wave propagation specific effects in biological tissue: Mapping of the photon transport equation to Maxwell’s equations,” Opt. Express 16, 17792–17807 (2008). [CrossRef] [PubMed]

10

10. C. C. Handapangoda, M. Premaratne, and D. M. Paganin, “Simulation of embedded photonic crystal structures for blood glucose measurement using Raman spectroscopy,” 2008 International conference on Nanoscience and Nanotechnology, (Melbourne, Australia, February 2008).

].

Many three-dimensional photon transport problems have been modeled by reducing the 3D propagation to equivalent 1D models. However, such mapping inevitably discard the inhomogeneities and anisotropic characteristics of such media. For example, realistic analysis of light propagation in biological tissue must be carried in a framework which accounts for both weak anisotropy and inhomogeneity. Moreover, if sensors are implanted in such scattering media, 1D models in general has difficulty accounting for both scattering and diffraction effects. Also, specific surface features of a model problem may be more accurately incorporated into threedimensional photon transport. Owing to these reasons, 3D models can provide quantitatively superior realistic results compared with its one or two dimensional counterparts.

Most of the recent developments on solving the three-dimensional transient PTE rely on approximate methods such as the diffusion approximation and spherical harmonics approximation, which are used to simplify the PTE [4

4. Z. M. Tan and P. F. Hsu, “Transient radiative transfer in three-dimensional homogeneous and non-homogeneous participating media,” J. Quant. Spectrosc. Radiat. Transfer 73, 181–194 (2002). [CrossRef]

]. These approximate methods have failed to predict the correct propagation speed within the medium and also the solution accuracy is not satisfactory except for specific cases that consider optically thick media at asymptotically longer time scales [4

4. Z. M. Tan and P. F. Hsu, “Transient radiative transfer in three-dimensional homogeneous and non-homogeneous participating media,” J. Quant. Spectrosc. Radiat. Transfer 73, 181–194 (2002). [CrossRef]

]. Tan et al. [4

4. Z. M. Tan and P. F. Hsu, “Transient radiative transfer in three-dimensional homogeneous and non-homogeneous participating media,” J. Quant. Spectrosc. Radiat. Transfer 73, 181–194 (2002). [CrossRef]

] have proposed an integral equation technique to model the transient radiative transfer in 3D homogeneous as well as non-homogeneous participating media. Guo et al. [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

] have formulated a complete transient three-dimensional discrete ordinates method to solve the transient PTE in a nonhomogeneous turbid medium with a rectangular enclosure. Guo et al. [11

11. Z. Guo and K. Kim, “Ultrafast-laser-radiation transfer in heterogeneous tissues with the discrete-ordinates method,” Appl. Opt. 42, 2897–2905 (2003). [CrossRef] [PubMed]

] have modeled heterogeneous biological tissues, using the discrete ordinates method, incorporating Fresnel specularly reflecting boundary condition. They treated the laser intensity as having a diffuse part and a specular part. The reflectivity at the tissue-air interface was calculated by using the Snell’s law and the Fresnel equations [11

11. Z. Guo and K. Kim, “Ultrafast-laser-radiation transfer in heterogeneous tissues with the discrete-ordinates method,” Appl. Opt. 42, 2897–2905 (2003). [CrossRef] [PubMed]

]. Chai et al. [12

12. J. C. Chai, P. F. Hsu, and Y. C. Lam, “Three-dimensional trnasient radiative transfer modeling using the finite-volume method,” J. Quant. Spectrosc. Radiat. Transfer 86, 299–313 (2004). [CrossRef]

] have proposed a finite-volume method to calculate transient radiative transfer in a three-dimensional enclosure. Sawetprawichkul et al. [13

13. A. Sawetprawichkul, P. F. Hsu, and K. Mitra, “Parallel computing of three-dimensional monte carlo simulation of transient radiative transfer in participating media,” in Proceedings of the 8th AIAA/ASME Joint Thermophysics and Heat Transfer Conferece, (American Institute of Aeronautics and Astronautics, St. Louse, Missouri, 2002), pp. 1–10.

] have implemented the Monte Carlo method for modeling the three-dimensional transient photon transport in a parallel computer system.

Oliveira et al. [14

14. J. V. P. de Oliveira, A. V. Cardona, M. T. Vilhena, and R. C. Barros, “A semi-analytical numerical method for time-dependent radiative transfer prolems in slab geometry with coherent isotropic scattering,” J. Quant. Spectrosc. Radiat. Transfer 73, 55–62 (2002). [CrossRef]

] have proposed a semi-analytical numerical method for solving the one-dimensional transient PTE in slab geometry. They represented the time dependence of the radiance and its time derivative using a truncated series of Laguerre polynomials. This results in a set of time-independent equations which are then solved recursively using a Laplace transform method in the space variable with analytical inversion.

In this paper, we propose a novel solution method for solving the three-dimensional transient PTE, which can also be easily adopted for solving one and two dimensional cases. This is an extension of the model we developed for the one-dimensional transient case in [15

15. C. C. Handapangoda, M. Premaratne, L. Yeo, and J. Friend, “Laguerre Runge-Kutta-Fehlberg method for simulating laser pulse propagation in biological tissue,” IEEE J. Sel. Top. Quantum Electron. 14, 105–112 (2008). [CrossRef]

]. Current work uses a transformation which eliminates the partial derivative term with respect to time. We then expand the radiance using a Laguerre basis which transforms the transient PTE to the steady state version, which is then solved using a finite volume approach with the discrete ordinates method (DOM) [16

16. M. F. Modest, “The method of discrete ordinates” in Radiative Heat Transfer, 2nd ed. (Academic Press, Boston, 2003), pp. 513–522.

]. With the use of a Laguerre basis we can represent any arbitrary pulse shape using only a few Laguerre polynomials, considerably reducing subsequent computational overhead. Also, the causality of the system is implicitly imposed by the Laguerre expansion. The use of the discrete ordinates method has several advantages. It requires only a single formulation to invoke higher-order approximations of discrete ordinate quadrature sets; and also, the DOM is applicable to the complete anisotropic scattering phase function and inhomogeneous media [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

].

This paper is organized as follows. In Section 2, we provide a detailed mathematical description about the proposed model. This section is further broken down into three subsections. The first subsection carries a step by step derivation of the proposed algorithm for the more general three-dimensional case. In the second and third subsections, descriptions on how the model reduces to two and one dimensional geometries are provided. In Section 3, we present results obtained by numerically implementing the proposed algorithm. First, a validation of the model is provided using the one-dimensional case, with special parameters for which the results can be intuitively derived. Then, a comparison of the results obtained for all the three dimensions is provided against the results obtained using the algorithm proposed by Guo et al. in [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

]. We conclude this paper in Section 4 after highlighting relative advantages of the proposed method.

2. Formulation

This section introduces the proposed solution technique. In the first sub-section, we present a detailed mathematical description of the method, as applied for the three-dimensional PTE. In the second sub-section, we describe how this method can be adopted to solve the two-dimensional PTE. In the last section, a description on how to adopt the proposed method for the one-dimensional case is presented.

2.1. For three-dimensional PTE

The three dimensional transient photon transport equation (PTE) is given by [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

]

1vtI(x,y,z,s,t)+ξxI(x,y,z,s,t)+ηyI(x,y,z,s,t)+μzI(x,y,z,s,t)
σs(x,y,z)4π4πP(s;s)I(x,y,s,t)ds+σt(x,y,z)I(x,y,z,s,t)=F(x,y,z,s,t),
(1)

where I(x,y, z, s, t) is the light intensity (radiance) with unitsW.m -2.sr -1.Hz -1, (x,y, z) are the cartesian coordinates, s is the local solid angle, µ, ξ and η are direction cosines such that µ=cosθ, ξ=sinθ cosϕ, η=sinθ sinϕ and t is time; σt (x,y,z) and σs(x,y,z) represent the position dependent attenuation and scattering coefficients, respectively. The speed of light in the medium is v, P(s′; s) is the phase function and F(x,y, z, s, t) is the source term. Since we consider a constant refractive index, the speed of light in the medium is a constant. In this paper we assume that there are no internal sources inside the medium and hence F=0. Figure 1 shows the incident direction, s, the scattered direction, s′, and the scattering angle. Figure 2 shows the relationship between the Cartesian and spherical coordinates, showing the zenith and azimuthal angles.

Fig. 1. Illustration of the incident direction, scattered direction and the scattering angle.
Fig. 2. Relationship between Cartesian and spherical coordinates.

Before proceeding to solve the PTE, the following substitution is used to map Eq. (1) to a moving reference frame with the pulse:

τ=tx3vξy3vηz3vμ.
(2)

Using Eq. (2) in Eq. (1), we are able to eliminate the partial derivative term with respect to time in the original PTE, as shown below.

1vτI(x,y,z,s,τ)ξ3vξτI(x,y,z,s,τ)+ξxI(x,y,z,s,τ)η3vητI(x,y,z,s,τ)+ηyI(x,y,s,τ)μ3vμτI(x,y,z,s,τ)+μzI(x,y,z,s,τ)σs(x,y,z)4π4πP(s;s)I(x,y,z,s,t)ds+σt(x,y,z)I(x,y,z,s,τ)=0.
(3)

Equation (3) can be simplified to get,

ξxI(x,y,z,s,τ)+ηyI(x,y,z,s,τ)+μzI(x,y,z,s,τ)σs(x,y,z)4π4πP(s;s)I(x,y,z,s,τ)ds+σt(x,y,z)I(x,y,z,s,τ)=0.
(4)

Then the time dependency is represented using Laguerre polynomials. We use Laguerre functions to expand I(x,y, z, s,τ) in a summation.

Laguerre polynomials are the canonical solutions of the differential equation

xy+(1x)y+ny=0,

which is also known as the Laguerre’s equation [17

17. M. Abramomitz and I. A. Stegun, “Orthogonal polynomials” in Handbook of mathematical functions with formulas, graphs and mathematical tables, (Dover, New York, 1965), pp. 773–784.

]. These solutions are a sequence of orthogonal polynomials which can be generated by Rodrigues’ formula [17

17. M. Abramomitz and I. A. Stegun, “Orthogonal polynomials” in Handbook of mathematical functions with formulas, graphs and mathematical tables, (Dover, New York, 1965), pp. 773–784.

]

Ln(x)=exn!dndxn(exxn).

These polynomials satisfy the orthogonality property

0Ln(x)Lm(x)exdx=δnm,

where δnm is the Kronecker delta [18

18. G. A. Korn and T. M. Korn, “Tensor analysis” in Mathematical handbook for scientists and engineers: definitions, theorems and formulas for reference and review, 2nd ed., (Dover, New York, 2000), pp. 544.

]. Laguerre polynomials are causal [19

19. C. J. Rivero-Moreno and S. Bres, “Video spatio-temporal signatures using polynomial transforms,” Lect. Notes Compt. Sci. 3736, 50–59 (2005). [CrossRef]

] and thus the causality of the system is implicitly imposed. Also, any causal function can be expanded using a Laguerre basis because it forms a complete orthogonal basis in real space [20

20. G. A. Korn and T. M. Korn, “Special functions” in Mathematical handbook for scientists and engineers: definitions, theorems and formulas for reference and review, 2nd ed., (Dover, New York, 2000), pp. 848–856.

]. Since such expansions are a linear superposition of causal Laguerre polynomials, the causality property of the original system is implicitly retained.

We expand the radiance, I, as a summation of N Laguerre polynomials;

I(x,y,z,s,τ)=k=0NBk(x,y,z,s)Lk(τ),
(5)

where ∫ 0 Ln (τ)Lm (τ)e =δmn. Using Eq. (5) in Eq. (4) and taking moments (i.e. multiplying by Ln (τ)e and integrating over [0,∞)), the time dependence can be removed resulting in N uncoupled equations, one for each Laguerre coefficient, Bn.

ξxBn(x,y,z,s)+ηyBn(x,y,z,s)+μzBn(x,y,z,s)σs(x,y,z)4π4πP(s;s)Bn(x,y,z,s)ds+σt(x,y,z)Bn(x,y,z,s)=0,
(6)

where n=1, …,N. Thus, the transient photon transport equation has been reduced to a set of uncoupled steady state photon transport equations. Equation (6) can be solved using the discrete ordinates method as detailed in [16

16. M. F. Modest, “The method of discrete ordinates” in Radiative Heat Transfer, 2nd ed. (Academic Press, Boston, 2003), pp. 513–522.

], and outlined below.

The solid angle dependence of Eq. (6) can be discretized using the discrete ordinates method [21

21. S. Chandrasekhar, “Quadrature formulae” in Radiative Transfer, (Dover, New York, 1960), pp. 54–60.

], resulting in L coupled equations for each Laguerre coefficient.

ξixBn(x,y,z,si)+ηiyBn(x,y,z,si)+μizBn(x,y,z,si)σs(x,y,z)4πj=1LwjP(sj;si)Bn(x,y,z,sj)+σt(x,y,z)Bn(x,y,z,si)=0,
(7)

where wj are quadrature weights, i=1, …,L and n=1, …,N. We divide the tissue specimen into infinitesimally small control volumes and apply the discrete ordinates method. We integrate Eq. (7) over a control volume:

VBnixdV=(Bni)xuAxu(Bni)xdAxd,
(8)

where VxΔyΔz is the volume of the control volume, (Bn)xu and (Bn)xd are average values of Bn over the faces Axu and Axd of the control volume, respectively. Similarly,

VBniydV=(Bni)yuAyu(Bni)ydAyd,
(9)
VBniydV=(Bni)zuAzu(Bni)zdAzd.
(10)

Here, Bin=Bn(x,y, z, si). Using Eq. (8), Eq. (9) and Eq. (10) in Eq. (7),

ξi((Bni)xuAxu(Bni)xdAxd)+ηi((Bni)yuAyu(Bni)ydAyd)+μi((Bni)zuAzu(Bni)zdAzd)
σs(x,y,z)4πVj=1LwjP(sj;si)(Bnj)p+σt(x,y,z)V(Bni)p=0,
(11)

(Bni)p=(γi)x(Bni)xd+(1(γi)x)(Bni)xu,
(12)
=(γi)y(Bni)yd+(1(γi)y)(Bni)yu,
(13)
=(γi)z(Bni)zd+(1(γi)z)(Bni)zu.
(14)

The weights γx, γy and γz can be determined using the scheme proposed by Lathrop [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

, 22

22. K. D. Lathrop, “Spatial differnecing of the transport equation: positive vs accuracy,” J. Comput. Phys. 4, 475–498 (1968). [CrossRef]

]:

(γi)x=max(0.5,1ξi/Δx2(ηi/Δy+μi/Δz)+σt),(γi)y=max(0.5,1ηi/Δy2(ξi/Δx+μi/Δz)+σt),(γi)z=max(0.5,1μi/Δz2(ξi/Δx+ηi/Δy)+σt),.

Using Eq. (12), Eq. (13) and Eq. (14) in Eq. (11) and eliminating (Bin)xd, (Bin)yd and (Bin)zd,

(Bni)p=V(Si)p+ξiAxu(Bni)xu/(γi)x+ηiAyu(Bni)yu/(γi)y+μiAzu(Bni)zu/(γi)zσtV+ξiAxu/(γi)x+ηiAyu/(γi)y+μiAzu/(γi)z,
(15)

where

(Si)p=σs(x,y,z)4πj=1LwjP(sj;sj)(Bnj)p.
(16)

Dividing the numerator and the denominator of the right hand side of Eq. (15) by V, we get

(Bni)p=(Si)p+(ξi/((γi)xΔx))(Bni)xu+(ηi/((γi)yΔy))(Bni)yu+(μi/((γi)zΔz))(Bni)zuσt+ξi/((γi)xΔx)+ηi/((γi)yΔy)+μi/((γi)zΔz).
(17)

In this work, the spatial differencing scheme suggested by Fiveland [23

23. W. A. Fiveland, “Three-dimensional radiative heat-transfer solutions by the discrete-ordinates method,” J. Thermophys. Heat Transfer 2, 309–316 (1988). [CrossRef]

] is used to determine the step sizes Δx, Δy and Δz. According to this scheme

Δx<ξiminσt(1(γi)x),Δy<ηiminσt(1(γi)y),Δz<μiminσt(1(γi)z)..

The algorithm of discrete ordinates for solving the PTE, which is reduced to the steady state version, given by Eq. (7) starts at corner control volumes. The (Bin)xu, (Bin)yu and (Bin)zu values for each corner control volume are given by the Laguerre coefficients of the boundary condition. Hence, the centre values, (Bin)p, can be determined from Eq. (17). Once (Bin)p is calculated (Bin)xd, (Bin)yd and (Bin)zd can be determined using Eq. (12) to Eq. (14). Also, (Si)p, which is initially assumed to be equal to zero, is updated using the calculated (Bin)p value from Eq. (16). (Bin)xd, (Bin)yd and (Bin)zd of the current control volume are equal to (Bin)xu, (Bin)yu and (Bin)zu of the adjacent control volumes, respectively. A more detailed description on this method can be found in [16

16. M. F. Modest, “The method of discrete ordinates” in Radiative Heat Transfer, 2nd ed. (Academic Press, Boston, 2003), pp. 513–522.

].

The overall algorithm proposed in this paper is summarized in ALGORITHM 1 below:

ALGORITHM 1

1. Use the transformation given by Eq. (2) in the transient PTE, Eq. (1).

2. Expand the radiance, I(x,y, z, s, t), using a Laguerre basis.

3. Expand the time-dependent boundary condition using a Laguerre basis.

4. Divide the specimen geometry into infinitesimally small control volumes (finite volume approach).

5. Apply the discrete ordinates method to Eq. (6) in order to solve for each Laguerre coefficient of the radiance, iteratively until convergence is achieved.

6. Combine the Laguerre coefficients to obtain the radiance values.

The main steps of the discrete ordinate method using the finite volume approach are listed in ALGORITHM 2 below:

1. Start at a corner control volume. (Bin)xu, (Bin)yu and (Bin)zu for this control volume are known from the boundary condition.

2. (Si)p is assumed to be zero initially for all control volumes.

3. Calculate (Bin)p of the current control volume from Eq. (17).

4. Determine (Bin)xd, (Bin)yd and (Bin)zd using Eq. (12) to Eq. (14).

5. The quantities calculated in step 4 become (Bin)xu, (Bin)yu and (Bin)zu for adjacent control volumes.

6. Repeat steps 3 to 5 until all control volumes are traversed.

8. Repeat steps 2 to 7 until all discrete ordinates corresponding to the current quadrant are covered.

10. Repeat steps 1 to 9 until all the corners are considered.

11. Update (Si)p for all discrete ordinates using Eq. (16).

12. Repeat steps 1 to 11 until (Si)p is converged.

13. Repeat steps 1 to 12 for all Laguerre coefficients of I.

2.2. For two-dimensional PTE

The two dimensional transient photon transport equation (PTE), for a medium without any internal sources, is given by

1vtI(x,z,s,t)+ξxI(x,z,s,t)+μxI(x,z,s,t)
σs(x,z)4π4πP(s;s)I(x,z,s,t)ds+σt(x,z)I(x,z,s,t)=0,
(18)

Instead of the transformation given by Eq. (2), we have to use the transformation

τ=tx2z2,
(19)

in order to make the partial derivative term with respect to time of Eq. (18) vanish. After using Eq. (2) in Eq. (18) to eliminate the time derivative term, the radiance, I, is expanded using a Laguerre basis using Eq. (5) and take moments as done for the three-dimensional case. This results in the two-dimensional steady state PTE:

ξxBn(x,z,s)+μzBn(x,z,s)σs(x,z)4π4πP(s:s)Bn(x,z,s)ds+σt(x,z)Bn(x,z,s)=0.
(20)

Then, ALGORITHM 2 can be used Eq. (20) to solve for Bn and using Eq. (5) the radiance, I, can be determined. However, when applying ALGORITHM 2, in the equations used for the three-dimensional case, the terms corresponding to the variable y vanish.

2.3. For one-dimensional PTE

The one dimensional transient photon transport equation (PTE), for a medium without any internal sources, is given by

1vtI(z,s,t)+μzI(z,s,t)σs(z)4π4πP(s;s)I(z,s,t)ds+σt(z)I(z,s,t)=0.
(21)

For this case, instead of the transformation given by Eq. (2), we have to use the transformation

τ=tz.
(22)

Use of Eq. (22) in Eq. (21) eliminates the time derivative term in Eq. (21). Then, expanding the transformed one-dimensional PTE using a Laguerre basis and taking moments we get,

μzBn(z,s)σs(z)4π4πP(s;s)Bn(z,s)ds+σt(z)Bn(z,s)=0,
(23)

where n=1, …,N. Then we discretize the solid angle, s, of Eq. (23) using the discrete ordinates method, which results in

μizBn(z,si)σs(z)4πj=1LwjP(sj;si)Bn(z,sj)+σt(z)Bn(z,si)=0,
(24)

where i=1, …,L and n=1, …,N. Thus, there are L coupled equations for each Laguerre coefficient, Bn and Eq. (24) corresponds to the i th discrete ordinate of s. We can write this set of equations in matrix form as;

zABn+σtBnσs4πPWBn=0,
(25)

where B n=[Bn(z, si] L,1, P=[P(s j;s i)] L,L, and, A is a L by L diagonal matrix with diagonal elements µi. The matrix W is also a L by L diagonal matrix with diagonal elements wj.

Rearranging, Eq. (25) can be written as

zBn=YBn,
(26)

where Y=A1[σs4πPWσtI].. Hence, the original transient PTE is reduced to a one-variable ordinary differential equation. The boundary condition should also be expanded using a Laguerre basis as described for the three-dimensional case. Thus, Eq. (26), subject to the given boundary condition, can be solved using the 4th order Runge-Kutta-Fehlberg(RKF) method [24

24. S. C. Chapra and R. P. Canale, “Runge-Kutta methods” in Numerical Methods for Engineers, 4th ed., (McGraw-Hill,New York,2002), pp. 719–720.

]. A more detailed description about solving the one-dimensional PTE using this technique (the Laguerre Runge-Kutta-Fehlberg method) can be found in our paper [15

15. C. C. Handapangoda, M. Premaratne, L. Yeo, and J. Friend, “Laguerre Runge-Kutta-Fehlberg method for simulating laser pulse propagation in biological tissue,” IEEE J. Sel. Top. Quantum Electron. 14, 105–112 (2008). [CrossRef]

]. In that paper, we have represented the solid angle using the zenith and azimuthal angles. However, the technique is essentially the same, the only difference being the matrix sizes.

3. Numerical results and discussion

We have simulated the proposed technique using Matlab. First, we validate the proposed technique for special cases, for which the results are intuitive and readily computable using the one-dimensional PTE. Then we provide a comparison for all the three cases; one, two and three dimensions, with results obtained from the technique proposed by Guo et al. [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

]. In [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

] they have provided a complete transient three-dimensional discrete ordinates method to solve the three-dimensional PTE in a rectangular enclosure. In that paper, they have extended the discrete ordinates algorithm described above in ALGORITHM 2, to include the transient term. We adopted their method for the two and one dimensional cases as well for the purpose of comparison. In the figures, we assigned the label “Laguerre DOM” for our method and “Transient DOM” for the method proposed in [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

]. In the simulations, all the units were normalized. A brief account on the normalization used is provided below.

Fig. 3. Incident Gaussian pulse (with I 0=1 in Eq. (27)).

3.1. Normalization of units

The PTE is linear in intensity, I, and thus arbitrary units for intensity can be used. Also, it is a first order, linear differential equation in time and thus time units can also be selected to enhance the numerical accuracy. It is very clear that input pulse parameters influence the numerical accuracy and efficiency of the algorithm. The input pulse was taken to be a Gaussian pulse described mathematically by

f(t)=I0e((tt0)T)2,
(27)

Fig. 4. (Color online) Short laser pulse incident on the biological tissue layer

In this paper, we have set Ts=T/1.5. We choose this normalization factor due to the fact that the Laguerre approximation of the Gaussian pulse is very accurate for pulses with T=1.5 or greater. This effect is shown in Fig. 5 and Fig. 6. Figure 5 presents a comparison of direct Laguerre approximation (i.e. without scaling with Ts) with the exact plot of the Gaussian pulse, for four different pulse widths. This shows that Laguerre approximation for pulses with T=1 and T=1.5 are very accurate. However, for the T=1 case, a very small bump appears for 6.5<x<7.5 (see Fig. 5b). Such inaccuracies can be always eliminated by choosing larger number of Laguerre coefficients. However, if the computational efficiency is also of concern, then better scaling can be used to reduce the required number of Laguerre coefficients. Therefore, we have chosen the T=1.5 and t 0=4 approximation. A pulse with different T and t 0 values should be scaled and shifted accordingly to obtain T=1.5 and t 0=4 before the Laguerre fit is carried out. The results obtained by the simulations should be scaled and shifted back (reversing this process). This idea is illustrated in Fig. 6. Four different Gaussian pulses with different T and t 0 values are scaled and shifted (up or down) to obtain the Gaussian pulse with T=1.5 and t 0=4 and a Laguerre fit is carried out. Then the scaling and shifting is reversed and the resulting fit is plotted with the exact plot of the corresponding pulse. It is clearly shown in Fig. 6 that this scaled Laguerre fit produces very accurate approximations.

Fig. 5. Comparison of Laguerre fitting, without scaling, with the exact plot of the corresponding Gaussian function.

We have set Zs=v×T̄. Here, T̄ can be chosen to suit the particular application. However, these scaling factors should be chosen carefully so that the matrices that are used remain well-conditioned. For the simulations presented in this paper, without loss of generality, we have chosen T/T̄=1.5 so that Zs=v×T1.5. In this paper we use a line over the letters to denote quantities with normalized units.

Fig. 6. Comparison of Laguerre fitting, after scaling with the factor Ts, with the exact plot of the corresponding Gaussian function.

3.2. Computational complexity

For the simulations we have considered an observation volume of -1<x̄<+1,-1<ȳ<+1 and 0<z̄<2. However, we have assumed refractive index-matched boundaries in all the three dimensions. Therefore, in order to minimize effects from internal reflections at the boundaries the simulation volume was taken to be larger than the observation volume in all the three dimensions (-5<x̄<+5,-5<ȳ<+5 and 0<z̄<5). The total number of voxels used in the finite volume computations was 108000. The number of voxels considered in the observation volume was 1728. The number of discrete ordinates taken for the simulations was 80. We have used the discrete ordinates for S 8 approximation given in Table 16.1 of reference [16

16. M. F. Modest, “The method of discrete ordinates” in Radiative Heat Transfer, 2nd ed. (Academic Press, Boston, 2003), pp. 513–522.

]. The two main drawbacks of the discrete ordinates method are false scattering and the ray effect [4

4. Z. M. Tan and P. F. Hsu, “Transient radiative transfer in three-dimensional homogeneous and non-homogeneous participating media,” J. Quant. Spectrosc. Radiat. Transfer 73, 181–194 (2002). [CrossRef]

, 16

16. M. F. Modest, “The method of discrete ordinates” in Radiative Heat Transfer, 2nd ed. (Academic Press, Boston, 2003), pp. 513–522.

]. False scattering is due to the spatial discretization errors. When a single collimated beam is traced through an enclosure by the discrete ordinates method, the beam will gradually widen as it moves farther away from its point of origin; and this unphysical effect, even in the absence of real scattering, is called false scattering [16

16. M. F. Modest, “The method of discrete ordinates” in Radiative Heat Transfer, 2nd ed. (Academic Press, Boston, 2003), pp. 513–522.

]. False scattering can be reduced by using finer meshes [16

16. M. F. Modest, “The method of discrete ordinates” in Radiative Heat Transfer, 2nd ed. (Academic Press, Boston, 2003), pp. 513–522.

]. The ray effect is due to the errors in angular discretization and can be reduced by increasing the sizes of the meshes [16

16. M. F. Modest, “The method of discrete ordinates” in Radiative Heat Transfer, 2nd ed. (Academic Press, Boston, 2003), pp. 513–522.

]. Therefore, if a finer spatial mesh is used to reduce the false scattering, a finer angular quadrature scheme should be used to reduce the ray effect [16]. Thus, if one uses the Laguerre DOM method proposed in this paper with a finer mesh, it should be accompanied by an increased number of discrete ordinates in order to minimize numerical errors.

The Gaussian pulse shape can be represented accurately using about 40 Laguerre polynomials. Increasing the number of polynomials up to 63 enables further refinement in the results. The incident Gaussian pulse in Fig. 3 can be approximated by 20 Laguerre polynomials with a normalized root mean squared deviation (NRMSD) of 8.07% in the observation window. The NRMSD can be reduced to 1.5% by using 40 polynomials. A further reduction up to 0.39% can be achieved by using 63 polynomials. However, increasing this number further does not improve the approximation because machine accuracy hinder further numerical improvements.

The simulation presented in this section was carried out in a computer with 2 Intel CPUs at 2.5 GHz. For the three-dimensional simulation of the proposed method (Laguerre DOM) shown in Fig. 13 to Fig. 16, with 63 Laguerre polynomials, it took 9 minutes and 1 second. For simulating the same problem using Transient DOM the same machine took 34 minutes and 6 seconds. In [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

] the authors state that their method may take several minutes to dozens of hours for simulations depending on the specified problem and the grid size. The efficiency of the proposed Laguerre DOM might be even better for more complex problems than the one considered here.

Figure 4 shows a short laser pulse of Gaussian shape incident on a biological tissue layer at (x̄=0,ȳ=0, z̄=0). This is the model used for the simulations presented in this paper. For all the simulations, we have used T̄=1.5 and t0¯=4 for the input pulse given by Eq. (27), the normalized velocity to be v̄=1 and the normalized thickness of the tissue layer, at which the results were obtained, to be z̄=2. The isotropic phase function (P(s′;s)=1) was used in the simulations. However, the proposed method is able to handle any other phase function.

Fig. 7. Variation of radiance with time, along the incident angle, at z̄=2 without any scattering or absorption.

Figures 7, 8 and 9 are used for validating the proposed technique for the one dimensional case. These figures were obtained for special cases for which the results are intuitive. For Fig. 7, the tissue medium was assumed to be lossless, without any scattering or absorption. That is, σs and σt in the PTE are set to zero. With this setup, there should be neither a decay in the intensity, nor scattering to other directions. Figure 7 shows the variation of radiance with time, along the direction of incidence (ie. µ=1), at z̄=2. As expected, at z̄=2, we have obtained the same Gaussian pulse, without any loss, but with the corresponding time delay.

Fig. 8. (Color online) Variation of radiance with time and µ, at z̄=2 with absorption but no scattering.

For Fig. 8, the medium was assumed to have a normalized absorption coefficient, σa¯, of 0.2 but without any scattering (hence, the normalized attenuation coefficient, σt¯=0.2 as well). The direction of incidence for this simulation was taken to be at µ=0.5. Figure 8 shows the variation of radiance with time and the cosine of the zenith angle, µ, at z=2. As expected, a decayed radiance profile exists along the incident direction, µ=0.5, but since there was no scattering no radiance values are seen along other directions.

Fig. 9. (Color online) Variation of radiance with time and µ, at z̄=2 with scattering and absorption.

For Fig. 9, the medium was assumed to have a normalized absorption coefficient, σa¯, of 0.02 and a normalized scattering coefficient, σs¯, of 0.98 (hence, σt¯=1). The direction of incidence for this simulation was taken to be at µ=0.17. Figure 9 shows the variation of radiance with time and µ. In this figure, the scattering to other directions is exhibited clearly.

Figures 10 to 18 show comparisons of the proposed method with the method proposed in [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

], for one, two and three dimensional cases, respectively. In figures 7 to 9, we have shown the radiance profile so that it is possible to explain the known results when there is no attenuation, when only absorption is present, and when both absorption and scattering is present. However, for the comparisons we have obtained irradiance profiles, which have more physical meaning, specially in practical measurements. For figures 10 to 18, σt¯=1, σs¯=0.98 and the direction of incidence is µ=1.

Fig. 10. (Color online) Comparison: Variation of irradiance with time for 1D PTE.

Figure 10 shows the variation of irradiance with time at z̄=2, for the one-dimensional PTE using the Laguerre DOM and the Transient DOM.

Fig. 11. (Color online) Variation of irradiance with time and x coordinate at z̄=2 for 2D PTE, using Laguerre DOM.

Figures 11 and 12 show the variation of irradiance with time and along x-axis at z̄=2, for the two-dimensional PTE using the two methods. In Fig. 11 one can observe that the time shift of the peak value of the irradiance profile increases as x̄ moves away from x̄=0. The reason for this is that light takes the shortest time to reach (x̄=0, z̄=2) point and the time increases as the point moves away from x̄=0 central axis. However, in the simulation obtained using Transient DOM shown in Fig. 12 this physical phenomenon is not clearly visible.

Figures 13 to 16 show the variation of intensity with time, along x and y axes, at z̄=2, for the three-dimensional PTE using Laguerre DOM. In these figures the value of the normalized irradiance (ranging from 0 to 16×10-4) is represented by the color scale. Figures 13 to 16 were obtained using the same data set. Figure 13 shows a slice plane at ȳ=0. Thus, in Fig. 13 the variation of irradiance with time along the x-axis on (ȳ=0, z̄=2) plane is shown. As expected, the light reaches (x̄=0,ȳ=0) point first as rays traveling at an angle to the central axis (x̄=0,ȳ=0) take longer to reach z̄=2 plane (because the diagonal distance is more than the central axis distance and the speed is constant). For the same reason the irradiance profile decays first at (x̄=0,ȳ=0) and later at other points. The more we move away from the central axis, the longer it takes the irradiance profile to decay.

Fig. 12. (Color online) Variation of irradiance with time and x coordinate at z̄=2 for 2D PTE, using Transient DOM.
Fig. 13. (Color online) Variation of irradiance with time, x coordinate and y coordinate at z̄=2 for 3D PTE, using Laguerre DOM (Slice plane at ȳ=0, z̄=2).

Figure 14 shows the variation of irradiance along x and y axes on z̄=2 plane at t̄=6. With isotropic scattering and normal incidence we have obtained a circular spatial distribution on the xy plane. From Fig. 13 it can be seen that at t̄=6 points around (x̄=0,ȳ=0) are receiving the second half of the gaussian pulse (i.e. the maximum irradiance has been reached earlier) and points further away from the central axis are receiving light corresponding to the maximum or the first half of the incident pulse. This phenomenon is reflected in Fig. 14. However, in the simulation result shown in Fig. 13 and Fig. 14, an uncharacteristically low number of photons can be seen along the central axis and the majority of diffuse photons appear off axis resulting in a well defined shape. The reason for this phenomenon is not very clear yet, and we intend to investigate this issue using Monte-Carlo techniques in the future. Figures 15 and 16 show that the irradiance profile at z̄=2 is symmetrical along x and y axes due to isotropic scattering and normal incidence.

Fig. 14. (Color online) Variation of irradiance with time, x coordinate and y coordinate at z̄=2 for 3D PTE, using Laguerre DOM (Slice plane at z̄=2, t̄=6).
Fig. 15. (Color online) Variation of irradiance with time, x coordinate and y coordinate at z̄=2 for 3D PTE, using Laguerre DOM.

Figures 17 and 18 show the variation of intensity with time, along x and y axes, at z̄=2, for the three-dimensional PTE using Transient DOM proposed in reference [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

]. The temporal profile in Fig. 17 obtained using this method is different to that in Fig. 13. Since we use a Gaussian shaped input, the irradiance profile should start decaying at (x̄=0,ȳ=0, z̄=2) before at other points. The further away the observation point moves from the central axis, the longer it should take for the decaying to start. This difference in Fig. 13 and Fig. 17 may be due to the strong numerical diffusion and false propagation in Transient DOM of Guo and Kumar [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

]. In [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

] the authors state that an obvious disadvantage of transient analysis using their method is that the numerical diffusion and false propagation are inevitable. In order to minimize the numerical diffusion, the spatial grid and the time step are required to be as fine as possible. The size of time step affects the transient behavior. Figure 3 and Fig. 4 in [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

] show influences of time step and spatial grid on the temporal behavior, with a comparison to Monte Carlo simulations. Figure 3 clearly shows that the Transient DOM cannot capture the abrupt rise of transmittance in the early time period. In Laguerre DOM, however, there is no discretization in timescale. We represent the entire time domain of the input pulse as a truncated series of Laguerre polynomials and then solve the PTE for each Laguerre coefficient separately. Also, using the substitution in Eq. (2) we map the PTE to a moving reference frame with the pulse.

Fig. 16. (Color online) Variation of irradiance with time, x coordinate and y coordinate at z̄=2 for 3D PTE, using Laguerre DOM.
Fig. 17. (Color online) Variation of irradiance with time, x coordinate and y coordinate at z̄=2 for 3D PTE, using Transient DOM (Slice plane at ȳ=0, z̄=2).
Fig. 18. (Color online) Variation of irradiance with time, x coordinate and y coordinate at z̄=2 for 3D PTE, using Transient DOM (Slice plane at z̄=2, t̄=6).

In [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

], the authors have used Duhamel’s superposition theorem, which can only be applied to linear systems, to incorporate the time-dependent boundary conditions. This theorem relates the solution of the problem subject to the time-dependent boundary condition to the solution subject to time-independent unit step input. The Laguerre DOM represents the time-dependency of the boundary condition using a Laguerre expansion. Then, a set of time-independent equations corresponding to each Laguerre coefficient is solved separately using DOM. Thus, one possible explanation for the difference in temporal and spatial irradiance profiles might be due to the strong numerical diffusion and false propagation in Transient DOM of Guo and Kumar [1

1. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

].

4. Conclusion

This paper introduced a novel technique for modeling multidimensional transient photon transport, for applications in bio-sensing and in short pulse propagation through turbid media. The proposed method uses a transformation to eliminate the transient term in the transient PTE. It then expands the radiance using a Laguerre basis to automatically account for causality while providing an efficient basis suitable for calculations. This reduces the original transient PTE to a steady state version. Hence, the discrete ordinates method, using a finite volume approach can be used to solve for the radiance.

Since the time dependence is expanded using a Laguerre basis all the sampling points in the time domain is obtained in a single execution, as opposed to time marching techniques used in existing solution methods. This makes the proposed algorithm much faster when one requires the intensity profile at a particular point or a plane over a time interval. In addition, the use of the Laguerre expansion to represent the time dependency enables modeling the system with any arbitrary input pulse shape, using only a few Laguerre polynomials. Specifically, the Gaussian pulse shape, which is used in many practical applications can be accurately represented using a few Laguerre polynomials, as opposed to the discrete sampling used in most of the existing models. Also, this expansion implicitly imposes the causality of the system.

References and links

1.

Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys. Heat Transfer 16, 289–296 (2002). [CrossRef]

2.

K. F. Evans, “The spherical harmonics discrete ordinate method for three-dimensional atmospheric radiative transfer,” J. Atmos. Sci. 55, 429–446 (1998). [CrossRef]

3.

G. L. Stephens, “Radiative transfer through arbitrarily shaped optical media. part I: A general method of solution,” -J. Atmos. Sci. 45, 1818–1836 (1988). [CrossRef]

4.

Z. M. Tan and P. F. Hsu, “Transient radiative transfer in three-dimensional homogeneous and non-homogeneous participating media,” J. Quant. Spectrosc. Radiat. Transfer 73, 181–194 (2002). [CrossRef]

5.

C. C. Handapangoda, M. Premaratne, D. M. Paganin, and P. R. D. S. Hendahewa, “Technique for handling wave propagation specific effects in biological tissue: Mapping of the photon transport equation to Maxwell’s equations,” Opt. Express 16, 17792–17807 (2008). [CrossRef] [PubMed]

6.

K. W. Johnson, J. J. Mastrototaro, D. C. Howey, R. L. Brunelle, P. L. Burden-Brady, N. A. Bryan, C. C. Andrew, H. M. Rowe, D. J. Allen, B. W. Noffke, W. C. McMahan, R. J. Morff, D. Lipson, and R. S. Nevin, “In vivo evaluation of an electroenzymatic glucose sensor implanted in subcutaneous tissue,” Biosens. Bioelectron. 7, 709–714 (1992). [CrossRef] [PubMed]

7.

D. Moatti-Sirat, F. Capron, V. Poitout, G. Reach, D. S. Bindra, Y. Zhang, G. S. Wilson, and D. R. Thevenot, “Towards continuous glucose monitoring: in vivo evaluation of a miniaturized glucose sensor implanted for several days in rat subcutaneous tissue,” Diabetologia 35, 224–230(1992). [CrossRef] [PubMed]

8.

V. Poitout, D. Moatti-Sirat, G. Reach, Y. Zhang, G. S. Wilson, F. Lemonnier, and J. C. Klein, “A glucose monitoring system for on line estimation in man of blood glucose concentration using a miniaturized glucose sensor implanted in the subcutaneous tissue and a wearable control unit,” Diabetologia 36, 658–663 (1993). [CrossRef] [PubMed]

9.

C. C. Handapangoda, M. Premaratne, and D. M. Paganin, “Simulation of a device concept for noninvasive sensing of blood glucose levels,” in Information and Automation for Sustainability, 2007. ICIAFS 2007. Third International Conference, (Melbourne, Australia, December2007), pp. 31–34.

10.

C. C. Handapangoda, M. Premaratne, and D. M. Paganin, “Simulation of embedded photonic crystal structures for blood glucose measurement using Raman spectroscopy,” 2008 International conference on Nanoscience and Nanotechnology, (Melbourne, Australia, February 2008).

11.

Z. Guo and K. Kim, “Ultrafast-laser-radiation transfer in heterogeneous tissues with the discrete-ordinates method,” Appl. Opt. 42, 2897–2905 (2003). [CrossRef] [PubMed]

12.

J. C. Chai, P. F. Hsu, and Y. C. Lam, “Three-dimensional trnasient radiative transfer modeling using the finite-volume method,” J. Quant. Spectrosc. Radiat. Transfer 86, 299–313 (2004). [CrossRef]

13.

A. Sawetprawichkul, P. F. Hsu, and K. Mitra, “Parallel computing of three-dimensional monte carlo simulation of transient radiative transfer in participating media,” in Proceedings of the 8th AIAA/ASME Joint Thermophysics and Heat Transfer Conferece, (American Institute of Aeronautics and Astronautics, St. Louse, Missouri, 2002), pp. 1–10.

14.

J. V. P. de Oliveira, A. V. Cardona, M. T. Vilhena, and R. C. Barros, “A semi-analytical numerical method for time-dependent radiative transfer prolems in slab geometry with coherent isotropic scattering,” J. Quant. Spectrosc. Radiat. Transfer 73, 55–62 (2002). [CrossRef]

15.

C. C. Handapangoda, M. Premaratne, L. Yeo, and J. Friend, “Laguerre Runge-Kutta-Fehlberg method for simulating laser pulse propagation in biological tissue,” IEEE J. Sel. Top. Quantum Electron. 14, 105–112 (2008). [CrossRef]

16.

M. F. Modest, “The method of discrete ordinates” in Radiative Heat Transfer, 2nd ed. (Academic Press, Boston, 2003), pp. 513–522.

17.

M. Abramomitz and I. A. Stegun, “Orthogonal polynomials” in Handbook of mathematical functions with formulas, graphs and mathematical tables, (Dover, New York, 1965), pp. 773–784.

18.

G. A. Korn and T. M. Korn, “Tensor analysis” in Mathematical handbook for scientists and engineers: definitions, theorems and formulas for reference and review, 2nd ed., (Dover, New York, 2000), pp. 544.

19.

C. J. Rivero-Moreno and S. Bres, “Video spatio-temporal signatures using polynomial transforms,” Lect. Notes Compt. Sci. 3736, 50–59 (2005). [CrossRef]

20.

G. A. Korn and T. M. Korn, “Special functions” in Mathematical handbook for scientists and engineers: definitions, theorems and formulas for reference and review, 2nd ed., (Dover, New York, 2000), pp. 848–856.

21.

S. Chandrasekhar, “Quadrature formulae” in Radiative Transfer, (Dover, New York, 1960), pp. 54–60.

22.

K. D. Lathrop, “Spatial differnecing of the transport equation: positive vs accuracy,” J. Comput. Phys. 4, 475–498 (1968). [CrossRef]

23.

W. A. Fiveland, “Three-dimensional radiative heat-transfer solutions by the discrete-ordinates method,” J. Thermophys. Heat Transfer 2, 309–316 (1988). [CrossRef]

24.

S. C. Chapra and R. P. Canale, “Runge-Kutta methods” in Numerical Methods for Engineers, 4th ed., (McGraw-Hill,New York,2002), pp. 719–720.

OCIS Codes
(170.4580) Medical optics and biotechnology : Optical diagnostics for medicine
(170.6930) Medical optics and biotechnology : Tissue

ToC Category:
Medical Optics and Biotechnology

History
Original Manuscript: September 30, 2009
Revised Manuscript: October 29, 2009
Manuscript Accepted: December 4, 2009
Published: December 7, 2009

Virtual Issues
Vol. 5, Iss. 1 Virtual Journal for Biomedical Optics

Citation
Chintha C. Handapangoda and Malin Premaratne, "Implicitly causality enforced solution of multidimensional transient photon transport equation," Opt. Express 17, 23423-23442 (2009)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-26-23423


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. Z. Guo, and S. Kumar, "Three-dimensional discrete ordinates method in transient radiative transfer," J. Thermophys. Heat Transfer 16, 289-296 (2002). [CrossRef]
  2. K. F. Evans, "The spherical harmonics discrete ordinate method for three-dimensional atmospheric radiative transfer," J. Atmos. Sci. 55, 429-446 (1998). [CrossRef]
  3. G. L. Stephens, "Radiative transfer through arbitrarily shaped optical media. part I: A general method of solution," J. Atmos. Sci. 45, 1818-1836 (1988). [CrossRef]
  4. Z. M. Tan, and P. F. Hsu, "Transient radiative transfer in three-dimensional homogeneous and non-homogeneous participating media," J. Quantum Spectrosc. Radiat. Transfer 73, 181-194 (2002). [CrossRef]
  5. C. C. Handapangoda, M. Premaratne, D. M. Paganin, and P. R. D. S. Hendahewa, "Technique for handling wave propagation specific effects in biological tissue: Mapping of the photon transport equation to Maxwell’s equations," Opt. Express 16, 17792-17807 (2008). [CrossRef] [PubMed]
  6. K. W. Johnson, J. J. Mastrototaro, D. C. Howey, R. L. Brunelle, P. L. Burden-Brady, N. A. Bryan, C. C. Andrew, H. M. Rowe, D. J. Allen, B. W. Noffke, W. C. McMahan, R. J. Morff, D. Lipson, and R. S. Nevin, "In vivo evaluation of an electroenzymatic glucose sensor implanted in subcutaneous tissue," Biosens. Bioelectron. 7, 709-714 (1992). [CrossRef] [PubMed]
  7. D. Moatti-Sirat, F. Capron, V. Poitout, G. Reach, D. S. Bindra, Y. Zhang, G. S. Wilson, and D. R. Thevenot, "Towards continuous glucose monitoring: in vivo evaluation of a miniaturized glucose sensor implanted for several days in rat subcutaneous tissue," Diabetologia 35, 224-230 (1992). [CrossRef] [PubMed]
  8. V. Poitout, D. Moatti-Sirat, G. Reach, Y. Zhang, G. S. Wilson, F. Lemonnier, and J. C. Klein, "A glucose monitoring system for on line estimation in man of blood glucose concentration using a miniaturized glucose sensor implanted in the subcutaneous tissue and a wearable control unit," Diabetologia 36, 658-663 (1993). [CrossRef] [PubMed]
  9. C. C. Handapangoda, M. Premaratne, and D. M. Paganin, "Simulation of a device concept for noninvasive sensing of blood glucose levels," in Information and Automation for Sustainability, 2007. ICIAFS 2007. Third International Conference, (Melbourne, Australia, December 2007), pp. 31-34.
  10. C. C. Handapangoda, M. Premaratne, and D. M. Paganin, "Simulation of embedded photonic crystal structures for blood glucose measurement using Raman spectroscopy," 2008 International conference on Nanoscience and Nanotechnology, (Melbourne, Australia, February 2008).
  11. Z. Guo, and K. Kim, "Ultrafast-laser-radiation transfer in heterogeneous tissues with the discrete-ordinates method," Appl. Opt. 42, 2897-2905 (2003). [CrossRef] [PubMed]
  12. J. C. Chai, P. F. Hsu, and Y. C. Lam, "Three-dimensional transient radiative transfer modeling using the finite volume method," J. Quantum Spectrosc. Radiat. Transfer 86, 299-313 (2004). [CrossRef]
  13. A. Sawetprawichkul, P. F. Hsu, and K. Mitra, "Parallel computing of three-dimensional monte carlo simulation of transient radiative transfer in participating media," in Proceedings of the 8th AIAA/ASME Joint Thermophysics and Heat Transfer Conferece, (American Institute of Aeronautics and Astronautics, St. Louse, Missouri, 2002), pp. 1-10.
  14. J. V. P. de Oliveira, A. V. Cardona, M. T. Vilhena, and R. C. Barros, "A semi-analytical numerical method for time dependent radiative transfer prolems in slab geometry with coherent isotropic scattering," J. Quantum Spectrosc. Radiat. Transfer 73, 55-62 (2002). [CrossRef]
  15. C. C. Handapangoda, M. Premaratne, L. Yeo, and J. Friend, "Laguerre Runge-Kutta-Fehlberg method for simulating laser pulse propagation in biological tissue," IEEE J. Sel. Top. Quantum Electron. 14, 105-112 (2008). [CrossRef]
  16. M. F. Modest, "The method of discrete ordinates" in Radiative Heat Transfer, 2nd ed. (Academic Press, Boston, 2003), pp. 513-522.
  17. M. Abramomitz, and I. A. Stegun, "Orthogonal polynomials" in Handbook of mathematical functions with formulas, graphs and mathematical tables, (Dover, New York, 1965), pp. 773-784.
  18. G. A. Korn, and T. M. Korn, "Tensor analysis" in Mathematical handbook for scientists and engineers: definitions, theorems and formulas for reference and review, 2nd ed., (Dover, New York, 2000), pp. 544.
  19. C. J. Rivero-Moreno, and S. Bres, "Video spatio-temporal signatures using polynomial transforms," Lect. Notes Compt. Sci. 3736, 50-59 (2005). [CrossRef]
  20. G. A. Korn, and T. M. Korn, "Special functions" in Mathematical handbook for scientists and engineers: definitions, theorems and formulas for reference and review, 2nd ed., (Dover, New York, 2000), pp. 848-856.
  21. S. Chandrasekhar, "Quadrature formulae" in Radiative Transfer, (Dover, New York, 1960), pp. 54-60.
  22. K. D. Lathrop, "Spatial differencing of the transport equation: positive vs accuracy," J. Comput. Phys. 4, 475-498 (1968). [CrossRef]
  23. W. A. Fiveland, "Three-dimensional radiative heat-transfer solutions by the discrete-ordinates method," J. Thermophys. Heat Transfer 2, 309-316 (1988). [CrossRef]
  24. S. C. Chapra, and R. P. Canale, "Runge-Kutta methods" in Numerical Methods for Engineers, 4th ed., (McGraw-Hill, New York,2002), pp. 719-720.

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited