Functions¶

Support¶

support.add_referencefield(B0_name, npath, dmeas, B0, B_si, benchmark)[source]¶

This function returns total computed magnetic field, along with its intensity, inclination and declination at all points on a specified path using the following steps:,

Rotating the vector from the model’s coordinate system to align with the paleomagnetism coordinate system.
Incorporating a (global) reference field, denoted as B0.
Applying standard paleomagnetic equations for calculation as referenced in [Tauxe, 2010].

Subsequently, the function outputs the individual (x,y,z) components of the total magnetic field, as well as the calculated intensity, inclination, and declination values. These results are saved into a file titled: measurements_path_refField{BO_name}.ascii

Parameters

B0_name (str) – name of reference field (used in writing file).
npath (scalar(int)) – amount of points on path.
dmeas (array_like(float)) – 1D array(npath) containing distance between two points of the path (from meyer), used for plotting.
B0 (array_like(float)) – 1D array(3) containing 3 components (0=x;1=y;2=z), components of reference field to be added.
B_si (array_like(float)) – 2D array(3,npath), containing components (0=x;1=y;2=z,npath) of the computed magnetic field due to the underlying magnetized matter (anomalous field) at each point of the path.
benchmark (str) – number associated with benchmark, see Benchmarks.

Returns

B_siB0 (array_like(float)) - 2D array(3,npath), containing components (0=x;1=y;2=z,npath) of the total magnetic field (anomalous + reference field) at each point of the path.
In_siB0 (array_like(float)) - 1D array(npath), containing intensity of the total magnetic field (anomalous + reference field) at each point of the path.
Ic_siB0 (array_like(float)) - 1D array(npath), containing inclination of the total magnetic field (anomalous + reference field) at each point of the path.
Dc_siB0 (array_like(float)) - 1D array(npath), containing declination of total magnetic field (anomalous + reference field) at each point of the path.

support.compute_analytical_solution(x, y, z, R, Mx, My, Mz, xcenter, ycenter, zcenter, benchmark)[source]¶

Returns analytical solution, see Benchmarks, of the components of the magnetic field for each benchmark.

Parameters

x (array_like(float)) – 1D array(npath) containing x coordinate of each computation point.
y (array_like(float)) – 1D array(npath) containing y coordinate of each computation point.
z (array_like(float)) – 1D array(npath) containing z coordinate of each computation point.
R (scalar(float)) – radius of the sphere (if applicable)
Mx (scalar(float)) – x component magnetization of magnetized body (constant).
My (scalar(float)) – y component magnetization of magnetized body (constant).
Mz (scalar(float)) – z component magnetization of magnetized body (constant).
xcenter (scalar(float)) – x coordinate of center of the sphere.
ycenter (scalar(float)) – y coordinate of center of the sphere.
zcenter (scalar(float)) – z coordinate of center of the sphere.
benchmark (str) – number associated with benchmark, see Benchmarks.

Returns

h_fs (array_like(float)) - 2D array(3,npath) containing components (0=x;1=y;2=z) produced by the analytical solution for the magnetic field for chosen benchmark (and model setup) and each point.

support.compute_magnetization(lat, Mint=7.5)[source]¶

Returns components of magnetization, assuming constant magnetization all over the globe: (not weakening towards the equator), magnetization intensity based on case study (7.5 A/m) is default.

Parameters

lat (scalar(float)) – lattitude in degrees (range: -90:90)
Mint (scalar(float)) – Intensity assigned to the magnetized matter in

Returns

M_lat (array_like(float)) - 1D array(3), containing components (0=x;1=y;2=z) of the magnetization in model configuration (not PMAG).

support.get_IGRF(lat, long=15, h=1.3, date=datetime.datetime(2018, 4, 1, 0, 0))[source]¶

Returns the components of the IGRF magnetic field at a given location, height and date. Default values are for Mt. Etna case study at various latitudes. Uses ppigrf package.

Parameters

lat (float) – Latitude in degrees (range: -90 to 90)
long (float) – Longitude in degrees
h (float) – Height above the Earth’s surface in km
date (datetime) – Date for which the IGRF field is evaluated

Returns

IGRF_E (float) - East component of the magnetic field [T]
IGRF_N (float) - North component of the magnetic field [T]
IGRF_D (float) - Downward component of the magnetic field [T]

Return type

tuple(float, float, float)

support.is_point_near_diagonal(x, y, Lx, Ly, eps)[source]¶

This function checks if a point (x, y) is near to either of the two diagonals of a rectangle with sides of length Lx and Ly.

Distance units are irrelevant but must be consistent.

Parameters

x (scalar(float)) – x-coordinate of the observation point.
y (scalar(float)) – y-coordinate of the observation point.
Lx (scalar(float)) – length of the edge of a rectangle along the x-axis.
Ly (scalar(float)) – length of the edge of a rectangle along the y-axis.
eps (scalar(float)) – tolerance value for proximity check.

Returns

bool (bool) - True if (x, y) is near (<eps) to either diagonal, False otherwise.

