An Approximation to the Heidler Function With an Analytical Integral

An analytical integral of a lightning current model is required in the analysis of lightning events including the induced effects and frequency analysis of lightning strikes. Previous work in this area has produced very specific forms of the Heidler function that are used to represent lightning current waveshapes. This work, however, focuses on a generic approximation, with parameters that can be modified to produce any required lightning current waveshape. The approximation function proposed has an analytical solution to the integral allowing for a true representation of Maxwell's equations for electromagnetic fields. Moreover, the approximation is designed in the Laplace domain allowing for a complete analytical frequency domain analysis. The function has parameters that can be changed to obtain different waveshapes (such as 10/350 and 0.25/100). The characteristics of the approximation are compared with those of the Heidler function to confirm that it is applicable for use with in IEC 62305-1 lightning protection standard applications. The maximum error in amplitude for the first and subsequent lightning stroke currents is less than 1.5%. The maximum error in the derivative is greater but still less than 8%. This is within the parameters defined in the standard.


An Approximation to the Heidler Function With an Analytical Integral
Brett R. Terespolsky, Hugh G. P. Hunt , and Ken J. Nixon , Member, IEEE Abstract-An analytical integral of a lightning current model is required in the analysis of lightning events including the induced effects and frequency analysis of lightning strikes.Previous work in this area has produced very specific forms of the Heidler function that are used to represent lightning current waveshapes.This work, however, focuses on a generic approximation, with parameters that can be modified to produce any required lightning current waveshape.The approximation function proposed has an analytical solution to the integral allowing for a true representation of Maxwell's equations for electromagnetic fields.Moreover, the approximation is designed in the Laplace domain allowing for a complete analytical frequency domain analysis.The function has parameters that can be changed to obtain different waveshapes (such as 10/350 and 0.25/100).The characteristics of the approximation are compared with those of the Heidler function to confirm that it is applicable for use with in IEC 62305-1 lightning protection standard applications.The maximum error in amplitude for the first and subsequent lightning stroke currents is less than 1.5%.The maximum error in the derivative is greater but still less than 8%.This is within the parameters defined in the standard.Index Terms-Current, frequency domain, Heidler function, lightning protection, waveshape.

I. INTRODUCTION
M ATHEMATICAL lightning current models are used in many areas of research ranging from the design of lightning protection systems (LPSs), to the understanding of electric and magnetic fields associated with lightning discharges [1], [2], [3].Lightning current models are typically used as design tools and for further research into the understanding of the effects of a lightning strike.A key use is as the channel base current for return stroke modeling [4], [5], [6], [7], [8], [9].The Heidler function, which is a standard lightning current model, cannot be integrated analytically making it impossible to utilize Maxwell's equations in analyzing lightning events.It is also not possible to represent the Heidler function in the frequency domain without using approximations [7].This research presents a suitable The authors are with the Johannesburg Lightning Research Laboratory (JLRL) School of Electrical and Information Engineering, University of the Witwatersrand, Johannesburg 2001, South Africa (e-mail: hugh.hunt@wits.ac.za.).
Color versions of one or more figures in this article are available at https://doi.org/10.1109/TEMC.2024.3372642.
Digital Object Identifier 10.1109/TEMC.2024.3372642replacement to the Heidler function for situations where an analytical integral or a continuous Fourier transform is required.This replacement falls under the category of an approximation to the Heidler function and as such, it can be tailored to any waveshape required, as with the Heidler, by varying the parameters.An investigation is carried out to determine the accuracy of the approximation.This is done by using computer simulations of the approximation and the results are compared to the Heidler to determine the viability of this approximation in the design of LPSs.All the results are based on the waveshapes defined in the IEC 62305-1 so that there is a known control.
The contribution of this study is to develop a function that approximates the Heidler function in the time domain.This approximation has the advantage of having an analytical integral and taking the same form as that of the Heidler function; only the parts of the Heidler function that cannot be integrated are replaced in producing this approximation.This approximation is easy to use and the parameters for creating different lightning current waveshapes are determined easily.In a system design or simulation, it would be trivial to replace the Heidler function with this approximation.Most importantly, it can be utilized in Maxwell's equations or finding the Fourier transform (frequency components).

II. BACKGROUND
The application of lightning return stroke models (RSMs), the waveshapes used, and previous approximations to the Heidler function are discussed in the following sections.

