1. Introduction

Ultraviolet (UV) communications allow non-light-of-sight (NLOS) transmission of information through scattering of UV light by the molecules and aerosols in the atmosphere [1–5]. They do not require accurate pointing and tracking found in typical outdoor optical communications [6]. In NLOS systems, transmit photons can arrive at the receiver via single or multiple scattering in the atmosphere. In the single scatter NLOS models, photons between the transmitter and the receiver are assumed to be scattered only once. This requires that the transmit and receive cones intersect in a common volume (CV) [7]. Although path loss evaluations in single scatter models are still computationally challenging, simplified models have also been studied, e.g., [8,9].

A multiple scattering NLOS channel model is more accurate, comprehensive, and essential in configurations where there is little or no CV. A well-recognized approach toward this is to use a full Monte Carlo simulation [1]. However, this is computationally unfeasible for high path losses, a typical scenario in UV applications. In [2], a novel Monte Carlo method (MCM) is described using a combined approach of random number generations and analytic calculations. An extension of this work is given in [10] with further insights along with an approximate model which is valid at short distances.

In general, Monte Carlo methods are computationally intensive, particularly when used with Mie theory phase functions [11]. An integral model for higher order scattering is given in [12]. This non-Monte Carlo approach simplifies only for the case of narrow beams. Thus, to our best knowledge, there is no non-Monte Carlo method that can significantly reduce the computational requirements for multiple scattering events in general transmitter/receiver configurations. Therefore, in this paper, we present a probability-based sampling method (PSM) that completely eliminates random number generations. PSM is focused on single and double scattering events, noting that these two events closely represent the total signal [2]. In fact, when there is a large CV, single scattering alone accurately represents the received signal. Due to its low complexity, PSM can find applications in analyzing networks of transmitters and receivers, e.g., [4,13,14], which is challenging, if not impossible, to study using MCM. PSM can be used to study hybrid systems [15]. PSM can also serve as a quick validation for the MCM results. It is also possible to extend the double scattering idea of PSM to higher order scattering at the cost of numerical complexity, but this is left as a potential future work.

This paper makes the following contributions. First, we present a mathematical framework to analyze single and double scattering events for NLOS UV channels. Since a photon can travel a long distance before undergoing a scattering event, a simple fixed uniform sampling of the photon path with a small step size is impractical. Similarly, the sampling of the scattering angles is also unclear. Our novel PSM overcomes these challenges using probability-based sampling. Note that PSM is different from importance sampling, e.g., [16], that still requires Monte-Carlo random number generations, which are completely eliminated in the proposed PSM. Second, we provide closed form expressions for the sample photon emission angles as well as the intersection points on the receiver cone for realistic non-planar settings. The computational requirements are reduced by a factor of hundreds over Monte Carlo simulation methods. In addition, the channel impulse responses obtained by PSM and MCM are shown to match well. Finally, PSM is applied to analyze the effects of dust and fog aerosols and useful insights are derived. At high aerosol densities, dust and fog provide different amounts of path losses, and for given transmitter/receiver inclinations, a near-planar set up shows more impact from the aerosols.

Notations: We use bold lower case letters for column vectors. The notations $[\cdot ]^{T}$, $\| \cdot \|$, $|\cdot |$, $\mbox {Re}(\cdot )$, $\ln (\cdot )$ and $\overrightarrow {AB}$ denote transpose, Euclidean norm, magnitude, real part, natural logarithm, and a vector from point $A$ to point $B$ respectively. A new spatial direction, denoted by an angle pair $(\theta , \phi )$ with respect to a vector $\textbf {u}$, means that the new direction is inclined at an angle of $\theta$ with respect to $\textbf {u}$, and $\phi$ is the angle that the new direction makes with respect to a reference direction when projected on a plane orthogonal to $\textbf {u}$. The components of a vector $\textbf {u}_{i}$ are denoted by extending its subscripts with $x$, $y$ and $z$ as $u_{i,x}$, $u_{i,y}$ and $u_{i,z}$.

2. UV channel analysis

We consider a transmitter/receiver configuration as shown in Fig. 1 in the Cartesian coordinate system. The location of the receiver (R) is $[0, 0, 0]$ while the transmitter (T) is located at $[0, r, 0]$ and so the range is $r$. The transmitter inclination and azimuth angles are denoted by $\theta _{T}$ and $\phi _{T}$ respectively. Similarly, $\theta _{R}$ and $\phi _{R}$ respectively denote the receiver inclination and azimuth angles. The inclination angle is measured with respect to the positive $z$-axis and the azimuth angle is measured with respect to the positive $x$-axis. The transmit beamwidth is $\beta _{T}$ and the receiver field of view is $\beta _{R}$. The transmitter and receiver axial directions are denoted by the unit vectors $\textbf {u}_{T}$ and $\textbf {u}_{R}$ respectively. We assume that the field of view of the receiver is uniform. Further, the detector area is assumed to be a disc perpendicular to the photon direction as in [2] and $S_{R}$ is the disc area.

 

Fig. 1. Transmitter/receiver configuration and the display of a single scattering event leading to detection at the receiver.

Download Full Size | PDF

When a photon leaves the transmitter, it travels an exponentially distributed random distance before interacting with air molecules or aerosols. Such an interaction results in either the photon being absorbed or scattered. The corresponding probabilities are $k_{a}/k_{e}$ and $k_{s}/k_{e}$ respectively, where $k_{a}$, $k_{s}$ and $k_{e}$ are the absorption, scattering and extinction coefficients respectively, and $k_{e}=k_{a}+k_{s}$. These coefficients for specific atmospheric conditions are discussed in Section 4.

2.1 Photon emission directions