support.read_header(topo_file)[source]¶

This function checks if a DEM file has a header, and if so, reads in the values.
These values are stored in dictionary {}. Header keys checked for: ncols, nrows, xllcorner, yllcorner, cellsize, and NODATA_value.
This function works optimally when the header is in a standardized ASCII format.

Parameters

topo_file (string) – the opened (with read statement) “topofile”, containing the possible header to be read in.

Returns

header (array_like(string)) - dictionary containing values read in from header.
statement (bool)) - bool set to True if the full header was read in completely, in case of an unexpected header term found or no header found, it is set to False.

support.shift_observation_points_edge(x, y, Lx, Ly, nelx, nely, nelz, xm, ym)[source]¶

The original subroutine of [Blakely, 1995], was initially designed exclusively for use with the faces of a polyhedron. In contrast, our model s subdivides the top and bottom of the hexahedron elements into four triangles, see magnetostatics.compute_B_surface_integral_wtopo(). This introduces additional singularities, and thus spatial issues, if the observation point p lies on the diagonal of the top or bottom face of an element.
This function adjusts observation points positioned near the diagonal of additional triangles on the top (or bottom) of a hexahedron cell.
This function uses support.is_point_near_diagonal(), passing the length of an element in the x- and y-direction, and the x- and y- coordinates of each observation point with the x-coordinate and y-coordinate of each node on the top surface (see nnz-1) subtracted. This ensures testing of proximity to the diagonals of all elements. The function is designed to stop after first modification, as within macroscopic purposes we do not expect an observation point to be shifted to another diagonal. Afterwards, it states any modification within the “message”.
Both x- and y-coordinates are shifted using a random value between -1 and 1 times the artificial distance (factor).
This spatial problem isn’t restricted to the domain bounds. While generating an artificial distance for singularities within the main calculation function (magnetostatics.facmag()) similar to the solution of [Blakely, 1995, Bott, 1963] for edge alignment, would be preferred. The problem emerges from our decision to utilize the magnetostatics.facmag() in this particular manner (by subdividing the top and bottom into additional triangles). Since this setup is invoked via an external function, and because it falls outside the scope of the function’s original intention (to handle planes as polyhedron sides), we have opted to create this external function as well. Furthermore, there are several setups that do not require subdivision, and subsequently use the magnetostatics.compute_B_surface_integral_cuboid(), hence shifting there could only potentially introduce inaccuracies. This is avoided by the introduction of this function, only called for in case function magnetostatics.compute_B_surface_integral_wtopo() is employed, see Benchmarks, Synthetic topopography: flank simulations, Case study: Mount Etna.

Benchmarks have established that the epsilon should be at least of a factor 1e-5 or larger. It is important to note that this value is not dynamically scaled and should be treated with careful consideration, particularly if there are any modifications to the core implementation of the model.

Parameters

x (array_like(float)) – 1D array(NV) containing the x-coordinate of each domain nodes.
y (array_like(float)) – 1D array(NV) containing the y-coordinate of each domain nodes.
Lx (scalar(float)) – length of domain along the x-axis.
Ly (scalar(float)) – length of domain along the y-axis.
nelx (scalar(int)) – amount of elements (domain) in x-direction.
nely (scalar(int)) – amount of elements (domain) in y-direction.
nelz (scalar(int)) – amount of elements (domain) in z-direction.
xm (scalar(float)) – x-coordinate of the observation point.
ym (scalar(float)) – y-coordinate of the observation point.

Returns

xm (scalar(float)) - adjusted x coordinate of the observation point.
ym (scalar(float)) - adjusted x coordinate of the observation point.
message (str)) - message indicating any adjustments made to the observation point.

support.topography(x, y, A, llambda, dir, slopex, slopey)[source]¶

Returns topography (height) value at each point defined by x,y coordinates passed.

Parameters

x (scalar(float)) – x coordinate of computation point.
y (scalar(float)) – y coordinate of computation point.
A (scalar(float)) – amplitude of sine function used to approximate wavy topography of ridges and gullies.
llambda (scalar(float)) – wavelength of sine function used to approximate wavy topography of ridges and gullies.
cos_dir (scalar(float)) – cosine of direction, that defines the direction of the wave function along the surface (default = perpendicular to slope).
sin_dir (scalar(float)) – sine of direction, that defines the direction of the wave function along the surface (default = perpendicular to slope).
slopex (scalar(float)) – the tangent of the angle (in rad) of the surface that defines the slope in x direction for each respective side.
slopey (scalar(float)) – the tangent of the angle (in rad) of the surface that defines the slope in y direction for each respective side.

Returns

h_fs (scalar(float)) - the height value for the point passed to function for the respective position on the flank (chosen by subbench).

Visualization¶

tools.export_line_measurements(N, x, y, z, filename, B_vi, B_si, B_th)[source]¶

