**3. Solving the radiative transfer equation: Monte Carlo simulation**

The Monte Carlo name was coined by Nicholas Metropolis [7]. The modern version of the method was invented in the late 1940s and found use early use in the studies of neutron transport for the design of nuclear weapons [8, 9]. More recently, Monte Carlo techniques evolved and are now used to solve problems in physical sciences, economics, engineering, and pure mathematics.

Monte Carlo methods may be used to solve any problem that have a probabilistic interpretation, meaning that if the probability of occurrence of each separate event in a sequence of events is known, one can determine the probability that the entire sequence of events will occur. By tracing the fate of millions of photons according to statistical probabilities, a solution for the RTE is built. Each simulated photon path is randomly distinct from the others, as determined by the probabilities of absorption and scattering in the underwater channel.

The Monte Carlo simulations of photon propagation offer a flexible yet rigorous approach toward solving numerically the RTE. In the last decades, the exponential improvement in computer speeds attenuated the problem of the high computational cost of Monte Carlo simulations.

In the following subsections, a detailed description for building a solution for the RTE, step by step, will be shown.

As seen in [10], there are four main Monte Carlo methods for photon migration in turbid media, that is, four different ways to build photon's trajectories. The four methods are referred as the albedo-weight method (AW), the albedo-rejection method (AR), the absorption-scattering path length rejection method (ASPR), and the microscopic Beer-Lambert law method (mBLL).

The AW method considers a source emitting packets of many photons and after each interaction a fraction of the photons in the packet are absorbed, i.e., lost. In this type of simulation, the packet is emitted with an initial weight of 1 and at each interaction the current weight is multiplied by the albedo to account for the loss of photons by absorption and the packet continues propagating in the scattered direction with a reduced weight. This method has in its favor that there are no wasted computations because all photon packets eventually reach the detector. This approach was favored by several research papers [11–16].

In this chapter, we present a simulation using the AR method; the main advantage of this method is that it mimics what happens in nature by tracking individual photons emitted by some source. Naturally for this case, the computational time of the absorbed photons are wasted.

It is seen in the investigation conducted by Sassaroli that the total probability of detecting a photon is equivalent in the four methods and they have statistical equivalence, with some differences regarding the convergence time. When comparing the AW to the AR method, as one may expect, the AR will require, in most scenarios, more launched photons for convergence. However, it is also noted that in the AR method, the photon is detected or lost similarly to what happens in a timeresolved experiment in which photons are detected using the time-correlated single photon counting technique and the noise on the simulated response reproduces the Poisson statistics that typically characterize the noise in experiments. Therefore, as we already stated, the AR method reproduces a more physical simulation of

*dL r*ð Þ *; θ; ϕ*

*Monte Carlo Radiative Transfer Modeling of Underwater Channel*

putting it in terms of the photon path length <sup>τ</sup> <sup>¼</sup> <sup>Ð</sup> *<sup>r</sup>*

of the beam between photon path lengths τ and τ þ *d*τ is:

*L r*ð Þ¼ *; θ; ϕ L*ð Þ 0*; θ; ϕ e*

*L*ð Þ¼ *τ; θ; ϕ L*ð Þ 0*; θ; ϕ e*

The exponential decay of radiance can be explained in terms of the fate of individual photons if the probability of any photon being absorbed or scattered out

*pT*ð Þτ *d*τ ¼ *e*

ð*<sup>x</sup>*<sup>2</sup> *x*1

�τ

Mobley [3] notes that *pT*ð Þτ satisfies the normalization condition for a probabil-

with *x*<sup>1</sup> ¼ 0 and *x*<sup>2</sup> ¼ ∞. The corresponding cumulative distribution function

selecting a random number *q* from a uniform distribution, 0≤*q*≤ 1, and solving

�τ

�τ

τ ¼ �ln 1ð Þ¼� � *q* ln ð Þ*q :* (38)

*PT*ð Þ¼ τ 1 � *e*

*q* ¼ *PT*ð Þ¼ τ 1 � *e*

For homogenous water, where *c*(*r*) does not depend on *r*, where τ ¼ *cr*, the

Recalling that the photon path length, or step size, is the geometric distance a photon will travel between a scattering event and an absorption event, after determining the distance the photon will travel from origin, one must decide if the next interaction will be an absorption or a scattering event. For this, we recall the definition of the single scattering albedo, *ωo*, that is the ratio between the scattering and the absorption coefficients. The single scattering albedo is referred sometimes as the probability of photon survival, because this ratio is the probability that the photon will undergo scattering instead of absorption. As we have mentioned before, in this Monte Carlo simulation, we used the AR method and to numerically

*<sup>r</sup>* ¼ � <sup>1</sup> *c*

implement this, we must select again a random number from a uniform

distribution*,* 0≤ *q*≤1, to determine which event will take place:

starting point. Integrating Eq. (31), we have:

*DOI: http://dx.doi.org/10.5772/intechopen.85961*

ity density function (PDF):

photon will travel a geometric distance, *r*:

**3.2 Absorption events and scattering angles**

(CDF) is:

for τ, gives:

**107**

where *c r*ð Þ is the beam attenuation coefficient and *r* is the distance from a given

� Ð *r*

*dr* ¼ �cð Þ*<sup>r</sup>* <sup>L</sup>ð Þ *<sup>r</sup>; θ ϕ ,* (31)

<sup>0</sup> *c r*<sup>0</sup> ð Þ*dr*<sup>0</sup> :

�*τ*

<sup>0</sup> *c r*<sup>0</sup> ð Þ*d r*<sup>0</sup> ð Þ*,* (32)

*:* (33)

*d*τ*:* (34)

*,* (36)

*,* (37)

ln ð Þ*q :* (39)

*px*ð Þ *x dx* ¼ 1*,* (35)

**Figure 2.** *Flow chart for the AR Monte Carlo simulation.*

propagation than the AW method, where the simulated photons are energy packets rather than physical photons. In **Figure 2**, we show a flow chart of the Monte Carlo algorithm used in this simulation, each step in this chart will be explained with further detail in the following sections.

### **3.1 Photon path length**

The photon path length, or step size, determines how far the photon will travel before encountering a water molecule or particle, meaning the distance a photon travels between an absorption event and a scattering event. The path length is directly related to the water type and, consequently, to the beam attenuation coefficient introduced in Section 2.1, and shown for several water types in **Table 1**.