A. Application of Lightning RSMs
RSMs can be characterized into four primary classes, namely gas dynamic models, electromagnetic (EM) models, distributed circuit models and engineering models [10], [11], [12].The engineering models are the primary models used when considering induced effects on power lines.Fig. 1 shows the physical [see Fig. 1(a)] and RSMs [see Fig. 1(b)] for a lightning strike in close proximity to a transmission line.The lightning channel is modeled as a monopole with a perfectly conductive ground and a vertical current path [10], [13], [14].
The induced voltages can be calculated by integrating the different components of the EM fields as outlined by Agrawal [15].The lightning electromagnetic pulse (LEMP) equations need to be determined before calculating the EM fields.The equations for the EM fields can be seen in ( 1) and (2), respectively.These are defined by Uman [16] as the LEMP equations.All quantities are as appear in [16] These equations utilize the RSM, i(z, t).This expression for the lightning current along a path can be defined by different models as in [2], [10], [17].All the models hold the form shown in RSMs where u(t) is the Heaviside function, P (z ) is the height-dependent current attenuation factor and i(0, t) is the lightning channel base current, which as the name suggests is the current as measured from the base of the object that is struck [2], [10], [18] In short, there are following two steps of integration required to calculate induced currents on transmission lines [10], [17], [19]. 1) Integrate some model of the lightning channel base current to obtain the EM fields.
2) Integrate the EM fields using the equations defined by Agrawal to obtain the induced currents.The induced effects on a nearby transmission line can be calculated by modeling the lightning event as in Fig. 1(b) and choosing an appropriate lightning channel base current.From the above-mentioned equations, it can be seen that at least the second integral of this channel base current is required.Using an accurate lightning channel base current leads to very complex models that require long computer simulations due to the numerical integration required [10].

B. Current Waveshape Models
The IEC 62305-1 lightning protection standard details the Heidler function with a specific configuration as a standardized waveshape for lightning current simulations.There are several lightning current models defined in the literature with the two most commonly used being the Heidler and the double exponential.There are also applications that use a combination of the two [17], [18], [20].This section gives some background into these two equations and their properties.
1) Double Exponential Model: The double exponential function is defined in the following: where I 0 is the peak current [A], A is the peak current correction, α is the decay time constant [1/s], and β is the rise time constant The double exponential model is often used in place of the Heidler function because it can be integrated.However, the waveshape produced by the double exponential equation is not physically realistic because the maximum current steepness occurs at time t = 0 [2], [21], [22], [23].Moreover, this function does not allow for waveshapes that comply with Table A.1 in the IEC 62305-1 standard.For instance, a subsequent negative stroke for lightning protection Level I would have a maximum current steepness of about 545 kA/μs, which is far greater than the maximum value outlined in the standard of less than 200 kA/μs [24].Another disadvantage to using this equation is that it is not easily adjustable, the parameters are not easily obtained for different waveshapes [18].
2) Heidler Function: To avoid the disadvantages of the double exponential equation, the Heidler function is used.This can be seen in the following: where I 0 is the peak current [A], η is the peak current correction, τ 1 is the rise time constant [s], τ 2 is the decay time constant [s], and n h is the Heidler steepness factor.This equation more realistically approximates a lightning return stroke and the peak current correction can be calculated using the following [2], [18], [21], [23], [25]: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
The first derivative of the Heidler function (i h ) is given by The IEC 62305-1 standard makes use of a specialised form of the Heidler function where n h = 10.Using this form of the equation, there are values for the other parameters in the standard that can be used to simulate both first and subsequent strokes [1], [24].The disadvantage of the Heidler function is that there is no analytical integral and hence no continuous analytical expression for the Heidler function in the frequency domain [2], [7], [21].This creates problems when performing analyses, such as those mentioned in Section II-A previously.

C. Heidler Function Approximations
Zhang et al. [2] developed the pulse function.This function produces a maximum error of 0.5%, however, this function is a modified form of the double exponential function and requires complex methods to determine the parameters used in the equation.Heidler et al. [22] approximated the Heidler function utilized in the IEC 62305-1 standard.This approximation has an analytical integral but the equation is specific to the subsequent negative stroke.There is no general form of this approximation and so this equation cannot be easily manipulated.Their study also details several other approximations to the Heidler function, but all with other specific applications.They also all have the disadvantage that they cannot be integrated analytically.Delfino et al. [23] conducted a study in which they develop a Prony series approximation to the Heidler function.The mathematics required for this approximation is complex, limiting its use in engineering applications.Each scenario is different as there is no truly generalised form and the number of terms used affects the error associated with the approximation.This would be difficult to use for engineering applications.
Javor et al. [6], [18] developed a new approximation called the NCBC.This function is a piecewise function which means that there are discontinuities when taking the analytical derivative with step functions.The approximation also contains incomplete gamma functions, which are defined as integrals and hence further complicate the mathematics [26].There is no analytical Fourier transform to this approximation [27].Vujevic et al. [28], [29] defined their version of the exponential approximation that is fixed for n h = 10.The mathematics presented in their study is particularly complex and the authors introduce several new unknown parameters.The approximation is not as intuitive as the Heidler function and the frequency response presented is not as expected at higher frequencies; it does not roll off as it does in the standard.
Mestriner et al. [7] utilized an approximation of the Fourier transform of the Heidler function in the development of a new channel-based lightning current formula but this is an approximation using numerical methods or DFT or similar, and is not continuous.
All of the approximations discussed in the previous paragraphs have one or more disadvantages or limitations.Generally, the mathematics are complex making it difficult to recommend them as substitutes for the Heidler function in engineering applications.The Terespolsky approximation presented here, provides an approximation with an analytical integral, for use in engineering applications that does not have these limitations.