Export 1D line measurements to vtu file (plot with paraview), used for exporting line measurements with magnetic field strength (B) computed by analytical solution, surface integral or volume integral as vectors. If this function is used, but the volume integral is not used (default = not used, see Computational approach), or a setup does not have an analytical solution (flank simulations and Case study: Mount Etna), the values from surface integral computation are duplicated.

Parameters

N (scalar(int)) – number of measurement points of the line or path.
x (array_like(float)) – 1D array(N) containing x coordinate of each measurement point.
y (array_like(float)) – 1D array(N) containing y coordinate of each measurement point.
z (array_like(float)) – 1D array(N) containing z coordinate of each measurement point.
filename (str) – the name of the output file (default: line_measurements.vtu)
B_vi (array_like(float)) – 2D array(3,NV), containing components (0=x;1=y;2=z,N) of the computed magnetic field due to the underlying magnetized matter (anomalous field) at each measurement point computed by volume integral.
B_si (array_like(float)) – 2D array(3,N), containing components (0=x;1=y;2=z,N) of the computed magnetic field due to the underlying magnetized matter (anomalous field) at each measurement point computed by surface integral.
B_th (array_like(float)) – 2D array(3,N), containing components (0=x;1=y;2=z,N) of the computed magnetic field due to the underlying magnetized matter (anomalous field) at each measurement point computed by analytical solution.

Returns

tools.export_mesh_1D(npath, xpath, ypath, zpath, filename)[source]¶

Export 1D path to vtu file (plot with paraview), can be used in main visualize measurement path.

Parameters

npath (scalar(int)) – number of measurement points of the line.
xpath (array_like(float)) – 1D array(npath) containing x coordinate of each measurement point.
ypath (array_like(float)) – 1D array(npath) containing y coordinate of each measurement point.
zpath (array_like(float)) – 1D array(npath) containing z coordinate of each measurement point.
filename (str) – the name of the output file (default: path.vtu)

Returns

tools.export_mesh_2D(NV, nel, x, y, z, icon, filename)[source]¶

This function creates a VTU file that represents a 2D mesh surface, which can be utilized for visualization in Paraview.

Parameters

NV (scalar(int)) – number of points.
nel (scalar(int)) – number of elements.
x (array_like(float)) – 1D array(NV) containing x coordinate of each point.
y (array_like(float)) – 1D array(NV) containing y coordinate of each point.
z (array_like(float)) – 1D array(NV) containing z coordinate of each point.
icon (array_like(int)) – 2D array(4,nel), containing connectivity, marking points that define an element, see code.
filename (str) – the name of the output file (default: mesh_plane_measurements.vtu)

Returns

tools.export_mesh_3D(NV, nel, x, y, z, icon, filename, Mx, My, Mz, nnx, nny, nnz)[source]¶

This function creates a VTU file that represents a 3D mesh with height data exported as point data and magnetization vectors exported as cell data. The VTU file can be visualized using Paraview software.

Parameters

NV (scalar(int)) – number of nodes.
nel (scalar(int)) – number of elements.
x (array_like(float)) – 1D array(NV) containing x coordinate of each node.
y (array_like(float)) – 1D array(NV) containing y coordinate of each node.
z (array_like(float)) – 1D array(NV) containing z coordinate of each node.
icon (array_like(int)) – 2D array(8,nel), containing connectivity, marking nodes that define an element, see code.
filename (str) – the name of the output file (default: mesh.vtu)
Mx (array_like(float)) – 1D array(nel) x component magnetization of each cell.
My (array_like(float)) – 1D array(nel) y component magnetization of each cell.
Mz (array_like(float)) – 1D array(nel) z component magnetization of each cell.
nnx (scalar(int)) – number of nodes, x direction.
nny (scalar(int)) – number of nodes, y direction.
nnz (scalar(int)) – number of nodes, z direction.

Returns

tools.export_plane_measurements(NV, nel, x, y, z, icon, filename, B_vi, B_si, B_th)[source]¶

Export 2D plane measurements to vtu file (plot with paraview), used for exporting plane measurements with magnetic field strength (B) computed by analytical solution, surface integral or volume integral as vectors. If this function is used, but the volume integral is not used (default = not used, see Computational approach), or a setup does not have an analytical solution (flank simulations and Case study: Mount Etna), the values from surface integral computation are duplicated.

Parameters

NV (scalar(int)) – number of measurement points of the plane.
nel (scalar(int)) – number of elements of the measurement plane.
x (array_like(float)) – 1D array(NV) containing x coordinate of each measurement point.
y (array_like(float)) – 1D array(NV) containing y coordinate of each measurement point.
z (array_like(float)) – 1D array(NV) containing z coordinate of each measurement point.
icon (array_like(int)) – 2D array(4,nel), containing connectivity, marking point that define an element, see code.
filename (str) – the name of the output file (default: plane_measurements.vtu)
B_vi (array_like(float)) – 2D array(3,NV), containing components (0=x;1=y;2=z,NV) of the computed magnetic field due to the underlying magnetized matter (anomalous field) at each measurement point computed by volume integral.
B_si (array_like(float)) – 2D array(3,NV), containing components (0=x;1=y;2=z,NV) of the computed magnetic field due to the underlying magnetized matter (anomalous field) at each measurement point computed by surface integral.
B_th (array_like(float)) – 2D array(3,NV), containing components (0=x;1=y;2=z,NV) of the computed magnetic field due to the underlying magnetized matter (anomalous field) at each measurement point computed by analytical solution.