We consider a UV source with uniformly distributed radiation through the transmit aperture. This can be easily extended to other distributions if needed. Let $(\theta _{0}, \phi _{0})$ denote a photon emission direction with respect to $\textbf {u}_{T}$, where $\theta _{0}$ denotes the angle of the photon direction with respect to $\textbf {u}_{T}$ and $\phi _{0}$ is the angle that this direction makes with respect to a reference direction when projected on a plane $P_{T}$ orthogonal to $\textbf {u}_{T}$. Define $\kappa =1/(1-\cos (\beta _{T}/2))$. For a uniformly distributed radiation source, the probability density function (PDF) is

Consider the orthogonal plane $P_{T}$. This plane cuts the surface of the transmit cone into the outermost circle as shown in Fig. 2(a). Let $A_{0}$ be the point of intersection between $P_{T}$ and the transmit axis. In the absence of any scattering, a photon emitted from the transmitter $T$ must pass through $P_{T}$. Therefore, the photon directions from the transmitter can be equivalently described by the corresponding points on plane $P_{T}$. From the infinite number of possible photon emission directions from $T$, we want to sample only a finite number of representative directions. Equivalently, select a finite number of sample points on $P_{T}$ that can accurately represent photon emissions. Each such point carries the probability of photon emission through its neighborhood.

 

Fig. 2. Illustration of sampling and grid circles for $N_{c}=2$: (a) intersection of transmitter cone with an orthogonal plane $P_{T}$, and (b) various circles on plane $P_{T}$. The dashed lines represent grid circles and the continuous lines represent sampling circles.

Download Full Size | PDF

To select the sample points on $P_{T}$, we consider a number $N_{c}$ of concentric sampling circles centered with $A_{0}$ on plane $P_{T}$. Let there be $N_{i}$ sample points spread equally on the circumference of the $i$-th sampling circle of radius $r_{i}$, $i=1, 2, \ldots , N_{c}$. We also set $N_{0}=1$ for the axial direction represented by $A_{0}$. Thus, there are a total of $N_{s}=\sum _{i=0}^{N_{c}}N_{i}$ sample directions of photon emission.

The sampling circles are arranged within $N_{c}+1$ grid circles of radii $r_{i}'$, so that $r_{i}' < r_{i} < r_{i+1}'$, $i=1, 2, \ldots , N_{c}$. Figure 2 displays an example of two sampling and three grid circles. All photon emissions through the annular region between the grid circles of radii $r_{i}'$ and $r_{i+1}'$ are represented by the $N_{i}$ sample points on the sampling circle of radius $r_{i}$. Each sample point represents a probability of photon emissions of $1/N_{s}$ through its neighborhood. Thus, the point $A_{0}$ represents photon emissions through the first circular region of radius $r_{1}'$ with a probability of $P_{0}=1/N_{s}$. Let $\theta _{i}'$ be the angle $\angle A_{0}T A_{i}'$, where $A_{i}'$ is an arbitrary point on the grid circle of radius $r_{i}'$. Then using Eq. (1), $P_{0}=\kappa (1-\cos \theta _{1}')$ giving $\theta _{1}'=\cos ^{-1}(1-1/\kappa N_{s})$. The probability of photon emission through the $i$-th annular region between grid circles of radii $r_{i}'$ and $r_{i+1}'$ is $P_{i}=\kappa (\cos \theta _{i}'-\cos \theta _{i+1}')$. Since there are $N_{i}$ uniformly distributed points on the corresponding sampling circle of radius $r_{i}$, the photon emission probability through each segment of the annular region represented by a sample point is $P_{i}/N_{i}=1/N_{s}$. This leads to $\cos \theta _{i+1}'$ $=\cos \theta _{i}'-N_{i}/(\kappa N_{s})$. Using this recursive relation, the grid circle angles can be obtained as

Each photon emission direction is identified as $(\theta _{i},\phi _{i,k})$ with reference to the transmit axial direction ($\textbf {u}_{T}$), where $\phi _{i,k}=2\pi k/N_{i}$, $k=0, 1, \ldots , N_{i}-1$. Since the values of $N_{i}$ are not known beforehand, we present an algorithm (Appendix I) to obtain the values of $\theta _{i}$ and $N_{i}$ through iterations. A photon direction unit vector $\textbf {u}_{i}$, with its Cartesian components $(u_{i,x},u_{i,y},u_{i,z})$, corresponding to point $A_{i}$ can be obtained using the equations in [2,17] for $|u_{T,z}| \neq 1$ as

When $|u_{T,z}| = 1$, they become $u_{i,x}=\sin \theta _{i} \cos \phi _{i,k}$, $u_{i,y}=\sin \theta _{i}\sin \phi _{i,k}$, and $u_{i,z}=u_{T,z}\cos \theta _{i}$.

2.2 Single scattering

We now consider the probability of detection of a photon under a single scattering event when the photon is emitted along a sample photon direction. Single scattering detection occurs only if the transmit and the receive beams have a CV as in Fig. 1. Therefore, a photon from the transmitter has to reach the CV before encountering any scattering. In general, a photon along a sample path originating from $T$ may or may not reach the CV. If the photon enters/leaves a CV without encountering any scattering then we have entry/exit points as shown in Fig. 1. Since the entry/exit point is on the surface of the receiver cone, the angle formed at the receiver by the entry/exit point and the vector $\textbf {u}_{R}$ must be $\beta _{R}/2$. Imagine a plane $P_{T}$ orthogonal to $\textbf {u}_{T}$ and passing through the entry or the exit point corresponding to a sample photon path from $T$. The entry/exit point can be referred to as $A_{i}$. Let $\textbf {p}$ denote the location vector of $A_{i}$. Therefore, $\mathbf {p}=\overrightarrow {RA_{i}}=[0 \hspace {2mm} r \hspace {2mm} 0]^{T}+\overrightarrow {TA_{0}}+\overrightarrow {A_{0}A_{i}}$. Writing $\overrightarrow {TA_{0}}=s\mathbf {u}_{T}$, where $s$ is the length from T along the transmit axis, and denoting the unit vector along $\overrightarrow {A_{0}A_{i}}$ as $\tilde {\textbf {u}}_{i}$, we can write $\textbf {p} = \textbf {r}+s(\mathbf {u}_{T}+\tilde {s}_{i} \tilde {\textbf {u}}_{i})$, where $\textbf {r}=[0 \hspace {2mm} r \hspace {2mm} 0]^{T}$, $\tilde {s}_{i}=\tan \theta _{i}$ and $\theta _{i}$ is the angle $\angle A_{0}TA_{i}$. Define the normalized vector $\textbf {u}_{i}=(\mathbf {u}_{T}+\tilde {s}_{i} \tilde {\textbf {u}}_{i})/\| \mathbf {u}_{T}+\tilde {s}_{i} \tilde {\textbf {u}}_{i}\|$. Then $\textbf {p}=\textbf {r}+s_{i} \textbf {u}_{i}$, where $s_{i}=s\| \mathbf {u}_{T}+\tilde {s}_{i} \tilde {\textbf {u}}_{i}\|$. Since the entry/exit point $\textbf {p}$ lies on the receiver cone surface, we get

Thus, in case of single scattering, the location of the ESP from the receiver for the given sample path number $i$ and ESP number $k$ is $\textbf {d}_{1}=\textbf {r}+\bar {s}_{i,k}\textbf {u}_{i}$ for $k=1, 2, \ldots N_{r}$. Define $E_{k}$ as the combined event of $E_{k}'$ and the event that the photon encounter results in scattering. Let us also define an event, $E_{R}$ as the event that the scattered photon reaches the receiver with no further scattering. The event $D_{1}$ occurs, given $R_{i}$, if one of the events, $E_{k}$, $k=1,2,\ldots N_{r}$, occurs and $E_{R}$ occurs. Since the events, $E_{k}$, $k=1, 2, \ldots , N_{r}$, are mutually exclusive, the probability $P[D_{1}|R_{i}]=$ $\sum _{k=1}^{N_{r}} P[E_{k}, E_{R}|R_{i}]$ $=\sum _{k=1}^{N_{r}} P[E_{R}|E_{k}, R_{i}]P[E_{k}|R_{i}]$. Obviously, $P[E_{k}|R_{i}]=$ $(k_{s}/k_{e})P[E_{k}'|R_{i}]$. The probability $P[E_{k}'|R_{i}]$ is $\int _{a_{k}}^{b_{k}} k_{e} \exp (-k_{e}x)dx$ $=\exp (-k_{e}a_{k})-\exp (-k_{e}b_{k})$, where $a_{k}= s_{i,k-1}'$ and $b_{k}=s_{i,k}'$. However, by our design, each segment has equal probability, and hence $P[E_{k}'|R_{i}]$ is $P(i) = (1/N_{r})(\exp (-k_{e}s_{i}^{(1)})-\exp (-k_{e}s_{i}^{(2)}))$ which saves numerical operations. From [2], the probability $P[E_{R}|E_{k}, R_{i}]$ is approximately obtained as