III. TERESPOLSKY APPROXIMATION
The Terespolsky approximation of the Heidler function is given in (8).The detailed derivation is described in Section IV where and n a is the approximation steepness factor.
Obtaining the integral of the approximation is trivial and can be seen in the following, where C is some arbitrary constant of integration because it is an indefinite integral: The following equation shows the analytical derivative of the approximation: The derivation and properties of the Terespolsky approximation, which requires the decomposition of the Heidler function, are detailed in the next section.

IV. DERIVATION OF THE APPROXIMATION
Insight into the solution is gained by decomposing the Heidler function and identifying the components that prevent it from being integrated analytically [30], [31].The approximation function has certain criteria: first, it must approximate the Heidler function in the time domain (within certain error limits); second, it must have an analytical solution to its integral; finally, it must take the same form as the Heidler function so that it can be used interchangeably for different applications.With these criteria in mind the approximation function is developed.

A. Decomposing the Heidler Function
An example of the Heidler function can be seen in Fig. 2.This function is defined by (5) in Section II-B2.
The shorthand version of the equation is Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.where and Equations ( 12) and ( 13) are the rise and fall components of the Heidler function respectively.The rise component cannot be analytically integrated, preventing an analytical transform of the Heidler function into the frequency domain.These equations will be analyzed in detail in Sections IV-A1 and IV-A2 and an approximation that overcomes this limitation of the Heidler function is developed.

1) Heidler Rise Function:
The rise time part of the Heidler function is defined by (12) and is shown in Fig. 3.This function takes the form of an S-curve, and can be readily modified to represent any lightning waveshape rise time.This can be achieved by varying n h (the steepness factor) and τ 1 (the rise time constant).Therefore this function meets the criteria that it can approximate any lightning waveform.However, this function cannot be integrated and hence cannot be transformed into the frequency domain.In order to solve this limitation another S-curve must be developed that approximates this one and which can also be integrated.
There are numerous forms of the S-curve, such as those given in (14a)-(14f) These equations do not allow for an easy manipulation of parameters to tailor the waveshape to a specific rise time and steepness.Many of these functions can also not be integrated, which would not solve the limitation of the Heidler function.Therefore, an alternative is required that can easily be modified.
2) Heidler Fall Function: The part of the Heidler function that controls the decay time and shape is given by ( 13) and shown in Fig. 4.This function meets all the criteria outlined previously-it can easily be altered to change the decay time and it can be trivially integrated.Moreover, because it is expressed Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
as e f (t), the product of it and another function can be integrated using integration by parts.This function is simply a complex shift of the signal in the frequency domain because of the rule of Laplace transforms shown in the following [32], [33]: Therefore, there is no need to redefine the decay part of the equation in any way and the approximation function can still be defined as follows: where I 0 , η, and y(t) as per the Heidler function ( 5) and x a (t) is the rise function of the approximation.The next section details the development of the approximation function.

B. Developing an Approximation to the Heidler Function
The development of the proposed approximation is different to those taken in the literature.The only part of the Heidler function that is approximated is the rise function (x h (t)) and the approximation is developed in the Laplace domain ensuring that the time domain equation can be integrated analytically and it is already represented continuously in the frequency domain [30].This implies that the rise equation that is developed must be transformed into the time domain.The decay equation and peak current can be used with this to obtain the overall approximation equation.
The step response of an nth order, real and negative pole creates an S-curve in the time domain.The following equation shows the start of the approximation in the Laplace domain, the step response of an nth order pole: where ω 0 is the rise time constant [rad/s] and n a is the approximation steepness factor.Taking the inverse Laplace transform of this equation, the rise function of the approximation is found to be The constants in this equation can be varied to obtain different steepness factors and rise times (see Section IV-C1).This function is shown in Fig. 5 and is comparable to Fig. 3, which shows the rise time function of the Heidler function.
By substituting the approximated rise function back into the Heidler function given by (11), the overall approximation can be obtained.The next section details the properties of the approximation and its frequency domain representation.