Returns

tools.export_spiral_measurements(N, x, y, z, filename, B_vi, B_si, B_th)[source]¶

Export 1D spiral measurements to vtu file (plot with paraview), used for exporting spiral measurements with magnetic field strength (B) computed by analytical solution, surface integral or volume integral as vectors. If this function is used, but the volume integral is not used (default = not used, see Computational approach), or a setup does not have an analytical solution (flank simulations and Case study: Mount Etna), the values from surface integral computation are duplicated.

Parameters

N (scalar(int)) – number of measurement points of the spiral.
x (array_like(float)) – 1D array(N) containing x coordinate of each measurement point.
y (array_like(float)) – 1D array(N) containing y coordinate of each measurement point.
z (array_like(float)) – 1D array(N) containing z coordinate of each measurement point.
filename (str) – the name of the output file (default: line_measurements.vtu)
B_vi (array_like(float)) – 2D array(3,N), containing components (0=x;1=y;2=z,N) of the computed magnetic field due to the underlying magnetized matter (anomalous field) at each measurement point computed by volume integral.
B_si (array_like(float)) – 2D array(3,N), containing components (0=x;1=y;2=z,N) of the computed magnetic field due to the underlying magnetized matter (anomalous field) at each measurement point computed by surface integral.
B_th (array_like(float)) – 2D array(3,N), containing components (0=x;1=y;2=z,N) of the computed magnetic field due to the underlying magnetized matter (anomalous field) at each measurement point computed by analytical solution.

Returns

Magnetostatics¶

magnetostatics.NNN(r, s, t)¶

Q1 basis functions inside the [-1:1]x[-1:1]x[-1:1] reference element

Parameters

r (scalar(float)) – Local coordinate in x-direction (range: -1 to 1)
s (scalar(float)) – Local coordinate in y-direction (range: -1 to 1)
t (scalar(float)) – Local coordinate in z-direction (range: -1 to 1)

Returns

N (array_like(float)) - 1D array(8) containing the values of the 8 Q1 basis functions at (r, s, t).

magnetostatics.compute_B_quadrature(xmeas, ymeas, zmeas, x, y, z, icon, Mx, My, Mz, nqdim)¶

Solves volume integral, numerical solution, as the volume integral is parameterized by the number of quadrature points per dimension (nqdim). Computes magnetic field components based on 2^3 quadrature point integration produced by a single hexahedron (cuboid) carrying a magnetization vector (Mx,My,Mz) assumed to be constant inside the element.

Parameters

xmeas (scalar(float)) – x coordinate of observation point.
ymeas (scalar(float)) – y coordinate of observation point.
zmeas (scalar(float)) – z coordinate of observation point.
x (array_like(float)) – 1D array(NV) containing x coordinate of each observation point.
y (array_like(float)) – 1D array(NV) containing y coordinate of each observation point.
z (array_like(float)) – 1D array(NV) containing z coordinate of each observation point.
icon (array_like(int)) – 1D array(8), containing connectivity, marking nodes that define an element, see code.
Mx (scalar(float)) – x component magnetization [A/m] of the element.
My (scalar(float)) – y component magnetization [A/m] of the element.
Mz (scalar(float)) – z component magnetization [A/m] of the element.
nqpts (scalar(int)) – the number of quadrature points per dimension.

Returns

B_vi (array_like(float)) - 1D array(3) containing 3 components (0=x;1=y;2=z) of the magnetic field strength [T] at a point computed by volume integral method.

magnetostatics.compute_B_surface_integral_cuboid(xmeas, ymeas, zmeas, x, y, z, icon, Mx, My, Mz)¶

This function computes the magnetic field at a point (defined my x,y,z-coordinate) produced by a cuboid element employing facmag() function on each face. Here again the magnetization vector (Mx,My,Mz) is assumed to be constant inside the element. It uses the icon array to identify the x-,y-,z-coordinates associated with each node, before calling facmag(), to compute the total magnetic field strength for each face.
Magnetic field strength \(\mathbf{B}\) is computed in Tesla [T]. Note that the negative of the computed values is passed, as the numbering of our nodes was counterclockwise, not clockwise as in [Blakely, 1995].
Distance units are irrelevant but must be consistent.  For numbering of the nodes, see comments in code.

Parameters

xmeas (scalar(float)) – x coordinate of observation point.
ymeas (scalar(float)) – y coordinate of observation point.
zmeas (scalar(float)) – z coordinate of observation point.
x (array_like(float)) – 1D array(NV) containing x coordinate of each node.
y (array_like(float)) – 1D array(NV) containing y coordinate of each node.
z (array_like(float)) – 1D array(NV) containing z coordinate of each node.
icon (array_like(int)) – 1D array(8), containing connectivity, marking nodes that define an element, see code.
Mx (scalar(float)) – x component magnetization [A/m] of the element.
My (scalar(float)) – y component magnetization [A/m] of the element.
Mz (scalar(float)) – z component magnetization [A/m] of the element.