It is important to observe that many of the probability calculations are eliminated by carefully using sampling operations based on equal probabilities.

2.3 Double scattering

Signal detection through double scattering requires that the first scattering occurs in the transmit cone and the second scattering occurs in the receive cone. This is shown in Fig. 3. We discretize the sample photon path, i.e., the vector direction $\textbf {u}_{i}$, into $N_{t}$ transmitter segments. The segment end points and the ESP for each segment can be obtained from Eqs. (9), (10) and (11) using $s_{i}^{(1)}=0$ and $s_{i}^{(2)}=\infty$. The ESP is then

 

Fig. 3. Example of a double scattering event leading to photon detection at the receiver.

Download Full Size | PDF

Consider a unit vector $\textbf {v}_{k}$ from the tip of the vector $\bar {s}_{i,n}\textbf {u}_{i}$ which is the location for the first scattering event. Let $\textbf {v}_{k}$ be oriented at an angle pair of $(\bar {\theta }_{k'}, \bar {\phi }_{k''})$ with respect to $\textbf {u}_{i}$, where $k'=1, 2, \ldots , N_{a}$, and $k''=1, 2, \ldots , N_{p}$. We will call this direction as the $k$-th angle pair direction, $k=1, 2, \ldots , N_{a}N_{p}$. Note that $\textbf {v}_{k}$ can be obtained using Eqs. (4)–(6) by replacing $\textbf {u}_{T}$ with $\textbf {u}_{i}$ and $\textbf {u}_{i}$ with $\textbf {v}_{k}$. Suppose the direction of $\textbf {v}_{k}$ meets the surface of the receiver cone at a distance $b_{i,n,k}$ from the tip of vector $\bar {s}_{i,n}\textbf {u}_{i}$. This provides the entry/exit points on the receiver cone. The location vector for this point is $\textbf {p}=\textbf {r}+\bar {s}_{i,n}\textbf {u}_{i}+b_{i,n,k} \textbf {v}_{k}$. Defining $\textbf {q}=\textbf {r}+\bar {s}_{i,n}\textbf {u}_{i}$, we get $\|\textbf {p}\|^{2}=\|\textbf {q}\|^{2}+2b_{i,n,k} \gamma +b_{i,n,k}^{2}$, where $\gamma =\textbf {q}^{T}\textbf {v}_{k}$. Similar to the single scatter case, the entry/exit point on the surface of the receiver cone requires $\textbf {u}_{R}^{T}\textbf {p}=\|\textbf {p}\| m_{R}$. Squaring this equation gives the following quadratic equation,

A) $b_{i,n,k}^{(m)}$ is positive for both $m$. Under this scenario, we have two cases: A1) $\textbf {q}$ is ICV. In this case, the vector $\textbf {v}_{k}$ meets the true and the virtual receiver cone surface. The virtual receiver cone is the mirror image of the true receiver cone with respect to the ground plane. Therefore, we set $b_{i,n,k}^{(1)}=0$, and $b_{i,n,k}^{(2)}$ to the smaller of the two solutions of the quadratic equation. The first solution is zero because the second scattering can occur immediately after the first scattering within the CV. A2) $\textbf {q}$ is OCV. In this case, accept both solutions as the entry/exit point solutions.

