Simulation study on the detection of size, shape and position of three different scatterers using Non-standard time domain time inverse Maxwell’s algorithm

Kisalaya Chakrabarti; James B. Cole

doi:10.1364/OE.18.004148

1. Introduction

In medical diagnosis, detection of malignant tumors of various size and shape from back scattered waves [1] with the application of the Inverse theorem is widespread. The challenge lies how well this method can detect the tumors when the refractive index of the malignant cells and the normal cells differs by a small quantity. In this paper we have investigated three different reconstructive spectrums of three different shapes of scatterer (circle, rectangle and triangle) using NSFD-TI algorithm based on Non-Standard Finite Difference [2,4] and also we have observed the effect of refractive index upon the size of the circular scatterer.

2. Model and Calculation Methods

Maxwell’s equations for a non-conducting dielectric are given by:

\partial_{t} B (x, t) = - \nabla \times E (x, t)

\partial_{t} D (x, t) = \nabla \times H (x, t) - J (x, t)

\nabla \cdot B (x, t) = 0

\nabla \cdot D (x, t) = ρ (x)

where, H is the magnetic field, B is the magnetic induction, E the electric field, D the electric displacement, J the electric current density, and ρ is the density of free charge. Materials, such as iron or glass, modify the electric and magnetic fields. B describes how the material affects the magnetic field, and D describes how it affects electric field. The vector position is x =

(x, y, z)

,

\partial_{t} = \partial / \partial t

, and

\nabla = (\partial_{x}, \partial_{y}, \partial_{z})

is a vector differential operator. The curl of a vector is defined by

\nabla \times V = (\partial_{y} V_{z} - \partial_{z} V_{y}, \partial_{z} V_{x} - \partial_{x} V_{z}, \partial_{x} V_{y} - \partial_{y} V_{z})

, and the divergence of V is

\nabla \cdot V = \partial_{x} V_{x} + \partial_{y} V_{y} + \partial_{z} V_{z}

, where

V = (V_{x}, V_{y}, V_{z})

.

In most optical materials there is no free charge and $ρ = 0$ . In a non-conducting dielectric, such as glass, no electric current flow so J = 0. Also to an excellent approximation

B (x, t) = μ H (x, t)

where, μ, the magnetic permeability, is a constant. In general, the electrical permittivity, ε depends upon the frequency, but when this dependence is weak, it is a good approximation to take

D (x, t) = ε (x) E (x, t)

The permittivity is different for different materials, however so it is a function of position. For example, in glass

ε ≅ 1.2

. Assuming that Eq. (5) and Eq. (6) hold, Maxwell’s equations for a non-magnetic, non-conducting medium with no free charges or electrical currents is

μ \partial_{t} H (x, t) = - \nabla \times E (x, t)

ε (x) \partial_{t} E (x, t) = \nabla \times H (x, t)

Now, let us introduce what is called a Finite Difference (FD) approximation. The first derivative is approximately given by

