#include <complex.h>
#include <stdio.h>
#include <stdbool.h>
#include <fftw3.h>

Macros
#define	PI 3.1415926535897932384626433832795029L

#define	TWO_PI (2.0 * PI)

#define	FOUR_PI (4.0 * PI)

#define	MU0 1.256637061435917e-6

#define	EPS0 8.854187817620e-12

#define	RHO_CU 1.689e-8

Functions
double *	linspace (double a, double b, size_t n, double u[])

double *	logspace (double a, double b, size_t n, double u[])

double	wave_length (double f, double sigma, double ep, double mur)

int	alipio_soil (_Complex double sigma, _Complex double epsr, double sigma0, _Complex double s, double h, double g, double eps_ratio)

int	smith_longmire_soil (_Complex double sigma, _Complex double epsr, double sigma0, _Complex double s, double erinf)

bool	equal_points (const double point_1, const double point_2)

bool	equal_points_tol (const double point_1, const double point_2, double tol)

double	vector_length (const double start_point[3], const double end_point[3])

double	snrm2_ (int n, double x, int *incx)

int	complex_matrix_file (size_t m, size_t n, const _Complex double a, int lda, FILE fp)

int	double_matrix_file (size_t m, size_t n, const double a, int lda, FILE fp)

int	zbesi_ (double zr, double zi, double fnu, int kode, int n, double cyr, double cyi, int nz, int *ierr)

void	print_matrix (char desc, int m, int n, const _Complex double a, int lda)

void	print_matrix_row (char desc, int m, int n, const _Complex double a, int lda)

int	matrix_copy (const _Complex double source, _Complex double target, size_t lds, size_t ldt, size_t nline, size_t ncol)

int	transpose_copy (const _Complex double source, _Complex double target, size_t lds, size_t ldt, size_t nline, size_t ncol)

int	pc_copy (const _Complex double source, _Complex double target, size_t lds, size_t ldt, size_t nline, size_t ncol)

int	pl_copy (const _Complex double source, _Complex double target, size_t lds, size_t ldt, size_t nline, size_t ncol)

int	pcl_copy (const _Complex double source, _Complex double target, size_t lds, size_t ldt, size_t nline, size_t ncol)

int	inv_laplace_trans (double f, _Complex double g, _Complex double *s, double tmax, size_t nt, int filter)

double	heidler (double t, double imax, double tau1, double tau2, int n)

Macro Definition Documentation

◆ EPS0

#define EPS0 8.854187817620e-12

Vacuum permittivity $ \varepsilon_0 $

Definition at line 25 of file auxiliary.h.

◆ FOUR_PI

#define FOUR_PI (4.0 * PI)

math constant $ 4\pi $

Definition at line 21 of file auxiliary.h.

◆ MU0

#define MU0 1.256637061435917e-6

Vacuum permeability $ \mu_0 $

Definition at line 23 of file auxiliary.h.

◆ PI

#define PI 3.1415926535897932384626433832795029L

High Performance Hybrid Electromagnetic Model calculations in C.

All parameters' units are in the SI base if omitted.

Constants, auxiliary functions and routines.math constant $ \pi $

Definition at line 17 of file auxiliary.h.

◆ RHO_CU

#define RHO_CU 1.689e-8

Copper resistivity $ \rho_{cu} $

Definition at line 27 of file auxiliary.h.

◆ TWO_PI

#define TWO_PI (2.0 * PI)

math constant $ 2\pi $

Definition at line 19 of file auxiliary.h.

Function Documentation

◆ alipio_soil()

int alipio_soil	(	_Complex double *	sigma,
		_Complex double *	epsr,
		double	sigma0,
		_Complex double	s,
		double	h,
		double	g,
		double	eps_ratio
	)

Calculates the soil parameters $ \sigma(s) $ and $ \varepsilon_r(s) $ based on the Alipio-Visacro model [1].

$ \sigma = \sigma_0 + \sigma_0 \times h(\sigma_0)\left( \frac{s}{\text{1 MHz}} \right)^g $