Returns

B_si (array_like(float)) - 1D array(3) containing 3 components (0=x;1=y;2=z) of the magnetic field strength in [T] at each point computed by surface integral method for a cuboid element.

magnetostatics.compute_B_surface_integral_wtopo(xmeas, ymeas, zmeas, x, y, z, icon, Mx, My, Mz)¶

This function computes the magnetic field at a point (defined my x,y,z-coordinate) produced by a hexahedron element which vertical sides are planar. Only the top and bottom faces can contain 4 nodes which are not co-planar.
In light thereof we subdivide the top and bottom faces into four triangles. This feature is needed in the case topography is prescribed at the top (or bottom) of the domain and the vertical position of the nodes are modified.
It uses the icon array to identify the x-,y-,z-coordinates associated with each node, before calling facmag(), to compute the total magnetic field strength for each face and triangle.
Magnetic field strength \(\mathbf{B}\) is computed in Tesla [T]. Note that the negative of the computed values is passed, as the numbering of our nodes was counterclockwise, not clockwise as in [Blakely, 1995].
Distance units are irrelevant but must be consistent. For numbering of the nodes, see comments in code.

Parameters

xmeas (scalar(float)) – x coordinate of observation point.
ymeas (scalar(float)) – y coordinate of observation point.
zmeas (scalar(float)) – z coordinate of observation point.
x (array_like(float)) – 1D array(NV) containing x coordinate of each node.
y (array_like(float)) – 1D array(NV) containing y coordinate of each node.
z (array_like(float)) – 1D array(NV) containing z coordinate of each node.
icon (array_like(int)) – 1D array(8), containing connectivity, marking nodes that define an element, see code.
Mx (scalar(float)) – x component magnetization [A/m] of the element.
My (scalar(float)) – y component magnetization [A/m] of the element.
Mz (scalar(float)) – z component magnetization [A/m] of the element.

Returns

B_si (array_like(float)) - 1D array(3) containing 3 components (0=x;1=y;2=z) of the magnetic field strength in tesla [T] at a point computed by surface integral method for a cuboid element with topography on top and/or bottom surface.

magnetostatics.compute_B_surface_integral_wtopo_noise(xmeas, ymeas, zmeas, x, y, z, icon, Mx, My, Mz, noise)¶

This function computes the magnetic field produced by a hexahedron element which vertical sides are planar. Only the top and bottom faces can contain 4 nodes which are not co-planar, and are therefor subdivided into four triangles.
Numbering of nodes is same as previous function. The noise is added to the extra node on the top/bottom surface.
Note that when computation are caried out and noise is added to this middle point, it could be that the observation point is then location within the domain, but we do not correct for it nor test it. So, this situation must be avoided during setup of any testing employing this function.
Distance units are irrelevant but must be consistent.

Parameters

xmeas (scalar(float)) – x coordinate of observation point.
ymeas (scalar(float)) – y coordinate of observation point.
zmeas (scalar(float)) – z coordinate of observation point.
x (array_like(float)) – 1D array(NV) containing x coordinate of each node.
y (array_like(float)) – 1D array(NV) containing y coordinate of each node.
z (array_like(float)) – 1D array(NV) containing z coordinate of each node.
icon (array_like(int)) – 1D array(8), containing connectivity, marking nodes that define an element, see code.
Mx (scalar(float)) – x component magnetization [A/m] of the element.
My (scalar(float)) – y component magnetization [A/m] of the element.
Mz (scalar(float)) – z component magnetization [A/m] of the element.
noise (scalar(float)) – the height substracted or added to the middle node on the top and/or bottom surface of each element representing noise potentially smoothed by the DEM.

Returns

B_si (array_like(float)) - 1D array(3) containing 3 components (0=x;1=y;2=z) of the magnetic field strength in [T] at a point computed by surface integral method for a cuboid element with topography on top and/or bottom surface and noise.

magnetostatics.cross(a, b)¶

Produces vector (cross) product of two vectors. \(\vec{a} \times \vec{b} = \vec{c}\)

Parameters

a (array_like(float)) – 1D array(3) containing 3 components (0=x;1=y;2=z) of vector a.
b (array_like(float)) – 1D array(3) containing 3 components (0=x;1=y;2=z) of vector b.

Returns

c (array_like(float)) - 1D array(3) containing 3 components (0=x;1=y;2=z) of vector c..

magnetostatics.dNNNdr(r, s, t)¶

Derivatives of the Q1 basis functions with respect to the local coordinate r.

Parameters

r (scalar(float)) – Local coordinate in x-direction (range: -1 to 1)
s (scalar(float)) – Local coordinate in y-direction (range: -1 to 1)
t (scalar(float)) – Local coordinate in z-direction (range: -1 to 1)

