Printer Friendly

An Insight into Creeping Electromagnetic Waves around the Human Body.

1. Introduction

As the inherent presence of the human body in body-centric communication systems [1] affects radiation properties of antennas and radio channel, there is considerable interest in characterization of electromagnetic wave propagation in vicinity of the human body with the goal to predict and enhance link budget and optimize the position of antennas on the body. Wireless links where the body surface takes the role in transmission channel are generally considered to allow for deterministic modelling to some extent, which gives rise to research in Norton and creeping waves (see, e.g., [27]). To provide an adequate canonical electromagnetic model, human body is typically modelled as cylindrical, spherical, or elliptical [6] structure representing either the torso or head, where link budget around the body is studied in various circumstances both numerically and experimentally.

Instead of using full geometrical optics calculations [4], the purpose of this paper is to gain an intuitive insight into the scattered field around the human body torso in terms of creeping waves. Thus we develop a simple model of propagation for on-body communication link by focusing on the most dominant effect which predicts general behaviour and attenuation of the electric field around the human body.

The problem is formulated as a canonical boundary-value problem in the spectral domain where the poles of Green's function for the two fundamental (T[M.sub.z] and T[E.sub.z]) field configurations are calculated and shown to correspond to creeping wave propagation coefficients around the human body giving rise to the general trend of field attenuation around the body. In the next part of the paper we exploit the example of practical application of the model by analyzing wearable PIFA antenna from [7] with respect to the developed theoretical model and discuss the limitations of the model by performing measurements on realistic on-body link.

2. Theoretical Model

To build a theoretical model we approximate human body as an infinite lossy dielectric cylinder with constitutive parameters corresponding to a muscle tissue [8], as shown in Figure 1 (i.e., we create an analytical muscle-equivalent phantom). We assume monochromatic field (with [e.sup.-i[omega]t] time dependence) and proceed with the analysis in circular cylindrical coordinate system to formulate Green's function of the problem and hence calculate the field distribution around the body in the plane of source. To simplify the treatment of general vector problem we first note that any vector field in source-free region can be derived using two scalar functions (so-called Debye potentials), using the prescription described, for example, in [9-11]. These functions may be regarded as the generating functions of two partial fields, while the superposition of these fields gives rise to total field. To find these functions one needs to solve inhomogeneous scalar equations (subject to proper boundary conditions [9]) of the type

[mathematical expression not reproducible], (1)

where [psi] = [psi] ([rho], [phi], z) represents the sought partial field function (i.e., the Debye potential). In cylindrical coordinates for the source terms for these two functions we choose infinitesimal magnetic and electric dipoles of amplitude A pointing into z-direction, respectively. These sources radiate in the presence of cylinder of radius a shown in Figure 1(a) representing the human body, which is the core of the proposed model. The source radial coordinate is [rho] = b while the other two coordinates are taken as zero. To further specialize (1) onto cylindrical problems one typically assumes that the axial (z-) coordinate is of infinite extent so the unknown field dependency in axial direction may be eliminated by applying Fourier transformation to wave number space [k.sub.z] (this is called guided wave representation [11]). In that case the Debye potentials are directly proportional to axial (z-) components of the electric ([E.sub.z]) and magnetic field ([H.sub.z]) [9], respectively, so instead of potentials one can resort to solving the scalar Helmholtz differential equations for the two axial field components which is arguably more intuitive approach in the given geometry. The total vector field is then obtained as a linear combination of partial fields arising from [E.sub.z] and [H.sub.z] field functions, which are termed as T[M.sub.z] and T[E.sub.z] modes, respectively [9,11]. In addition, since the cross section of the observed dielectric object is circular and hence periodic in angular ([phi]-) direction, we may apply Fourier series in angular direction and thereby further reduce the problem to one-dimensional in radial ([rho]-) direction. In that way the sought field function can be expressed as

[mathematical expression not reproducible] (2)

where [k.sub.z] is the wave number in axial direction and m is the angular wave number, respectively. The remaining radial dependency follows by inserting spectral domain field function [mathematical expression not reproducible] into (1) which is then reduced to Bessel's differential equation describing radial transmission line [9, 10]. The general solution in spectral domain then follows as

[mathematical expression not reproducible] (3)

