**1. Introduction**

As a prerequisite for evaluating lightning-induced voltages in distribution networks, it is of grea<sup>t</sup> significance to calculate the electromagnetic field generated by a lightning return stroke over a lossy ground accurately and e fficiently. At present, there are two kinds of methods for the calculation: Finite-di fference time-domain method (FDTD) and Cooray–Rubinstein (CR) formula. The FDTD method [1–3] has grea<sup>t</sup> applicability when considering of the e ffect of lossy ground, such as the lightning electromagnetic field over layered ground. However, the space and time are required to be discretized in FDTD, which results in a huge amount of computation time and relatively complicated programming, especially for the absorbing boundary. In order to reduce the computation time, analytical formulae are an alternative choice. The Sommerfeld integral [4] can be used to rigorously evaluate the lightning electromagnetic field generated by a lightning return stroke over a lossy ground; however, the highly oscillatory and slow convergen<sup>t</sup> integrand make it di fficult to evaluate the integral efficiently. The Cooray–Rubinstein formula [5–8] is a good choice to overcome this di fficulty, and has become a widely used method. In the past years, most studies have mainly focused on e fficient evaluation of the CR formula in the time domain [9–18]. Essentially, the main task of these methods is how to describe the integral kernel function of the CR formula in the time domain and accelerate the calculation of the convolution.

In practice, the implementation of the CR formula to evaluate the lightning electromagnetic field over a lossy ground requires a two-step procedure. The first step involves calculating the lightning electromagnetic field over a perfectly conducting ground, and the second step is to evaluate the CR formula. Actually, compared with the second step, the calculation of the ideal lightning

electromagnetic field (*Er*,*ideal* and *<sup>H</sup>*ϕ,*ideal*) accounted for most of the computation time. Therefore, increasing the e fficiency of the first step is more conducive to improving the overall computational efficiency. The relevant studies about this tend to be neglected. The formulae for evaluating *Er*,*ideal* and *<sup>H</sup>*ϕ,*ideal* have been established based on the dipole method, and they are composed by integrals with respect to the lightning channel. The common method to evaluate these integrals are mainly based on the numerical integration, by means of a discretization of the lightning channel. However, only a su fficiently small discretization step is essential to ge<sup>t</sup> an accurate result, which leads to a relatively large number of calculations and results in a lengthy computation time. Besides, the programming is relatively complicated because the propagation of the lightning current along the channel must be considered.

In order to increase the e fficiency and simplify the programming, an improved method is proposed in this paper. Firstly, regarding the return stroke current in the lightning channel as a series of current sources distributed along the channel, the original time-varying system can be transformed into a time-invariant system. Then, the calculation formulae are rewritten in the frequency domain by means of the Fourier transform, which simplifies the integral and di fferential operations about time in the original formula, so that the formulae are reduced from having two variables (*z*- , *t*) in the time domain to having only one variable (*z*-) in the frequency domain. Finally, a series of analytical formulae according to the integrals are derived, which e ffectively simplify the calculation procedure. Additionally, compared with the conventional method, the e fficiency of the proposed method can be increased, which can be examined by the simulation example and discussion.

#### **2. Method for Calculating the Lightning Electromagnetic Field over Perfectly Conducting Ground**

### *2.1. Review of the Lightning Return Stroke Model*

In order to calculate the lightning electromagnetic field, the lightning return stroke model must be established firstly. Basically, there are four kinds of lightning return stroke models: The gas dynamic or physical mode, the electromagnetic model, the distributed-circuit model, and the engineering model [19–22]. The engineering model is adopted here due to its wide utilization. There are several engineering models, such as the Bruce–Golde model (BG), transmission line model (TL), traveling current source model (TCS), modified transmission line exponential decay model (MTLE), and modified transmission line linear model (MTLL) [23]. In this paper, the MTLE model is adopted, which can be illustrated as Figure 1 and described by:

$$\begin{cases} \ i(z',t) = i(0, t - \frac{z'}{v}) \exp(-\frac{z'}{a}), & t \ge \frac{z'}{v} \\\ i(z',t) = 0, & t < \frac{\overline{z'}}{v} \end{cases} \tag{1}$$