Returns

dNdr (array_like(float)) - 1D array(8) containing derivatives of the basis functions with respect to r at (r, s, t).

magnetostatics.dNNNds(r, s, t)¶

Derivatives of the Q1 basis functions with respect to the local coordinate s.

Parameters

r (scalar(float)) – Local coordinate in x-direction (range: -1 to 1)
s (scalar(float)) – Local coordinate in y-direction (range: -1 to 1)
t (scalar(float)) – Local coordinate in z-direction (range: -1 to 1)

Returns

dNds (array_like(float)) - 1D array(8) containing derivatives of the basis functions with respect to s at (r, s, t).

magnetostatics.dNNNdt(r, s, t)¶

Derivatives of the Q1 basis functions with respect to the local coordinate t.

Parameters

r (scalar(float)) – Local coordinate in x-direction (range: -1 to 1)
s (scalar(float)) – Local coordinate in y-direction (range: -1 to 1)
t (scalar(float)) – Local coordinate in z-direction (range: -1 to 1)

Returns

dNdt (array_like(float)) - 1D array(8) containing derivatives of the basis functions with respect to t at (r, s, t).

magnetostatics.facmag(Mx, My, Mz, x0, y0, z0, x, y, z, n)¶

Function FACMAG computes the magnetic field strength B due to surface charge on a polygonal face. Repeated calls on each face of a polyhedron builds the field of said arbitrary polyhedron, algorithm from [Bott, 1963].
The polygon is limited to 10 corners (n). Requires functions rot(), line(), plane() as define above. Distance units are irrelevant but must be consistent.
A spatial problem, see Computational approach, is mitigated due to the introduction of a small (1e-20) artificial distance for p, in case p lies on a plane aligning with an edge of an element.
Directly translated from ANSI-standard FORTRAN 77 subroutines from [Blakely, 1995].

Parameters

Mx (scalar(float)) – x component magnetization [A/m] of the element.
My (scalar(float)) – y component magnetization [A/m] of the element.
Mz (scalar(float)) – z component magnetization [A/m] of the element.
x0 (scalar(float)) – x coordinate of observation point (xmeas).
y0 (scalar(float)) – y coordinate of observation point (ymeas).
z0 (scalar(float)) – z coordinate of observation point (zmeas).
x (array_like(float)) – 1D array(n+1) containing x coordinate for each polygon corner (and wrapping around).
y (array_like(float)) – 1D array(n+1) containing y coordinate for each polygon corner (and wrapping around).
z (array_like(float)) – 1D array(n+1) containing z coordinate for each polygon corner (and wrapping around).
n (scalar(int)) – amount of corners defining the polygon surface, max 10.

Returns

B (array_like(float)) - 1D array(3) containing 3 components (0=x;1=y;2=z) of the magnetic field strength [T] at a point due to the surface charges of a polygon surface.

magnetostatics.line(x0, y0, z0, x1, y1, z1, x2, y2, z2)¶

Function LINE determines the intersection (x,y,z) of two lines.  First line is defined by points (x1,y1,z1) and (x2,y2,z2).
Second line is perpendicular to the first and passes through point (x0,y0,z0).Distance between (x,y,z) and (x0,y0,z0) is returned as r. Computation is done by a transformation of coordinate systems, and defining two vectors T0 and T2 in the process. T0 = vector from point 1 (x1,y1,z1) to point 0 (x0,y0,z0). T2 = vector from point 1 (x1,y1,z1) to point 2 (x2,y2,z2).
Directly translated from ANSI-standard FORTRAN 77 subroutines from [Blakely, 1995].

Parameters

x0 (scalar(float)) – x coordinate of point 0.
y0 (scalar(float)) – y coordinate of point 0.
z0 (scalar(float)) – z coordinate of point 0.
x1 (array_like(float)) – x coordinate of point 1.
y1 (array_like(float)) – y coordinate of point 1.
z1 (array_like(float)) – z coordinate of point 1.
x2 (array_like(float)) – x coordinate of point 2.
y2 (array_like(float)) – y coordinate of point 2.
z2 (array_like(float)) – z coordinate of point 2.

Returns

x (scalar(float)) - x coordinate of intersection point between the lines
y (scalar(float)) - y coordinate of intersection point between the lines
z (scalar(float)) - z coordinate of intersection point between the lines
v1 (scalar(float)) - the negative of the projection of vector T0 on T2 (dot product of normalized T2 and T0).
v2 (scalar(float)) - the length of vector T2 minus the projection of vector T0 on T2.
r (scalar(float)) - distance between intersection point and point 0 (x0,y0,z0)

magnetostatics.plane(x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3)¶