As stated in [4], from the derivation of the RTE, the radiance in a particular direction ð Þ *θ; ϕ* decays due to absorption and scattering according to:

*Monte Carlo Radiative Transfer Modeling of Underwater Channel DOI: http://dx.doi.org/10.5772/intechopen.85961*

$$\frac{d\mathcal{L}(r,\theta,\phi)}{dr} = -\mathbf{c}(r)\mathcal{L}(r,\theta\,\phi),\tag{31}$$

where *c r*ð Þ is the beam attenuation coefficient and *r* is the distance from a given starting point. Integrating Eq. (31), we have:

$$L(r, \theta, \phi) = L(\mathbf{0}, \theta, \phi) e^{-\int\_0^r c(r') d(r')},\tag{32}$$

putting it in terms of the photon path length <sup>τ</sup> <sup>¼</sup> <sup>Ð</sup> *<sup>r</sup>* <sup>0</sup> *c r*<sup>0</sup> ð Þ*dr*<sup>0</sup> :

$$L(\mathfrak{r}, \theta, \phi) = L(\mathbf{0}, \theta, \phi)e^{-\mathfrak{r}}.\tag{33}$$

The exponential decay of radiance can be explained in terms of the fate of individual photons if the probability of any photon being absorbed or scattered out of the beam between photon path lengths τ and τ þ *d*τ is:

$$p\_T(\mathfrak{r})d\mathfrak{r} = e^{-\mathfrak{r}}d\mathfrak{r}.\tag{34}$$

Mobley [3] notes that *pT*ð Þτ satisfies the normalization condition for a probability density function (PDF):

$$\int\_{\varkappa\_1}^{\varkappa\_2} p\_x(\varkappa) d\varkappa = 1,\tag{35}$$

with *x*<sup>1</sup> ¼ 0 and *x*<sup>2</sup> ¼ ∞. The corresponding cumulative distribution function (CDF) is:

$$P\_T(\mathbf{r}) = \mathbf{1} - e^{-\mathbf{r}},\tag{36}$$

selecting a random number *q* from a uniform distribution, 0≤*q*≤ 1, and solving for τ, gives:

$$q = P\_T(\mathfrak{r}) = \mathbf{1} - \mathfrak{e}^{-\mathfrak{r}},\tag{37}$$

$$
\pi = -\ln\left(1 - q\right) = -\ln\left(q\right). \tag{38}
$$

For homogenous water, where *c*(*r*) does not depend on *r*, where τ ¼ *cr*, the photon will travel a geometric distance, *r*:

$$r = -\frac{1}{c} \ln \left( q \right). \tag{39}$$

### **3.2 Absorption events and scattering angles**

Recalling that the photon path length, or step size, is the geometric distance a photon will travel between a scattering event and an absorption event, after determining the distance the photon will travel from origin, one must decide if the next interaction will be an absorption or a scattering event. For this, we recall the definition of the single scattering albedo, *ωo*, that is the ratio between the scattering and the absorption coefficients. The single scattering albedo is referred sometimes as the probability of photon survival, because this ratio is the probability that the photon will undergo scattering instead of absorption. As we have mentioned before, in this Monte Carlo simulation, we used the AR method and to numerically implement this, we must select again a random number from a uniform distribution*,* 0≤ *q*≤1, to determine which event will take place:

propagation than the AW method, where the simulated photons are energy packets rather than physical photons. In **Figure 2**, we show a flow chart of the Monte Carlo algorithm used in this simulation, each step in this chart will be explained with

The photon path length, or step size, determines how far the photon will travel before encountering a water molecule or particle, meaning the distance a photon travels between an absorption event and a scattering event. The path length is directly related to the water type and, consequently, to the beam attenuation coefficient introduced in Section 2.1, and shown for several water types in **Table 1**. As stated in [4], from the derivation of the RTE, the radiance in a particular

direction ð Þ *θ; ϕ* decays due to absorption and scattering according to:

further detail in the following sections.

*Flow chart for the AR Monte Carlo simulation.*

*Wireless Mesh Networks - Security, Architectures and Protocols*

**3.1 Photon path length**

**Figure 2.**

**106**

$$\text{event} = \begin{cases} \text{absorption, for } q \le \alpha\_o \\ \text{scattering, for } q > \alpha\_o \end{cases} \tag{40}$$

From this description, one may clearly see where the name "Albedo-Rejection" comes from: if the photon is absorbed, it must be terminated, and its position is recorded for further analysis; if the photon was scattered, the scattering angles must be calculated, and its direction of propagation is updated.

### *3.2.1 The Henyen-Greenstein phase function*

Unlike the FSO channel, in an underwater environment, the photons will encounter a much larger number of particles, and thus multiple scattering events will play a more important role in the received optical power, when compared to other optical wireless communication channels.

One of the models proposed to describe the SPF is the Henyen-Greenstein phase function (HG) [12, 17, 18]. The HG was originally proposed by Henyen and Greenstein to describe the scattering of light by interstellar dust [19]:

$$\widetilde{\mathfrak{F}}(\mathbf{g},\theta) = p\_{HG} = \frac{1}{4\pi} \frac{1 - \mathbf{g}^2}{(1 + \mathbf{g}^2 - 2\mathbf{g}\cos\theta)^\frac{3}{2}},\tag{41}$$

where *g* is referred to as the anisotropy factor and is also the mean cosine of the scattering angle cos *θ*, which takes values between �1 and 1. When *g* ¼ �1, there is complete backscattering; for *g* ¼ 0, the scattering is isotropic; and *g* ¼ 1 gives complete forward scattering. Using *g* ¼ 0*:*924, the HG is a good approximation for photon propagation in water and a good fit for the Petzold's measurements.

To find the scattering angle, Eq. (41) must be solved for *θ*:

$$\cos \theta = \frac{1}{2g} \left\{ 1 + \text{g}^2 - \left( \frac{1 - \text{g}^2}{1 - \text{g} + 2 \text{g} \, q} \right)^2 \right\},\tag{42}$$

where *q* is a uniformly distributed random number between 0 and 1. The complete method for obtaining Eq. (42) may be found elsewhere [4].