B) $b_{i,n,k}^{(m)}$ is negative for one or both $m$. Again, two cases arise: B1) For the ICV case, if both solutions are negative, then one solution is set to zero and the other solution is set to infinity because the photon path is moving away from the receiver. On the other hand, if the solutions have different signs, then $b_{i,n,k}^{(1)}$ is set to zero and the positive solution from the quadratic equation is set as $b_{i,n,k}^{(2)}$. This is because $b_{i,n,k}^{(2)}$ is the exit point and the second scattering can occur anywhere between the first scattering point and the exit point for the receiver to be able to detect it. B2) For the OCV case, we set $b_{i,n,k}^{(2)}$ to infinity and the positive solution from the quadratic equation is set to $b_{i,n,k}^{(1)}$ when the solutions have different signs. In this case, $b_{i,n,k}^{(1)}$ is the entry point, and the exit point is at infinity. When the solutions are both negative for the OCV case, they are not used.

Finally, we discretize the photon path inside the receiver cone into $N_{r}$ receiver segments. The photon path inside the receiver cone and a second scattering event example is shown in Fig. 3. Each segment has an equal photon absorption/scattering probability. The $l$-th segment ESP is obtained using Eq. (11) as

3. Discussion on computational complexity

We compare PSM and MCM in terms of numerical operations needed for path loss calculations. Since many computations in the MCM are decided based on random numbers, it is challenging to provide an exact number of operations. To simplify the presentation, functions, such as logarithms, trigonometric functions, square roots, exponential functions, are simply referred to as nonlinear functions (NLF), and we list only multiplications (include divisions) and NLF computations. Some operations in PSM, e.g., $\cos (\beta _{R}/2)$, $u_{T,x}u_{T,z}$, can be computed a priori and are not included in the computations. The computations within phase functions are also not included assuming similar impact on both methods. The sample photon direction selections involve negligible operations and can also be performed a priori.

Consider single scattering first. The main tasks are listed in Table 1. Each task includes all related operations. For example, solving Eq. (8) also includes computing $\textbf {u}_{i}$ and $m_{i}$, $i=1, \ldots , N_{s}$. Computing $\textbf {d}_{1}$ also includes evaluating $\bar {s}_{i,k}$. The total number of NLF evaluations is about $4N_{s}N_{r}+7N_{s}$ and the number of multiplications is about $12N_{s}N_{r}+28N_{s}$ when all the transmit paths enter the CV. In contrast, lower bounds on the total number of NLF evaluations and multiplications in MCM can be given as $10\mathcal {N}$ and $40 \mathcal {N}$ respectively, where $\mathcal {N}$ is the number of Monte Carlo photon generations. Additionally, MCM also requires operations for random number generations. If $\mathcal {N}=10^{7}$, the NLF calculations and multiplications are at least $10^{8}$ and $4\times 10^{8}$ respectively. For PSM, using $N_{s}=10$ and $N_{r}=10$ for single scattering, the NLF evaluations and multiplications would be about $470$ and $1480$ respectively. Thus, for single scattering, PSM provides computational advantage by many orders of magnitude.

Table 1. PSM operations. A row in a category assumes completion of the previous task.

View Table

Let us next consider the operations needed for double scattering. Solving Eq. (15) also includes finding all the relevant quantities, e.g., $\textbf {v}_{k}$, $\rho m_{b}$, $\gamma$ and others, associated with Eq. (15). Similarly, $\textbf {d}_{2}$ computation also requires finding $\bar {b}_{i,n,k,l}$, and $P[E_{R}|S_{l},A_{k},T_{n},R_{i}]$ requires $\textbf {u}_{R}^{T}\textbf {d}_{2}$, $\|\textbf {d}_{2}\|$ and all the quantities associated with it. Some of these quantities can be simplified and expressed in terms that have already been calculated in relation to Eq. (15). For example, using $\textbf {d}_{2}=\textbf {q}+\bar {b}_{i,n,k,l}\textbf {v}_{k}$, we can get $\textbf {u}_{R}^{T}\textbf {d}_{2}$ and $\|\textbf {d}_{2}\|$ in terms of $\rho$, $\gamma$ and $m_{b}$. The total number of NLF operations and multiplications are found to be about $4N_{A}+3N_{A}/N_{r}+N_{F}$ and $13N_{A}+18N_{A}/N_{r}+N_{M}$ respectively, where $N_{F}$ $=2N_{a}+$ $2N_{p}+$ $N_{s}$ and $N_{M}=14N_{s}N_{p}N_{a}$ $+13N_{s}N_{t}+$ $3N_{s}N_{a}+$ $2N_{s}+$ $2N_{t}$. It is important to note that most computations are caused due to second and third rows for double scattering in the table. However, these are triggered only when the hypothesized direction $\textbf {v}_{k}$ meets the receiver cone, i.e., when Eq. (15) has real valid solutions. Since most $\textbf {v}_{k}$ vectors do not meet the receiver cone in general, these row operations are bypassed. Considering $N_{s}=N_{a}=N_{p}=N_{r}=10$ and $N_{t}=50$, we get $N_{A}=5\times 10^{5}$, thus producing about $2.2\times 10^{6}$ NLF operations ignoring the bypass savings. On the other hand, MCM requires simulations of several millions of photons (e.g., at least 10 million) to get reliable results for double scattering, thereby involving NLF operations in the order of $10^{8}$. Although this shows an operational gain of about 45, the actual computer time found in our numerical results section shows significantly larger gain due to reasons such as bypass savings.

