Instruction FACTA UNIVERSITATIS Series: Electronics and Energetics Vol. 32, N o 4, December 2019, pp. 555-569 https://doi.org/10.2298/FUEE1904555V © 2019 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND THE INFLUENCE OF CONDUCTIVE PASSIVE PARTS ON THE MAGNETIC FLUX DENSITY PRODUCED BY OVERHEAD POWER LINES  Slavko Vujević, Tonći Modrić University of Split, Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture, Split, Croatia Abstract. There has been apprehension about the possible adverse health effects resulting from exposure to power frequency magnetic field, especially in the overhead power lines vicinity. Research work on the biological effects of magnetic field has been substantial in recent decades. Various international regulations and safety guidelines, aimed at the protection of human beings, have been issued. Numerous measurements are performed and different numerical algorithms for computation of the magnetic field, based on the Biot- Savart law, are developed. In this paper, a previously developed 3D quasistatic numerical algorithm for computation of the magnetic field (i.e. magnetic flux density) produced by overhead power lines has been improved in such a way that cylindrical segments of passive conductors are also taken into account. These segments of passive conductors form the conductive passive contours, which can be natural or equivalent, and they substitute conductive passive parts of the overhead power lines and towers. Although, their influence on the magnetic flux density distribution and on the total effective values of magnetic flux density is small, it is quantified in a numerical example, based on a theoretical background that was developed and presented in this paper. Key words: cylindrical segments of passive conductors, magnetic flux density, self and mutual potential coefficients, self and mutual impedances 1. INTRODUCTION Extremely high-voltage overhead power lines are among the most significant sources of extremely low-frequency (ELF) electric and magnetic fields. These fields change very slowly over time and are therefore considered quasistatic [1, 2]. Hence, electric and magnetic fields are computed separately. Potential long-term health effects of exposure arising from the power distribution system including overhead power lines, due to their proximity to residential areas and field levels they emit, have been extensively studied over Received February 22, 2019; received in revised form April 23, 2019 Corresponding author: Tonći Modrić University of Split, Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture, Rudjera Boškovica 32,21000 Split, Croatia (E-mail: tmodric@fesb.hr)  556 S. VUJEVIĆ, T. MODRIĆ the last few decades [3–6]. Epidemiological studies have suggested that long-term low-level exposure to ELF (50-60 Hz) magnetic field might be associated with an increased risk of childhood leukaemia. However, there is no firm evidence of such adverse health effects related to prolonged exposure to ELF magnetic field. The International Commission on Non-Ionizing Radiation Protection (ICNIRP) recommends the exposure limits related to short-effects, called basic restrictions in their guidelines [7]. The exposure limits outside the body, called reference levels, are set to 1000 μT for occupational exposure and to 200 μT for general public exposure for 50 Hz magnetic flux density. The assessment of human exposure to the magnetic field originating from overhead power lines is based on measurements or computations. Measurements are performed in accordance with international standards [8– 10]. Various methods for computation of ELF magnetic field, such as method of moments (MoM), finite-difference time-domain method (FDTD), finite element method (FEM), charge simulation method (CSM) and surface charge simulation method (SCSM), are well- known [11]. The magnetic field of the overhead power lines is computed using the Biot-Savart law. In simplified two-dimensional (2D) numerical algorithms [12–14] overhead power line conductors are infinitely long straight thin-wire horizontal lines parallel to the flat Earth’s surface. The number of line sources equals the number of overhead power line phase conductors and shield wires, and the contribution of each of them is taken into account. In three-dimensional (3D) numerical algorithms [15–19] the catenary form of the overhead power line conductors can be taken into account and therefore, more accurate computation results can be obtained. The basis of this paper is a previously developed 3D algorithm for computation of the magnetic field (i.e. magnetic flux density) produced by overhead power lines. The catenary conductors are approximated by a set of straight thin-wire cylindrical conductor segments. Moreover, the cylindrical segments of passive conductors are also taken into account using closed current contours, which substitute conductive passive parts of the overhead power lines and towers. Because of the currents induced in them, they may have the influence on the magnetic flux density distribution. As anticipated, the influence is not as pronounced as in computation of the electric field intensity, especially in close vicinity of the towers and other passive parts that strongly distort the electric field. Therefore, many researchers ignore their effect, but nevertheless, a theoretical background for taking into account conductive passive parts and their effect on the magnetic flux density distribution is developed herein. Moreover, this influence has been quantified for the first time so far. Expressions for self and mutual potential coefficients of cylindrical conductor segments are given. Equations for self and mutual impedances per unit length of the conductive passive contours are derived and included in the system of linear equations for computation of currents in natural and equivalent conductive passive contours. Finally, the sum of the contributions of cylindrical segments of active and passive conductors is taken into account for computation of the magnetic flux density. In the numerical example, two different cases are observed, the first where conductive passive parts (CPPs) are neglected and the second where they are taken into account. The obtained results of the magnetic flux density distribution are shown and compared with available results from the literature. The Influence of Conductive Passive Parts on the Magnetic Flux Density 557 2. CONDUCTIVE PASSIVE PARTS OF THE OVERHEAD POWER LINES The magnetic flux density distribution at the arbitrary field point T (x, y, z) in the air of a two-layer medium can be computed using the well-known Biot-Savart law. One of the advanced 3D numerical algorithms for computation, sufficiently accurate as computation module HIFREQ of the CDEGS software package, is presented in detail in [19]. In addition to the cylindrical segments of active conductors with known currents, the cylindrical segments of passive conductors can also be taken into account. In these cylindrical segments of passive conductors (i.e. current contours), the currents are induced and therefore, they have influence on the magnetic flux density distribution. Conductive passive parts of the overhead power lines and towers can be described using closed current contours, approximated by a set of cylindrical conductor segments. The contours can be natural (Fig. 1) or equivalent, that substitute parts of conductive passive surfaces. Examples of conductive passive surfaces are overhead power line towers, fences or any other conductive passive parts in high-voltage substations, which can be modelled using straight thin-wire cylindrical conductor segments or using subparametric spatial 2D finite elements, as in [20]. Hence, a network model of conductive passive surfaces is used herein. The cylindrical conductor segments, that form the conductive passive contours, are oriented from the start to the end point of the segment. The unit vector 0s  is assigned to each cylindrical segment. Fig. 1 Closed passive contour approximated with 5 cylindrical conductor segments 2.1. Currents of the Conductive Passive Contours The system of linear equations for computation of currents in conductive passive contours, written in matrix form, can be expressed as follows:                                                  S NS S KS NSNK KS NK KS NS KS K NK K KK NKNK KK NK KK NK KK I I ZZ ZZ I I ZZ ZZ         1 ,1, ,11,11 ,1, ,11,1 (1) where NK is the total number of conductive passive contours, NS is the total number of active cylindrical conductor segments, K ik I is the phasor of the ik-th conductive passive 558 S. VUJEVIĆ, T. MODRIĆ contour current, S k I is the phasor of the k-th active cylindrical conductor segment current, KK jkik Z , is the mutual impedance of the ik-th and jk-th conductive passive contour, KS kik Z , is the mutual impedance of the ik-th conductive passive contour and k-th cylindrical segment of active conductor. Self impedance of the ik-th conductive passive contour is described by:          1 1 1 , , 1 , , ,1 , 2 ik NS is ik NS jk NSjs ikik jsis ik NS is ikik isis ik is iku is KK ikik LjLjZZ  (2) where ik NS is the total number of segments of the ik-th conductive passive contour, iku is Z ,1 is the internal impedance per unit length of the is-th cylindrical conductor segment of the ik-th conductive passive contour, ik is  is the length of the is-th cylindrical conductor segment of the ik-th conductive passive contour, ikik isis L , , is the external inductance of the is-th cylindrical conductor segment of the ik-th conductive passive contour, ikik jsis L , , is the mutual inductance of the is-th and js-th cylindrical conductor segments of the ik-th conductive passive contour,  is the circular frequency and j is the imaginary unit. The internal impedance per unit length of the cylindrical conductor segment of the natural conductive passive contour is described by following expression [21, 22]: 1 0 0 0 1 0 ( ) 2 ( ) u J k rk Z r J k r         (3) where k is the complex wave number, 0 r is the radius of the cylindrical conductor segment, v  is the electrical conductivity of the cylindrical conductor segment, 0 0 ( )J k r is the complex Bessel function of the first kind of order zero, 01( )J k r is the complex Bessel function of the first kind of order one. The complex wave number k is defined by the following equation: )4/(exp2  jfk vv (4) where  is the magnetic permeability of the cylindrical conductor segment, f is the time- harmonic current frequency. The complex Bessel function of the first kind of order Nn  can be written as [23, 24]: 2 0 2 ( ) ( 1) ! ( )! n m m n m k r J k r m n m                 (5) Conductive passive surface can be replaced by a contour formed by a set of equivalent cylindrical conductor segments. Radius of these segments is equal to [23]: vv v f d r   2 1 2 (6) where d is the skin depth of the wave into the conductive surface. The Influence of Conductive Passive Parts on the Magnetic Flux Density 559 Internal impedance per unit length of these conductor segments can be described using the following expression: 1 (1 ) 2 22 pu v v vv v E Z Z j rr H                (7) where v Z is the wave impedance of the medium from which the conductive passive surface is made, vE is the phasor of electric field intensity on the surface of the conductor, vH is the phasor of magnetic field intensity on the surface of the conductor, p is the magnetic permeability of the surface. Mutual impedance of the ik-th and the jk-th conductive passive contour is described by the following equation: KK ikjk ik NS is jk NS js jkik jsis KK jkik ZLjZ , 1 1 , ,,      (8) where jkNS is the total number of segments of the jk-th conductive passive contour, jkik jsis L , , is the mutual inductance of the is-th cylindrical conductor segment of the ik-th conductive passive contour and js-th cylindrical conductor segment of the jk-th conductive passive contour. Mutual impedance of the ik-th conductive passive contour and k-th cylindrical segment of active conductor is defined by the following expression:    ik NS is ik kiskik LjZ KS 1 ,, (9) where ik kis L , is the mutual inductance of the is-th cylindrical conductor segment of the ik- th conductive passive contour and k-th cylindrical segment of active conductor. 2.2. Self and Mutual Inductances of the Cylindrical Conductor Segments Self inductance of the cylindrical conductor segment is described by the following expression: isis un isis PSSεL ,00,  (10) where isis un PSS , is the self potential coefficient of the is-th cylindrical conductor segment in homogeneous unbounded dielectric medium with permittivity 0  , which can be computed as described in detail in chapter 3. Mutual inductance of the is-th and js-th cylindrical conductor segment, which can be cylindrical segment of active conductor or part of conductive passive contour, is described using the following expressions: ,, 0 0 0 0 ,( ) un is jsis js is js js isL s s PSS L       (11) where is s 0  is the unit vector of the is-th cylindrical conductor segment, jss0  is the unit vector of the js-th cylindrical conductor segment, jsis un PSS , are the mutual potential 560 S. VUJEVIĆ, T. MODRIĆ coefficients of the is-th and js-th cylindrical conductor segments in homogeneous unbounded dielectric medium with permittivity 0  , which can be computed as described in detail in chapter 3. 3. SELF AND MUTUAL POTENTIAL COEFFICIENTS OF THE CYLINDRICAL CONDUCTOR SEGMENTS Self potential coefficients of the cylindrical conductor segments in homogeneous unbounded dielectric medium are described by the following expression: 0 0 ( , ) 4 π ε un i ii P r PSS    (12) where auxiliary function P is defined [25] as:             vvv v vP 22 222 ln2),(    (13) According to [26], two cylindrical conductor segments can be parallel or nonparallel. Two parallel cylindrical conductor segments, i-th and j-th, in homogeneous unbounded dielectric medium are shown in Fig. 2. Further, i-th cylindrical conductor segment with endpoints T1 (u1, vi) and T2 (u2, vi) is observed in the local coordinate system (u, v) of the j-th cylindrical conductor segment. Fig. 2 Two parallel cylindrical conductor segments in homogeneous unbounded dielectric medium Mutual potential coefficients of two parallel cylindrical conductor segments, i-th and j-th segment, in homogeneous unbounded dielectric medium can be obtained from the following expression: 1 2 3 4 0( ) / 4 un ijPSS C C C C ε      (14) where auxiliary function Ck (k = 1, 2, 3, 4) is described by: 2222 ln ikkikkk vwwvwwC       (15) The Influence of Conductive Passive Parts on the Magnetic Flux Density 561 2/21 juw  (16) 2/12 juw  (17) 2/13 juw  (18) 2/24 juw  (19) In a case of two nonparallel cylindrical conductor segments, there is always one and only one pair of parallel planes, 1 and 2 in which these segments lie (Fig. 3). In a limiting case, these two planes overlap and then, nonparallel segments lie on intersecting straight lines. Fig. 3 Two nonparallel cylindrical conductor segments in homogeneous unbounded dielectric medium Mutual potential coefficients of two nonparallel cylindrical conductor segments in homogeneous unbounded dielectric medium, defined by using the Galerkin-Bubnov method, are described by:         2 1 2 1 0 ε4 1 ij ij un R dd PSS (20) where the mutual distance between points on the axes of the segments is equal to:  cos2 222 DR ij (21) where D is the distance between the parallel planes on which segments lie,  is the angle between lines on which the segments lie,  is the distance of the observed point on the axis of the i-th segment and the start point O1,  is the distance of the observed point on the axis of the j-th segment and the start point O2. 562 S. VUJEVIĆ, T. MODRIĆ After the double integration in (20) is carried out, the following expression, known as Cejtlin’s formula [25], can be obtained: 1 1 2 2 1 2 2 1 0 0 ( , ) ( , ) ( , ) ( , ) 4 4 un ij A A A A PSS                     (22) where 1 ( , ) ln ( cos ) ln( cos ) 2 tan tan sin 2 ij ij ij A R R RD D                                   (23) Self potential coefficients of the i-th cylindrical conductor segment, with linear charge density i λ approximated by a constant, in the air of the two-layer medium (Fig. 4) can be computed, using the well-known image method: sii un rii un ii PSSkPSSPSS  (24) where sii un PSS is the mutual potential coefficient of the i-th cylindrical conductor segment and its image in homogeneous unbounded dielectric medium with permittivity 0, whereas rk is the reflection coefficient derived for a point current source [23, 27, 28] which can be approximated to high accuracy by 1 r k for power line frequencies as a consequence of assumption that the Earth’s conductivity is infinite. Fig. 4 Cylindrical conductor segment in the air of the two-layer medium and its image If i-th and j-th cylindrical conductor segments are in the air of the observed two-layer medium (Fig. 5), their mutual potential coefficient can be expressed by image method: sij un rij un ij PSSkPSSPSS  (25) where sij un PSS is the mutual potential coefficient of the i-th cylindrical conductor segment and image of the j-th cylindrical conductor segment in homogeneous unbounded dielectric medium with permittivity . 0  The Influence of Conductive Passive Parts on the Magnetic Flux Density 563 Fig. 5 Two cylindrical conductor segments and the image from one of them 4. NUMERICAL EXAMPLE In order to estimate the influence of conductive passive parts on the magnetic flux density distribution, a computer program was developed, in which these passive parts can be considered on the basis of the presented theory. In the numerical example, two spans between three identical towers of a typical 400 kV overhead power line, each carrying three phases with two conductors in the bundle per phase and two shield wires are observed (Fig. 6). Detailed input data concerning the tower geometry, the maximum and minimum heights of all conductors and sags, radii of all phase conductors and shield wires, the length of the overhead power line span, as well as electrical input data are given in [20]. The maximum allowed symmetrical currents for cross section of the chosen phase conductors and symmetrical operating conditions have been assumed. Two different cases are observed. In the first case, only phase conductors and shield wires (16 catenaries, each approximated using 60 thin-wire cylindrical segments of active and passive conductors) are taken into account, whereas conductive passive parts (CPPs) are neglected. In the second case, in addition to aforementioned catenaries, a central tower is approximated using 68 thin-wire cylindrical segments of passive conductors and 40 conductive passive contours are taken into account. The electrical conductivity of the cylindrical conductor segments  is equal to 7 MS/m, while the magnetic permeability of the cylindrical conductor segments  is equal to 500. Computation of the magnetic flux density distribution is carried out at height of 1 m above the Earth’s surface in the close vicinity of a central tower along observed x- and y-axes in a total of 500 points. Fig. 6 Simplified representation of the overhead power line 564 S. VUJEVIĆ, T. MODRIĆ Figures 7–10 present computed effective (rms) values of the total magnetic flux density and its components along x-axis, whereas Figures 11–14 present computed effective (rms) values of the total magnetic flux density and its components along y-axis for the aforementioned two cases. Maximum absolute deviations of the computed total magnetic flux density distribution along x- and y-axes for two cases in the chosen example are equal to 0.15 % and 0.89 %, respectively. As expected, according to well-known parameters affecting the magnetic flux density distribution, these absolute deviations due to conductive passive parts are small. Nevertheless, they have not been quantified so far. The maximum computed value of the magnetic flux density, obtained in this example, in the close vicinity of a central tower (Fig. 7) is equal to 16.84 μT, as well as the maximum computed value obtained under the midspan of the overhead power lines, which is equal to 31.83 μT, are substantially less than the exposure limits given in [7]. Fig. 7 Distribution of the total magnetic flux density along x-axis Fig. 8 Magnetic flux density x-component along x-axis The Influence of Conductive Passive Parts on the Magnetic Flux Density 565 Fig. 9 Magnetic flux density y-component along x-axis Fig. 10 Magnetic flux density z-component along x-axis Fig. 11 Distribution of the total magnetic flux density along y-axis 566 S. VUJEVIĆ, T. MODRIĆ Fig. 12 Magnetic flux density x-component along y-axis Fig. 13 Magnetic flux density y-component along y-axis Fig. 14 Magnetic flux density z-component along y-axis The Influence of Conductive Passive Parts on the Magnetic Flux Density 567 In order to verify the accuracy of the presented algorithm, the magnetic flux density results computed by EFC-400EP software [29] are shown in several points are compared to computed results obtained by numerical algorithm given herein (Fig. 15) and a very good agreement can be seen. Detailed input data of a 400 kV overhead power line are given in [29]. Fig. 15 Comparison of computed magnetic flux density results obtained by presented algorithm with results computed by EFC-400EP software Table 1 shows percent errors (p.e.) of magnetic flux density results obtained by presented algorithm with respect to results computed by EFC-400EP software, in chosen points, given in Fig. 15, along observed x-axis. Table 1 Percent errors of magnetic flux density results obtained by presented algorithm with respect to results computed by EFC-400EP software x (m) p.e. (%) 0 3.276 2.5 3.033 5.0 2.727 7.5 0.995 10.0 0.389 12.5 2.021 15.0 2.541 17.5 3.317 20.0 1.588 22.5 3.503 25.0 4.807 5. CONCLUSION In this paper, a 3D quasistatic numerical model for taking into account conductive passive parts of the overhead power lines and their effect on the computation of the magnetic flux density distribution is presented. The catenary conductors of the overhead 568 S. VUJEVIĆ, T. MODRIĆ power line span are approximated by a set of straight thin-wire cylindrical conductor segments. Besides cylindrical segments of active conductors, the cylindrical segments of passive conductors are also taken into account using closed current contours, which can be natural or equivalent. These conductive passive parts have small influence on the magnetic flux density distribution, which has been quantified herein. Primarily, it is due to extremely low-frequency of the magnetic flux density produced by overhead power lines. An originally developed theoretical background is described in detail, including expressions for self and mutual potential coefficients of cylindrical conductor segments and expressions for self and mutual impedances per unit length of the conductive passive contours. An influence of conductive passive parts on the magnetic flux density is shown and quantified in the chosen numerical example of a typical 400 kV overhead power line. REFERENCES [1] R. G. Olsen and P. S. Wong, "Characteristics of low frequency electric and magnetic fields in the vicinity of electric power lines", IEEE Transactions on Power Delivery, vol. 7, no. 4, pp. 2046–2055. [2] R. Fitzpatrick, Maxwell’s Equations and the Principles of Electromagnetism. Infinity Science Press LLC, Hingham, 2008. [3] N. Wertheimer and E. Leeper, "Electrical wiring configurations and childhood cancer", American Journal of Epidemiology, vol. 109, no. 3, pp. 273–284, 1979. [4] J. C. Teepen and J. A. van Dijck, "Impact of high electromagnetic field levels on childhood leukemia incidence", International Journal of Cancer, vol. 131, no. 4, pp. 769–778, 2012. [5] International Agency for Research on Cancer, Monographs on the Evaluation of Carcinogenic Risks to Humans, Non-Ionizing Radiation, Part 1: Static and Extremely Low-Frequency (ELF) Electric and Magnetic Fields, vol. 80, IARCPress, Lyon, France, 2002. [6] World Health Organization, Extremely low frequency fields, Environmental Health Criteria Monograph No. 238, WHO Press, Geneva, 2007. [7] International Commission on Non-Ionizing Radiation Protection, Guidelines for limiting exposure to time-varying electric and magnetic fields (1 Hz to 100 kHz), Health Physics, vol. 99, no. 6, pp. 818–836, 2010. [8] IEEE Standard Procedures for Measurement of Power Frequency Electric and Magnetic Fields from AC Power Lines, IEEE Standard 644-1994. DOI: 10.1109/IEEESTD.1995.122621. [9] Measurement of DC magnetic, AC magnetic and AC electric fields from 1 Hz to 100 kHz with regard to exposure of human beings – Part 1: Requirements for measuring instruments, IEC Std 61786–1, 2013. [10] Measurement of DC magnetic, AC magnetic and AC electric fields from 1 Hz to 100 kHz with regard to exposure of human beings – Part 2: Basic standard for measurements, IEC Standard 61786–2, 2014. [11] P. Zhou, Numerical Analysis of Electromagnetic Fields. Springer-Verlag, Berlin Heidelberg, 1993. [12] C. Garrido, A. F. Otero and J. Cidrás, "Low-frequency magnetic fields from electrical appliances and power lines", IEEE Transactions on Power Delivery, vol. 18, no. 4, pp. 1310–1319, 2003. [13] G. Filippopoulos, D. Tsanakas, "Analytical calculation of the magnetic field produced by electric power lines", IEEE Transactions on Power Delivery, vol. 20, no. 2, pp. 1474–1482, 2005. [14] F. Moro and R. Turri, "Fast analytical computation of power-line magnetic fields by complex vector method", IEEE Transactions on Power Delivery, vol. 23, no. 2, pp. 1042–1048, 2008. [15] G. Lucca, "Magnetic field produced by power lines with complex geometry", European Transactions on Electrical Power, vol. 21, no. 1, pp. 52–58, 2011. [16] A. Z. El Dein, "Magnetic-field calculation under EHV transmission lines for more realistic cases", IEEE Transactions on Power Delivery, vol. 24, no. 4, pp. 2214–2222, 2009. [17] J. C. Salari, A. Mpalantinos and J. I. Silva, "Comparative analysis of 2- and 3-D methods for computing electric and magnetic fields generated by overhead transmission lines", IEEE Transactions on Power Delivery, vol. 24, no. 1, pp. 338–344, 2009. [18] B. Zemljaric, "Calculation of the connected magnetic and electric fields around an overhead-line tower for an estimation of their influence on maintenance personnel", IEEE Transactions on Power Delivery, vol. 26, no. 1, pp. 467–474, 2011. The Influence of Conductive Passive Parts on the Magnetic Flux Density 569 [19] T. Modrić, S. Vujević and D. Lovrić, "3D computation of the power lines magnetic field", Progress in Electromagnetics Research M, vol. 41, pp. 1–9, 2015. [20] T. Modrić and S. Vujević, "Computation of the electric field in the vicinity of overhead power line towers", Electric Power Systems Research, vol. 135, pp. 68–76, 2016. [21] S. Vujević, V. Boras and P. Sarajčev, "A novel algorithm for internal impedance computation of solid and tubular cylindrical conductors", International Review of Electrical Engineering, vol. 4, no. 6, pp. 1418–1425, 2009. [22] D. Lovrić, V. Boras and S. Vujević, "Accuracy of approximate formulas for internal impedance of tubular cylindrical conductors for large parameters", Progress in Electromagnetics Research M, vol. 16, pp. 177– 184, 2011. [23] J. A. Stratton, Electromagnetic Theory. McGraw-Hill Book Company, New York and London, 1941. [24] M. R. Spiegel, S. Lipschutz and J. Liu, Mathematical Handbook of Formulas and Tables, fourth ed., McGraw-Hill Education, New York, 2012. [25] L. R. Neiman and P. L. Kalantarov, Theoretical Fundamentals of Electrical Engineering, Part 3: Theory of Electromagnetic Field, Gosenergoizdat, Moscow, Leningrad, 1959 (in Russian). [26] W. H. McCrea, Analytical Geometry of Three Dimensions. Dover Publications, New York, 2006. [27] S. Vujević and P. Sarajčev, "Potential distribution for a harmonic current point source in horizontally stratified multilayer medium", COMPEL: The International Journal for Computation and Mathematics in Electrical and Electronic Engineering, vol. 27, no. 3, pp. 624–637, 2008. [28] T. Takashima, T. Nakae and R. Ishibashi, "High frequency characteristics of impedances to ground and field distributions of ground electrodes", IEEE Transactions on Power Apparatus and Systems, vol. PAS- 100, no. 4, pp. 1893–1900, 1981. [29] S. Carsimamovic, Z. Bajramovic, M. Rascic, M. Veledar, E. Aganovic and A. Carsimamovic, "Experimental results of ELF electric and magnetic fields of electric power systems in Bosnia and Herzegovina", In Proceedings of EUROCON 2011, International Conference on Computer as a Tool, Lisbon, Portugal, 2011, pp. 1–4.