f^{'} (x) ≅ \frac{f (x + \frac{h}{2}) - f (x - \frac{h}{2})}{h}

For convenience let us define the difference operator

d_{x} f (x) = f (x + \frac{h}{2}) - f (x - \frac{h}{2})

with this definition the FD approximation to

f^{'}

is

f^{'} (x) ≅ \frac{d_{x}}{h} f (x)

Defining the vector difference operator:

d = (d_{x}, d_{y}, d_{z})

and,

\nabla f (x) ≅ \frac{d}{h} f (x)

where,

Δ x = Δ y = Δ z = h

. The second derivative

f^{″}

can be approximated by applying Eq. (11) twice:

f^{″} (x) ≅ \frac{d_{x}^{2}}{h^{2}} f (x) .

and from the definition of

d_{x}

it is easy to show that

d_{x}^{2} = d_{x} d_{x}

is given by:

d_{x}^{2} f (x) = f (x + h) + f (x - h) - 2 f (x)

Defining,

d^{2} = d_{x}^{2} + d_{y}^{2} + d_{z}^{2}

We now have,

\nabla^{2} f (x) ≅ \frac{d^{2}}{h^{2}} f (x)

and also replacing the derivatives of Maxwell’s Eq. (7) by FD approximations we obtain

d_{t} H (x, t) = - \frac{1}{μ} \frac{Δ t}{h} d \times E (x, t)

Here the electromagnetic field components are arranged on the grid in such a way those central FDs can be taken. Expanding

d_{t} H (x, t)

, and solving for

H (x, t + \frac{Δ t}{2})

gives

H (x, t + \frac{Δ t}{2}) = H (x, t - \frac{Δ t}{2}) - \frac{1}{μ} \frac{Δ t}{h} d \times E (x, t)

Doing the same for Eq. (8) we have

d_{t} E (x, t + \frac{Δ t}{2}) = - \frac{1}{ε (x)} \frac{Δ t}{h} d \times H (x, t + \frac{Δ t}{2})

and we compute

d_{t} E (x, t + \frac{Δ t}{2})

, rather than

d_{t} E (x, t)

in order to use the updated value of H. Expanding

d_{t} E (x, t + \frac{Δ t}{2})

, and solving for

E (x, t + Δ t)

gives

E (x, t + Δ t) = E (x, t) - \frac{1}{ε (x)} \frac{Δ t}{h} d \times H (x, t + \frac{Δ t}{2})

Equation (19) and Eq. (21) are called the Yee algorithm [3]. The notation does not show it, but each component of H and E is evaluated at a different position. For example to compute

d_{t} H_{z} \sim d_{x} E_{y} - d_{y} E_{x}

, if

H_{z}

lies on

(x, y, z)

then

E_{y}

must be known at

(x \pm \frac{h}{2}, y, z)

, and

E_{x}

at

(x, y \pm \frac{h}{2}, z)

.

To develop the Nonstandard-FDTD (NS-FDTD) model or in short NSFD we first construct an alternative Standard-FDTD (S-FDTD) model or SFD of Maxwell’s equations. Defining $H^{'} = μ h H / Δ t$ , we obtain

d_{t} H^{'} (x, t) = - d_{1} \times E (x, t),

d_{t} E (x, t + \frac{Δ t}{2}) = \frac{1}{ε μ} \frac{Δ t^{2}}{h^{2}} d_{1} \times H^{'} (x, t + \frac{Δ t}{2})

We introduced in earlier publication [4,7], the vector difference operator

d_{0}

which is also highly isotropic in nature with the property that

d_{1} \cdot d_{0} = d_{0} \cdot d_{1} = d_{0}^{2} .

whereas,

d_{1} \cdot d_{1} = d_{0}^{2}

,

d_{0} \cdot d_{0} \neq d_{0}^{2}

. In two dimensions

d_{0} = (d_{x}^{(0)}, d_{y}^{(0)})

, where

d_{x}^{(0)} = d_{x} + (\frac{1 - γ}{4}) d_{x} d_{y}^{2}

d_{y}^{(0)} = d_{y} + (\frac{1 - γ}{4}) d_{y} d_{x}^{2}

and d_{0} = d_{1} (1 + \frac{1 - γ}{4} d_{x} d_{y})

Here, γ which is also a function of

k = k_{0} \sqrt{ε}

(assuming uniform μ), given by the expression

γ = \frac{2}{3} - \frac{1}{90} {(k h)}^{2} - \frac{1}{15120} {(k h)}^{4} (11 - 5 \sqrt{2}) - \dots

is chosen to minimize the error of the NSFD approximation [4].

Making the substitutions $d_{1} \to d_{0}$ and $Δ t^{2} / (h^{2} ε_{i} μ) \to u_{0} {(\sqrt{ε})}^{2}$ in Eq. (23) we obtain a new NSFD model of Maxwell’s equations,

d_{t} H^{'} (x, t) = - d_{1} \times E (x, t)

d_{t} E (x, t + \frac{Δ t}{2}) = \frac{u_{0}^{2}}{ε} d_{0} (\sqrt{ε}) \times H^{'} (x, t + \frac{Δ t}{2})

Solving Eq. (28) and Eq. (29) for

H^{'} (x, t + \frac{Δ t}{2})

and

E (x, t + Δ t)

, we obtain the NS-Yee algorithm

H^{'} (x, t + \frac{Δ t}{2}) = H^{'} (x, t - \frac{Δ t}{2}) - d_{1} \times E (x, t)

E (x, t + Δ t) = E (x, t) + \frac{u_{0}^{2}}{ε} d_{0} (\sqrt{ε}) \times H^{'} (x, t + \frac{Δ t}{2})

This NSFD algorithm is somewhat different from the one introduced in [4]. We find that it gives better results when μ is constant. To reverse the time direction we need to reverse the sequence of the above computations; doing so we get,

H^{'} (x, t - \frac{Δ t}{2}) = H^{'} (x, t + \frac{Δ t}{2}) + d_{1} \times E (x, t + Δ t)

and

E (x, t) = E (x, t + Δ t) - \frac{u_{0}^{2}}{ε} d_{0} (\sqrt{ε}) \times H^{'} (x, t - \frac{Δ t}{2})

Here, we observed that we basically get the same equations as we already derived for forward propagating wave. This depicts that even if we iterate the wave in time reverse direction the nature of the wave remain same. For this purpose we have applied the periodic boundary conditions instead of Mur boundary conditions because the later causes field explosion when we iterate in time reverse condition [5]. Let, the computational domain be

x = 0, h, 2 h, \dots N_{x} h

. In periodic boundary we define,

ψ (- h, t) = ψ (N_{x} h, t)

where, waves that exit the computational domain at

x = 0

re-enter it from

x = N_{x} h

.

Algorithm Stability

When any algorithm is iterated is can be become unstable. Let A be an algorithm to solve for an unknown vector ψ. Let $φ_{0}$ be the initial solution vector. Then $φ_{1} = A φ_{0}$ , $φ_{2} = A φ_{1}$ , ⋯, and thus $φ_{n} = A^{n} φ_{0}$ . If $\underset{n \to \infty}{\lim | φ_{n} |} \to \infty$ when $| ψ |$ is finite, the algorithm is said to be unstable.

It can be shown that the D-dimensional FDTD algorithm is stable if $\frac{v Δ t}{h} \leq \frac{\sqrt{D}}{D}$ . Taking finite-difference algorithm of the form:

ψ (t + Δ t) = p ψ (t - Δ t) + 2 q ψ (t)

Defining

t = n Δ t

and

ψ^{n} = ψ (n Δ t)

we get,

ψ^{n + 1} = p ψ^{n - 1} + 2 q ψ^{n}

Now we postulate a solution:Inserting

ψ^{n} = λ^{n}

and dividing by

λ^{n - 1}

we get,

λ^{2} - 2 q λ - p = 0

Solving for λ we get,

λ_{\pm} = q \pm \sqrt{q^{2} + p}

Thus the most general solution is:

ψ^{n} = c_{+} λ_{+}^{n} + c_{-} λ_{-}^{n}

The condition for the stability will be

| ψ^{n} |

being finite as

n \to \infty

is

| λ_{\pm} | \leq 1

.

We want finite oscillatory solutions of the difference Eq. (36), for the oscillation the solutions must have an imaginary part. If $p < 0 \Rightarrow p = - | p |$ and $p + q^{2} < 0$ and $λ_{\pm} = q \pm \sqrt{q^{2} - | p |} = q \pm i \sqrt{| p | - q^{2}} \Rightarrow {| λ_{\pm} |}^{2} = q^{2} - | p | + q^{2} = | p |$ . Thus $| p | \leq 1$ and $q^{2} < | p |$ are the conditions for stability for an oscillating solution.

Applying these results to the wave equation given below:

(\partial_{t}^{2} - v^{2} \nabla^{2}) ψ (x, t) = 0,

Generalized Finite Difference model is:

(d_{t}^{2} - u_{0}^{2} d^{2}) ψ (x, t) = 0

and the FDTD algorithm is

ψ (x, t + Δ t) = - ψ (x, t - Δ t) + (2 + u_{0}^{2} d^{2}) ψ (x, t)

While in S-FDTD algorithm,

d^{2} = d_{1}^{2}

and

u_{0}^{2} = \frac{v^{2} Δ t^{2}}{h^{2}}

, whereas in NS-FDTD

d^{2} = d_{0}^{2}

and

u_{0}^{2} = \frac{\sin^{2} (ω Δ t / 2)}{\sin^{2} (k h / 2)}

. To direct the FDTD algorithm in the form of (34) we need to compute

d^{2} ψ

. Substituting

ψ = φ_{k} = e^{i (k . x - ω t)}

we can write,

d^{2} φ_{k} = - 2 d^{2} φ_{k} (x, t)

where the values of

d^{2}

vary case by case [4]:

For SFD:

d^{2} = d_{1}^{2} = 2 \sin^{2} (k h / 2) : 1Dimension

d^{2} = d_{1}^{2} = 2 [\sin^{2} (k_{x} h / 2) + \sin^{2} (k_{y} h / 2)] : 2Dimensions

Max

(d_{1}^{2}) = 2 D

where, D = Dimension

For NSFD:

d_{0}^{2} = d_{1}^{2} = d_{x}^{2} and Max (d_{0}^{2}) = 2 : 1Dimension

also,

d_{2}^{2} = 1 - \cos (k_{x} h / 2) \cos (k_{x} h / 2) : 2Dimensions

Till now we have not specified the value ofγ. For the stability we optimize the value of γ by

γ_{0}

(also termed as weight coefficient for

d_{1}^{2}

and

d_{2}^{2}

)Where

γ_{0} = 1 - \frac{\sin^{2} ({k^{'}}_{x} h / 2) + \sin^{2} ({k^{'}}_{y} h / 2) - \sin^{2} (k h / 2)}{2 \sin^{2} ({k^{'}}_{x} h / 2) \sin^{2} ({k^{'}}_{y} h / 2)}

and

({k^{'}}_{x}, {k^{'}}_{y}) = k (2^{- 1 / 4}, \sqrt{1 - 2^{- 1 / 2}})

.Thus for the NSFD algorithm we have,

d_{0}^{2} = γ_{0} d_{1}^{2} + (1 - γ_{0}) d_{2}^{2}

From the expanded value of γ given above, we get Max (

γ_{0}

) = 2/3. The above equation can be rewritten as:

d_{0}^{2} (k) = γ_{0} d_{1}^{2} (k) + (1 - γ_{0}) d_{2}^{2} (k)

where,

k = k (\cos θ, \sin θ)

and Max

(d_{0}^{2}) = 8 / 3

in 2 Dimensions.If we write,

ψ (x, t) = ψ_{x}^{t}

then the FDTD algorithm (41) becomes:

φ_{x}^{t + 1} = - φ_{x}^{t - 1} + 2 (1 - d^{2} u_{0}^{2}) φ_{x}^{t}

This equation is in the form of (36) and the solutions are:

λ_{\pm} = 1 - u_{0}^{2} d^{2} \pm \sqrt{{(1 - u_{0}^{2} d^{2})}^{2} - 1}

To accomplish the stability condition

| λ_{\pm} | \leq 1

, we have to have

{(1 - u_{0}^{2} d^{2})}^{2} \leq 1

.

If $d^{2} > 0$ and $u_{0}^{2} > 0$ , for non-vanishing monochromatic plane waves, the stability condition turns out to be [4]:

u_{0}^{2} \leq \frac{2}{d^{2}}

While

d^{2}

is a function of the magnitude and direction of k, stability is guaranteed by inserting the maximum value of

d^{2}

, Max

(d^{2})

instead of

d^{2}

.

u_{0}^{2} \leq \frac{2}{Max (d^{2})}

Considering NSFD algorithm, we have

d_{0}^{2} = d_{1}^{2} = d^{2}

and Max

(d_{0}^{2}) = 2, 8 / 3

for dimensions, D = 1, 2 respectively. Defining earlier,

u_{0}^{2} = \frac{\sin^{2} (ω Δ t / 2)}{\sin^{2} (k h / 2)}

the stability condition (48) becomes,

\frac{\sin (ω Δ t / 2)}{\sin (k h / 2)} \leq \sqrt{\frac{2}{Max(d_{0}^{2})}}

Replacing,

ϖ = ω Δ t, \bar{k} = k h

and

ϖ = c \bar{k}

, where the value of c is to be found out so that (49) does not violate. In D dimensions expressing

Max(d_{0}^{2}) = m_{d}^{n s f d}

the Eq. (49) turns out to be:

\sin (c \bar{k} / 2) \leq \sqrt{\frac{2}{m_{d}^{n s f d}}} \sin (\bar{k} / 2)

In the range

0 \leq θ \leq π / 2

,

\sin (θ)

rises monotonically with increasing θ, so it is obvious that

\sin θ_{1} \leq \sin θ_{2} \Rightarrow θ_{1} \leq θ_{2}

. The Nyquist limit can be articulated in the form

k < π / \sqrt{D}

. Hence we shall never figure out a FDTD algorithm which disobeys

\frac{λ}{h} > 2 \sqrt{D}

condition, and presuming that,

c \leq 1

we conclude from (50) that

c < \frac{2 \sqrt{d}}{π} \arcsin [\sqrt{\frac{2}{m_{d}^{n s f d}}} \sin (\frac{π}{2 \sqrt{D}})]

Putting the numerical values of

m_{d}^{n s f d}

from above, we have

c_{1}^{nsfd} = 1 : In One Dimension

c_{1}^{nsfd} ≃ \frac{2 \sqrt{2}}{π} \arcsin [\sqrt{\frac{3}{2}} \sin (\frac{π}{2 \sqrt{2}})] ≅ 0.84 : In Two Dimensions

3. Verifications and Practical Tests

In Fig. 1 we compare FDTD calculations of Mie scattering [6] for an infinite dielectric cylinder of radius $0.65 λ_{0}$ and refractive index $n = 1.8$ . The cylinder is parallel to the z-axis. An infinite plane wave propagating in the $+ x$ direction impinges upon the cylinder. The electric field is polarized parallel to the cylinder axis, and the scattered intensity is visualized in shades of red for the analytic solution (left) the NS-FDTD algorithm (middle) and the S-FDTD algorithm (right) with a discretization of $λ_{0} / h = 8$ outside the cylinder, and $λ_{0} / n h ≅ 4.4$ inside. Similar result was published by one of the author in his previous publication [7].

Fig. 1 Mie scattering from an infinite dielectric cylinder. Scattered intensity is visualized in shades of red. E is polarized parallel cylinder axis. Infinite plane wave impinges from the left. Cylinder radius = $0.65 λ_{0}$ , $λ_{0}$ = vacuum wavelength, refractive index $n = 1.8$ .

Download Full Size | PDF

In Fig. 2(a) the scattered fields shown in Fig. 1 are plotted as a function of angle the (θ) from the $+ x$ -axis on a circular contour of radius $λ_{0}$ centered on the axis of the cylinder. At discretization $λ_{0} / h = 8$ the root mean square error for the NS-FDTD and S-FDTD algorithms is $ε_{NS-8} = 0.04$ , and $ε_{S-8} = 0.20$ , respectively. In Fig. 2(b) at discretization $λ_{0} / h = 24$ we compare the S-FDTD calculation with the analytic solution. The root mean square error is $ε_{S-24} = 0.04$ the same as the NS-FDTD algorithm at $λ_{0} / h = 8$ . At $λ_{0} / h = 8$ the NS-FDTD algorithm delivers the same accuracy as the S-FDTD one does at $λ_{0} / h = 24$ . Because of the lower coarser discretization, the computational cost of NSFDTD algorithm is only $1 / 27$ th that of the Standard FDTD one.

Fig. 2 Angular distribution of scattered intensity about a circular contour of radius $λ_{0}$ centered on the axis of the cylinder. (a) Analytic solution (black) compared with the NSFD calculation (NSFD-8, red) and the SFD one (SFD-8, blue) at $λ_{0} / h = 8$ . (b) Analytic solution (black) compared with the SFD (SFD-24, blue) one at $λ_{0} / h = 24$ .

Download Full Size | PDF

Understanding the fact that NSFD algorithm is better than S-FDTD one, we extended our study for the Time Inverse algorithm on the detection of size, shape of different scatterers. for this purpose we have generated a weighted source array continuous wave Gaussian beam that turns on at the time t = 1. The beam width of the Gaussian pulse taken for the simulation purpose is 1.25 $λ_{0}$ , where $λ_{0}$ is the wavelength of the central frequency of the pulse. For the simulation we have chosen three types of the scatterer; Circular, Rectangular and triangular. Now we have allowed the light to be passed through the scatterer and it further moves to a distance which is very close to the periodic boundary of the computational space ( $N_{x} =$ 140h), then we remove the scatterer and iterate the forward moving wave in the time reverse in order to trace the location and shape of the scatterer and examine how well it will match with the actual location and size. For this experiment, we have confined our observation on TE mode only, i.e. we will compute x-intensity and y-intensity as depicted in Fig. 3 .

Fig. 3 Time Inverse Computed Intensities using NS-FDTD for: (a)Circular scatterer (diameter = 2 $λ_{0}$ )(b) Rectangular scatterer (side = 2 $λ_{0}$ ) (c)Triangular scatterer(base = 2 $λ_{0}$ ,height = 1 $λ_{0}$ ),where, $λ_{0} / h = 10$

Download Full Size | PDF

In the next experiment, our objective is to investigate the effect of refractive index on the size of the scatterer. The minimum detectable size in case is called distinguishable size. The refractive index of the scatterer is fixed and in all the cases is taken to 1.8. In this experiment, we have considered circular scatterer because it is homogeneous in shape and no steep edges and corners. Now, gradually we have changed the refractive index of the circular scatterer from 1.1 to 1.8 and compute the distinguishable size of the circular scatterer for a particular refractive index and plot the graph shown in Fig. 4 for TM Mode of polarization.

Fig. 4 Distinguishable Size vs. Refractive Index

Download Full Size | PDF

4. Results and discussions

From the Fig. 3 we noticed that the Time Inverse algorithm (TE mode) has two intensities, x-intensity and y intensity in which left one is the computed x-intensity and right one is the computed y-intensity. It is clearly seen from the figure that x-intensities gives the clear indication of the shape, which means the first one is circular, the middle one is the rectangular and the last one is the triangular scatterer. The y- intensity is perpendicularly oriented to the x-intensity and it does not differentiate the edges well as it does for the x-intensity. But considering both the intensities we can distinguish the exact shape, size and location of the scatterer by Non standard Time Domain Time Inverse algorithm. It is very helpful in medical imaging because a small cancerous tumor of arbitrary shape is perfectly detected and chance of misdetection or false detection is eliminated.

In Fig. 4, shown below applying the TM Mode of polarization of incoming light we have noticed that as the refractive index of the circular scatterer gradually increases the distinguishable size of the scatterer gradually decreases. One of the probable reasons is that, the scatterer having the refractive index close to vacuum, it size must be large enough in order to get some considerable contrast. We have also noticed that our Non-Standard Time Domain Inverse algorithm is good enough to detect the shape, size and location of the object in comparison to Standard Time Domain Inverse algorithm even if the refractive index is very close to unity. In the later case more time steps are required and the picture contrast is not good as that of NSFD algorithm.

Acknowledgements

We would like to acknowledge the Ministry of Education, Government of Japan for granting one of the authors (KC) with a scholarship during the period of this work.

References and links

1. L. W. Schmerr, Jr., in Fundamentals of Ultrasonic Nondestructive Evaluation: A Modeling Approach, (Plenum Press, New York, 1998).

2. R. Mickens, E, in Nonstandard finite difference models of differential equations, (World Scientific Publishing Co., Inc., River Edge, NJ, 1994).

3. K. S. Kunz, and R. J. Luebbers, in The Finite Difference Time Domain Method for Electromagnetics, (CRC Press, New York, 1993).

4. J. B. Cole, S. Banerjee, and M. Haftel, in Advances in the Applications of Nonstandard Finite Difference Schemes, ed. R. E. Mickens (World Scientific Singapore, 2007), Chap.4.

5. R. Sorrentino, L. Roselli, and P. Mezzanotte, “Time reversal in finite difference time domain method,” IEEE Microw. Guid. Wave Lett. 3(11), 402–404 (1993). [CrossRef]

6. P. W. Barber, and S. C. Hill, in Light Scattering by Particles: Computational Methods, (World Scientific, Singapore, 1990).

7. J. B. Cole, “High-Accuracy Yee algorithm Based on Nonstandard Finite Difference: New Developments and Verifications,” IEEE Trans. Antenn. Propag. 50(9), 1185–1191 (2002). [CrossRef]

8. H. Kudo, et al., “Numerical Dispersion and Stability Condition of the Nonstandard FDTD Method” Electronics and Communications in Japan, 85, 22–30(2002), http://www3.interscience.wiley.com/cgi-bin/fulltext/93514073/PDFSTART.

Simulation study on the detection of size, shape and position of three different scatterers using Non-standard time domain time inverse Maxwell’s algorithm

Abstract

1. Introduction

2. Model and Calculation Methods

Algorithm Stability

3. Verifications and Practical Tests

4. Results and discussions

Acknowledgements

References and links

Cited By

Figures (4)

Equations (57)

Optics Express