4. Atmospheric scattering parameters

The particles much smaller than the wavelength, $\lambda$, cause molecular scattering and can be modeled using the Rayleigh phase function,

The function $p_{mie}(\theta , \phi )$ can be approximated using the Henyey-Greenstein (HG) function

5. Numerical results

All figures use $\beta _{T}=17^{o}$, $\beta _{R}=30^{o}$ and $S_{R}=1.77\times 10^{-4} \hspace {1mm} \mbox {m}^{2}$. Figures 4, 5 and 6 use the HG model for Mie scattering and we use the same parameters as in [2] corresponding to $\lambda =260$ nm. Thus, $k_{s,ray}=0.266\times 10^{-3}$ $\mbox {m}^{-1}$, $k_{s,mie}=0.284\times 10^{-3}$ $\mbox {m}^{-1}$, $k_{a}=0.802\times 10^{-3}$ $\mbox {m}^{-1}$, $\gamma _{R}=0.017$, $f_{HG}=0.5$ and $g_{HG}=0.72$. The inclination angles for these figures are set at $\theta _{T}=70^{o}$ and $\theta _{R}=60^{o}$. All MCM results are obtained by simulating 10 million photons. In Figs. 7–9, we use the Mie theory model described in Appendix III for $\lambda = 250$ nm. We study two types of aerosols: fog and dust. From [18], the refractive index for dust at 250 nm is $\mu =1.53+j0.03$ and for fog [19], it is $\mu =1.362$. We use $k_{a,ray}=1.0926\times 10^{-3}$ $\mbox {m}^{-1}$ and $k_{s,ray}=3.2117\times 10^{-4}$ $\mbox {m}^{-1}$ for standard atmosphere [19]. We have chosen $\theta _{T}=70^{o}$, $\phi _{T}= - 80^{o}$, $\theta _{R}=60^{o}$ and $\phi _{R}=90^{o}$ for these figures. In the following, we refer to $N_{s}$, $N_{t}$, $N_{a}$, $N_{p}$ and $N_{r}$ as the PSM model parameters.

 

Fig. 4. Impact of PSM parameter variation when one parameter is varied keeping all other parameters fixed at 10.

Download Full Size | PDF

Figure 4 shows the root mean squared error (RMSE) between MCM and PSM when one of the PSM parameters is increased while keeping all the other parameters fixed at 10. We use $\phi _{T}=-90^{o}$. The RMSE is the positive square root of $(1/\nu )\sum _{i=1}^{\nu } (L_{MCM}^{(i)}-L_{PSM}^{(i)})^{2}$, where $L_{MCM}^{(i)}$ and $L_{PSM}^{(i)}$ are the path losses in dB obtained by MCM and PSM respectively for the $i$-th transmitter receiver configuration, and $\nu$ is the total number of configurations. We use $\nu =9$ configurations corresponding to the $(\phi _{R}, r)$ pair values of $(60^{o}, 20~\mbox {m})$, $(60^{o}, 90~\mbox {m})$, $(60^{o}, 160~\mbox {m})$, $(90^{o}, 20~\mbox {m})$, $(90^{o}, 90~\mbox {m})$, $(90^{o}, 160~\mbox {m})$, $(-90^{o}, 20~\mbox {m})$, $(-90^{o}, 90~\mbox {m})$, and $(-90^{o}, 160~\mbox {m})$. For single scattering, the gap between MCM and PSM is less than 1 dB when parameters are set at 10 or above. For double scattering, a large $N_{t}$ is needed, while other parameters can still be kept low. Since the transmitter segments have to cover a large (theoretically infinite) path in the transmit beam, a smaller $N_{t}$ adversely affects performance. Note that even for MCM, double scattering results are found to fluctuate within a couple of dB. Therefore, PSM results are quite compatible with MCM results.

 

Fig. 5. Path loss versus receiver azimuth angle at $r=50$ m.

Download Full Size | PDF

Figure 5 shows the path loss variation when the receiver azimuth angle ($\phi _{R}$) is changed keeping $\phi _{T}=-30^{o}$ at $r=50$ m. In PSM, $N_{t}=50$ with all other parameters set at 10. MCM results show scattering up to the fourth order. The single and the double scattering results of PSM and MCM match throughout the whole range. Finally, the overall scattering results from PSM include only single and double scattering, whereas the MCM overall results include scattering up to the fourth order. Nevertheless, PSM and MCM overall results agree too, as the impact of higher order scattering is negligible. With a 2.2 GHz processor and 8 GB RAM, the time taken to produce both single and double scattering path losses together for a typical receiver azimuth angle in the figure by MCM is about 212 seconds while it takes about 1.3 seconds for PSM. This gain in computing time by PSM can be raised significantly more by lowering the values of the PSM parameters at the cost of some decrease in performance.

 

Fig. 6. Channel impulse response at (a) $r=20$ m and (b) $r=60$ m.

Download Full Size | PDF

Figure 6 shows the impulse response at distances of 20 m and 60 m for $\phi _{T}=-30^{o}$ and $\phi _{R}=30^{o}$. Fixed photon paths of PSM from the transmitter to the receiver are used instead of random paths of MCM to obtain the time delays. The time delay axis is divided into multiple bins and the added probabilities in each bin are divided by the bin width and $S_{R}$. We use $N_{s}=N_{t}=50$, and also $N_{r}=50$ for single scattering while $N_{r}=10$ is used for double scattering with the other PSM parameters fixed at 10. The PSM and MCM impulse responses match well. At longer distances, the peak value of PSM is found to underestimate the MCM peak.