Function plane computes the intersection (x,y,z) of a plane and a perpendicular line.
The plane is defined by three points (x1,y1,z1), (x2,y2,z2), and (x3,y3,z3), in facmag the (first 3) polygon corners are passed. The line passes through (x0,y0,z0), in facmag the observation point is passed.
Computation is done by a transformation and inverse transformation of coordinates systems. It defines three vectors N0, N2 and N3 in the process. N0 = vector from point 1 (x1,y1,z1) to point 0 (x0,y0,z0). N2 = vector from point 1 (x1,y1,z1) to point 2 (x2,y2,z2). N3 = vector from point 1 (x1,y1,z1) to point 3 (x3,y3,z3).
Directly translated from ANSI-standard FORTRAN 77 subroutines from [Blakely, 1995].

Parameters

x0 (scalar(float)) – x coordinate of point 0.
y0 (scalar(float)) – y coordinate of point 0.
z0 (scalar(float)) – z coordinate of point 0.
x1 (scalar(float)) – x coordinate of point 1.
y1 (scalar(float)) – y coordinate of point 1.
z1 (scalar(float)) – z coordinate of point 1.
x2 (scalar(float)) – x coordinate of point 2.
y2 (scalar(float)) – y coordinate of point 2.
z2 (scalar(float)) – z coordinate of point 2.
x3 (scalar(float)) – x coordinate of point 3.
y3 (scalar(float)) – y coordinate of point 3.
z3 (scalar(float)) – z coordinate of point 3.

Returns

x (scalar(float)) - x coordinate of the intersection point (plane and line)
y (scalar(float)) - y coordinate of the intersection point (plane and line)
z (scalar(float)) - z coordinate of the intersection point (plane and line)
r (scalar(float)) - distance between intersection point and point 0 (x0,y0,z0)

magnetostatics.qcoords_1D(nqpts)¶

Determines the coordinates of the quadrature points depending on the amount chosen.

Parameters

nqpts (scalar(int)) – the number of quadrature points per dimension.

Returns

coordsq (array_like(float)) - 1D array(nqpts) containing coordinates of quadrature points.

magnetostatics.qweights_1D(nqpts)¶

Determines the weights of the quadrature points depending on the amount chosen.

Parameters

nqpts (scalar(int)) – the number of quadrature points per dimension.

Returns

weightsq (array_like(float)) - 1D array(nqpts) containing weights of quadrature points.

magnetostatics.rot(ax, ay, az, bx, by, bz, nx, ny, nz, px, py, pz)¶

Function ROT finds the sense of rotation of the vector from (ax,ay,az) to (bx,by,bz) with respect to a second vector through point (px,py,pz). The second vector has components given by (nx,ny,nz).
Returned parameter s is 1 if anticlockwise, -1 if clockwise, or 0 if colinear.
Directly translated from ANSI-standard FORTRAN 77 subroutines from [Blakely, 1995].

Parameters

ax (scalar(float)) – x coordinate of point a, defines start of the first vector.
ay (scalar(float)) – y coordinate of point a, defines start of the first vector.
az (scalar(float)) – y coordinate of point a, defines start of the first vector.
bx (scalar(float)) – y coordinate of point b, defines end of the first vector.
by (scalar(float)) – y coordinate of point b, defines end of the first vector.
bz (scalar(float)) – z coordinate of point b, defines end of the first vector.
nx (scalar(float)) – x component of second vector.
ny (scalar(float)) – y component of second vector.
nz (scalar(float)) – z component of second vector.
px (scalar(float)) – x coordinate of point p, second vector passes through.
py (scalar(float)) – y coordinate of point p, second vector passes through.
pz (scalar(float)) – z coordinate of point p, second vector passes through.

Returns

s (scalar(int)) - parameter that defines sense of rotation of a vector w.r.t. a second vector

Set measurement parameters¶

set_measurement_parameters.set_measurement_parameters(rDEM, sDEM, site, path, ho)[source]¶

Returns measurement parameters from [Meyer and de Groot, 2024] of the chosen (passed) parameters.

Parameters

rDEM (scalar(float)) – number associated with resolution of the DEM, see Case study: Mount Etna.
sDEM (scalar(float)) – number associated with size of the DEM, see Case study: Mount Etna.
site (scalar(float)) – number associated with site (FLUX), see Case study: Mount Etna.
path (scalar(float)) – number associated with specific paths of a site, see Case study: Mount Etna.
ho (scalar(float)) – number associated with height option for the path, see Case study: Mount Etna.

Range

rDEM - for all sites except site 3: either 2 meter or 5 meter, for site 3, only 2 meter available.
sDEM - different for most sites/rDEM combinations.
site - [1:6].
path - for all sites except site 6: [1:3], for site 6, only 1 available.
ho - for all sites except site 6: [1:2], where 1 relates to 1 meter above the surface and 2 relates to 1.8 meter above the surface. For site 6, only 1 available.

After DEM choice, for each site one can choose the path and height above the surface, options for each site are given in the last two columns, irrespective of DEM choice, all can be selected.

site	rDEM	sDEM	~size [m] 1	path	ho (height)
1	2m	1	2100	1-3	1 (1m) 2 (1.8m)
1	2m	2	1500
1	2m	3	1100
1	2m	4	700
1	2m	5	500
1	2m	6	50 (0m around)
1	2m	7	4000
1	2m	8	6000