$ \varepsilon_r = \frac{\varepsilon_\infty'}{\varepsilon_0} + \frac{\tan(\pi g / 2) \times 10^{-3}}{2\pi\varepsilon_0(\text{1 MHz})^g} \sigma_0 \times h(\sigma_0) s^{g - 1} $

Recommended values of $ h(\sigma_0) $, $ g $ and $ \varepsilon_\infty' / \varepsilon_0 $ are given in Fig. 8 of [1]:

Results	$ h(\sigma_0) $	$ g $	$ \varepsilon_\infty' / \varepsilon_0 $
mean	$ 1.26 \times (1000 \sigma_0)^{−0.73} $	0.54	12
relatively conservative	$ 0.95 \times (1000 \sigma_0)^{−0.73} $	0.58	8
conservative	$ 0.70 \times (1000 \sigma_0)^{−0.73} $	0.62	4

[1] R. Alipio and S. Visacro, "Modeling the Frequency Dependence of Electrical Parameters of Soil," in IEEE Transactions on Electromagnetic Compatibility, vol. 56, no. 5, pp. 1163-1171, Oct. 2014, doi: 10.1109/TEMC.2014.2313977.

Parameters

sigma	pointer to where the conductivity $f \sigma(s) $f is written in S/m
epsr	pointer to where the relative permitivitty $ \varepsilon_r(s) $ is written
sigma0	value of the soil conductivity in low frequency $ \sigma_0 $ in S/m
s	complex frequency $ s = c + j\omega $ of interest in rad/s
h	parameter $ h(\sigma_0) $
g	parameter $ g $
eps_ratio	parameter $ \varepsilon_\infty' / \varepsilon_0 $

Returns: 0 on success

◆ complex_matrix_file()

int complex_matrix_file	(	size_t	m,
		size_t	n,
		const _Complex double *	a,
		int	lda,
		FILE *	fp
	)

Prints a complex (double) matrix to a file.

Parameters

m	number of rows
n	number of columns
a	pointer to array where matrix is stored
lda	leading dimension of $ a $
fp	pointer to file stream

Returns: 0 on success

◆ double_matrix_file()

int double_matrix_file	(	size_t	m,
		size_t	n,
		const double *	a,
		int	lda,
		FILE *	fp
	)

Prints a real (double) matrix to a file.

Parameters

m	number of rows
n	number of columns
a	pointer to array where matrix is stored
lda	leading dimension of $ a $
fp	pointer to file stream

Returns: 0 on success

◆ equal_points()

bool equal_points	(	const double *	point_1,
		const double *	point_2
	)

Check if two points are the same using DBL_EPSILON as tolerance.

Parameters

point1	first point in $ \mathbf{R}^3 $
point2	second point in $ \mathbf{R}^3 $

Returns: true or false

See also: equal_points_tol

◆ equal_points_tol()

bool equal_points_tol	(	const double *	point_1,
		const double *	point_2,
		double	tol
	)

Check if two points are the same within a tolerance.

Parameters

point1	first point in $ \mathbf{R}^3 $
point2	second point in $ \mathbf{R}^3 $
tol	the tolerance $ \epsilon $ applied to each coordinate.

Returns: false if the difference between any coordinates is greater than tol (e.g. $ |x_1 - x_2| > \epsilon $).

See also: equal_points

◆ heidler()

double heidler	(	double	t,
		double	imax,
		double	tau1,
		double	tau2,
		int	n
	)

Heidler function to create lightning current waveforms [1]. For parameters' values, see e.g. [2]. Calculates $ i(t) = \frac{I_0}{\xi} \frac{(t / \tau_1)^n}{1 + (t / \tau_1)^n} e^{-t / \tau_2} $ where $ \xi = e^{-(\tau_1 / \tau_2) (n\,\tau_2 / \tau_1)^{1 / n}} $

[1] HEIDLER, Fridolin; CVETIĆ, J. A class of analytical functions to study the lightning effects associated with the current front. European transactions on electrical power, v. 12, n. 2, p. 141-150, 2002. doi: 10.1002/etep.4450120209

[2] A. De Conti and S. Visacro, "Analytical Representation of Single- and Double-Peaked Lightning Current Waveforms," in IEEE Transactions on Electromagnetic Compatibility, vol. 49, no. 2, pp. 448-451, May 2007, doi: 10.1109/TEMC.2007.897153.

Parameters

t	time $t$ in seconds
imax	current peak $ I_0 $ in Amperes
tau1	rise time $ \tau_1 $ in seconds
tau2	decay time $ \tau_2 $ in seconds
n	steepness expoent

Returns: current $ i(t) $ in Amperes

◆ inv_laplace_trans()

int inv_laplace_trans	(	double *	f,
		_Complex double *	g,
		_Complex double *	s,
		double	tmax,
		size_t	nt,
		int	filter
	)

$ with a real output. Uses FFTW, therefore is not thread-safe! Based on: GÓMEZ, Pablo; URIBE, Felipe A. The numerical Laplace transform: An accurate technique for analyzing electromagnetic transients on power system devices. International Journal of Electrical Power & Energy Systems, v. 31, n. 2-3, p. 116-123, 2009. @param f inverse transfomed function $ f(t) $ of size $ n $ from $ t = 0 $ to $ t = t_\max $ with uniform spacing $ \Delta t $ @param g sampled function $ g(s) $ of size $ \left \lfloor{n/2}\right \rfloor + 1 $ @param s ``frequencies'' $ s_k = c + j\omega_k $ of size $ \left \lfloor{n/2}\right \rfloor + 1 $, the value of $ c $ is assumed constant @param tmax maximum time $ t_\max $f of the inverse transformed function $ f(t) \$f @param nt size $ n $ of the inverse transformed function $ f(t) $f

Parameters

filter enum INLT_Filter that defines the $ \sigma(\omega) $ filter (data window) to apply to $ g(s) := \sigma(\omega) g(s)$

Returns: 0 on success

See also: laplace_trans; INLT_Filter; http://www.fftw.org/fftw3_doc/; planned_inv_laplace

◆ linspace()

double* linspace	(	double	a,
		double	b,
		size_t	n,
		double	u[]
	)

Fill array $ u $ with $ n $ linearly spaced numbers between $ a $ and $ b $ (included).

Parameters

a	first number
b	last number
n	size of array
u	array to be filled

Returns: pointer to filled array

See also: Source code taken from https://github.com/ntessore/algo

◆ logspace()

double* logspace	(	double	a,
		double	b,
		size_t	n,
		double	u[]
	)

Fill array $ u $ with $ n $ logarithmically spaced numbers between $ 10^a $ and $ 10^b $ (included).

Parameters

a	first power
b	last power
n	size of array
u	array to be filled

Returns: pointer to filled array

See also: Source code taken from https://github.com/ntessore/algo

◆ matrix_copy()

int matrix_copy	(	const _Complex double *	source,
		_Complex double *	target,
		size_t	lds,
		size_t	ldt,
		size_t	nline,
		size_t	ncol
	)

Copies the source matrix into target matrix considering COLUMN MAJOR storage. $ B := A $

Parameters

source	pointer to be copied
target	pointer where the copy will be stored
lds	leading dimension of the source array
ldt	leading dimension of the target array
nline	number of lines in source array to copy
ncol	number of columns in source array to copy

Returns: 0 on success

◆ pc_copy()

int pc_copy	(	const _Complex double *	source,
		_Complex double *	target,
		size_t	lds,
		size_t	ldt,
		size_t	nline,
		size_t	ncol
	)

pc_copy Copies the column permutation of source matrix into target matrix considering COLUMN MAJOR storage.

Parameters

source	pointer to be copied
target	pointer where the copy will be stored
lds	leading dimension of the source array
ldt	leading dimension of the target array
nline	number of lines in source array to copy
ncol	number of columns in source array to copy

Returns: 0 on success

◆ pcl_copy()

int pcl_copy	(	const _Complex double *	source,
		_Complex double *	target,
		size_t	lds,
		size_t	ldt,
		size_t	nline,
		size_t	ncol
	)

Copies the column and line permutation of source matrix into target matrix considering COLUMN MAJOR storage. $ B := A^{Pcl} $

Parameters

source	pointer to be copied
target	pointer where the copy will be stored
lds	leading dimension of the source array
ldt	leading dimension of the target array
nline	number of lines in source array to copy
ncol	number of columns in source array to copy

Returns: 0 on success

◆ pl_copy()

int pl_copy	(	const _Complex double *	source,
		_Complex double *	target,
		size_t	lds,
		size_t	ldt,
		size_t	nline,
		size_t	ncol
	)

Copies the line permutation of source matrix into target matrix considering COLUMN MAJOR storage. $ B := A^{Pl} $

Parameters

source	pointer to be copied
target	pointer where the copy will be stored
lds	leading dimension of the source array
ldt	leading dimension of the target array
nline	number of lines in source array to copy
ncol	number of columns in source array to copy

Returns: 0 on success

◆ print_matrix()

void print_matrix	(	char *	desc,
		int	m,
		int	n,
		const _Complex double *	a,
		int	lda
	)

Prints a COLUMN MAJOR matrix to stdio.

Parameters

desc	description to print before the matrix
m	number of rows
n	number of columns
a	column major matrix to be printed
lda	leading dimension of a

◆ print_matrix_row()

void print_matrix_row	(	char *	desc,
		int	m,
		int	n,
		const _Complex double *	a,
		int	lda
	)

Prints a ROW MAJOR matrix to stdio.

Parameters

desc	description to print before the matrix
m	number of rows
n	number of columns
a	row major matrix to be printed
lda	leading dimension of a

◆ smith_longmire_soil()

int smith_longmire_soil	(	_Complex double *	sigma,
		_Complex double *	epsr,
		double	sigma0,
		_Complex double	s,
		double	erinf
	)

◆ snrm2_()

double snrm2_	(	int *	n,
		double *	x,
		int *	incx
	)

LAPACK routine to compute the Euclidean norm of a vector $ \|x\|_2 $.

Parameters

n	Specifies the number of elements in vector $x$.
x	Array of size at least $(1 + (n-1) \times i)$
incx	Specifies the increment $ i $ for the elements of $x$

Returns: $ \|x\|_2 $

◆ transpose_copy()

int transpose_copy	(	const _Complex double *	source,
		_Complex double *	target,
		size_t	lds,
		size_t	ldt,
		size_t	nline,
		size_t	ncol
	)

Copies the transpose of source matrix into target matrix considering COLUMN MAJOR storage. $ B := A^{T} $

Parameters

source	pointer to be copied
target	pointer where the copy will be stored
lds	leading dimension of the source array
ldt	leading dimension of the target array
nline	number of lines in source array to copy
ncol	number of columns in source array to copy

Returns: 0 on success

◆ vector_length()

double vector_length	(	const double	start_point[3],
		const double	end_point[3]
	)

Computes the $ \|\cdot\|_2 $ (euclidian length) of the vector which starts at start_point and ends at end_point.

Parameters

start_point	array $(x,y,z)_1$ that defines the starting point
end_point	array $(x,y,z)_0$ that defines the ending point

Returns: $ \| (x,y,z)_1 - (x,y,z)_0 \|_2 $

◆ wave_length()

double wave_length	(	double	f,
		double	sigma,
		double	ep,
		double	mur
	)

Computes the harmonic electromagnetic wave length $ \lambda $ in a lossy medium.

Parameters

f	frequency (Hz)
sigma	medium (real) conductivity $ \sigma $ in $ \Omega \cdot m $
ep	medium electric permittivity $ \varepsilon = \varepsilon_r \, \varepsilon_0$
mur	medium relative magnetic permeability $ \mu_r $

Returns: $ \lambda $ in $ m $

◆ zbesi_()

int zbesi_	(	double *	zr,
		double *	zi,
		double *	fnu,
		int *	kode,
		int *	n,
		double *	cyr,
		double *	cyi,
		int *	nz,
		int *	ierr
	)

FORTRAN subroutine (from SLATEC) to calculate I-Bessel function, i.e., modified Bessel function of the first kind, with complex argument. $ cy = I_{fnu} (z) $

Parameters

zr	argument's real part $\Re(z)$
zi	argument's imaginary part $\Im(z)$
fnu	order of initial I function
kode	scaling 1: no scaling 2: $e^{-\|x\|}$
n	number of members of the sequence
cyr	result's real part $\Re(cy)$
cyi	result's imaginary part $\Im(cy)$
nz	number of components set to zero due to underflow
ierr	error flag

See also: http://netlib.org/amos/zbesi.f

sigma	pointer to where the conductivity $f \sigma(s) $f is written in S/m
epsr	pointer to where the relative permitivitty \( \varepsilon_r(s) \) is written
sigma0	value of the soil conductivity in low frequency \( \sigma_0 \) in S/m
s	complex frequency \( s = c + j\omega \) of interest in rad/s
h	parameter \( h(\sigma_0) \)
g	parameter \( g \)
eps_ratio	parameter \( \varepsilon_\infty' / \varepsilon_0 \)

n	Specifies the number of elements in vector \(x\).
x	Array of size at least \((1 + (n-1) \times i)\)
incx	Specifies the increment \( i \) for the elements of \(x\)

Results	\( h(\sigma_0) \)	\( g \)	\( \varepsilon_\infty' / \varepsilon_0 \)
mean	\( 1.26 \times (1000 \sigma_0)^{−0.73} \)	0.54	12
relatively conservative	\( 0.95 \times (1000 \sigma_0)^{−0.73} \)	0.58	8
conservative	\( 0.70 \times (1000 \sigma_0)^{−0.73} \)	0.62	4

point1	first point in \( \mathbf{R}^3 \)
point2	second point in \( \mathbf{R}^3 \)

t	time $t$ in seconds
imax	current peak \( I_0 \) in Amperes
tau1	rise time \( \tau_1 \) in seconds
tau2	decay time \( \tau_2 \) in seconds
n	steepness expoent

start_point	array \((x,y,z)_1\) that defines the starting point
end_point	array \((x,y,z)_0\) that defines the ending point

f	frequency (Hz)
sigma	medium (real) conductivity \( \sigma \) in \( \Omega \cdot m \)
ep	medium electric permittivity \( \varepsilon = \varepsilon_r \, \varepsilon_0\)
mur	medium relative magnetic permeability \( \mu_r \)

zr	argument's real part \(\Re(z)\)
zi	argument's imaginary part \(\Im(z)\)
fnu	order of initial I function
kode	scaling 1: no scaling 2: \(e^{-\|x\|}\)
n	number of members of the sequence
cyr	result's real part \(\Re(cy)\)
cyi	result's imaginary part \(\Im(cy)\)
nz	number of components set to zero due to underflow
ierr	error flag

Macros

Functions

Macro Definition Documentation

◆ EPS0

◆ FOUR_PI

◆ MU0

◆ PI

◆ RHO_CU

◆ TWO_PI

Function Documentation

◆ alipio_soil()

◆ complex_matrix_file()

◆ double_matrix_file()

◆ equal_points()

◆ equal_points_tol()

◆ heidler()

◆ inv_laplace_trans()

◆ linspace()

◆ logspace()

◆ matrix_copy()

◆ pc_copy()

◆ pcl_copy()

◆ pl_copy()

◆ print_matrix()

◆ print_matrix_row()

◆ smith_longmire_soil()

◆ snrm2_()

◆ transpose_copy()

◆ vector_length()

◆ wave_length()

◆ zbesi_()