Figure 7 displays path loss versus particle density ($M$) for $r=20$ m, $100$ m and $180$ m with particle radius $a=0.5$ $\mu$m [20]. Although haze and fog are studied in [19], our focus is on the impact of PSM parameters in a non-planar setting and development of further insights. The path loss is found to decrease with particle density in the given range, and higher fog density gives about 5 dB benefit as more particles help in getting the signal scattered toward the receiver. A significant increase in particle density will eventually raise the path loss significantly [19]. Interestingly, higher dust density is not helpful as much as fog density. We use three sets of PSM model parameters. In the low model, $N_{t}=50$ with other parameters fixed at 10. In the medium model, $N_{t}$ is same as in the low model while all the other parameters are raised to 12. In the high model, $N_{t}=60$ with the other parameters set at 18. We note that the path loss gap between the low model and the high model is 0.5 dB or less and hence the low PSM represents path loss results quite accurately, and a higher complex PSM model is not necessary. Similar results are also observed in Fig. 8 on path loss versus particle size keeping $M=10^{8}$ $\mbox {m}^{-3}$. For larger particle sizes, higher PSM parameters will be needed for more accurate results. However, in the ranges shown, the low PSM is close to the high PSM. Finally, we also studied the root mean squared (rms) delay spread [5]. In the given set up, the rms delay spread is found to decrease with particle density or size. The rms delay spread increases with $r$, but the proportional change in the delay spread is not affected by $r$. For example, in the case of fog at $r=20$ m, the rms delay spread decreases from 1.77 ns to 1.45 ns when $M$ increases from $10^{7}$ $\mbox {m}^{-3}$ to $10^{9}$ $\mbox {m}^{-3}$. This is an $18\%$ decrease. Even for $r=180$ m, a similar decrease is observed for fog as it decreases from 15.9 ns to 12.9 ns. For dust, the delay spread is found to be more robust to particle density changes.

 

Fig. 7. Effect of particle density on path loss for different PSM parameters.

Download Full Size | PDF

 

Fig. 8. Effect of particle size on path loss for different PSM parameters.

Download Full Size | PDF

 

Fig. 9. Path loss versus distance between the transmitter and the receiver.

Download Full Size | PDF

Finally, Fig. 9 studies path loss for changing distances between the transmitter and the receiver. We distinguish between near-planar and off-planar set ups. The particle size is $a=0.5$ $\mu$m. The low and high densities refer to $10^{7}$ $\mbox {m}^{-3}$ and $10^{9}$ $\mbox {m}^{-3}$ respectively. We use $\theta _{T}=70^{o}$ and $\theta _{R}=60^{o}$ with two azimuth pairs for ($\phi _{T}, \phi _{R})$: 1) near-planar $(-80^{o}, 90^{o})$ and 2) off-planar $(-40^{o}, 50^{o})$. Dust and fog give similar path losses for low densities. At high densities, however, fog gives lower path loss and the gap between fog and dust is larger for near-planar set up for the given inclination angles. As the paths get closer to the LOS, the impact of particle scattering becomes more important due to the nature of the Mie phase scattering function.

6. Conclusions

This paper proposes and analyzes a non-Monte Carlo method for calculating single and double scattering event characteristics in NLOS UV communications. The novel idea of the proposed method is the deterministic probability-based sampling of photon emissions and scattering events. The method improves upon the complexity of Monte-Carlo simulations by several orders of magnitude. Many of the computations can be performed off-line and stored. For example, when the transmitter is fixed and we want to study the effect of various receiver configurations, orientations or distance variations, many of the stored computations can be re-used. A complexity comparison for the method is presented. The impact of the model parameters on the accuracy of the method for Rayleigh and Mie scattering is studied. Accurate results within a dB of path loss are obtained using low PSM parameter values that can be raised to obtain better accuracy, if needed. Future work can include computationally simplified modeling of longer links with turbulence effects and their overall impact on system performance such as the bit error rates.

Appendix I

Sampling points (photon paths) on the transmit beam

An algorithm is proposed to distribute the $N_{s}$ sample points uniformly over the available angular space. We have to find the values of $N_{i}$ and $\theta _{i}$ for $i=1, 2, \ldots , N_{c}$. The steps are as follows.

• 1. Find the angle $\theta _{1}'$ for the first grid circle using $\kappa (1-\cos \theta _{1}')=1/N_{s}$.
• 2. Obtain an initial estimate on the number of circles $N_{c}$ using $(2N_{c}+1)\theta _{1}'=\beta _{T}/2$, so that the available angle space is equally divided. Use the ceiling function to get an integer $N_{c}$.
• 3. Obtain an initial estimate $\theta _{i}$ for each circle based on $\theta _{i}=2i \times \theta _{1}'$, $i=1,2, \ldots , N_{c}$.
• 4. To find the number of samples $N_{i}$ on the $i$-th sampling circle, we want to spread them equally over each circle. If we imagine a sphere of radius $\mathcal {R}$ with center at T, the circumference of the $i$-th sampling circle with the center at $A_{0}$ is $2\pi \mathcal {R} \sin \theta _{i}$. Therefore, to spread the sample points equally on each sampling circle, we require

  $$\frac{N_{1}}{\sin \theta_{1}}=\frac{N_{2}}{\sin \theta_{2}}=\cdots=\frac{N_{N_{c}}}{\sin \theta_{N_{c}}}.$$

  This will give estimates

  $$N_{i}=(N_{s}-1)\frac{\sin \theta_{i}}{\sum_{k=1}^{N_{c}} \sin \theta_{k}}.$$

  The actual number is obtained by rounding of $N_{i}$. If any $N_{i}$ is zero, then we reduce the number of circles $N_{c}$ by 1, and repeat Steps 3 and 4. If due to rounding, $\sum _{i=1}^{N_{c}}N_{i}$ becomes less than $N_{s}-1$, place the extra points $N_{s}-\sum _{i=1}^{N_{c}}N_{i}-1$ in the outermost sampling circle.