C. Function Definition and Properties
The proposed approximation to the Heidler function (8) has the advantage that it has an analytical integral and an analytical solution in the frequency domain.Moreover, it can still be  tailored to any waveshape, meaning that the steepness of the graph, the rise time, the fall time and the peak current can all be modified.This allows for analysis using a 10/350, 0.25/100, or any other lightning waveshapes required.Modifying I 0 , η, ω 0 , τ 2 , and n a in (8) will give the desired lightning current waveform.An example plot of this function is shown in Fig. 6.
As the approximation was designed as a modification of the Heidler function, it still takes the same form as in (11) with x h (t) replaced with x a (t).The difference is that x h (t) cannot be integrated but x a (t) can be integrated and transformed into the frequency domain analytically.The followingsections show how the steepness factor, rise time constant and fall time constant affect the shape and properties of the approximation function.
1) Steepness Factor: The steepness factor of the approximation, n a , changes the shape of the approximation function.Fig. 7 shows several plots of the approximation function with different    I where the values of I 0 , η, and τ 2 are obtained from the IEC 62305-1 standard and ω 0 is determined empirically for a 10/350 waveshape.As expected from (18), the steepness factor only affects the rise function (x a (t)) of the entire function.From the figure, it is evident that an increased steepness factor decreases the steepness of the rise time.The initial knee in the curve (upward bend) is delayed by a higher value of n a and hence the maximum instantaneous change in current occurs later in time.The figure also shows that the rise time of the waveshape changes with a change in n a but this is not the defining change.
Table II shows the peak current and rise time errors as the steepness factor changes.These errors are calculated as a percentage of a 10/350, 200 kA waveshape.The error in rise time is more pronounced than the peak current error.This table gives   an indication of the sensitivity of the function with a change in the steepness factor (n a ).
2) Rise Time Constant: The rise time constant, ω 0 , is used primarily to change the rise time (the number before the "/") of the waveshape.Fig. 8 shows the effect of changing the rise time constant in the approximation while keeping the rest of the parameters constant.The constant values used in this plot are tabulated in Table III where once again, I 0 , η, and τ 2 are obtained from the IEC 62305-1 standard and n a is determined empirically for a 10/350 waveshape.
This plot shows how the rise time changes with a change in the rise time constant.A greater rise time constant results in a faster rise time.The change in ω 0 also has an effect on the steepness of the graph, but the defining feature is the change in rise time.
Table IV shows the peak current and rise time errors as the rise time constant changes.Again, these errors are calculated as a percentage of a 10/350, 200 kA waveshape.As with the steepness factor, the error in rise time is more pronounced than Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.9 the error in peak current.This table gives an indication to the sensitivity of the function with a change in the rise-time constant (ω 0 ).
3) Fall Time Constant: The fall time constant, τ 2 , changes the decay time (the number after the "/") of the waveshape.Fig. 9 shows the effect of changing the fall time constant in the approximation function with other parameters all fixed.These values are seen in Table V where once again, I 0 and η are obtained from the IEC 62305-1 and ω 0 and n a are determined empirically for a 10/350 waveshape.
Several observations are made in this graph, namely the following.
1) This is the same effect as that seen in the Heidler function, which is as expected.2) As with the Heidler function, the decay time (the time to 50% of the peak current) is not the same as the decay time constant (τ 2 ).
3) The peak current varies slightly with the change in fall time constant as expected from (6).Table VI shows the peak current and rise time errors as the fall time constant changes.Again, these errors are calculated as a percentage of a 10/350, 200 kA waveshape.This table is different to those of the steepness factor and rise time constant (see Table II) because the fall time constant has a negligible effect on the rise part of the waveshape.This table gives an indication to the sensitivity of the function with a change in the fall time constant (τ 2 ).

D. Time Domain Properties
This section discusses the time domain properties of the approximation, namely the time derivative and the integral.Unlike the Heidler function, the approximation has an analytical integral.
1) Derivative: Equation (19) shows the analytical derivative of the approximation.This shows the rate of change of the current in the lightning stroke current model A graph of the above-mentioned function with the same parameters as those utilized in creating the graph in Fig. 6, can be seen in Fig. 10.
As is expected, the greatest rate of change of current occurs during the rise of the graph and the negative rate of change is comparatively small.
2) Integral: A primary reason for developing an approximation to the Heidler function is to have the ability to integrate such a function.Obtaining the integral of the approximation is trivial and can be seen in ( 9), where C is some arbitrary constant of integration because it is an indefinite integral.