where [J.sub.m] and [H.sup.(l).sub.m] denote Bessel functions and Hankel functions of first kind of order m, respectively, while [k.sub.0] is the free-space wave number. For each m and [k.sub.z] one basically only needs to find coefficients [C.sub.1] and [C.sub.2], which are readily obtained by applying the boundary conditions [9] at the interfaces of respective radial regions (Figure 1(b)). That way we obtain the sought field components (which give rise to T[M.sub.z] and T[E.sub.z] modes) in region around the body (denoted in Figure 1(b) as a < [rho] < b) for some (m, [k.sub.z]) as

[mathematical expression not reproducible] (4)

where [mathematical expression not reproducible] are spectral current density amplitudes of electric and magnetic sources, respectively; [k.sub.[rho]] = [square root of ([k.sup.2.sub.0] - [k.sup.2.sub.z])] is the radial wave number; b is the source radial distance.

Note that in (4) the term containing Bessel function can be interpreted as incident (primary) wave while the term containing Hankel function is interpreted as reflected (scattered) wave, which is the property of radial transmission lines [11]. The pertinent reflection coefficients [R.sub.TM] and [R.sub.TE] are calculated as

[mathematical expression not reproducible] (5)

where [k.sub.1] is the wave number of the dielectric of permittivity [[epsilon].sub.r1] representing human muscle tissue (nonmagnetic material is assumed). Furthermore, [J'.sub.m] and [H'.sup.(l).sub.m] are derivatives of Bessel and Hankel functions taken in radial direction, respectively.

Although the incident (primary) field configuration arising from the z-directed elementary source is either pure T[M.sub.z] or T[E.sub.z] depending on the type of excitation (electric or magnetic, resp.), the scattered field from some cylindrical obstacle (such as human body) in general needs to contain both types of fields to satisfy the boundary conditions at the cylinder surface [11]. The resulting total field hence suffers from mode coupling that blurs the physical insight.

Nevertheless, for practical applications we consider it interesting to study the behaviour of each field mode around the body separately since the total field is still linear superposition of the two principal modes, while the actual amplitude of each mode would depend on some particular excitation and the character of the obstacle. Thus to simplify the computation process and obtain an insight into dominant effects of the field behaviour for the two types of sources we implement a two-step approximation as explained below. Note that in the rest of the paper with the incident field we will denote the field excited by the z-directed elementary source without presence of the dielectric cylinder.

2.1. The Two-Step Approximation. As a first step in the approximation we focus on the case of [k.sub.z] = 0 which represents the dominant spectral contribution if the source and observation points lie in the xy plane and if the distance between them is large enough. Apart from simplifying the calculations, the assumption [k.sub.z] = 0 in addition enables uncoupling of T[M.sub.z] and T[E.sub.z] modes in scattered (and total) field and observing the behaviour of each mode around the human body separately, which is the primary goal of the analysis. Expression (2) is therefore reduced to

[mathematical expression not reproducible]. (6)

Since series (6) is slowly convergent we develop an alternative fast converging series by applying Watson's transformation technique [10] which is based on interpreting summation (6) as a residue summation of the integral of the type

[mathematical expression not reproducible] (7)

where v is continuous complex variable.

For such integral representation we are allowed to change the contour of integration so that instead of summing the residues at poles of the sine function (v = m) we perform summation of residues at poles of spectral domain functions [mathematical expression not reproducible] given in (4) for the two considered cases (T[M.sub.z] and T[E.sub.z]). The new poles arise from the nulls of denominator of reflection coefficients given in (5) and occur for some complex numbers [v.sub.k] which need to be determined. For each observed field configuration this eventually leads to the spatial field functional dependence given as [10]:

[mathematical expression not reproducible], (8)

where [A.sub.k] is the value of residue at the kth pole. Note that in this expression the complex poles [v.sub.k] can be interpreted as the angular propagation coefficients for the two waves traveling in the opposite directions around the cylinder, which are termed creeping waves [10]. The imaginary parts of [v.sub.k] represent the losses arising from radiation into space, so the field of each creeping wave is exponentially attenuated around the body. Compared to representation (6), the summation in terms of creeping waves (8) is fast converging due to presence of the imaginary parts of [v.sub.k].