• 5. Using Eqs. (2) and (3), obtain the true $\theta _{i}$ based on the $N_{i}$ estimates obtained from Step 4.
• 6. Iterate between Steps 4 and 5 until convergence.

In our studies, convergence has been obtained in 3 or 4 iteration steps. The computational load of the algorithm is negligible in the overall context of PSM implementation.

Appendix II

PSM implementation summary

Figure 10 summarizes the operations. In Block 1, the transmitter and receiver parameters are $\theta _{T}$, $\theta _{R}$, $\phi _{T}$, $\phi _{R}$, $\beta _{T}$, $\beta _{R}$, $r$, and $S_{R}$. The atmospheric parameters are given in Section 4 and the PSM parameters are given in the first paragraph in Section 5. For single scattering, solving Eq. (8) in Block 3A requires computing $\textbf {u}_{i}$ using Eqs. (4)–(6) and the other terms needed in Eq. (8), e.g., $m_{i}$. The solution adjustment there refers to the discussion immediately after Eq. (8). Block 4A refers to the calculation of $\bar {s}_{i,k}$ in Eq. (11) performed on the paths that pass through the CV. Finally, Block 5A requires $P(i)$ calculation and $P[E_{R}|E_{k},R_{i}]$ using Eq. (12). In the case of double scattering, Block 3B involves Eq. (14). The procedure to obtain $\textbf {v}_{k}$ in Block 4B is given immediately before Eq. (15) and the solution adjustment is described immediately after Eq. (15). The receiver ESPs specified in Block 5B require Eq. (16). The computational procedure to obtain $P[E_{R}|S_{l},A_{k},T_{n},R_{i}]$ in Eq. (18) is described immediately after the equation. The path losses in dB for single and double scattering are given by $- 10\log _{10} P[D_{i}]$, $i=1,2$, and the overall path loss is $- 10\log _{10} (P[D_{1}]+P[D_{2}])$. Finally, the time delays required for obtaining the channel impulse response can be obtained for all the valid photon paths using $(\bar {s}_{i,k}+\|\textbf {d}_{1}\|)/c$ and $(\bar {s}_{i,n}+ \bar {b}_{i,n,k,l}+\|\textbf {d}_{2}\|)/c$ for the single and double scattering events respectively, where $c$ is the speed of light in air.

 

Fig. 10. Block diagram of the PSM operations for path loss calculation.

Download Full Size | PDF

Appendix III

Calculation steps in Mie theory

We assume the refractive index of the atmosphere to be one. Suppose $\mu$ is the particle refractive index. We calculate the coefficients $a_{n}$ and $b_{n}$ for the scattered field as [11]

(22)$$a_{n}=\frac{\tilde{\mu} \mu^{2} j_{n}(\mu x)[xj_{n}(x)]'-\tilde{\mu}_{1}j_{n}(x)[\mu xj_{n}(\mu x)]'}{\tilde{\mu}\mu^{2} j_{n}(\mu x)[xh_{n}^{(1)}(x)]' - \tilde{\mu}_{1}h_{n}^{(1)}(x)[\mu x j_{n}(\mu x)]'},$$

where $\tilde {\mu }$ is the permeability of the atmosphere, $\tilde {\mu }_{1}$ is the permeability of the particle, $j_{n}(x)=\sqrt {\pi /2x}J_{n+0.5}(x)$, $J_{n}(x)$ is the Bessel function of first kind, $h_{n}^{(1)}(x)=j_{n}(x)+jy_{n}(x)$, $y_{n}(x)=\sqrt {\pi /2x} Y_{n+0.5}(x)$, $Y_{n}(x)$ is the Bessel function of second kind, and a prime mark denotes the derivative with respect to the argument. Similarly,

(23)$$b_{n}=\frac{\tilde{\mu}_{1} j_{n}(\mu x)[xj_{n}(x)]'-\tilde{\mu} j_{n}(x)[\mu xj_{n}(\mu x)]'}{\tilde{\mu}_{1} j_{n}(\mu x)[xh_{n}^{(1)}(x)]' - \tilde{\mu} h_{n}^{(1)}(x)[\mu x j_{n}(\mu x)]'}.$$

Using $\omega =\cos \theta$ for the scattering angle $\theta$, $S_{1}(\omega )$ and $S_{2}(\omega )$ are obtained as

(24)$$S_{1}(\omega)=\sum_{n=1}^{n_{max}} \frac{2n+1}{n(n+1)}(a_{n}\pi_{n}(\omega)+b_{n}\tau_{n}(\omega)),$$

where $n_{max}$ is the nearest integer to $(2+x+4x^{1/3})$, and the functions $\pi _{n}(\omega )$ and $\tau _{n}(\omega )$ can be obtained recursively as

$$ \pi_{n}(\omega)=\frac{2n-1}{n-1}\omega \pi_{n-1}(\omega)-\frac{n}{n-1}\pi_{n-2}(\omega) \hspace{2mm} \mbox{and} \hspace{2mm} \tau_{n}(\omega)=n\omega \pi_{n}(\omega)-(n+1)\pi_{n-1}(\omega) $$

with the initial conditions: $\pi _{0}(\omega )=0$, $\pi _{1}(\omega )=1$, $\tau _{0}(\omega )=0$, and $\tau _{1}(\omega )=\omega$. Similarly,