E. Frequency Domain Properties
Another reason for using an approximation to the Heidler function is to obtain a Fourier transform and hence a current density or frequency response.The approximation given by (8) makes this process even more trivial.Rather than using the fact that the approximation can be integrated (see Section IV-D2), the approximation is developed in the Laplace domain.Hence the Fourier transform is found by using the complex shifting property of Laplace [see (15)] and replacing s with jω.The result of this can be seen in the following: By using dimensional analysis, it can be seen that in (20), I a (jω) has the units of A/Hz.This implies a current density across the angular frequency spectrum.In order to obtain current density as a function of frequency, ω should be replaced by 2πf Equations ( 20) and ( 21) are complex functions and the modulus is required to plot the current density of the approximation.The modulus is given by ( 22) and a plot of this equation, based on the waveshape shown in Fig. 6, can be seen in Fig. 11 |I a (jf )| = I 0 η V. VERIFICATION METHODOLOGY To verify the accuracy of the proposed approximation to the Heidler function, mathematical modeling software, such as MATLAB [34], Mathematica [35], or Maxima [36], can be used to simulate both functions and assess the level of accuracy using the parameters detailed in Table B.1 of the IEC 62305-1 standard [1] as a control.These parameters represent the 10/350 and the 0.25/100 waveshapes (first positive and subsequent negative impulses, respectively).
Initial values are estimated for the parameters of the approximation from the Heidler function parameters.These values are then empirically optimized in order to minimize the error between the approximation and the Heidler waveshapes.These tabulated and calculated parameters are utilized in ( 5) and ( 8), respectively, to evaluate the accuracy of the approximation.An n a of 33 is found to be appropriate for the waveshapes defined in the standard which have an n h of 10 (see Section VII for more).
The evaluation includes three simulations per waveshape.These are the current waveshape, the change in current (first derivative) and the current density (Fourier transform).The absolute value of the difference between the two functions is determined and the maximum error is defined from this as a percentage of the Heidler function.The first derivative is evaluated in the same manner as the current waveshape.
The error is quantified by determining the error as a function of time as in (23).The maximum error is found by equating the first time derivative of the error function to zero as in (24) and solving for t.This time is substituted back into the error function to obtain the maximum error in kA.This value is found as a percentage of the peak value of the Heidler function where e

(t) is the error function [A] and e (t) is the derivative of error function [A/s].
The current density is evaluated by using the parameters calculated for the approximation in (21).As there is no analytical Fourier transform of the Heidler function, this is only an indication and no quantification of error can be obtained.

VI. COMPARISON WITH THE HEIDLER FUNCTION
The simulations used in this study are those of two of the short lightning current waveshapes detailed in the IEC 62305-1 standard [1].These are the first positive and subsequent negative impulses (10/350 and 0.25/100, respectively).

A. First Positive Impulse (10/350)
The first positive impulse as defined by the IEC 62305-1 has a front time of 10 μs and a decay time of 350 μs (parameters shown Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.  in Table VII) [1].A graph depicting both the Heidler function (solid line) and the approximation (dashed line) waveshapes can be seen in Fig. 12.
As can be seen in the figure, the approximation closely follows the waveshape produced by the Heidler function.In the case of the 10/350 waveshape, with the parameters defined in Table VII, the maximum error is calculated to be as 1.38%.This error is seen to occur during the rise part of the waveshape, which is expected because the decay functions are identical.
The next comparison made is between the first derivatives of both the Heidler function and the approximation [( 7) and (10), respectively].This shows the difference in the instantaneous change in current of the two waveshapes.The same parameters are used, i.e., those in Table VII.The graph showing both of these waveshapes can be seen in Fig. 13.The error is more pronounced in the derivative.Using the same method as previous to obtain the maximum error as a percentage of the Heidler function maximum, the maximum error is calculated to be 7.91%.
The maximum dI/dt occurs during the rise time of the function.Another characteristic of the plot is that the exponential decay is much longer than the rise.This causes the negative component of the derivative to be much smaller but longer than the rise time component.In both Figs. 12 and 13, the time scale goes up to 200 μs.This is because the tails of the two waveshapes are the same.Therefore, the resolution needs to be shown on the  rise time of the waveshapes.The amplitudes shown in the two graphs are large enough to show the maximum values (200 kA in Fig. 12 and 27.5 kA/μs in Fig. 13).
For completeness, the current density of the approximated first positive impulse is plotted in Fig. 14.This is produced using (21) and the parameters in Table VII