In the second step of the approximation we observe the contribution of the poles with the smallest imaginary part for both T[M.sub.z] and T[E.sub.z] case. If the concerned cylinder is electrically large enough, the azimuthal field distribution around the cylinder can be adequately approximated with only one pole [10]. Apart from simplifying the calculations, this way we can express the most contributing portion of the field around the human body in terms of only two waves propagating from the source in two opposite angular directions, portion of which is radiated into free space when traveling around the curvature. This means that the used approach provides the clear and intuitive insight into dominant effect of the interaction between human body and field source (antenna) and enables one to predict the properties of the communication channel around the human body. Note however that such representation is most accurate for shadow region in the vicinity of cylinder, that is, where no direct field from the source is present [10,11] and contribution of higher order poles is attenuated.

2.2. Application to Human Body Model. To illustrate the accuracy of the developed approach we calculate the poles [v.sub.k] for the human body model we developed for our research. We assume that the cylinder radius a is 14 cm, which closely corresponds with torso circumference of human volunteer who participated in our research (it can be also regarded to correspond to typical adult male torso). The frequency of interest is 2.45 GHz (ISM band), while the relative permittivity and conductivity of the analyzed cylinder are taken as 53 and 1.7 S/m, respectively, which corresponds to the muscle tissue at the given frequency [8].

By applying Olver's uniform representation for Bessel and Hankel functions of complex order [12], the sought poles [v.sub.k] for the two analyzed cases are found by plotting the denominators of (5) as functions of complex variable and searching for their zeros. The calculated loci of poles in complex v-plane are illustrated in Figure 2, while the exact values of the poles with smallest imaginary part (i.e., the propagation coefficients of the dominant creeping waves) are summarized in Table 1. Note also that there are two sets of poles in Figure 2 which arise from the zeros of Hankel (slanted line) and Bessel functions (approximately horizontal line), respectively. The latter poles are associated with the field inside the body [13]; however they are of less importance in this analysis due to large attenuation inside the body.

In Figure 3 we show normalized magnitudes of the electric field around the used human body model for the T[M.sub.z] and T[E.sub.z] cases. The solid line shows the angular field distribution of dominant creeping waves calculated using (8) and the results from Table 1. The source is assumed to be placed at [phi] = 0[degrees] while it can also be seen that the interference of the two opposite creeping waves occurs around the region opposite to the source (i.e., around [phi] = 180[degrees]). To demonstrate the relevance of the used approximations from previous section we have also added the full-wave solution calculated directly from expression (2), which is shown with diamonds.

By comparison of relevant cases shown in Figure 3 it can be seen that the developed simplified model (with [k.sub.z] = 0 and one-pole approximations) adequately represents the general trend of field decay for both the T[M.sub.z] and T[E.sub.z] modes. The minor discrepancies which nevertheless occur in the lit regions were actually expected due to already mentioned limitations of the model [10, 14]. To further illustrate the influence of the used two-step approximation, in Figure 4 we show an intermediate step where we compare the field distribution with assumption of [k.sub.z] = 0 calculated exactly as a full Fourier summation in m (which is equivalent to considering all poles in calculation) and compare it with the used one-pole approximation. By comparing with Figure 3 it can be seen that the discrepancy in results principally arises from assumption [k.sub.z] = 0. Note that in both cases the agreement between full summation and assumption of only one pole is very good, meaning that the approximation with one pole simplifies the calculation process without introducing practically significant errors. This justifies the use of the proposed simple and intuitive model for real-life applications, for example, as a guideline in wearable antenna design and communication link budget around human body.

From Figures 3 and 4 it can also be seen that T[E.sub.z] modes exhibit considerably smaller slope of field decay around the body compared to T[M.sub.z] modes so we expect they would generally prevail in communication when present in significant amount. Thus to enhance the communication around the body one principally needs to consider the sources that excite the T[E.sub.z] modes, that is, which efficiently couple energy into the axial magnetic field ([H.sub.z]).

2.3. Small Dipole Radiation. According to [9-11], the incident radiated field is pure T[M.sub.z] or T[E.sub.z] if the sources are axially (z-) oriented while the type of field mode depends only on whether such source is of electric or magnetic type, respectively. On the other hand, the angular ([phi]-) and radial ([rho]-) sources, both electric and magnetic, in general couple the energy into both types of modes [9, 10]. Note that in reality true magnetic sources which would hence be the best choice for communication around the body do not exist, so one needs to reside to equivalent sources such as in [3] where normally oriented dipole with respect to the human body was found to exhibit smaller attenuation rate around human body torso compared to the axially oriented one. Using the proposed model based on T[M.sub.z] or T[E.sub.z] field decomposition this finding [3] gets firmer theoretical ground.