1	5m	1	2100
1	5m	2	300
1	5m	3	200
					2	5
					2	5
2/5	2m	n/a	2100	1-3	1 (1m) 2 (1.8m)	1 (25cm) 2 (75cm) 3 (125cm) 4 (175cm)

2/5	5m	1	2100
2/5	5m	2	300
2/5	5m	3	200


3	2m	1	2100	1-3	1 (1m) 2 (1.8m)
3	2m	2	300
3	2m	3	200


4/6	2m	n/a	2100

4/6	5m	1	2100

4	5m	2	200	1-3	1 (1m) 2 (1.8m)
4	5m	3	100	1-3	1 (1m) 2 (1.8m)

6	2m	2	1100	1	1 (1m)
6	5m	3	500
6	5m	4	300
6	5m	5	200 (50m arnd)
6	5m	6	100 (20m arnd)

Returns

Lx (scalar(float)) - length of domain in x-direction in meters.
Ly (scalar(float)) - length of domain in y-direction in meters.
Lz (scalar(float)) - length of domain in x-direction in meters.
nelx (scalar(int)) - amount of elements in x-direction of domain.
nely (scalar(int)) - amount of elements in y-direction of domain.
nelz (scalar(int)) - amount of elements in z-direction of domain.
xllcorner (scalar(float)) - x-coordinate of the DEM cut corner, used to align domain with field measurement coordinates in (WGS84) UTM Easting (meters), zone 33N.
yllcorner (scalar(float)) - y-coordinate of the DEM cut corner, used to align domain with field measurement coordinates in (WGS84) UTM Northing (meters), zone 33N.
npath (scalar(int)) - amount of measurement points for the field path from [Meyer and de Groot, 2024].
zpath_height (scalar(float)) - height of the field measurement path above the surface from [Meyer and de Groot, 2024] in meters.
pathfile (string) - name of the file associated with the field path from [Meyer and de Groot, 2024], containing all coordinates of measurement points and magnetic field strength measurement data (converted to UTM).
pathfile (string) - name of the file associated with the DEM cut chosen.
IGRFx (scalar(float)) - x-component of the IGRF (global reference field) in Tesla.
IGRFy (scalar(float)) - y-component of the IGRF (global reference field) in Tesla.
IGRFz (scalar(float)) - z-component of the IGRF (global reference field) in Tesla.

Artificial DEM¶

art_DEM.generate_fractal_map(size, roughness, sigma, seed=None)[source]¶

Generates a square fractal map using the diamond-square algorithm, followed by Gaussian filtering.

The diamond-square algorithm is a fractal method used to generate random landscapes or textures. The algorithm works by dividing the map into squares and diamonds, then randomly displacing the midpoint of each to create roughness. Gaussian filtering is applied at the end to smooth out the map for a more natural look.

Parameters

size (scalar(int)) – Determines the size of the map which will be (2^size) + 1.
roughness (scalar(float)) – The initial roughness factor for the diamond-square algorithm.
sigma (scalar(float)) – Standard deviation for Gaussian filter, affecting the smoothness.
seed ((optional(scalar(int)))) – Optional seed for the random number generator.

Returns

dem (array_like(float)) - A 2D numpy array representing the fractal map.

art_DEM.generate_pathfile(xllcorner, yllcorner, Lx, rough_length, npath, seed=None, filename='art_path.txt')[source]¶

Generate a path file for artifical data sampling paths in an ASCII format.

Parameters

xllcorner (scalar(float)) – The western x-coordinate of the lower left corner of the grid.
yllcorner (scalar(float)) – The southern y-coordinate of the lower left corner of the grid.
Lx (scalar(float)) – The length of the domain along the x-axis.
rough_length (scalar(float)) – The rough length of the path, minor changes will happen through this function, so length is eventually different to passed length.
npath (scalar(int)) – Number of points in the path.
seed (scalar(int)) – Optional seed for the random number generator.
filename (string) – The name of the text file to write the path coordinates.

art_DEM.write_dem_to_ascii(dem, ncol, nrow, cellsize, xllcorner, yllcorner, filename='art_dem.ascii')[source]¶

Writes the digital elevation model (DEM) data to an ASCII file in grid format.

The ASCII format is a simple text-based format for storing raster data which is widely used in GIS software. Each line contains the data for a single row of the grid, with values separated by spaces.

Parameters

dem (array_like(float)) – 2D array containing elevation data.
ncol (scalar(int)) – Number of columns in the DEM grid.
nrow – Number of rows in the DEM grid.
cellsize (scalar(float)) – Size of each cell in the grid.
xllcorner (scalar(float)) – The western x-coordinate of the lower left corner of the DEM grid.
yllcorner (scalar(float)) – The southern y-coordinate of the lower left corner of the DEM grid.
filename (string) – The name of the ASCII file to write.

Footnotes

1: These rough sizes are rounded off values, in case the DEM cut was not square, the value denotes the rounded length of the shortest side of the DEM.