Fourier Transform Application in the Computation of Lightning Electromagnetic Field

Atmospheric discharge is one of the most interesting and powerful natural phenomenon hiding from men its undiscovered features and secrets for centuries. Lightning discharges have been studied in many theoretical and experimental ways. However, application of Fourier transform has been introduced in this research only in recent few decades. It proved very useful in solving many lightning research problems.


Introduction
Atmospheric discharge is one of the most interesting and powerful natural phenomenon hiding from men its undiscovered features and secrets for centuries. Lightning discharges have been studied in many theoretical and experimental ways. However, application of Fourier transform has been introduced in this research only in recent few decades. It proved very useful in solving many lightning research problems.
There are two main groups of problems in lightning studies that involve using Fourier transform. One deals with determining how the energy is distributed over a continuous frequency spectrum for the quantity of interest. Channel-base currents, induced voltages and currents, electric and magnetic field components, so as integrals and derivatives of the same functions, distribute components over the entire frequency range in different ways. These are non-periodic functions scattering their energy throughout the frequency spectra.
Another important group of problems to be solved by Fourier transform deals with the calculation of lightning induced effects at different distances from the lightning discharge and risk assessment for buildings, various structures, people and property in an external impulse electromagnetic field. These calculations are based on experimental results for the measured lightning electromagnetic field (LEMF) and channel-base currents, which are used in different types of lightning stroke models including lossy ground effects. Lightning discharge channel is usually modeled by a thin vertical antenna at a lossy ground of known electrical parameters. Both ground and air are treated as linear, isotropic and homogeneous half-spaces. Even for such a simple approximation of the lightning channel -calculations can not be easily done in time domain, but transformation to frequency domain is used instead. Once the calculations are done in frequency domain, way back to time domain is made by Inverse Fourier transform applied to the obtained results.
In both groups of problems an impulse function which has the analytical derivative, integral and integral transformations is very useful. New functions proposed by the author for representing lightning currents are presented in this Chapter, Section 2. These can be also used in other high voltage technique calculations. The main problem for a user of the impulse functions already given in literature to approximate some quantity is the choice of parameters so to obtain desired waveshape characteristics or values adequate to experimentally measured. Parameters of the new channel-base current (NCBC) function (Javor, 2008;Javor & Rancic, 2011) are calculated according to IEC 62305 standard lightning currents (International Electrotechnical Commission, Technical Committee 81 [IEC, TC 81], 2006), and the procedure to choose function parameters is also explained in (Javor & Rancic, 2011). Functions to represent other typical lightning currents are proposed in (Javor, 2011b(Javor, , 2011c, such as long stroke current (LSC), and two-rise function (TRF) as a multi-peaked current. These functions can be used to obtain the desired peak value, rise-time, decaying time to half of the peak value, current steepness, integral of the function (representing also impulse charge), or integral of the square of the function (representing also specific energy), etc. Fourier transform is obtained analytically and the results are presented in Section 3. Application of the Fourier transform in LEMF computation is shown in Section 4. Based on these results, some conclusions are given in Section 5.

NCBC function
Any mathematical function capable to approximate the impulse lightning channel-base current in electromagnetic field calculations can be also used to represent far lightning electromagnetic field as an external excitation inducing currents and voltages in objects and structures inside such a field. Such functions are necessary in lightning stroke modeling.
There is an overview of lightning return-stroke models in (Rakov & Uman, 1998, 2006 where they are classified into four classes according to the type of governing equations. An adequate return-stroke model would be the one that provides simultaneously an approximation to the experimentally measured channel-base current, to the lightning electric and magnetic field waveshapes and intensities at various distances, and to the observed return-stroke speeds. Functions implied in these models are based either on the double-exponential (DEXP) function (Bewley, 1929;Bruce & Golde, 1941) or Heidler's function (Heidler, 1985). DEXP first order derivative at t=0 + is too large, thus causing numerical problems in LEMF calculations. It also has physically non-realistic convex waveshape in the rising part. Heidler's function reproduces concave rising part and its derivative is equal to zero at t=0+. It is also used for representing lightning currents in the International standard IEC 62305 (IEC/TC 81, 2006). In order to obtain better agreement with experimental results (Berger et al., 1975) one function was proposed (Nucci et al., 1990) as a linear combination of the DEXP and Heidler's function. Another pulse function was also proposed in (Feizhou & Shange, 2004). All these functions need peak correction factors, and their parameters cannot be easily chosen according to the desired waveshape. However, NCBC function parameters in the rising part can be chosen independently from parameters in the decaying part. The exact rising time to the peak value can be chosen in advance, and the decaying time to half of the peak value can be selected in the approximation procedure. The maximum value of the function can be chosen without the peak correction factor. Maximum current steepness can be adjusted (Javor & Rancic, 2011) using analytical expression for the first derivative. Experimentally measured impulse charge and specific energy values (Berger et al., 1975;Anderson & Eriksson, 1980) can be achieved using analytically obtained integral of NCBC function and integral of the square of the function (Javor, 2011a).
NCBC function (Fig. 1) where a and b i are parameters, c i coefficients, n the chosen number of expressions in the decaying part, so that the total sum of n weighting coefficients c i is equal to unit, and t m is the rise-time to the maximum current value I m . For n=1, c 1 =1 and b 1 =b, NCBC function reduces to CBC function (Javor & Rancic, 2006) with four parameters (I m , t m , a and b). In the special case, for n=1, a=4 and b=0.0312596735, CBC function reduces to High-Voltage Pulse (HVP) function 1.2/50 s (Velickovic & Aleksic, 1986). Impulse duration time is defined as t i =t k -t a' , for t k the time in which the current decreased to half of its peak value (Fig. 1). The rising part of the function and its front rise-time are given in Fig. 2. The front rise-time is defined as t c =t b' -t c' , for t b' and t c' determined as the time values corresponding to the points B' and A', obtained from intersecting horizontal lines for the maximum I m and the zero function value (time axis) with the line drawn through the points A, for i(t)=0.3I m (Fig. 2) or for i(t)=0.1I m in some other definitions, and B, for i(t)=0.9I m .

Fig. 1. Normalized NCBC function
NCBC function is an analytically prolonged mathematical function (but still continuous, so as its first derivative, whereas higher order derivatives are not continuous at the point of function maximum I m ), the parameter a in the rising part can be chosen to approximate the front of the waveshape independently from parameters b i and weighting coefficients c i in the decaying part, which facilitates the approximation procedure. NCBC function belongs to C 1 differentiability class. Parameters of NCBC function can be chosen so that it represents waveshape of the often used DEXP function with parameters given in (Bruce & Golde, 1941). DEXP function i(t)=I m [exp(-αt)-exp(t)] for I m =11kA, t m =0.5826 s, α=3 . 10 4 s -1 , and =10 7 s -1 , has the decreasing time to half of the peak value of approximately 23 s. It can be www.intechopen.com Fourier Transform Applications 60 approximated with NCBC function for n=1 and I m =11kA, t m =0.5826 s, a=0.5 and b=0.02, having also the decreasing time to half of the peak value of about 23 s. If using lightning stroke models and electromagnetic theory relations, lightning electric and magnetic field components above perfectly conducting ground in general have three terms, depending on the integral of the channel-base current function, on the function itself, and on the function derivative. Consequently, the DEXP function having just a convex rising part and great values of the first derivative at t=0 + makes numerical problems in LEMF calculations. Fig. 2. The rising part of the normalized CBC function representing current t c /t i =1.2/50 s NCBC function, for a > 1, satisfies the demand of having the first derivative equal to zero at t=0 + and the concave to convex rising part. The value of parameter a has to be chosen greater than 1 in order to obtain one saddle point in the rising part (Javor & Rancic, 2006). Heidler's function also has the first derivative equal to zero at t=0 + , but it needs peak correction factor, and its rise-time to the maximum current value cannot be chosen in advance. Heidler's function doesn't have analytical integral and its Fourier transform cannot be obtained analytically, but just numerically (Heidler & Cvetic, 2002).
NCBC function can approximate also other channel-base currents used in lightning returnstroke modeling, as the one proposed in (Nucci et al., 1990). This function is used in many papers and approximates well experimental results for subsequent negative strokes. In order to represent the same current two terms are needed in the decaying part of NCBC function (n=2) and its parameters are calculated as I m =11kA, t m =0.472 s, a=1.1, b 1 =0.16, c 1 =0.34, b 2 =0.0047, and c 2 =0.66.

Parameters of NCBC function
The normalized NCBC function for different values of parameter a, in the first 1 s is presented in Fig. 3, for t m =0.5826 s. The parameter a determines the rising part of the www.intechopen.com function. NCBC function is presented in Fig. 4 for different values of parameter b, for t m =1.9 s, but in a longer time period (300 s), as parameter b determines the decaying part of the function. If the rising time to the maximum value t m has values e.g. 0.5, 1, 2, 5, or 10 s, and other parameters are a=1.5 and b=0.01, the changes in the function waveshape are presented in Fig. 5.
The first derivative of CBC function (NCBC for n=1) is presented in Fig. 6 up to 1 s, for b = 0.03 (but note that parameter b is irrelevant for the rising part), and different values of a, t m and I m .

Integral of NCBC function
Integral of NCBC function is calculated as: Impulse charge is defined in IEC 62305 standard (IEC/TC 81, 2006) as the integral of the channel-base current function. For standard lightning currents impulse charge is calculated in (Javor, 2011a), so as specific energy as the integral of the square of NCBC function.

Fourier transform of NCBC function
Whether a function is given analytically, graphically or numerically, its Fourier transform can be obtained analytically, numerically or using some of the commercial programs. Faster or slower rising/decaying of the function and its waveshape at the end of the covered impulse duration time determine the number N of needed points for Fast Fourier Transform.
For NCBC function Fourier and Laplace transforms are calculated analytically. It should be noted that we obtain Fourier transform for any Fourier-transformable function which is zero for t<0 if we substitute variable s with j2πf in its Laplace transform.
The analytical expression for the unilateral Laplace transform of NCBC function is: as defined in (Abramowitz, Stegun, 1970).
Fourier transform of NCBC function is obtained from (4) for s= j2πf, so: It can be obtained also numerically, or by using some computer program such as FAS (Walker, 1996) with the application of FFT (Fast Fourier Transform). Any frequency domain calculations can be done just for the positive frequencies, as for FFT of real functions the following relations are valid: This feature makes the time for computing FFT twice shorter for real functions, which is very useful if the number of FFT points is great so that calculations are time-consuming. For smaller number of FFT points the computing time is shortened. FFT results given for NCBC function representing lightning current indicate that the major part of its power is in the lower frequency range. The modulus of FFT decays fast as the function of frequency, so it is better to use some other than linear scaling of frequencies in order to cover the lower part of the frequency range. For very high frequencies, the values of real and imaginary parts of FFT of the pulse functions near the end of the frequency interval can be taken as zeros, as they are relatively very small comparing to the values at lower frequencies. Fourier transform of Heidler's function is calculated approximately (Andreotti et al., 2005). Some FFT results for Heidler's function are also given in (Vujevic & Lovric, 2010). NCBC/CBC have analytical solutions for Fourier transform which enables obtaining analytical solution for LEMF results in the case of perfectly conducting ground.
The modulus of Fourier transform of HVP function 1.2/50 s is presented in Fig. 7, for I(f) normalized to I(0), as the function of log 10 f. For f 2 ≈6.4kHz≈10 3.8 Hz the normalized modulus of FFT is |I(f 2 )/I(0)|≈0.3657. The results are obtained for FFT calculated in N=8192 points, for the time step ∆T≈19.064ns, frequency step ∆f≈6.4kHz, and limit frequency f g ≈52.455 MHz. From Fig. 7 is obvious that using frequencies higher than a few MHz is not necessary.
For some frequencies up to the chosen limit frequency the results are presented in Table 1 for FFT of HVP function 1.2/50 s calculated in N=8192 points, for ∆T≈19.064ns, ∆f≈6.4kHz, and f g ≈52.455MHz. For the chosen number of points for FFT, the corresponding time interval is T=N∆T≈156.17 s and the frequency step is ∆f=(N∆T) -1 = T -1 . For such values there are 8192 frequencies in the chosen frequency interval [-0.5f g ,+0.5 f g ], and the highest positive frequency is f 4097 ≈26.2275MHz. FFT results for |I(f)|≈10 -n |I(0)| for n=1, 2, 3, …, 8, are also given in Table 1, corresponding approximately to frequencies f 5 , f 35 , f 108 , f 352 , f 2478 , f 3912 , f 4078 , f 4095 , and they also point out to the fast decaying of FFT modulus. As (6) and (7) are valid for the real functions, this makes the time for computing twice shorter. This is very useful if the number of FFT points is great and relevant calculations in frequency domain are timeconsuming. For smaller number of FFT points the total computing time is shortened, but the sampling interval in time domain is critical in that case.  9 shows imaginary part of FFT for HVP function in the same considered frequency range [-0.5f g ,+0.5f g ], and enlarged rectangle area for fє [-200∆f, +200∆f]. It can be seen that Eq. (7) is valid.
FFT results given in Table 1 for this pulse function indicate that the major part of its power is in the lower frequency range. The modulus of FFT decays rapid as the function of frequency which can be noticed both in Table 1 and Fig. 7. Fig. 9. Imaginary part of FFT for HVP in the range [-0.5f g ,+0.5f g ]=[-4096∆f,+4096∆f] and enlarged rectangle area for frequencies fє [-200∆f ,+200∆f] As the value |I(f 2 )/I(0)|≈0.3657 is obtained already for the frequency f 2 ≈6.4kHz (Fig. 7), it is better to use some other then linear scaling of frequencies in order to cover better the lower part of the frequency range. For very high frequencies, the values of real and imaginary part of FFT near the end of the frequency interval can be taken as zeros, as they are relatively very small comparing to their values at low frequencies. That doesn't make much influence on the computation results if IFFT (Inverse Fast Fourier Transform) is done for obtaining time domain results based on frequency domain calculations. Both FFT and IFFT can be done using program FAS (Walker, 1996(Walker, ) in 128, 256, 512, 1024(Walker, , 2048 or 8192 points. The input data can be given as analytical functions (as in the case of NCBC or other channel-base current functions) or discrete values in the corresponding set of values for that quantity (as in the case of experimentally measured data).
In order to analyze lightning electromagnetic field in a far zone, as an external excitation of some receiving antenna structure, vertical electric field and azimuthal magnetic field can be represented with the waveshape of NCBC function according to the measured data and its Fourier transform should be known. Fourier transform of the channel-base current is also necessary if calculating equivalent voltage source as a product of the input impedance and the channel-base current in frequency domain (Moini et al. 2000, Shoory et al. 2005, instead of the current source itself at the channel-base. A few computer programs for antenna analysis are used for calculations (Richmond, 1974(Richmond, , 1992Bewensee, 1978;Burke & Poggio, 1980Djordjevic et al., 2002), and most of these use the voltage generator as excitation (Grcev et al., 2003).  Real and imaginary parts of FFT in the frequency range fє [-200∆f ,+200∆f] are presented in Figs. 11 and 12.
Log-log dependence of FFT modulus on frequency is presented in Fig. 13   The results show that for frequencies higher than 10MHz FFT modulus of these channelbase current functions is less than 1 o / oo of its value at 10kHz, so for frequencies higher than 10MHz calculations of LEMF are not needed in frequency domain.

Electromagnetic field calculation using antenna model of the lightning channel at a lossy ground
Based on experimentally measured characteristics of natural lightning (Berger et al., 1975;Lin et al. 1979;Anderson & Eriksson 1980) and also artificially initiated lightning discharges, LEMF can be estimated using some of the models from literature (Rakov & Uman, 1998, 2006, and lightning currents are assumed to propagate with attenuation and distortion while distributing charge along the channel. In the case of perfectly conducting ground calculations are simple in time domain, but full-wave approach in the case of a lossy ground is complex to implement, so Fourier transform is applied. After calculations are done in frequency domain, the conversion to time domain is performed by Inverse Fast Fourier transform (IFFT). A solution of the Sommerfeld's problem is required for the lossy ground. A few alternatives are proposed in literature.

Alternatives to full-wave time domain computation
The problem of LEMF calculation can be easily solved directly in time domain if the ground is treated as perfectly conductive. In such a case there exist vertical electric and azimuthal magnetic field components in the observed point at the ground surface. Horizontal component of electric field is zero at the perfectly conducting ground surface, but non-zero above the ground surface. However, it exists above, under, and at the surface of the lossy ground. Vertical component of the lightning electric and azimuthal component of the magnetic field can be easily (but in that case approximately) determined at the distances greater than a kilometer under the assumption of perfectly conducting ground. For smaller distances, propagation above ground of finite conductivity results in the noticeable attenuation of the high-frequency components of electric and magnetic field, and thus in appearance of the horizontal electric field at the surface. Finite conductivity has greater impact on horizontal than on vertical electric field, so calculation of horizontal component requires rigorous computation or, at least, acceptable approximations.
Approximate formulas in frequency domain are often used for calculation of the horizontal electric field in air, up to heights of tens of meters above the ground surface. These formulas can be integrated in the calculation of LEMF in time domain, but the obtained expressions are much more complex. There are simple approximations: the assumption of perfectly conducting ground, "wavetilt" formula, Cooray's approach and Rubinstein's approach. Cooray, 1992, proposed the calculation of horizontal electric field at the surface of finitely conductive ground using azimuthal magnetic induction and the expression for ground surface impedance. He showed that this simple formula provides very accurate results at the distances of about 200m. Rubinstein, 1996, proposed expression for the horizontal electric field with two terms: 1) horizontal electric field calculated under the assumption of perfectly conducting ground, and 2) the correction factor, given as a product of the magnetic field calculated under the same assumption and the function similar as in "wavetilt" formula which represents the effect of finite conductivity. The basic assumption of Rubinstein's approximation is σ 1 >> 0 r1 , and that finite ground conductivity does not affect the horizontal magnetic field at the surface. If this is not the case, then more general formula can be written, known in literature as the Cooray-Rubinstein's formula (Cooray, 2002). Wait gave generalization of Cooray-Rubinstein's formula and the exact evaluation of horizontal electric field, showing under which circumstances this general expression reduces to Cooray-Rubinstein's formula (Wait, 1997). Lindquist, 1983, andCooray, 1987, using the attenuation function in time domain proposed by Wait, 1956, included effects of the finite conductivity, and obtained results for the electric field that are in better agreement with experiments. Terms for approximate formulas in time domain are complex, so the approach in frequency domain is preferred.

Method of images
Method of images gives an approximate solution to the Sommerfeld's integral. Complex image technique often uses one or more terms of exponential series to approximate plane wave reflection coefficient. Thus, multiple discrete and/or continuous image sources are introduced, and this technique also proved not to be limited to a quasi-static range alone. Using this technique, different authors (Shubair & Chow, 1993;Yang & Zhou, 2004;Popovic & Petrovic, 1993;and Petrovic, 2005) obtained results for Sommerfeld's integral kernel for vertical dipoles above a lossy half-space which were used for the comparison with the new Two-image approximation (TIA). Approximate formulas are often valid for a limited range of ground electrical parameters, field point distances, or heights of dipoles above ground, but TIA approximation of Sommerfeld's integral kernel proposed in (Rancic & Javor, 2006 has the advantage of being valid in a wide range of lossy ground electrical parameters, for various heights of vertical electric dipoles above the ground, and for possibility to calculate electromagnetic field in both near and far zone. LEMF can be calculated using thin wire antenna modeling of a lightning channel assumed as vertical, without branches and reflections from the end. If using electromagnetic model, boundary conditions are fulfilled at the wire antenna surface, and a voltage or current source is assumed at the channel-base. Unknown current distribution along the antenna is determined by solving some system of integro-differential equations of EFIE or MFIE type satisfying boundary condition for electric or magnetic field components, respectively. One of those is e.g. System of integral equations of two potentials (SIE-TP), (Rancic, 1995), which is of EFIE type. Current distribution along the antenna can be approximated e.g. by polynomials (Popovic, 1970) with unknown complex coefficients. These can be determined by using Method of Moments (MoM) (Harrington, 1968) i.e. by matching in a sufficient number of points along the antenna. Based on that current distribution LEMF is calculated using electromagnetic theory relations. Thus, the problem is solved in full (without mentioned approximations of some LEMF components using others) on the basis of approximate solution of Sommerfeld's problem of a dipole radiating at the arbitrary height above the lossy half-space, which is a classical problem in electromagnetics.
The problem of a dipole radiation in the presence of a lossy medium was noticed yet in 1909, by Sommerfeld, who determined the solution in the form of superposition of inhomogeneous plane waves. His work had such a great influence on later theoretical research in this area that the solution of more complex problems which shows the influence of a lossy medium on the properties of linear antennas is marked as Sommerfeld's solution. Sommerfeld was the first who treated all four cases of elementary Hertz's dipoles in the air applying his formulation through the cylindrical waves. Van der Pol, 1931, proposed approximate method for solving Sommerfeld's integrals. Class of integrals obtained by Sommerfeld's formulation which strictly defines boundary conditions on the surface of discontinuity is known in literature as "integrals of Sommerfeld's type". These integrals have limits of integration 0 and infinity, and the integrand is a product of Bessel's function (J 0 or J 1 ), exponential, and one more function, so it is of very complex shape. Depending on the complexity of the model, this function has different forms so that different problems appear in its numerical integration. Since these integrals are highly oscillatory and slowly decreasing functions they don't have solutions in the closed form. In the analysis of antennas above lossy media, there is an inevitable and major problem to determine Sommerfeld's integral as accurate as possible in a wide range of values of electrical parameters of the medium, for different positions of current elements, different frequencies and positions of the observed point in the field. There are many different ways to approximate Sommerfeld's integral which can be found in literature. The model that includes complex mathematical functions in the case of lossy ground gives result equivalent to the field of infinite number of current images. Simpler approximations have limitations as for the values of electrical parameters of the medium, the position of the point in the field (near or far zone) and positions of matching points near or at the interface of two half-spaces. In order to include the influence of electrical parameters of a real ground on antenna characteristics, direct Sommerfeld's approach to the problem can be also used, but it is very complex for solving from both the analytical and numerical point of view. In addition to the direct approach to this problem, it is possible to modify the path of integration of Sommerfeld's integral, or to use interpolation and speed up the calculation, or to use tables of Sommerfeld's integrals, but there were also some attempts to solve the problem by obtaining approximations in closed form for some special cases of Sommerfeld's integrals. A suitable solution would be the one to determine efficiently Sommerfeld's integrals for arbitrary positions and orientations of dipoles, arbitrary positions of the points in LEMF and a wide range of frequencies.
One approximate solution is obtained by use of approximation of the spectral reflection coefficient of a plane wave with exponential series, having one or more terms which is named in literature as the method of images. For a horizontal electric dipole one complex image is the simplest approximation, which was first introduced by Wait, 1955, and later justified in (Wait & Spies, 1969). An important contribution was given by Bannister, (Bannister, 1966), who also showed that techniques of the theory of images in the case of finitely conductive ground are not restricted just to a quasi-static extent of the problem (Bannister, 1978). He later published a review paper "Summary of image theory expressions for the quasi-static fields of antennas at or above the earth's surface" (Bannister, 1979). Bannister extended the validity of the approximation to high-frequency problems by introducing exponential function exp(-u 0 d) where u 0 =(α 2 + 0 2 ) 1/2 and 0 is the complex propagation constant in the air. Thus the effects of propagation are included, the same solution obtained, and for both approximations there is the assumption that the refraction coefficient of ground-to-air is much greater than 1 (n 10 = 1 / 0 >>1). Mahmoud & Metwally, 1981, used discrete, so as the combination of discrete and continuous images. Mahmoud, 1984, extended the theory to multiple images and to a layered ground. Distributed images were also used and named "the exact theory of continuous images" in (Lindell & Alanen, 1984). By using Prony's method and nonlinear optimization technique Chow et al., 1991, analytically expressed Sommerfeld's integrals in closed form of spatial complex images. Shubair & Chow, 1993, discussed the impact of complex images on the problem of vertical antennas in the presence of lossy media and took advantage of the spatial Green's function in closed form in order to obtain a superposition of the sources influence, quasi-dynamic image and three complex images. Arand et al., 2003, used the method of discrete complex images and Generalized pencil of function technique to obtain the locations and impact of current sources. A procedure to approximate Sommerfeld's integrals was carried out by Popovic & Petrovic, 1993, named "a simple near-exact solution" which uses a few images determined so to approximately satisfy the boundary conditions in a limited area of the interface and to yield accurate field values in the domain of the antenna. Electromagnetic models including influence of the lossy ground on LEMF calculation often use approximations of Sommerfeld's integral kernel on the basis of theory of images in frequency domain. Takashima et al., 1980, proposed a modified theory of images for smaller distances from the point of observation to the vertical electric dipole. As in lightning research is necessary to treat also larger distances from the point of observation to the vertical electric dipole, this approximation can not satisfy the needs of computation in the range of distances of interest. For larger distances from the point of observation to the vertical electric dipole complete expressions, derived by King, 1990, were used for EM field of a vertical electric dipole over the lossy medium, if the condition | 2 | 2 >>| 1 | 2 or | 2 |>>3| 1 | is satisfied. King & Sandler, 1994, confirmed with their results the validity of these expressions over certain types of ground.
The new approximation TIA can be classified into two-image approximations. The basic idea for obtaining a new approximation of the spectral reflection coefficient is matching the approximation of reflection coefficient as a two-terms function (with three unknown constants) at two points, but also matching its first derivative at one point, which resulted in a very efficient approximation of Sommerfeld's integral. TIA gives good results for modeling in both near and far zone, as well as physically conceivable arrangement of images and representation of the problem (as the distances of the images are real values and not complex as in other approximations). TIA is similar to two-image approximations of Sommerfeld's integral, and therefore is their abbreviation retained, but the way of deriving the corresponding expressions is different. Sommerfeld's integral kernel results are compared to the results from literature for different values of lossy ground electrical parameters and various heights of vertical dipoles above the ground (Javor & Rancic, 2009). For the spectral reflection coefficient the following approximation is used 00 0 where A and B are unknown complex coefficients, d 0 is the distance from the source to the second image, and u 0c is the characteristic value chosen so that u 0c = 0 . Unknown constants are determined from equations obtained by matching the value of spectral reflection coefficient in two points (u 0 →∞ and u 0 = 0 ), and its first derivative at u 0 = 0 . This results in B=R ∞ , A=R 0 -R ∞ and d 0 = 0 -1 (1+n 10 -2 ), for R ∞ and R 0 the quasi-stationary reflection coefficients Instead of a complex distance of the second image the value is selected as d im =| 0 -1 (1+n 10 -2 )|, based on numerical experiments.

Antenna model of a lightning discharge at a lossy ground
The simplest approximation of a LD channel is a vertical transmitting antenna, with an excitation at its base, positioned at a lossy ground surface as in Fig. 15. For the application of an electromagnetic model and equations of antenna theory, it is necessary to divide this vertical antenna into segments of the length not greater than approximately half of the wavelength corresponding to the frequency for which the analysis is done. Vertical rod antenna with an excitation by an ideal Dirac's voltage source in the channel-base, having voltage u(t)=U (t) and frequency f, is presented in Fig. 15. The total height h of the antenna models the height of a LD channel, which may be as high as several thousands of meters in natural conditions. In such a model is assumed that the antenna is of a circular cross section with the radius a, so that a<< 0 , where 0 is the wave-length in the air corresponding to the frequency f. The corresponding angular frequency is =2πf. The antenna is divided into N segments, so that h=l 1 + l 2 +...+l N , and the length of each segment is l k >>a k , whereas radius a k << 0 . It can be simply taken that a k =a for all segments, for k=1, 2,..., N. The division into segments is necessary because of the wide range of frequencies of interest and for some of those the antenna should be divided into hundreds of segments. The lossy ground is treated as homogeneous, linear and isotropic half-space of electrical parameters: specific conductivity σ 1 , electric permittivity 1 = 0 r1 , and magnetic permeability 1 = 0 . Other parameters of two half-spaces are defined as: the complex conductivity σ i =σ i +j i , complex propagation constant i =(j i σ i ) 1/2 , for i=0 denoting air and i=1 denoting ground of the relative complex permittivity r1 = r1 -j i1 = r1 -j60σ 1 0 , and ground to the air refraction index n 10 = 1 / 0 =( r1 ) 1/2 .
Polynomial approximation is used to represent the current distribution along k-th segment with axis s k ', so that in (13) and (14), φ(s)=-div Π 0 =-∂П/∂z 0 represents the scalar potential at the antenna surface, П sn (s)= П z0 (s) Hertz vector tangential component, whereas Z n '(s) is the impedance per unit length, s n the matching point, and s the local coordinate along the n-th conductor, so that 0≤s≤ l n . П z0 (s) and φ 0 (s) are potentials in the upper half-space to be calculated from: for K 0 (r 1k )=exp(-0 r 1k )/r 1k the standard potential kernel, and S 00 ( r 2k ) the Sommerfeld's integral kernel defined as LEMF calculations were also done for the antenna models of cage structures (Javor, 2003) and in the case of lightning protection rods at a lossy ground (Javor & Rancic, 2009).

Current distribution along the vertical antenna at the lossy ground
The results obtained with TIA approximation are shown in Fig. 16 for real and imaginary part of the current along the vertical antenna having height h=300m, radius a=0.05m, and resistivity per unit length R'=0.1Ω/m, at the lossy ground of relative dielectric constant r1 =2 and r1 =10, specific conductivity σ 1 =10 -1 S/m or 10 -5 S/m, for frequency f=3MHz, the number of antenna segments N=30, and the degree of polynomial approximation n k =3, for k=1,...,30. Fig. 16. Real and imaginary part of the current along the antenna for h=300m and radius a=0.05m, for different lossy ground parameters and f=3MHz

Input impedance of a vertical antenna modeling the lightning discharge channel
Input impedance, Z ul , is important for the antenna analysis in frequency domain and presents integral characteristic of the antenna structure, which also allows checking the accuracy of TIA for Sommerfeld's integral kernel. Satisfactory agreement for the polynomial degree n>2 was confirmed. It is enough to choose a polynomial degree n=2 if the length of the antenna is not greater than 0.6 0 , for N=1 segment of the antenna. The polynomial degree should not take values n>8 as the polynomial approximation of the antenna current is not appropriate for those. For σ 1 0 <10 -1 the obtained values for the input impedance/admittance are not dependent on the normalized conductivity, but approximately equal to the values of the input impedance/admittance in the case of a perfect dielectric of relative dielectric constant r1 .

Frequency domain results for lightning electric field
Results for vertical and radial electric field as functions of the normalized radial distance 0 ρ for the normalized antenna height h/ 0 =0.25, normalized radius of the antenna a/ 0 =0.0005 at the ground of specific conductivity σ 1 =0.01S/m and for different values of electric permittivity r1 as parameter, are presented in Fig. 19, and for r1 =10 and for different values of specific conductivity σ 1 as parameter in Fig. 20. Vertical and radial component of the electric field at radial distances r from 50 to 250m (0.5 0 ≤r≤2.5 0 ) from Fig. 19. Vertical and radial electric field of the quarter-wavelength monopole antenna at the lossy ground, for σ 1 =0.01S/m and r1 as parameter, as functions of normalized distance Fig. 20. Vertical and radial electric field of the quarter-wavelength monopole antenna at the lossy ground, for r1 =10 and σ 1 as parameter, as the functions of normalized distance the base of the antenna having height h=300m, circular cross section of radius a=5cm at the ground surface of parameters r1 =10 and σ 1 =0.01S/m, for the frequency f=3MHz, the polynomial degree of the current distribution approximation along each of the segments n k =3, and the number of segments 20, 30 and 50, is presented in Fig. 21. The results obtained for the electric field in the points at a height z=1.5m above the ground surface, for the radial distances 0.5 0 ≤r≤2.5 0 , differ a little from results for the field at the ground surface (z=0).

Time domain results for lightning electromagnetic field
Results for vertical electric and azimuthal magnetic field components are presented in Figs. 22-24 for different distances from the channel-base: r=500m, 5km and 100km. For CBC function with parameters I m =11kA, t m =0.5826 s, a=1.5 and b=0.02, and for NCBC function with parameters I m =11kA, t m =0.472 s, a=1.1, b 1 =0.16, c 1 =0.34, b 2 =0.0047, and c 2 =0.66, and two different decaying constants =2000m and =4500m, the results are compared to the results from (Nucci, 1990) calculated for perfectly conducting ground using the same Modified Transmission Line Model with Exponential Decay (MTLE) with the decaying constant =2000m (Nucci et al., 1990), the same return-stroke speed v=1.3 . 10 8 m/s and the channel having height H=2600m and radius a=0.05m. For the same v, H and a, for the distributed resistance R'=0.1Ω/m along the antenna, driven by a Dirac delta current source connected across a 3.25m gap, the results are presented also for the perfectly conducting ground (Shoory et al., 2005). Shoory et al., 2005, used the electromagnetic model and a similar procedure to here presented: EFIE type equation, Method of Moments (MoM), method of images for the approximate solution of Sommerfeld's integrals (another approximation), concept of current source, NEC2 computer program (Burke & Poggio, 1981) and FFT/IFFT in 8192 points, i.e. calculations for 4092 positive frequencies. It is interesting that results very similar to (Shoory et al., 2005) would be obtained with NCBC & MTLE, but for =6000m. Fig. 23. Vertical electric and azimuthal magnetic field at the ground surface (z=0) for r=5km Fig. 24. Vertical electric and azimuthal magnetic field at the ground surface (z=0) for r=100km

Conclusion
Fourier transform proved to be very successful in lightning research. It enables calculations in frequency domain which are more suitable for including lossy ground effects than in time domain. It also provides information about frequency spectra of the quantities of interest in lightning research. For antenna modeling of a lightning stroke in frequency domain, Sommerfeld's integral is calculated efficiently using Two-image approximation.
Results are in good agreement with the results from literature for various ground electrical parameters, heights of vertical dipoles above ground, in the near and far field. Based on these results, it can be concluded that the effects of lossy ground are greater on horizontal than on vertical electric field, and that specific conductivity influences more than electrical permittivity. Fourier transform application has to be further investigated in terms of optimal choice of FFT parameters in order to reduce computing time which can be important in antenna modeling of lightning discharge channels and in analysis of lightning electromagnetic fields.

Acknowledgment
This work is supported by Project Grant III44004 (2011-2014) financed by Ministry of Education and Science, Republic of Serbia.