To clarify this matter we perform simulations in CST Microwave Studio to calculate the electric field around the body generated by small electric dipole (8 mm length, 0.5 mm thickness) at the distance of 1 mm from the body, which represents an elementary electric source. We consider the cases when the dipole is axially (z-), angularly ([phi]-), and radially ([rho]-) oriented with respect to the body (Figure 5). The body is modelled in the same way as in previous calculations, while the cylinder height is 65 cm which can be considered electrically large enough to conform to the proposed model. The applied voltage on all the antenna terminals is 1 volt. The calculated electric field around the body (angular distance from the source is given in degrees) is given in Figure 6.

The axially oriented electric dipole excites the T[M.sub.z] mode only and thus exhibits the larger field decay rate around the body compared to two other dipole orientations, as predicted theoretically. By comparison of Figure 6 with Figure 3 it can be seen that the slope around the body indeed corresponds to the one predicted for T[M.sub.z] modes. On the other hand, the angular and radial dipoles both exhibit the T[E.sub.z]-like behaviour of field decay around the body. This occurs since angular and radial dipoles excite both types of modes [10]; however the T[E.sub.z] component predominates in the total field due to smaller attenuation (see Figure 3). In Figure 6 it can also be seen that the total field amplitude of radial dipole is around 10 dB higher than that of the angular one, merely due to its smaller coupling to the lossy dielectric cylinder (i.e., human body). The illustrations of total and instantaneous fields are given in Figures 7 and 8 for the cases of axial and radial dipole, respectively, by which one can visualize the difference between T[M.sub.z] and T[E.sub.z] field behaviour around the body. It can be seen that when dipole is radially oriented, larger primary field is radiated into the air giving also rise to larger field amplitude around the body. In addition, by observing field magnitudes and snapshots (Figures 7 and 8, resp.), one notices that there is virtually no direct penetration of field through the body. This is expected since the skin depth of the muscle tissue at ISM frequency range is only about 2.3cm [8].

3. Practical Verification

To verify the proposed model in practice we perform an analysis of the on-body PIFA antenna which has already been used for investigation of various aspects of on-body communications [7] as well as for research of conductive textiles [15]. This type of antenna was found suitable for applications in wearable systems since it is low profile (4 mm height), integrable in clothes, straightforward to manufacture, and with good radiation characteristics in the presence of the human body. The prototype is shown in Figure 9, together with the scheme and dimensions of its patch in mm (more details can be found in [7, 15]). In addition, Figure 9(c) illustrates the dominant radiation mechanism. The radiation occurs mainly due to the presence of the electric field at the open slot opposite to feeding point. This field can be regarded as an equivalent magnetic source [M.sub.y] [9-11] which approximately acts as an angular ([phi]-directed) magnetic source when the antenna is vertically oriented to the body in a way shown in Figure 10 (i.e., parallel to z-axis).

According to theory discussed in Section 2, the angular magnetic source would contribute to both T[E.sub.z] and T[M.sub.z] modes. This is checked using CST by calculating all the field components around the cylinder (Figure 11) the same way as it was previously done for dipoles. The [E.sub.z] field component (i.e., T[M.sub.z] modes) decays fast around the body as predicted. In the total electric field the radial component predominates and by comparing the results to Figures 3 and 6, it can be seen that the T[E.sub.z] field behaviour around the body prevails, confirming the theoretical predictions. This means that using appropriate magnetic sources it is possible to have similar field decay around the body as with radial electric dipole, while at the same time keeping the antenna low profile which is suitable for wearing on body.