B. Subsequent Negative Impulse (0.25/100)
The subsequent negative impulse has a front duration of 0.25 μs and a decay time of 100 μs [1].This implies a much sharper rise time than that of the first positive impulse.The comparison is shown in Fig. 15 using the parameters provided in Table VIII.
The maximum error is calculated to be 1.36% (as before, the error is seen during the rise part of the waveshape).The decay Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.time is significantly greater than the front duration (around 400 times greater), so there is almost no decay in the graph shown.This is so that the variation in the rise part of the function can be seen.The first time derivative of both the Heidler function (solid line) and the approximation (dashed line) are shown in Fig. 16.Again, the error in the first time derivative is more pronounced than in the actual waveshape and this error is calculated to be 7.86% (as before, the error is seen during the rise part of the waveshape).The decay time is so large in comparison  to the rise time, that the negative dI/dt is negligible and in most engineering applications can be assumed to be zero.The maximum change in current is ten times greater than that of the first return stroke.The current density of the approximated subsequent negative impulse can be seen in Fig. 17.This is produced using (21) and the parameters in Table VIII.As expected, there are higher frequency components in the subsequent negative impulse than in the first positive impulse, but the amplitude is lower.

VII. DISCUSSION
The following sections discuss the results obtained in Sections VI-A and VI-B.This is followed by observations with respect to the frequency domain analysis.

A. Lightning Stroke Current
According to the IEC 62305-1 (see Table C.3), the tolerance on peak current is ±10% for both first and subsequent impulses [1].Similarly, the tolerance on the rise time is ±20% for both the first and subsequent impulses.In order for an approximation to be a suitable replacement for the Heidler function, it must fit within these boundaries.Table IX shows the calculated peak current and rise time errors for both functions (Heidler and approximation) and both waveshapes (first and subsequent impulses).The errors for the first positive impulse are calculated as a percentage of a 10/350, 200 kA waveshape and the errors for the subsequent negative impulse are calculated as a percentage of a 0.25/100, 50 kA waveshape.As seen in the table, the approximation is also well within the allowed tolerance and is therefore suitable as a replacement to the Heidler function.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.I 0 (peak current), η (peak current correction), and τ 2 (decay time constant) are the same for both the Heidler function and the approximation.This is expected, as only x h (t) is replaced in the Heidler function with x a (t) in the approximation.As these are the respective rise time functions, it follows that the amplitude and decay time parts of the equation are unaffected.This implies that the approximation is easily interchangeable with the Heidler function.
Fig. 18(a) and (b) shows the effects of changing the steepness factor, n h , on the Heidler function and its first time derivative, respectively.A steeper graph implies a faster rise time and hence the steepness factor affects the rise time of the waveshape, however, the defining feature change is the steepness of the rise time.The half peak value remains at the same point in time for different steepness factors.The maximum instantaneous change in current occurs later in time and with a greater amplitude for a greater n h .This occurs shortly after the first knee (the upward bend) in the curve.Fig. 18(c) and (d) shows the effects of changing the steepness factor, n a , on the approximation and its first time derivative, respectively.The change in the rise time of the waveshape is not as significant as with the Heidler function but the upward trend does still begin later in time.As expected, the maximum instantaneous change in current still occurs later in time but the peak value is decreased with an increase in n a .
Fig. 19(a) and (b) shows the effects of changing the rise time constant, τ 1 , on the Heidler function and its first time derivative, respectively.A smaller rise time constant produces a faster rise time of the waveshape and hence the steepness is also increased.The defining feature is the change in rise time.Again the peak instantaneous change in current occurs later in time with an increase in τ 1 but converse to what is seen with the steepness factor, the peak amplitude decreases.
Fig. 19(c) and (d) shows the effects of changing the rise time constant, ω 0 , on the approximation and its first time derivative, respectively.A greater ω 0 results in a faster rise time.Again this affects the steepness of the rise time graph but the defining feature is the change in rise time.The derivative shows that the increase in ω 0 increases the amplitude of the peak instantaneous change in current, but this occurs earlier in time unlike in the other figures.ω 0 is proportional to the inverse of time (ω 0 ∝ 1 t ) and therefore the opposite of what is seen with the Heidler function (and τ 1 ) is expected with the approximation (and ω 0 ).
The steepness factors and the rise time constants of the Heidler function and the approximation do not have the same effect on their respective waveshapes and hence these values for the approximation are calculated empirically.