Because scattering is a three-dimensional process, besides the polar angle (*θ*Þ defined above, the azimuthal ð Þ *ϕ* angle must also be defined. For an unpolarized beam, the azimuthal angle has an equal probability to take any value between 0 and 2π [3]. Selecting again a random number, the azimuthal angle *ϕ* can be found by:

$$
\phi = 2\pi q.\tag{43}
$$

To overcome this limitation, Haltrin [21] proposed a two-term Heyen-Greenstein (TTHG) phase function that matched Petzold's results better at small scattering angles and provided the elongation seen at larger angles. The TTHG is

*Comparison between the HG and Petzold's average phase function (left) and large (right) scattering angles.*

<sup>þ</sup> ð Þ <sup>1</sup> � <sup>α</sup> *pHG <sup>θ</sup>; <sup>g</sup>*<sup>2</sup>

where *pHG* is the HG given by Eq. (41), α is a weighting factor, and *g*<sup>1</sup> is given a value near 1, which will make the TTHG increase more strongly at small angles. The parameter *g*<sup>2</sup> takes a negative value that will make the TTHG increase, as the angle *θ*

. The values of α and *g*<sup>2</sup> may be calculated by:

*,* (44)

*pTTHG*ð Þ¼ *θ* α*pHG θ; g*<sup>1</sup>

*Comparison of the HG (blue) and the Petzold's average phase function (red).*

*Monte Carlo Radiative Transfer Modeling of Underwater Channel*

*DOI: http://dx.doi.org/10.5772/intechopen.85961*

given by:

**Figure 4.**

**Figure 3.**

approaches 180°

**109**

Comparing the HG phase function using Eq. (41) with the Petzold's avg. Phase Function extracted from [3], we can see in **Figure 3** that the HG is indeed an adequate approximation for the SPF.

The HG phase function found great use in ocean optics simulations and is heavily used in Monte Carlo simulations and other relevant channel modeling works by several researchers [20].

### *3.2.2 The two-term Henyey-Greenstein phase function*

Despite being an approximate fit, it can be seen in **Figure 4** that the HG phase function shows noticeable differences from the Petzold's measurements for small angles, *θ* < 20° , and large angles *θ* > 150° .

*Monte Carlo Radiative Transfer Modeling of Underwater Channel DOI: http://dx.doi.org/10.5772/intechopen.85961*

event <sup>¼</sup> absorption*,*for *<sup>q</sup>* <sup>≤</sup>*ω<sup>o</sup>*

Unlike the FSO channel, in an underwater environment, the photons will encounter a much larger number of particles, and thus multiple scattering events will play a more important role in the received optical power, when compared to

function (HG) [12, 17, 18]. The HG was originally proposed by Henyen and Greenstein to describe the scattering of light by interstellar dust [19]:

4π

where *g* is referred to as the anisotropy factor and is also the mean cosine of the scattering angle cos *θ*, which takes values between �1 and 1. When *g* ¼ �1, there is complete backscattering; for *g* ¼ 0, the scattering is isotropic; and *g* ¼ 1 gives complete forward scattering. Using *g* ¼ 0*:*924, the HG is a good approximation for photon propagation in water and a good fit for the Petzold's measurements.

<sup>1</sup> <sup>þ</sup> *<sup>g</sup>*<sup>2</sup> � <sup>1</sup> � *<sup>g</sup>*<sup>2</sup>

where *q* is a uniformly distributed random number between 0 and 1. The com-

Because scattering is a three-dimensional process, besides the polar angle (*θ*Þ defined above, the azimuthal ð Þ *ϕ* angle must also be defined. For an unpolarized beam, the azimuthal angle has an equal probability to take any value between 0 and 2π [3]. Selecting again a random number, the azimuthal angle *ϕ* can be found by:

Comparing the HG phase function using Eq. (41) with the Petzold's avg. Phase

Despite being an approximate fit, it can be seen in **Figure 4** that the HG phase function shows noticeable differences from the Petzold's measurements for small

.

Function extracted from [3], we can see in **Figure 3** that the HG is indeed an

The HG phase function found great use in ocean optics simulations and is heavily used in Monte Carlo simulations and other relevant channel modeling works

� �<sup>2</sup> ( )

<sup>e</sup>βð Þ¼ *<sup>g</sup>; <sup>θ</sup> pHG* <sup>¼</sup> <sup>1</sup>

To find the scattering angle, Eq. (41) must be solved for *θ*:

plete method for obtaining Eq. (42) may be found elsewhere [4].

cos *<sup>θ</sup>* <sup>¼</sup> <sup>1</sup> 2*g*

adequate approximation for the SPF.

*3.2.2 The two-term Henyey-Greenstein phase function*

, and large angles *θ* > 150°

by several researchers [20].

angles, *θ* < 20°

**108**

One of the models proposed to describe the SPF is the Henyen-Greenstein phase

<sup>1</sup> � *<sup>g</sup>*<sup>2</sup> <sup>1</sup> <sup>þ</sup> *<sup>g</sup>* ð Þ <sup>2</sup> � <sup>2</sup>*gcos<sup>θ</sup>* <sup>3</sup>

1 � *g* þ 2*gq*

2

*ϕ* ¼ 2π*q:* (43)

*,* (41)

*,* (42)

From this description, one may clearly see where the name "Albedo-Rejection" comes from: if the photon is absorbed, it must be terminated, and its position is recorded for further analysis; if the photon was scattered, the scattering angles must

�

be calculated, and its direction of propagation is updated.

*Wireless Mesh Networks - Security, Architectures and Protocols*

*3.2.1 The Henyen-Greenstein phase function*

other optical wireless communication channels.

scattering*,*for *q*>*ω<sup>o</sup>*

*:*

(40)

**Figure 3.** *Comparison of the HG (blue) and the Petzold's average phase function (red).*

**Figure 4.** *Comparison between the HG and Petzold's average phase function (left) and large (right) scattering angles.*

To overcome this limitation, Haltrin [21] proposed a two-term Heyen-Greenstein (TTHG) phase function that matched Petzold's results better at small scattering angles and provided the elongation seen at larger angles. The TTHG is given by:

$$p\_{TTHG}(\theta) = \mathfrak{a}p\_{HG}(\theta, \mathfrak{g}\_1) + (\mathfrak{1} - \mathfrak{a})p\_{HG}(\theta, \mathfrak{g}\_2),\tag{44}$$

where *pHG* is the HG given by Eq. (41), α is a weighting factor, and *g*<sup>1</sup> is given a value near 1, which will make the TTHG increase more strongly at small angles. The parameter *g*<sup>2</sup> takes a negative value that will make the TTHG increase, as the angle *θ* approaches 180° . The values of α and *g*<sup>2</sup> may be calculated by:

*Wireless Mesh Networks - Security, Architectures and Protocols*

$$\mathbf{g}\_2 = -0.3061446 + 1.000568 \mathbf{g}\_1 - 0.0182633 \mathbf{g}\_1^2 + 0.03643748 \mathbf{g}\_1^3,\tag{45}$$

$$\alpha = \frac{\mathbf{g}\_2(\mathbf{1} + \mathbf{g}\_2)}{(\mathbf{g}\_1 + \mathbf{g}\_2)(\mathbf{1} + \mathbf{g}\_2 - \mathbf{g}\_1)}. \tag{46}$$

values of *g*<sup>1</sup> against the Petzold's average phase function, we can see that *g*<sup>1</sup> ¼ 0*:*9809, as proposed in [17], provides a good fit at small angles, where the phase function has its largest values, and maintains the elongation at larger angles

Despite the attempts from several researchers of developing analytical phase functions for natural waters, we can see that the ones described in the previous sections do not provide a good fit over the total range of *θ*. In [22], Kirk computed

seen before, a lot of scattering happens. In the results presented later in Section 4, we computed a CDF of Petzold's data over the whole range of scattering angles and used linear interpolation to compute *θ* for each random value selected. The accuracy obtained by this method comes with the drawback of increased computation time. In [23], Haltrin developed an empirical fit using strong regression for all Petzold's measurements. One of the strongest points of this model is that the Phase Function may be obtained only by using the scattering coefficient and the single

> 5 *n*¼1

> > 5 *n*¼1

ð Þ �<sup>1</sup> *<sup>n</sup>*

ð Þ �<sup>1</sup> *<sup>n</sup>*

<sup>p</sup> � <sup>16</sup>*:*722b <sup>þ</sup> <sup>5</sup>*:*932b ffiffiffi

*knθ <sup>n</sup>* ð Þ<sup>2</sup> � � � � *,* (48)

*knθ <sup>n</sup>* ð Þ<sup>2</sup> � � � � *,* (49)

*p*ð Þ*θ* sin ð Þ*θ dθ* ¼ 1*,* (50)

*k*<sup>1</sup> ¼ 1*:*188 � 0*:*688ω0*,* (52) *k*<sup>2</sup> ¼ 0*:*1 3ð Þ *:*07 � 1*:*90ω<sup>0</sup> *,* (53) *k*<sup>3</sup> ¼ 0*:*01 4ð Þ *:*58 � 3*:*02ω<sup>0</sup> *,* (54) *k*<sup>4</sup> ¼ 0*:*001 3ð Þ *:*24 � 2*:*25ω<sup>0</sup> *,* (55) *k*<sup>5</sup> ¼ 0*:*0001 0ð Þ *:*84 � 0*:*61ω<sup>0</sup> *:* (56)

*b*

<sup>p</sup> *,* (51)

. Unfor-

, angles where, as

a CDF from Petzold's data for *θ* over the range of 2.5–180° in intervals of 5°

tunately, his CDF did not include data for angles smaller than 2*:*5°

βð Þ¼ *θ* exp ψ 1 þ ∑

*<sup>b</sup>* exp <sup>ψ</sup> <sup>1</sup> <sup>þ</sup> <sup>∑</sup>

where *b* is the scattering coefficient and the remaining parameters are:

*b*

The strong regression given by these equations can be used for an empirical model for the phase functions. As seen from **Figure 7**, this model represents a better fit to the Petzold's measurements than the phase functions presented before,

In our simulations, we found that when comparing the CDF computed by directly interpolating the Petzold's measurements, the CDF computed from the empirical phase function showed some divergence. Our CDF showed that for

especially in the regions where the phase function is the strongest.

characteristic of the Petzold's phase function.

*DOI: http://dx.doi.org/10.5772/intechopen.85961*

*Monte Carlo Radiative Transfer Modeling of Underwater Channel*

*3.2.3 Haltrin empirical phase function*

scattering albedo as inputs:

with:

**111**

and the scattering phase function:

pð Þ¼ *θ*

4π

<sup>ψ</sup> <sup>¼</sup> <sup>2</sup>*:*<sup>598</sup> <sup>þ</sup> <sup>17</sup>*:*<sup>748</sup> ffiffiffi

Haltrin also derived a relationship that relates the scattering coefficient, *b*, with the anisotropy factor, *g*1:

$$\mathbf{g\_1} = \mathbf{1} - \frac{\mathbf{0}.001247}{b}.\tag{47}$$

After calculating the parameters using Eqs. (44)–(47) and computing the phase function for clear ocean waters with *<sup>b</sup>* <sup>¼</sup> <sup>0</sup>*:*<sup>037</sup> *<sup>m</sup>*�<sup>1</sup> using Eq. (44), we now compare the resulting phase function with the Petzold's phase function for clear ocean water. We can see from **Figure 5** that the TTHG represents a better approximation for small angles and presents the characteristic elongation at large angles. Despite that, in **Figure 6**, it can also be seen that the TTHG with *g*<sup>1</sup> ¼ 0*:*9943 calculated using Eq. (47), differs greatly from the Petzold's average phase function data at intermediate angles, 20–110° . Comparing the results for the TTHG for different

**Figure 6.**

*Comparison of scattering angles for two phase functions with Petzold's data for different values of the anisotropy factor g*1*.*

### *Monte Carlo Radiative Transfer Modeling of Underwater Channel DOI: http://dx.doi.org/10.5772/intechopen.85961*

values of *g*<sup>1</sup> against the Petzold's average phase function, we can see that *g*<sup>1</sup> ¼ 0*:*9809, as proposed in [17], provides a good fit at small angles, where the phase function has its largest values, and maintains the elongation at larger angles characteristic of the Petzold's phase function.

## *3.2.3 Haltrin empirical phase function*

*<sup>g</sup>*<sup>2</sup> ¼ �0*:*<sup>3061446</sup> <sup>þ</sup> <sup>1</sup>*:*000568*g*<sup>1</sup> � <sup>0</sup>*:*01826332*g*<sup>2</sup>

*Wireless Mesh Networks - Security, Architectures and Protocols*

the anisotropy factor, *g*1:

intermediate angles, 20–110°

**Figure 5.**

**Figure 6.**

*factor g*1*.*

**110**

<sup>α</sup> <sup>¼</sup> *<sup>g</sup>*<sup>2</sup> <sup>1</sup> <sup>þ</sup> *<sup>g</sup>*<sup>2</sup>

*g*<sup>1</sup> þ *g*<sup>2</sup>

*Comparison of small and large angles for two phase functions,TTHG and HG, with Petzold's data.*

*Comparison of scattering angles for two phase functions with Petzold's data for different values of the anisotropy*

<sup>1</sup> <sup>þ</sup> *<sup>g</sup>*<sup>2</sup> � *<sup>g</sup>*<sup>1</sup>

Haltrin also derived a relationship that relates the scattering coefficient, *b*, with

After calculating the parameters using Eqs. (44)–(47) and computing the phase function for clear ocean waters with *<sup>b</sup>* <sup>¼</sup> <sup>0</sup>*:*<sup>037</sup> *<sup>m</sup>*�<sup>1</sup> using Eq. (44), we now compare the resulting phase function with the Petzold's phase function for clear ocean water. We can see from **Figure 5** that the TTHG represents a better approximation for small angles and presents the characteristic elongation at large angles. Despite that, in **Figure 6**, it can also be seen that the TTHG with *g*<sup>1</sup> ¼ 0*:*9943 calculated using Eq. (47), differs greatly from the Petzold's average phase function data at

*<sup>g</sup>*<sup>1</sup> <sup>¼</sup> <sup>1</sup> � <sup>0</sup>*:*<sup>001247</sup>

<sup>1</sup> þ 0*:*03643748*g*<sup>1</sup>

*:* (46)

*<sup>b</sup> :* (47)

. Comparing the results for the TTHG for different

3

*,* (45)

Despite the attempts from several researchers of developing analytical phase functions for natural waters, we can see that the ones described in the previous sections do not provide a good fit over the total range of *θ*. In [22], Kirk computed a CDF from Petzold's data for *θ* over the range of 2.5–180° in intervals of 5° . Unfortunately, his CDF did not include data for angles smaller than 2*:*5° , angles where, as seen before, a lot of scattering happens. In the results presented later in Section 4, we computed a CDF of Petzold's data over the whole range of scattering angles and used linear interpolation to compute *θ* for each random value selected. The accuracy obtained by this method comes with the drawback of increased computation time.

In [23], Haltrin developed an empirical fit using strong regression for all Petzold's measurements. One of the strongest points of this model is that the Phase Function may be obtained only by using the scattering coefficient and the single scattering albedo as inputs:

$$\mathfrak{F}(\theta) = \exp\left[\mathfrak{w}\left(\mathbf{1} + \sum\_{n=1}^{5} (-\mathbf{1})^{n} k\_{n} \theta^{\left(\frac{\mathbf{1}}{2}\right)}\right)\right],\tag{48}$$

and the scattering phase function:

$$\mathbf{p}(\theta) = \frac{4\pi}{b} \exp\left[\psi \left(\mathbf{1} + \sum\_{n=1}^{5} (-\mathbf{1})^{n} k\_{n} \theta^{\left(\frac{n}{2}\right)}\right)\right],\tag{49}$$

with:

$$\frac{1}{2}\int\_{0}^{\pi}p(\theta)\sin\left(\theta\right)d\theta=\mathbf{1},\tag{50}$$

where *b* is the scattering coefficient and the remaining parameters are:

$$
\Psi = 2.598 + 17.748\sqrt{b} - 16.722\mathbf{b} + 5.932\mathbf{b}\sqrt{b},\tag{51}
$$

$$k\_1 = 1.188 - 0.688a\_0,\tag{52}$$

$$k\_2 = 0.1(3.07 - 1.90a\_0),\tag{53}$$

$$k\_3 = 0.01(4.58 - 3.02a\_0),\tag{54}$$

$$k\_4 = 0.001(3.24 - 2.25\alpha\_0),\tag{55}$$

$$k\_5 = 0.0001(0.84 - 0.61a\_0).\tag{56}$$

The strong regression given by these equations can be used for an empirical model for the phase functions. As seen from **Figure 7**, this model represents a better fit to the Petzold's measurements than the phase functions presented before, especially in the regions where the phase function is the strongest.

In our simulations, we found that when comparing the CDF computed by directly interpolating the Petzold's measurements, the CDF computed from the empirical phase function showed some divergence. Our CDF showed that for

It can be seen in [25] that with *n* ¼ 1*:*16 and *μ* ¼ 3*:*4586, the FF phase function gives a good fit to the average Petzold phase function over all scattering angles. This

When comparing the CDF computed from the FF phase function and the CDF

In [25], a semi-analytical model for a phase function was proposed by Sahu and Shanmugam (SS). The SS phase function predicts the light scattering properties for all forward angles, namely 0° ≤*θ* ≤90° as function of the real index of refraction and the slope parameter of the hyperbolic distribution. The phase function was divided

. The authors state that

*,*

�*<sup>x</sup>* <sup>þ</sup> *bm*ð Þþ *<sup>x</sup> cm,* (61)

*<sup>m</sup>y* þ *km,* (63)

*<sup>y</sup>*<sup>2</sup> <sup>þ</sup> *em* sin 4ð Þþ <sup>π</sup>*<sup>y</sup> <sup>f</sup> <sup>m</sup><sup>y</sup>* <sup>þ</sup> *gm,* (62)

*<sup>y</sup>*<sup>2</sup> <sup>þ</sup> *om* sin 4ð Þþ <sup>π</sup>*<sup>y</sup> pmy* <sup>þ</sup> *qm,* (64)

*x* ¼ μ � 3*,* (65) *y* ¼ *n* � 1*,* (66)

(60)

from Petzold's data for harbor waters, we get very close probabilities over all

*3.2.5 Sahu and Shanmugam (SS) phase function with flat back approximation*

this division was necessary, since the phase function is highly nonlinear with a

� � <sup>¼</sup> *<sup>P</sup>*1ð Þ ln ð Þ*<sup>θ</sup>* <sup>2</sup> <sup>þ</sup> *<sup>P</sup>*2ð Þþ ln ð Þ*<sup>θ</sup> <sup>P</sup>*<sup>3</sup> <sup>0</sup>*:*1<sup>∘</sup> <sup>≤</sup> *<sup>θ</sup>* <sup>≤</sup>5<sup>∘</sup>

*<sup>y</sup>*<sup>2</sup> <sup>þ</sup> *im* sin 4ð Þþ <sup>π</sup>*<sup>y</sup> <sup>j</sup>*

*Comparison of the FF and the Petzold average phase function over all scattering angles and for small scattering*

*Pm* ¼ *ame*

*am* <sup>¼</sup> *dm*

*bm* <sup>¼</sup> *hm*

*cm* <sup>¼</sup> *lm*

. This phase function is given by:

*<sup>P</sup>*1ð Þ ln ð Þ*<sup>θ</sup>* <sup>3</sup> <sup>þ</sup> *<sup>P</sup>*2ð Þ ln ð Þ*<sup>θ</sup>* <sup>2</sup> <sup>þ</sup> *<sup>P</sup>*3ð Þþ ln ð Þ*<sup>θ</sup> <sup>P</sup>*<sup>4</sup> <sup>5</sup><sup>∘</sup> <sup>≤</sup>*<sup>θ</sup>* <sup>≤</sup> <sup>90</sup><sup>∘</sup>

in two parts, one from 0.1 to 5° and the other from 5 to 90°

*Monte Carlo Radiative Transfer Modeling of Underwater Channel*

comparison is shown in **Figure 8**.

*DOI: http://dx.doi.org/10.5772/intechopen.85961*

scattering angles.

change of slope around 5°

8 < :

log eβð Þ*θ*

where

**Figure 8.**

*angles.*

**113**

### **Figure 7.**

*Comparison of the empirical phase function, the TTHG and HG, against the Petzold's measurements for small and large angles.*

harbor waters with *<sup>b</sup>* <sup>¼</sup> <sup>1</sup>*:*8241 and *<sup>ω</sup>*<sup>0</sup> <sup>¼</sup> <sup>0</sup>*:*83, *<sup>P</sup> <sup>θ</sup>* <sup>≤</sup> 10° � � <sup>¼</sup> <sup>0</sup>*:*64; while in the empirical phase function, we found *<sup>P</sup> <sup>θ</sup>* <sup>≤</sup>10° � � <sup>¼</sup> <sup>0</sup>*:*84.

The above presented divergence resulted in a small difference in the final received power; however, the empirical phase function largely overestimated the channel bandwidth, especially in harbor waters. This is due to the smaller scattering angles generated by this empirical phase function, causing the photons to travel in a straighter line when compared to the case of scattering with the CDF computed from Petzold's data.

### *3.2.4 The Fournier-Forand phase function*

As seen the previous sections, several analytical and empirical formulae were used as empirical fits for ocean waters but none, specially the analytical ones provide a perfect fit for all scattering angles.

In [24], Fournier and Forand derived an analytic expression for the phase function of a set of particles that have a hyperbolic particle size distribution, with each particle scattering according to the anomalous diffraction theory. The Fournier-Forand (FF) phase function is given by:

$$\begin{split} \tilde{\beta}\_{FF}(\theta) &= \frac{1}{4\pi(1-\delta)^{2}\delta^{\mathrm{v}}} \Bigg\{ [\mathbf{v}(1-\delta)-(1-\delta^{\mathrm{v}})]+[\delta(1-\delta^{\mathrm{v}})-\mathbf{v}(1-\delta)]\sin^{-2}\left(\frac{\theta}{2}\right) \Bigg\}, \\ &+\frac{1-\delta\_{180}^{\mathrm{v}}}{16\pi(\delta\_{180}-1)\delta\_{180}^{\mathrm{v}}} (3\cos^{2}\theta-1), \end{split} \tag{57}$$

where

$$\mathbf{v} = \frac{\mathfrak{Z} - \mu}{2},\tag{58}$$

$$\delta = \frac{4}{\Im(n-1)} \sin^2 \left(\frac{\theta}{2}\right). \tag{59}$$

Being *n* in Eq. (59) the real index of refraction of the particles, μ in Eq. (58) the slope parameter of the hyperbolic distribution, and δ<sup>180</sup> in Eq. (57) is δ evaluated at *<sup>θ</sup>* <sup>¼</sup> 180° .

*Monte Carlo Radiative Transfer Modeling of Underwater Channel DOI: http://dx.doi.org/10.5772/intechopen.85961*

It can be seen in [25] that with *n* ¼ 1*:*16 and *μ* ¼ 3*:*4586, the FF phase function gives a good fit to the average Petzold phase function over all scattering angles. This comparison is shown in **Figure 8**.

When comparing the CDF computed from the FF phase function and the CDF from Petzold's data for harbor waters, we get very close probabilities over all scattering angles.

### *3.2.5 Sahu and Shanmugam (SS) phase function with flat back approximation*

In [25], a semi-analytical model for a phase function was proposed by Sahu and Shanmugam (SS). The SS phase function predicts the light scattering properties for all forward angles, namely 0° ≤*θ* ≤90° as function of the real index of refraction and the slope parameter of the hyperbolic distribution. The phase function was divided in two parts, one from 0.1 to 5° and the other from 5 to 90° . The authors state that this division was necessary, since the phase function is highly nonlinear with a change of slope around 5° . This phase function is given by:

$$\log\left(\widetilde{\boldsymbol{\beta}}(\boldsymbol{\theta})\right) = \begin{cases} \boldsymbol{P}\_{1}\left(\ln\left(\boldsymbol{\theta}\right)\right)^{2} + P\_{2}\left(\ln\left(\boldsymbol{\theta}\right)\right) + P\_{3}\left0.1^{\*} \leq \boldsymbol{\theta} \leq \boldsymbol{\mathfrak{S}}^{\*} \\ \boldsymbol{P}\_{1}\left(\ln\left(\boldsymbol{\theta}\right)\right)^{3} + P\_{2}\left(\ln\left(\boldsymbol{\theta}\right)\right)^{2} + P\_{3}\left(\ln\left(\boldsymbol{\theta}\right)\right) + P\_{4}\left\mathfrak{S}^{\*} \leq \boldsymbol{\theta} \leq \boldsymbol{\mathfrak{M}}^{\*} \end{cases},\tag{60}$$

where

harbor waters with *<sup>b</sup>* <sup>¼</sup> <sup>1</sup>*:*8241 and *<sup>ω</sup>*<sup>0</sup> <sup>¼</sup> <sup>0</sup>*:*83, *<sup>P</sup> <sup>θ</sup>* <sup>≤</sup> 10° � � <sup>¼</sup> <sup>0</sup>*:*64; while in the

The above presented divergence resulted in a small difference in the final received power; however, the empirical phase function largely overestimated the channel bandwidth, especially in harbor waters. This is due to the smaller scattering angles generated by this empirical phase function, causing the photons to travel in a straighter line when compared to the case of scattering with the CDF computed

*Comparison of the empirical phase function, the TTHG and HG, against the Petzold's measurements for small*

As seen the previous sections, several analytical and empirical formulae were used as empirical fits for ocean waters but none, specially the analytical ones

In [24], Fournier and Forand derived an analytic expression for the phase function of a set of particles that have a hyperbolic particle size distribution, with each particle scattering according to the anomalous diffraction theory. The Fournier-

<sup>v</sup> <sup>¼</sup> <sup>3</sup> � <sup>μ</sup>

<sup>3</sup>ð Þ *<sup>n</sup>* � <sup>1</sup> <sup>2</sup> *sin*<sup>2</sup> *<sup>θ</sup>*

Being *n* in Eq. (59) the real index of refraction of the particles, μ in Eq. (58) the slope parameter of the hyperbolic distribution, and δ<sup>180</sup> in Eq. (57) is δ evaluated at

<sup>δ</sup>*<sup>v</sup>* v 1ð Þ� � <sup>δ</sup> <sup>1</sup> � <sup>δ</sup><sup>v</sup> <sup>½</sup> ð Þ� þ <sup>δ</sup> <sup>1</sup> � <sup>δ</sup><sup>v</sup> ½ � ð Þ� v 1ð Þ � <sup>δ</sup> sin�<sup>2</sup> *<sup>θ</sup>*

2 � �

� � � �

<sup>2</sup> *,* (58)

*:* (59)

2

(57)

empirical phase function, we found *<sup>P</sup> <sup>θ</sup>* <sup>≤</sup>10° � � <sup>¼</sup> <sup>0</sup>*:*84.

*Wireless Mesh Networks - Security, Architectures and Protocols*

from Petzold's data.

**Figure 7.**

*and large angles.*

<sup>e</sup>β*FF*ð Þ¼ *<sup>θ</sup>* <sup>1</sup>

where

*<sup>θ</sup>* <sup>¼</sup> 180°

**112**

.

þ

*3.2.4 The Fournier-Forand phase function*

provide a perfect fit for all scattering angles.

Forand (FF) phase function is given by:

<sup>1</sup> � <sup>δ</sup><sup>v</sup> 180 <sup>16</sup>π δð Þ <sup>180</sup> � <sup>1</sup> <sup>δ</sup>*<sup>v</sup>*

180

3cos2 *<sup>θ</sup>* � <sup>1</sup> � �*,*

<sup>δ</sup> <sup>¼</sup> <sup>4</sup>

<sup>4</sup>πð Þ <sup>1</sup> � <sup>δ</sup> <sup>2</sup>

$$P\_m = a\_m e^{-\chi} + b\_m(\chi) + c\_m,\tag{61}$$

$$a\_m = \frac{d\_m}{\mathcal{Y}^2} + e\_m \sin\left(4\pi\mathcal{Y}\right) + f\_m \mathcal{Y} + \mathcal{g}\_{m'} \tag{62}$$

$$b\_m = \frac{h\_m}{\mathcal{y}^2} + i\_m \sin\left(4\pi\mathcal{y}\right) + j\_m\mathcal{y} + k\_m,\tag{63}$$

$$c\_m = \frac{l\_m}{\mathcal{Y}^2} + o\_m \sin \left(4\pi\mathcal{Y}\right) + p\_m \mathcal{Y} + q\_m. \tag{64}$$

$$
\mathfrak{x} = \mathfrak{y} - \mathfrak{Z},
\tag{65}
$$

$$
\mathbf{y} = \mathbf{n} - \mathbf{1},
\tag{66}
$$

*Comparison of the FF and the Petzold average phase function over all scattering angles and for small scattering angles.*

with *<sup>m</sup>* <sup>¼</sup> <sup>1</sup>*,* <sup>2</sup>*,* 3 for 0*:*1° <sup>≤</sup>*<sup>θ</sup>* <sup>≤</sup> 5° and *<sup>m</sup>* <sup>¼</sup> <sup>1</sup>*,* <sup>2</sup>*,* <sup>3</sup>*,* 4 for 5° <sup>≤</sup>*<sup>θ</sup>* <sup>≤</sup>90° . The values for the regression coefficients *dm, em, f <sup>m</sup>, gm, hm, im, jm, km, lm, pm*, and *qm* may be found in [26] and are not shown here for brevity.

Sahu found, in his most recent works [14, 26], that the set of values, *n* ¼ 1*:*16 and μ ¼ 3*:*4319 gives the best fit to Petzold average phase function. Using Eqs. (60)–(66) and with the values mentioned above for *n* and μ, we compared the obtained values with the Petzold average phase function and found that it indeed provides an excellent fit at all scattering angles, but specially at small angles where most of the scattering will happen. Results for both parts of the phase function, 0.1–5° and 5–90° , are presented in **Figure 9**.

### *3.2.6 Computing phase functions for Monte Carlo simulations*

As seen in [27], phase functions satisfy the normalization condition:

$$2\pi \int\_0^\pi \tilde{\beta}(\theta) \sin \theta d\theta = 1,\tag{67}$$

and to determine the polar scattering angle:

$$q = 2\pi \int\_0^{\theta} \tilde{\mathbb{\hat{B}}}(\theta) \sin \theta \mathrm{d}\theta,\tag{68}$$

In **Table 2** and **Figures 10** and **11**, we compare the CDF obtained for three phase functions, the SS, FF, and the TTHG with *g* ¼ 0*:*9809 against the CDF for Petzold's average phase function. In **Figure 11**, we show the abscissa in logarithmic scale for a better visualization of the differences between phase functions for small values of *θ*. In this table, we can see that both the FF and SS phase functions are a great approximation for the Petzold's average phase function, while the TTHG heavily overestimates the amount of scattering up to 10° and underestimates the amount of

*θ* **(°)** *P***(***θ***PETZ)** *P***(***θ***FF)** *P***(***θ***SS) P(***θ***TTHG)** *:*126 0.009094 0.012355 0.004059 0.002382 *:*2 0.031941 0.041314 0.025992 0.011982 *:*251 0.046044 0.057930 0.039860 0.020878 *:*316 0.062621 0.076605 0.056320 0.034634 *:*398 0.081934 0.097344 0.075422 0.055321 *:*501 0.104274 0.120254 0.097285 0.085451 *:*631 0.129946 0.145647 0.122155 0.127856 *:*794 0.158838 0.173561 0.149991 0.184165 0.191217 0.204415 0.181067 0.254854 *:*512 0.352290 0.358677 0.335353 0.596932 *:*012 0.498215 0.505639 0.472517 0.788187 0.659192 0.665537 0.630018 0.897428 0.891501 0.882151 0.892093 0.974829 0.958125 0.951993 0.960280 0.989702 0.986110 0.983167 0.987188 0.996168 0.994130 0.992771 0.994681 0.998674 0.997280 0.996631 0.997562 0.999169 90 1 1 1 1

*Monte Carlo Radiative Transfer Modeling of Underwater Channel*

*DOI: http://dx.doi.org/10.5772/intechopen.85961*

effects of those deviations will be larger for waters where scattering is the main

One of the most useful information regarding the underwater channel that may be provided by the Monte Carlo simulation, which is unavailable in most simple models like the Beer law, is the channel temporal dispersion. As seen before in this chapter, seawater is a dispersive medium in which light will suffer the effect of multiple scattering events. Scattering will cause the photons to arrive at the receiver at different time intervals causing temporal dispersion and a penalty in the channel bandwidth. In recent years, several researchers employed the double gamma functions to model the impulse response in underwater channels [15, 28]. Despite the accuracy of the double gamma function models, we find that the method described

. As one may expect, the

scattering in intermediate angles, i.e., between 30° and 80°

**3.3 Temporal dispersion and the underwater channel bandwidth**

attenuation factor.

**115**

**Table 2.**

*3.3.1 Temporal dispersion*

*Computed CDF for selected phase functions.*

where *q* is a uniformly distributed random number between 0 and 1.

However, given the shape of most phase functions, it is more efficient [27] to use the equation for the phase function to build up a table of *θ* vs. eβð Þ*θ* , compute a CDF making use of the normalization condition, and then interpolate into the CDF the value of *θ* for each random number, *q*, drawn. As seen in Section 3.2.3, Kirk [22] computed a CDF for Petzold's measurements in San Diego Harbor for angles between 2*:*5 and 177*:*5° in intervals of 5° .

For the results presented in **Table 2**, we computed a CDF directly from Petzold average phase function extracted from [3] and used the same method described in [27] to obtain values for all angles between 0 and 90° . We chose to proceed with numerical integration for this range, and not from 0 to 180° , because as mentioned in Section 3.2.4, the SS phase function provides values only for forward scattering angles.

**Figure 9.** *Comparison between the SS and the Petzold average phase function over the forward scattering angles.*


*Monte Carlo Radiative Transfer Modeling of Underwater Channel DOI: http://dx.doi.org/10.5772/intechopen.85961*

with *<sup>m</sup>* <sup>¼</sup> <sup>1</sup>*,* <sup>2</sup>*,* 3 for 0*:*1° <sup>≤</sup>*<sup>θ</sup>* <sup>≤</sup> 5° and *<sup>m</sup>* <sup>¼</sup> <sup>1</sup>*,* <sup>2</sup>*,* <sup>3</sup>*,* 4 for 5° <sup>≤</sup>*<sup>θ</sup>* <sup>≤</sup>90°

and μ ¼ 3*:*4319 gives the best fit to Petzold average phase function. Using

As seen in [27], phase functions satisfy the normalization condition:

eβð Þ*θ* sin*θ*d*θ* ¼ 1*,* (67)

eβð Þ*θ* sin*θ*d*θ,* (68)

. We chose to proceed with

, because as mentioned

, are presented in **Figure 9**.

*3.2.6 Computing phase functions for Monte Carlo simulations*

and to determine the polar scattering angle:

between 2*:*5 and 177*:*5° in intervals of 5°

[27] to obtain values for all angles between 0 and 90°

numerical integration for this range, and not from 0 to 180°

2π ðπ 0

*q* ¼ 2π

ð*θ* 0

where *q* is a uniformly distributed random number between 0 and 1.

. For the results presented in **Table 2**, we computed a CDF directly from Petzold average phase function extracted from [3] and used the same method described in

in Section 3.2.4, the SS phase function provides values only for forward scattering

*Comparison between the SS and the Petzold average phase function over the forward scattering angles.*

However, given the shape of most phase functions, it is more efficient [27] to use the equation for the phase function to build up a table of *θ* vs. eβð Þ*θ* , compute a CDF making use of the normalization condition, and then interpolate into the CDF the value of *θ* for each random number, *q*, drawn. As seen in Section 3.2.3, Kirk [22] computed a CDF for Petzold's measurements in San Diego Harbor for angles

in [26] and are not shown here for brevity.

*Wireless Mesh Networks - Security, Architectures and Protocols*

0.1–5° and 5–90°

angles.

**Figure 9.**

**114**

the regression coefficients *dm, em, f <sup>m</sup>, gm, hm, im, jm, km, lm, pm*, and *qm* may be found

Sahu found, in his most recent works [14, 26], that the set of values, *n* ¼ 1*:*16

Eqs. (60)–(66) and with the values mentioned above for *n* and μ, we compared the obtained values with the Petzold average phase function and found that it indeed provides an excellent fit at all scattering angles, but specially at small angles where most of the scattering will happen. Results for both parts of the phase function,

. The values for

### **Table 2.** *Computed CDF for selected phase functions.*

In **Table 2** and **Figures 10** and **11**, we compare the CDF obtained for three phase functions, the SS, FF, and the TTHG with *g* ¼ 0*:*9809 against the CDF for Petzold's average phase function. In **Figure 11**, we show the abscissa in logarithmic scale for a better visualization of the differences between phase functions for small values of *θ*.

In this table, we can see that both the FF and SS phase functions are a great approximation for the Petzold's average phase function, while the TTHG heavily overestimates the amount of scattering up to 10° and underestimates the amount of scattering in intermediate angles, i.e., between 30° and 80° . As one may expect, the effects of those deviations will be larger for waters where scattering is the main attenuation factor.

### **3.3 Temporal dispersion and the underwater channel bandwidth**