where: *v*—the return stroke velocity along the lightning channel, α—the decaying constant, *<sup>i</sup>*(0,*t* − *z*- /*v*)—the delay of the channel base current, and exp(−*z*- /α)—the attenuation of the channel base current.

**Figure 1.** Lightning return stroke model.

*Appl. Sci.* **2020**, *10*, 4263

As for the channel base current, several functions are available, such as Rubinstein and Uman function, Bruce and Golde function, and Heidler function [24]. In this paper, the channel base current is represented as a sum of two Heidler functions, as given in Equation (2), and the typical parameters are listed in Table 1:

$$i(0,t) = \left[\frac{i\_{01}}{\eta\_1} \frac{\left(t/\tau\_{11}\right)^{n\_1}}{1 + \left(t/\tau\_{11}\right)^{n\_1}} \exp\left(-\frac{t}{\tau\_{12}}\right) + \frac{i\_{02}}{\eta\_2} \frac{\left(t/\tau\_{21}\right)^{n\_2}}{1 + \left(t/\tau\_{21}\right)^{n\_2}} \exp\left(-\frac{t}{\tau\_{22}}\right)\right].\tag{2}$$

where:

 $\eta\_1 = \exp\left(-\frac{\tau\_{11}}{\tau\_{12}} \left(n\_1 \frac{\tau\_{12}}{\tau\_{11}}\right)^{\frac{1}{n\_1}}\right)$  $\text{and } \eta\_2 = \exp\left(-\frac{\tau\_{21}}{\tau\_{22}} \left(n\_1 \frac{\tau\_{22}}{\tau\_{21}}\right)^{\frac{1}{n\_2}}\right).$ 


#### **Table 1.** The basic parameters for the channel base current function.

#### *2.2. Fundamental Formulation of Lightning Electromagnetic Field over Perfectly Conducting Ground*

The widely used model of the lightning electromagnetic field over a perfect conducting ground can be illustrated in Figure 2 [25]. The lightning channel can be regarded as a combination of an above-ground lightning channel and its mirror in the free space. Regarding the source as a vertical electric dipole of current *<sup>i</sup>*(*z*-,*<sup>t</sup>*)*dz*- located at *z*-, the formulas of the electromagnetic field at an observation point generated by the lightning channel can be obtained and shown as:

$$\begin{cases} \begin{aligned} &E\_{\text{ideal}} = \frac{1}{4\pi\varepsilon\_{0}} \Bigg{ \left\{ \begin{aligned} &\frac{3r(z-z')}{R\_{\star}^{3}} \int\_{0}^{h} i(z',\tau-R\_{\star}\prime/c) d\tau \\ &+\frac{3r(z-z')}{cR\_{\star}^{3}} i(z',t-R\_{\star}\prime/c) \\ &+\frac{r(z-z')}{c^{2}R\_{\star}^{3}} \frac{\partial i(z',\tau-R\_{\star}\prime/c)}{\partial t} \end{aligned} \right\} dz' + \begin{bmatrix} &\frac{3r(z+z')}{cR\_{\star}^{3}} \int\_{0}^{h} i(z',\tau-R\_{-\star}\prime/c) d\tau \\ &+\frac{r(z+z')}{cR\_{\star}^{3}} i(z',t-R\_{-\star}\prime/c) \\ &+\frac{r(z+z')}{c^{2}R\_{\star}^{3}} \frac{\partial i(z',t-R\_{-\star}\prime/c)}{\partial t} \end{aligned} \end{cases} \\\begin{aligned} &H\_{\text{ideal}} = \frac{1}{4\pi\varepsilon} \Bigg{pmatrix} \left\{ \begin{aligned} &i(0,0) = \frac{r}{R\_{\star}} \int\_{0}^{h} i(z',0-R\_{-\star}\prime/c) \\ &+\frac{r(z+z')}{cR\_{\star}^{3}} \frac{\partial i(z',t-R\_{-\star}\prime/c)}{\partial t} \end{aligned} \right\} \end{aligned}$$

where: *Eideal*(ϕ, *r*, *z*, *t*)—ideal horizontal electric field at observation point *P*(*<sup>r</sup>*, *z*), *Hideal*(ϕ, *r*, *z*, *t*)—ideal tangential magnetic field at observation point *P*(*<sup>r</sup>*, *z*), *h*—the height of the lightning channel, ε0—vacuum permittivity, *r*—horizontal distance from the lightning channel to point *P*(*<sup>r</sup>*, *z*), *Rz*- and *R*−*z*-—distance from the dipole to point *P*(*<sup>r</sup>*, *z*), *z* and *z*-—the z-coordinate of point *P*(*<sup>r</sup>*, *z*) and the dipole, respectively, and *c*—speed of light in vacuum.

**Figure 2.** Model for calculating the lightning electromagnetic field.

The calculation of Equation (3) is commonly based on a numerical integration due to the lack of the analytical formulae of Equation (3), which is called the conventional method in this paper. As can be seen, the integrands in Equation (3) are functions of two variables (*t* and *z*-), and the propagation of the current along the channel must be considered when implementing the numerical integration. Therefore, it is somewhat complicated and not easily programmed. Experience indicates that a su fficiently small computational step *dz*- is essential to ge<sup>t</sup> an accurate result, but this leads to a large number of calculations, resulting in a lengthy calculation time.

Therefore, we doubt whether or not some of the integral can be implemented in an analytical way; if so, the e fficiency can be improved, especially for the evaluation of the lighting-induced voltages of a distributed network, in which lightning electromagnetic fields with a grea<sup>t</sup> amount of observation points are required to be calculated. In the following sections, an e fficient method is achieved, which can calculate most of the integrals in Equation (3) analytically.

#### *2.3. Proposed Method for Calculating the Lightning Electromagnetic Field over Perfectly Conducting Ground*

In order to simplify Equation (3), firstly, we transform the system from time varying to time invariant. The process of the channel base current propagating along the lightning channel is time varying, but if regarded it as a current model with *i*(*z*- , *t*) = *ibase*(*t-z*- /*v*) distribution at all points along the channel direction, the system can be transformed into a time-invariant system. Then, we can transform the time-invariant system into the frequency domain based on the Fourier transform. Finally, the analytical formulae for most of the integral in Equation (3) can be obtained by means of some derivations.

Substituting Equations (1) and (2) into (3) and representing the equations in the frequency domain allows the formulas to be rearranged as:

$$\begin{split} E\_{\text{ideal}} &= \frac{1}{4\pi\varepsilon\_{0}} \Bigg{{}^{0}} \Bigg{{}^{0}} \Bigg{{}^{0} \Bigg{{}^{1}} \Bigg{{}^{3r(z-\omega\_{\star}^{\prime})}I(0,j\omega)} \Bigg{{}^{-2r(\frac{\omega}{\overline{\sigma}}+\frac{\omega}{\overline{\omega}})}e^{-\omega^{\prime}\left(\frac{\omega}{\overline{\sigma}}+\frac{\omega}{\overline{\omega}}\right)}e^{-\beta\overline{R}\_{z'}}}{{}^{3r(z-\omega\_{\star}^{\prime})}I(0,j\omega)} \Bigg{{}^{-2r(\frac{\omega}{\overline{\sigma}}+\frac{\omega}{\overline{\omega}})}e^{-\beta\overline{R}\_{z'}}} \\ & \underbrace{\frac{1}{-j\omega\frac{r(z-\omega^{\prime})}{c^{2}R\_{y'}^{2}}I(0,j\omega)e^{-\omega^{\prime}\left(\frac{\omega}{\overline{\sigma}}+\frac{\omega}{\overline{\omega}}\right)}e^{-\beta\overline{R}\_{z'}}}}\_{E\_{\text{shear}}} \underbrace{\Bigg{{}^{-2r(\frac{\omega}{\overline{\sigma}}+\frac{\omega}{\overline{\omega}})}I(0,j\omega)}\_{-j\omega\frac{r(z+\omega^{\prime})}{c^{2}R\_{z'}^{2}}I(0,j\omega)}\Bigg{{}^{-2r(\frac{\omega}{\overline{\sigma}}+\frac{\omega}{\overline{\omega}})}e^{-\beta\overline{R}\_{z'}}}\_{E\_{\text{manif}}}}\_{E\_{\text{manif}}} \\\end{split}$$

As can be seen in Equation (4), *Eabove* and *Eunder* have the same form, so only the derivation for *Eabove* is provided. The above-ground lightning channel-generated horizontal electric field can be divided into three parts, i.e., *EA*, *EB*, and *EC*, as shown in Equation (5):

$$E\_{\text{above}} = \frac{\mathbb{I}^{(0,j\_0)\_T}}{4\pi\epsilon\_0} \left\{ \underbrace{\frac{3}{4\nu} \int\_0^{\mathbb{R}} \underbrace{\frac{(z-z')}{R\_{z'}^{3}} \cdot e^{-z'\left(\frac{\omega}{\overline{\sigma}} + \frac{1}{\overline{\sigma}}\right) - \beta R\_{z'}} dz'}\_{\begin{subarray}{c} \mathbb{E}\_A \\ + \ j\omega \le \frac{1}{\overline{\sigma}} \end{subarray}} + \underbrace{\frac{3}{c} \int\_0^{\mathbb{R}} \underbrace{\frac{(z-z')}{R\_{z'}^{3}} \cdot e^{-z'\left(\frac{\omega}{\overline{\sigma}} + \frac{1}{\overline{\sigma}}\right) - e^{-\beta R\_{z'}} dz'} dz'}\_{\begin{subarray}{c} \mathbb{E}\_B \\ \mathbb{E}\_C \end{subarray}} \right\}, \begin{subarray}{c} \left\{ k = \frac{\omega}{\overline{\sigma}} \right\} \end{subarray} \right\}, \begin{subarray}{c} \left\{ k = \frac{\omega}{\overline{\sigma}} \right\} \end{subarray} \right\}, \begin{subarray}{c} \left\{ k = \frac{\omega}{\overline{\sigma}} \right\} \end{subarray} $$

By means of (*<sup>z</sup>*−*z*-) *R*<sup>5</sup> *z*- *dz*- = 1 3*R*<sup>3</sup> *z*- and *uvdz* = *uv* − *uvdz*, *EA* can be rewritten as:

$$E\_A = \underbrace{\left| \frac{1}{j\omega} \cdot \frac{1}{R\_{z'}^3} e^{-z'(\frac{j\omega}{\omega} + \frac{1}{\omega})} e^{-j\beta R\_{z'}} \right|^k}\_{E\_{A1}} + \underbrace{\frac{1}{j\omega} \cdot \left( \frac{j\omega}{\omega} + \frac{1}{\pi} \right) \int\_0^h \frac{1}{R\_{z'}^3} e^{-z'(\frac{j\omega}{\omega} + \frac{1}{\omega})} e^{-j\beta R\_{z'}} dz' + \underbrace{\left( -\frac{1}{\varepsilon} \int\_0^h \frac{(z - z')}{R\_{z'}^4} e^{-z'(\frac{j\omega}{\omega} + \frac{1}{\omega})} e^{-j\beta R\_{z'}} dz' \right)}\_{E\_{A2}} \tag{6}$$

$$E\_{A2} = \left\{ \underbrace{\frac{1}{j\omega} \cdot \left(\frac{j\omega}{\upsilon} + \frac{1}{a}\right) \frac{1}{R\_{x'}} \frac{1}{z - z'} \varepsilon^{-z'\left(\frac{j\omega}{\upsilon} + \frac{1}{a}\right)} e^{-j k R\_{x'}}}\_{E\_{A21}} \right\}\_{E\_{A21}} + \underbrace{\left(-\frac{1}{j\omega} \cdot \left(\frac{j\omega}{\upsilon} + \frac{1}{a}\right)\right) \int\_{0}^{k\_1} \frac{1}{R\_{x'}} \frac{1}{(z - z')^2} \varepsilon^{-z'\left(\frac{j\omega}{\upsilon} + \frac{1}{a}\right)} e^{-j k R\_{x'}} dz'}\_{E\_{A22}} \right\}\_{E\_{A23}} \tag{7}$$
 
$$\left\{ \underbrace{-\frac{1}{j\omega} \cdot \left(\frac{j\omega}{\upsilon} + \frac{1}{a}\right) \left(\frac{j\omega}{\upsilon} + \frac{1}{a}\right) \int\_{0}^{k\_1} \frac{1}{R\_{x'}} \frac{1}{z - z'} \varepsilon^{-z'\left(\frac{j\omega}{\upsilon} + \frac{1}{a}\right)} e^{-j k R\_{x'}} dz'}\_{E\_{A23}} + \underbrace{\left(-\frac{1}{\varepsilon} \cdot \left(\frac{j\omega}{\upsilon} + \frac{1}{a}\right)\right) \int\_{0}^{k\_1} \frac{1}{R\_{x'}^2} \varepsilon^{-z'\left(\frac{j\omega}{\upsilon} + \frac{1}{a}\right)} e^{-j k R\_{x'}} dz'}\_{E\_{A24}}\right) \right\}\_{E\_{A24}} \tag{7}$$

$$E\_{A22} = \left\{ \underbrace{-\frac{1}{j\omega} \cdot \left(\frac{j\omega}{\upsilon} + \frac{1}{\alpha}\right) \frac{R\_{z'}}{r^2(z-z')} e^{-z'\left(\frac{j\omega}{\upsilon} + \frac{1}{\alpha}\right)} e^{-j\beta R\_{z'}}}\_{+} + \underbrace{\frac{1}{j\omega} \cdot \left(\frac{j\omega}{\upsilon} + \frac{1}{\alpha}\right) \hbar \frac{1}{r^2} \int\_{0}^{\hbar} e^{-z'\left(\frac{j\omega}{\upsilon} + \frac{1}{\alpha}\right)} e^{-j\beta R\_{z'}} dz'}\_{E\_{A22}} \right. \tag{8}$$
 
$$\left\{ -\underbrace{\frac{1}{j\omega} \cdot \left(\frac{j\omega}{\upsilon} + \frac{1}{\alpha}\right) \left(\frac{j\omega}{\upsilon} + \frac{1}{\alpha}\right)}\_{E\_{A22}} \int\_{0}^{\hbar} \frac{R\_{z'}}{r^2(z-z')} e^{-z'\left(\frac{j\omega}{\upsilon} + \frac{1}{\alpha}\right)} e^{-j\beta R\_{z'}} dz'}\_{E\_{A22}} \right\}\_{(8)} \tag{9}$$

$$E\_{A3} + E\_B = \left\{ \underbrace{\frac{2}{c} \frac{1}{2R\_{z'}^2} e^{-z'(\frac{j\omega}{v} + \frac{1}{a})} e^{-j k R\_{z'}}}\_{=\left\{ \begin{array}{c} \frac{1}{c} \left( \frac{j\omega}{v} + \frac{1}{a} \right) \end{array} \right\}\_0^h + \underbrace{\frac{1}{c} \left( \frac{j\omega}{v} + \frac{1}{a} \right) \int\_0^h \frac{1}{R\_{z'}^2} e^{-z'(\frac{j\omega}{v} + \frac{1}{a})} e^{-j k R\_{z'}} dz'}\_{E\_{B2}} \right\}\_0^h. \tag{9}$$

$$E\_{A222} + E\_{A223} + E\_{A23} = \left\{ \underbrace{-\frac{1}{j\omega r^2} \left( \frac{j\omega}{\upsilon} + \frac{1}{\alpha} \right) \left( \frac{j\alpha}{\upsilon} + \frac{1}{\alpha} \right) \frac{e^{-jkR\_{\varphi}}}{jk} e^{-z'\left(\frac{j\omega}{\upsilon} + \frac{1}{\alpha}\right)}}\_{E\_{A231}} \right\}\_{E\_{A231}} \tag{10}$$
 
$$\left\{ \underbrace{-\frac{1}{j\omega r^2} \left( \frac{j\omega}{\upsilon} + \frac{1}{\alpha} \right) \left[ jk - \frac{1}{jk} \left( \frac{j\omega}{\upsilon} + \frac{1}{\alpha} \right)^2 \right]}\_{E\_{A241}} \right\}\_{E\_{A241}} \cdot \tag{11}$$

Considering *EB2* + *EA24* = 0 and *EB3* + *EC* = 0, Equation (5) can be rewritten as:

*Eabove* = *<sup>I</sup>*(0,*j*ω)*r* 4πε0 ⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ 1*j*ω · 1*R*<sup>3</sup>*z*- *e*<sup>−</sup>*z*-( *j*ω*v* + 1α )*e*<sup>−</sup>*jkRz h*0 -------------------------------------------------------- *EA*1 + 2*c* 1 2*R*<sup>2</sup>*z*- *e*<sup>−</sup>*z*-( *j*ω*v* + 1α )*e*<sup>−</sup>*jkRz h*0 ---------------------------------------------------- *EB*1 + 1*j*ω · *j*ω*v* + 1α 1*Rz*- 1 *z* − *ze*<sup>−</sup>*z*-( *j*ω*v* + 1α )*e*<sup>−</sup>*jkRz h*0 ---------------------------------------------------------------------------------------------------- *EA*21 +(− 1*j*ω · *j*ω*v* + 1α *Rz*- *r*<sup>2</sup>(*z* − *z*-)*e*−*z*-( *j*ω*v* + 1α )*e*<sup>−</sup>*jkRz h*0) -------------------------------------------------------------------------------------------------------------- *EA*221 +(− 1*j*ω*r*<sup>2</sup> · *j*ω*v* + 1α *j*ω*v* + 1α *e*<sup>−</sup>*jkRz*- *jk e*<sup>−</sup>*z*-( *j*ω*v* + 1α ) *h*0) ---------------------------------------------------------------------------------------------------------------------- *EA*231 <sup>+</sup>⎡⎢⎢⎢⎢⎣*<sup>j</sup>*ω1*v*<sup>1</sup>*c* − *cv*2 + 1α1*c* − 3*cv*2 − 1*j*ω 3*c* α2*v* − 1 (*j*ω)<sup>2</sup> *c*α<sup>3</sup> ⎤⎥⎥⎥⎥⎦ · 1*r*2 *h*0 *e*<sup>−</sup>*z*-( *j*ω*v* + 1α )*e*<sup>−</sup>*jkRz*- *dz*- ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- *EA*241 ⎫⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭ (11)

and it can be transformed into the time domain, i.e.:

$$E\_{\text{above}} = E\_{\text{arithmetic}} + E\_{\text{integral}}.\tag{12}$$

where:

*Earithmetic* = 1 4πε0 ⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ *rc* 1*R*2*h* + *r vRh* 1*<sup>z</sup>*−*h* − *Rh vr*(*<sup>z</sup>*−*h*) − *crv*2 *i* 0, *t* − *hv* − *Rhc e*<sup>−</sup> *h*α −*rc* 1*R*20 + *r vR*0 1*z* − *R*0 *vrz* − *crv*2 *i* 0, *t* − *R*0*c* + *rR*3*h* + 1α *rRh* 1*<sup>z</sup>*−*h* − 1α *Rh <sup>r</sup>*(*<sup>z</sup>*−*z*-) − 2*c* α*vr <sup>i</sup>*−<sup>1</sup> 0, *t* − *hv* − *Rhc e*<sup>−</sup> *h*α − *rR*30 + 1α *rR*0 1*z* − 1α *R*0 *rz* − 2*c* α*vr<sup>i</sup>*−<sup>1</sup> 0, *t* − *R*0*c* − *c r*α<sup>2</sup> *<sup>i</sup>*−<sup>2</sup> 0, *t* − *hv* − *Rhc e*<sup>−</sup> *h*α + *c r*α<sup>2</sup> *<sup>i</sup>*−<sup>2</sup> 0, *t* − *R*0*c* ⎫⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭ (13) 1 1− ·1·1 1− ·1·

$$E\_{\rm integral} = \begin{pmatrix} \frac{1}{v} \left( \frac{1}{c} - \frac{c}{v^2} \right) \cdot \frac{1}{r} \cdot \frac{d}{dt} i\_{\rm integral}(t) + \frac{1}{a} \left( \frac{1}{c} - \frac{3c}{v^2} \right) \cdot \frac{1}{r} \cdot i\_{\rm integral}(t) \\\ -\frac{3c}{a^2 v} \cdot \frac{1}{r} \cdot i\_{\rm integral}^{-1}(t) - \frac{c}{a^3} \cdot \frac{1}{r} \cdot i\_{\rm integral}^{-2}(t) \end{pmatrix} \tag{14}$$

where: *i*(0, *t*)—the channel base current, *i* 0, *t* − *hv* − *Rhc* and *i* 0, *t* − *R*0*c* —the channel base current with delay, *iintegral*(*t*) = 1 4πε0 *h*0 *i*(0, *t* − *zv* − *Rz*- *c* )*e*<sup>−</sup>*z*-( 1α )*dz*-—the integral component associated with the channel base current, *Rh*—the distance between the highest point of the lightning channel and point *P*(*<sup>r</sup>*, *z*), and *R*0—the distance between the lightning strike point and point *P*(*<sup>r</sup>*, *z*).

By means of the similar derivation, the formula for *Eunder* can also be achieved, only by replacing *<sup>z</sup>*–*z*- and *Rz*- by *z* + *z*- and *R*–*z*-, respectively. Finally, the horizontal electric field *Eideal* can be obtained by adding them together.

As can be seen from Equations (12)–(14), almost all of the integral can be evaluated in an analytical way, except for *Eintegral*. Moreover, the numerical calculation of *Eintegral* is very simple, because it only contains the integral of the current along the channel and its derivatives. Generally, comparing Equations (12)–(14) with Equation (3), it can be seen that the proposed method equates the integral operation with a simple arithmetic operation and a simple integral. The arithmetic operation involves fewer steps, requires less data, and produces more accurate results. The integral is much simpler than that of the conventional method. Thus, doing this makes the proposed method more efficient and easily programmed.