B. Lightning Stroke Current Derivative
A few observations can be made from Figs. 13 and 16.The first is that the peak change in current is seen at 50% of peak current value during the rise of the waveshape.This holds true for both the Heidler function and the approximation.This is expected as the rise time function is an S-curve which is most steep in the middle of the rise.The maximum change in current Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. in the subsequent stroke is ten times greater than that of the first stroke.This is critical when designing LPSs as the change in current is directly proportional to the voltage produced across an inductor (any wire).Therefore, the subsequent impulses are more likely to have an effect on more sensitive systems due to even the smallest inductances.The difference in the errors in the first impulse and the subsequent impulse are almost the equivalent-the maximum absolute error in the first impulse is 7.91% and 7.86% in the subsequent impulse.(21).The expected amplitude densities of the first positive and subsequent negative strokes are both clearly marked in the figure .If the values are multiplied by their respective frequencies, i.e., by using ( 21), the amplitude spectrum can be found.For the first positive stroke, the peak current components are between 500 Hz and 40 kHz with a rapid decay above 40 kHz.The subsequent negative stroke has lower amplitude current components with a wider frequency range as expected.The peak current components here are between about 1 kHz and 1 MHz.The subsequent negative stroke is more important when analyzing wideband or high-frequency systems/effects.

D. Further Work
As noted previously, there is an error associated with the current waveshape and its derivative when compared with those of the Heidler function.This error is within a tolerable range and is quantified.It can therefore be taken into account in very sensitive system designs.It is posited previously that this error could be consistent across all variations of the waveshape, however there are only two cases dealt with in this study (two of the waveshapes detailed in the IEC 62305-1) and therefore this hypothesis requires further testing.
Another area for further research is again related to the errors that are quantified previously.It is possible that the errors are as "high" as they are (no approximation will ever have zero error) because the numbers used for n a (approximation steepness factor) and ω 0 (approximation rise time constant) are determined empirically.These could be better calculated by using some form of error optimization algorithm.With these calculated numbers the error could be reduced.There may be some relationship between ω 0 and the parameters used in the Heidler function, particularly τ 1 (Heidler rise time constant).By the same notion there may be some relationship between n a and the parameters used in the Heidler function.These relationships would make plotting different waveshapes even easier than using the ratio method described previously.Further work is required to find such relationships.This however does not affect the validity of the approximation as a suitable replacement for the Heidler function with an analytical integral.

VIII. CONCLUSION
An approximation to the Heidler function that has an analytical integral has been developed and discussed.This is useful in situations where the EM fields produced by a lightning stroke need to be calculated.It is also useful in any scenario where the frequency components of a lightning stroke are required for evaluation.The properties of the approximation are discussed with reference to the IEC 62305-1 standard.The viability of the approximation as a suitable replacement for the Heidler function in the standard has been evaluated through the investigation of the simulations produced.
It can be seen that the waveshapes produced by the approximation are very similar to those produced by the Heidler function.The maximum error in amplitude for the first and subsequent lightning stroke currents is less than 1.5%.The maximum error in the derivative is greater but still less than 8%.This is within the parameters defined in the standard.The simulations detail the frequency response of the approximation for both waveshapes.There is no way of quantifying the error in this result because the Heidler function has no analytical integral and hence no Fourier transform.However the results are similar to what is expected in the lightning protection standard.All of this provides evidence that this approximation is a suitable replacement for the Heidler function when integrals and frequency spectra of lightning strokes are required.
Additional work is required to optimize the approximation.This would further reduce the error and make working with the approximation even easier.This work includes proving or disproving the hypothesis that the error remains constant for any waveshape produced by the approximation.There is only an empirical optimization of the parameters used in this study.A full optimization algorithm should be run in order to reduce the 1.5% inaccuracy between the approximation and the Heidler function.It is posited that there may be some relationship between the parameters used in the Heidler rise time function and the approximation rise time function.Finding this relationship would further simplify the use of the approximation in studies and system design.
The approximation that is developed is found to be a suitable replacement for the Heidler function with the errors quantified and an analytical integral.Evidence is provided to show that the approximation to the Heidler function can be used for the first and subsequent short strokes mentioned in the IEC 62305-1 with less than 1.5% error.

Manuscript received 9
October 2023; revised 26 January 2024; accepted 24 February 2024.Date of publication 19 March 2024; date of current version 15 August 2024.This work was supported in part by the National Research Foundation (NRF) through the Thuthuka programme under Grant TTK23030380641 and in part by DEHNAFRICA through their Johannesburg Lightning Research Laboratory.(Corresponding author: Hugh G. P. Hunt.)

Fig. 1 .
Fig. 1.Physical (a) and transmission line (b) models of a lightning strike close to a transmission line that induces overvoltages on the line.

Fig. 5 .
Fig. 5. Graph depicting the rise function of the approximation function (an S-curve).

Fig. 6 .
Fig. 6.Graph of an example Heidler function approximation lightning current waveshape in the form of a 10/350 with a 200 kA peak current.

Fig. 7 .
Fig. 7. Graph showing the effect of changing the steepness factor (n a ) in the approximation function while keeping all the other variables constant.