Finally, we perform measurements with the antennas mounted on human torso (Figure 10), by which we check the practical applicability of the proposed model. For this purpose two prototypes of PIFA antennas are used to create the on-body wireless link around torso, while the interest is focused to measurement of the antenna coupling around the body (i.e., the transmission parameter [S.sub.21], which is proportional to electric field). The measurements are performed at 21 points around the body (Figure 10(c)), repeated three times, and averaged to mitigate other influences (e.g., accidental body movement, breathing). The relevant results are given in Figure 12 together with comparison with analytic model for T[E.sub.z] case (which is normalized to the measurement starting point) and accompanying simulation in CST. The slope of the averaged measurements basically corresponds to T[E.sub.z] mode in the region of deeper shadow ([phi] > 90[degrees]) while some discrepancies occur when the antennas are visible, since the developed model is best suitable for shadow region [10,11,14], in accordance with Figure 3.

4. Conclusion

In this paper the creeping wave propagation around human torso is modelled by considering radiation of elementary sources over an infinite muscle-equivalent cylinder. The inverse Fourier transformation is calculated by two-step approximation (Watson's transformation and taking into account the most contributing spectral variable) in order to obtain a simple expression for the dominant wave propagation part. It was found that T[E.sub.z] portion of the total field (generated by axial magnetic field [H.sub.z]) prevails in communication via creeping waves when excited by the source, since it possesses smaller field decay (i.e., smaller imaginary part of dominant pole). Some efficient sources of T[E.sub.z] modes are radial electric dipole or microstrip-like radiating slot, the latter being of more interest for applications in body-centric sensor systems due to its low profile. Further analysis of the practical wearable antenna prototype (i.e., PIFA antenna) and corresponding measurements are in accordance with the theoretical predictions.

The used theoretical model can be readily further extended, for example, by taking into account more poles of Green's function or considering more layers of human body. In reality however, some outer influences (multipath, body movement, etc.) would eventually blur any deterministic model. Nevertheless, the presented results show that the proposed model provides an intuitive insight into the characteristics of propagation around the human torso and can thus serve as a practical guideline in preferable choice and design of wearable antennas and communication systems which operate in the vicinity of human body.

https://doi.org/10.1155/2017/2510196

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work has been fully supported by Croatian Science Foundation (HRZZ) under Projects IP-11-2013-6198 and IP 11-2013-3425.

References

[1] P. S. Hall and Y. Hao, Antennas and Propagation for Body-Centric Wireless Communications, Artech House, 2nd edition, 2012.

[2] C. Y. Seo, K. Saito, M. Takahashi, and K. Ito, "Quantitative estimation of scattering waves in cylinder-body model for body area networks: Comparison of analyses with unifrom cylinder-and slab-body models," Progress in Electromagnetics Research B, no. 22, pp. 145-170, 2010.

[3] T. Alves, B. Poussot, and J.-M. Laheurte, "Analytical propagation modeling of BAN channels based on the creeping-wave theory," IEEE Transactions on Antennas and Propagation, vol. 59, no. 4, pp. 1269-1274, 2011.

[4] T. Mavridis, L. Petrillo, J. Sarrazin, D. Lautru, A. Benlarbi-Delai, and P. De Doncker, "Theoretical and experimental investigation of a 60-GHz off-body propagation model," Institute of Electrical and Electronics Engineers. Transactions on Antennas and Propagation, vol. 62, no. 1, pp. 393-402, 2014.

[5] L. Petrillo, T. Mavridis, J. Sarrazin, D. Lautru, A. Benlarbi-Delai, and P. De Doncker, "Theoretical and experimental investigation of a 60-GHz off-body propagation model," Institute of Electrical and Electronics Engineers. Transactions on Antennas and Propagation, vol. 62, no. 1, pp. 393-402, 2014.

[6] R. Chandra and A. J. Johansson, "An analytical link loss model for on-body propagation around the body based on elliptical approximation of the torso with arms' influence included," IEEE Antennas and Wireless Propagation Letters, vol. 12, pp. 528-531, 2013.

[7] B. Ivsic, J. Bartolic, D. Bonefacic, and Z. Sipus, "Design and performance of miniaturized wearable antennas in UHF band," in Proceedings of the 8th European Conference on Antennas and Propagation, (EuCAP '14), pp. 1556-1560, The Hague, Netherlands, April 2014.

[8] C. Gabriel, "An internet resource for calculation of dielectric parameters of body tissues," http://niremf.ifac.cnr.it/tissprop.

[9] G. Tyras, Radiation and Propagation of Electromagnetic Waves, Academic Press, Cambridge, Mass, USA, 1969.

[10] A. Ishimaru, Electromagnetic Wave Propagation, Radiation, and Scattering, Prentice-Hall, NJ, USA, 1991.

[11] W. C. Chew, Waves and Fields in Inhomogeneous Media, vol. 522, IEEE Press, 1995.

[12] R. Paknys, "Evaluation of Hankel Functions with Complex Argument and Complex Order," IEEE Transactions on Antennas and Propagation, vol. 40, no. 5, pp. 569-578, 1992.

[13] W. Franz and P. Beckmann, "Creeping waves for objects of finite conductivity," IRE Transactions on Antennas and Propagation, vol. 4, no. 3, pp. 203-208, 1956.

[14] S. Sensiper, "Cylindrical Radio Waves," IRE Transactions on Antennas and Propagation, vol. 5, no. 1, pp. 56-70, 1957

[15] B. Ivsic, D. Bonefacic, and J. Bartolic, "Considerations on embroidered textile antennas for wearable applications," IEEE Antennas and Wireless Propagation Letters, vol. 12, pp. 1708-1711, 2013.

Branimir Ivsic, Davor Bonefacic, Zvonimir Sipus, and Juraj Bartolic

Faculty of Electrical Engineering and Computing, University of Zagreb, Zagreb, Croatia

Correspondence should be addressed to Branimir Ivsic; branimir.ivsic@fer.hr

Received 26 July 2017; Accepted 16 November 2017; Published 13 December 2017

Academic Editor: David Naranjo-Hernandez

Caption: FIGURE 1: Sketch of the analytical model for propagation around cylinder (torso). (a) Small dipole; (b) top view.

Caption: FIGURE 2: The scheme of pole loci in complex v-plane for the two considered cases (normalized to cylinder electric dimension in the air [k.sub.0]a.).

Caption: FIGURE 3: Calculated electric field around cylinder for two considered field modes using obtained angular propagation coefficient. Solid line: two-step approximation with one pole and [k.sub.z] = 0 (8); diamonds: full-wave solution (2).

Caption: FIGURE 4: Accuracy of the first step approximation ([k.sub.z] = 0). Diamonds: solution with only the first step approximation, that is, complete Fourier series expansion (6); solid line: solution with additional second step (one-pole) approximation (8).

Caption: FIGURE 5: Setup for simulations of small dipole in vicinity of human body. (a) Axial (z-) orientation; (b) angular ([phi]-) orientation; (c) radial ([rho]-) orientation.

Caption: FIGURE 6: Simulated total electric field around cylinder (human body) for three respective dipole orientations.

Caption: FIGURE 7: Plot of field magnitudes around the human body model together with field scale. (a) Axial dipole excitation; (b) Radial dipole excitation.

Caption: FIGURE 8: Snapshots of the field distribution. (a) Axial dipole excitation; (b) radial dipole excitation. The relevant field scale is given in Figure 7.

Caption: FIGURE 9: The on-body PIFA antenna. (a) The antenna geometry; (b) the manufactured prototype; (c) close view of E-field at the radiating slot and the equivalent magnetic current.

Caption: FIGURE 10: The realistic on-body link around torso. (a) Simulation in CST; (b) measurement setup; (c) cross section of volunteer torso with measurement points.

Caption: FIGURE 11: Simulated electric field around the body model for excitation with vertically oriented PIFA antenna, magnitudes of total field and field components.

Caption: FIGURE 12: Calculated and measured magnitudes of transmission parameter of wireless on-body link around torso.

TABLE 1: Dominant normalized poles of Green's function.

                  Real part   Imaginary part

T[M.sub.z] case     1.225         0.446
T[E.sub.z] case     1.08          0.222
COPYRIGHT 2017 Hindawi Limited
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2017 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Title Annotation:Research Article
Author:Ivsic, Branimir; Bonefacic, Davor; Sipus, Zvonimir; Bartolic, Juraj
Publication:Wireless Communications and Mobile Computing
Article Type:Report
Date:Jan 1, 2017
Words:4680
Previous Article:Investigation of Sphere Decoder and Channel Tracking Algorithms for Media-Based Modulation over Time-Selective Channels.
Next Article:Privacy-Preserving Meter Report Protocol of Isolated Smart Grid Devices.
Topics:

Terms of use | Privacy policy | Copyright © 2020 Farlex, Inc. | Feedback | For webmasters