(25)$$S_{2}(\omega)=\sum_{n=1}^{n_{max}} \frac{2n+1}{n(n+1)}(a_{n}\tau_{n}(\omega)+b_{n}\pi_{n}(\omega)).$$

Finally, the phase function is obtained by

(26)$$p_{mie}(\theta,\phi)=\frac{|S_{1}(\omega)|^{2}+|S_{2}(\omega)|^{2}}{4\pi \sum_{n=1}^{n_{max}}(2n+1)(|a_{n}|^{2}+|b_{n}|^{2})}.$$

Funding

U.S. Army Combat Capabilities Development Command Data & Analysis Center at White Sands Missile Range, New Mexico.

Acknowledgment

The authors would like to thank Alejandro Gongora for helpful discussions on system aspects and Prof. Charles Bruce for scattering theory clarifications. This work was supported by the U.S. Army CCDC Data & Analysis Center at White Sands Missile Range, New Mexico, through a grant with the NMSU Physical Science Laboratory. The first author also acknowledges support from the International Foundation for Telemetering (IFT).

Disclosures

The authors declare no conflicts of interest.

References

1. H. Ding, G. Chen, A. K. Majumdar, B. M. Sadler, and Z. Xu, “Modeling of non-line-of-sight ultraviolet scattering channels for communication,” IEEE J. Select. Areas Commun. 27(9), 1535–1544 (2009). [CrossRef]  

2. R. J. Drost, T. J. Moore, and B. M. Sadler, “UV communications channel modeling incorporating multiple scattering interactions,” J. Opt. Soc. Am. A 28(4), 686–695 (2011). [CrossRef]  

3. M. Noshad, M. Brandt-Pearce, and S. G. Wilson, “NLOS UV communications using M-ary spectral-amplitude-coding,” IEEE Trans. Commun. 61(4), 1544–1553 (2013). [CrossRef]  

4. M. J. Weisman, F. T. Dagefu, T. J. Moore, C. H. Arslan, and R. J. Drost, “Analysis of the low-probability-of-detection characteristics of ultraviolet communications,” Opt. Express 28(16), 23640–23651 (2020). [CrossRef]  

5. M. A. Elshimy and S. Hranilovic, “Non-line-of-sight single-scatter propagation model for noncoplanar geometries,” J. Opt. Soc. Am. A 28(3), 420–428 (2011). [CrossRef]  

6. D. K. Borah and D. G. Voelz, “Pointing error effects on free-space optical communication links in the presence of atmospheric turbulence,” J. Lightwave Technol. 27(18), 3965–3973 (2009). [CrossRef]  

7. M. R. Luettgen, J. H. Shapiro, and D. M. Reilly, “Non-line-of-sight single-scatter propagation model,” J. Opt. Soc. Am. A 8(12), 1964–1972 (1991). [CrossRef]  

8. Z. Xu, H. Ding, B. M. Sadler, and G. Chen, “Analytical performance study of solar blind non-line-of-sight ultraviolet short-range communication links,” Opt. Lett. 33(16), 1860–1862 (2008). [CrossRef]  

9. L. Wang, Z. Xu, and B. M. Sadler, “Non-line-of-sight ultraviolet link loss in noncoplanar geometry,” Opt. Lett. 35(8), 1263–1265 (2010). [CrossRef]  

10. R. J. Drost, T. J. Moore, and B. M. Sadler, “Ultraviolet scattering propagation modeling: analysis of path loss versus range,” J. Opt. Soc. Am. A 30(11), 2259–2265 (2013). [CrossRef]  

11. G. F. Bohren and D. R. Huffman, Absorption and Scattering by a Sphere, Absorption and Scattering of Light by Small Particles (Wiley, 1983).

12. R. Yuan, J. Ma, P. Su, and Z. He, “An integral model of two-order and three-order scattering for non-line-of-sight ultraviolet communication in a narrow beam case,” IEEE Commun. Lett. 20(12), 2366–2369 (2016). [CrossRef]  

13. R. J. Drost, M. J. Weisman, F. T. Dagefu, and C. H. Arslan, “A low-complexity approach to UV communication channel model evaluation,” in 2019 International Conference on Military Communications and Information Systems (ICMCIS), (2019), pp. 1–7.

14. F. Li, Y. Zuo, A. Li, Z. Du, and J. Wu, “Spatially correlated MIMO for exploiting the capacity of NLOS ultraviolet turbulent channels,” Opt. Express 27(21), 30639–30652 (2019). [CrossRef]  

15. K. Kumar and D. K. Borah, “Hybrid FSO/RF symbol mappings: Merging high speed FSO with low speed RF through BICM-ID,” in 2012 IEEE Global Communications Conference (GLOBECOM), (2012), pp. 2941–2946.

16. R. Yuan, J. Ma, P. Su, Y. Dong, and J. Cheng, “Monte-Carlo integration models for multiple scattering based optical wireless communication,” IEEE Trans. Commun. 68(1), 334–348 (2020). [CrossRef]  

17. D. E. Peplow, “Direction cosines and ploarization vectors for Monte Carlo photon scattering,” Nucl. Sci. Eng. 131(1), 132–136 (1999). [CrossRef]  

18. E. P. Shettle and R. W. Fenn, “Models for the aerosols of the lower atmosphere and the effects of humidity variations on their optical properties,” Tech. rep., Air Force Geophysics Lab, Massachusetts (1979).

19. C. Xu, H. Zhang, and J. Cheng, “Effects of haze particles and fog droplets on NLOS ultraviolet communication channels,” Opt. Express 23(18), 23259–23269 (2015). [CrossRef]  

20. J. N. Sanders-Reed and S. J. Fenley, “Visibility in degraded visual environments (DVE),” Proc. SPIE 10642, 106420S (2018). [CrossRef]