Fig. 8 .
Fig. 8. Graph showing the effect of changing the rise time constant (ω 0 ) in the approximation function while keeping all the other variables constant.

Fig. 9 .
Fig. 9. Graph showing the effect of changing the fall time constant (τ 2 ) in the approximation function while keeping all other variables constant.

Fig. 11 .
Fig. 11.Amplitude density of the approximation model produced from the waveshape plotted in Fig.

Fig. 12 .
Fig. 12. Graph of the first positive impulse (10/350) current model using both the Heidler function and the approximation.The time scale is up to 200 μs and the amplitude is as high as 210 kA.

Fig. 13 .
Fig. 13.Graph of the first time derivative of the first positive impulse (10/350) current model using both the Heidler function and the approximation.The time scale is up to 200 μs and the amplitude is as high as 30 kA/μs.

Fig. 14 .
Fig. 14.Current density of the approximation model produced from the waveshape plotted in Fig. 12.
. A current density is plotted because this is what is shown in Figure B.5 in the IEC 62305-1 standard.

Fig. 15 .
Fig.15.Graph of a subsequent negative impulse (0.25/100) current model using both the Heidler function and the approximation.The time scale is up to 5 μs and the amplitude is as high as 52 kA.

Fig. 16 .
Fig.16.Graph of the first time derivative of a subsequent negative impulse current model using both the Heidler function and the approximation.The time scale is up to 5 μs and the amplitude is as high as 300 kA/μs.

Fig. 17 .
Fig. 17.Current density of the approximation model produced from the waveshape plotted in Fig. 15.

Fig. 18 .
Fig. 18.Effects of varying the steepness factors, n h and n a , in both the Heidler function and the approximation respectively.The (a) Heidler function and the (b) Heidler function derivative are compared to the (c) approximation and its (d) derivative.

Fig. 19 .
Fig. 19.Effects of varying the rise time constants, τ 1 and ω 0 , in both the Heidler function and the approximation respectively.The (a) Heidler function and the (b) Heidler function derivative are compared to the (c) approximation and its (d) derivative.

Fig. 20 .
Fig. 20.Amplitude densities of the different lightning currents according to LPL-I as stipulated by the IEC 62305-1 with the results of the Fourier transform of the approximation (adapted from Figure B.7 in [1]).The solid lines are obtained from the standard and clearly marked.

Fig. 20
Fig. 20 shows an adaptation of Figure B.7 in the IEC 62305-1 standard with the results obtained from the Fourier transform of the approximation(21).The expected amplitude densities of the first positive and subsequent negative strokes are both clearly marked in the figure.If the values are multiplied by their respective frequencies, i.e., by using(21), the amplitude spectrum can be found.For the first positive stroke, the peak current components are between 500 Hz and 40 kHz with a rapid decay above 40 kHz.The subsequent negative stroke has lower amplitude current components with a wider frequency range as expected.The peak current components here are between about 1 kHz and 1 MHz.The subsequent negative stroke is more important when analyzing wideband or high-frequency systems/effects.

TABLE I CONSTANT
(8)UES USED IN(8)TO OBTAIN THE GRAPHS IN FIG.7

TABLE II VARIATION
OF RISE TIME AND PEAK CURRENT OF THE APPROXIMATION WITH A CHANGE IN THE STEEPNESS FACTOR (n a ) AS A PERCENTAGE OF A 200 KA 10/350 WAVESHAPE steepness factors but constant rise time constant (ω 0 ), fall time constant (τ 2 ), peak current (I 0 ), and correction factor (η).These values are tabulated in Table

TABLE III CONSTANT
(8)UES USED IN(8)TO OBTAIN THE GRAPHS PLOTTED IN FIG.8

TABLE IV VARIATION
OF RISE TIME AND PEAK CURRENT OF THE APPROXIMATION WITH A CHANGE IN THE RISE-TIME CONSTANT (ω 0 ) AS A PERCENTAGE OF A 200 KA 10/350 WAVESHAPE

TABLE V CONSTANT
VALUES USED IN (8) TO OBTAIN THE GRAPHS IN FIG.

TABLE VI VARIATION
OF RISE TIME AND PEAK CURRENT OF THE APPROXIMATION WITH A CHANGE IN THE FALL TIME CONSTANT (τ 2 ) AS A PERCENTAGE OF A 200 KA 10/350 WAVESHAPE

TABLE VIII PARAMETERS
USED IN (5) AND (8) TO PLOT THE WAVESHAPES SHOWN IN FIG. 15

TABLE IX COMPARISON
OF THE ERRORS FOR THE APPROXIMATION AND THE HEIDLER FUNCTION FOR BOTH THE FIRST AND SUBSEQUENT IMPULSES