seidart.routines.materials module¶

seidart.routines.materials.anisotropic_boolean(im: ndarray | list, matbool: ndarray | list, angvect: ndarray | list) → Tuple[ndarray, ndarray]¶

Determines if materials identified in an image are anisotropic and provides corresponding angular file names if applicable.

Parameters:

im (Union[np.ndarray, list]) – An array representing material IDs in an image.
matbool (Union[np.ndarray, list]) – Array indicating whether a material is anisotropic.
angvect (Union[np.ndarray, list]) – Array of angular file names for anisotropic materials.

Returns:

A tuple of two arrays; the first indicates anisotropy (boolean), the second contains angular file names.

Return type:

Tuple[np.ndarray, np.ndarray]

seidart.routines.materials.anisotropy_parameters(tensor: ndarray[tuple[Any, ...], dtype[floating]]) → Tuple[float, float, float]¶

Compute simple anisotropy measures from a symmetric 3x3 tensor.

Eigenvalues of the input tensor are used to derive a ratio, absolute birefringence, and normalized spread as scalar measures of anisotropy.

Parameters:: tensor (numpy.typing.NDArray[numpy.floating]) – Symmetric 3x3 tensor (e.g., permittivity or stiffness in a principal coordinate basis).
Returns:: A tuple (ratio, birefringence, spread) where: ratio is max(eigvals) / min(eigvals), birefringence is max(eigvals) - min(eigvals), and spread is birefringence / mean(eigvals).
Return type:: tuple[float, float, float]
Raises:: ValueError – If tensor is not 3x3 or if its smallest eigenvalue is effectively zero.

Note

These parameters provide a compact summary of anisotropy strength. A ratio close to 1 and small birefringence indicate near-isotropy, while larger values reflect stronger anisotropy.

# Example: anisotropy of a diagonal tensor
T = np.diag([1.0, 1.2, 1.5])
ratio, biref, spread = anisotropy_parameters(T)

seidart.routines.materials.astiffness2moduli(C: ndarray[tuple[Any, ...], dtype[floating]], transversely_isotropic: bool = True) → Tuple[float, float, ndarray[tuple[Any, ...], dtype[floating]]]¶

Compute isotropic bulk and shear moduli from an anisotropic stiffness tensor.

Hill (Voigt–Reuss) averages are used to derive effective isotropic bulk and shear moduli from a 6x6 anisotropic stiffness tensor. For transversely isotropic media, specialized approximate formulas are used; otherwise, general isotropic Voigt–Reuss expressions are applied.

Parameters:

C (numpy.typing.NDArray[numpy.floating]) – Anisotropic stiffness tensor in 6x6 Voigt notation (Pa).
transversely_isotropic (bool) – If True, use formulas tailored to transversely isotropic (TI) media. If False, use general isotropic Voigt–Reuss expressions.

Returns:

A tuple (K_iso, mu_iso, C_isotropic) where: K_iso is the effective isotropic bulk modulus (Pa), mu_iso is the effective isotropic shear modulus (Pa), and C_isotropic is the corresponding 6x6 isotropic stiffness tensor (Pa).

Return type:

tuple[float, float, numpy.typing.NDArray[numpy.floating]]

Raises:

ValueError – If C is not a 6x6 array or is singular.

Note

The Voigt bulk modulus is computed directly from stiffness components, while the Reuss bulk modulus is derived from the compliance matrix (inverse of C). The Hill bulk modulus is the arithmetic mean of these two bounds, and similarly for the shear modulus.

The resulting isotropic stiffness tensor is constructed from K_iso and mu_iso using standard isotropic relationships, with Lamé parameter \(\lambda = K_\mathrm{iso} - 2 \mu_\mathrm{iso} / 3\).

# Example: isotropic moduli from a TI stiffness tensor
C_ti = np.zeros((6, 6))
# ... fill C_ti with TI stiffness components ...
K_iso, mu_iso, C_iso = astiffness2moduli(C_ti, transversely_isotropic=True)

seidart.routines.materials.bond(R: ndarray) → ndarray¶

Calculates the 6x6 Bond transformation matrix from a 3x3 rotation matrix, useful for transforming stiffness or compliance matrices in crystallography and materials science.

Parameters:: R (np.ndarray) – The 3x3 rotation matrix.
Return M:: The 6x6 Bond transformation matrix.
Rtype M:: np.ndarray

seidart.routines.materials.bulk_modulus_air(T, P=101325)¶

Computes the bulk modulus of air (Pa) as a function of temperature (°C) and pressure (Pa), while accounting for variable air density.

Parameters:

T_C – Temperature in Celsius
P – Pressure in Pascals (default is atmospheric pressure, 101325 Pa)

Returns:

Bulk modulus in Pascals (Pa)

seidart.routines.materials.bulk_modulus_water(T: float) → Tuple[float, ndarray[tuple[Any, ...], dtype[floating]]]¶

Compute the bulk modulus of water as a function of temperature.

The bulk modulus is evaluated from empirical polynomial fits for liquid and supercooled water and returned in pascals. An equivalent isotropic stiffness matrix in 6x6 Voigt notation is also constructed with the bulk modulus in the normal components.

Parameters:: T (float) – Temperature in degrees Celsius. Valid range is approximately [-30, 100].
Returns:: A tuple (K, C) where: K is the bulk modulus in pascals and C is a 6x6 isotropic stiffness matrix in Voigt notation with normal components equal to K.
Return type:: tuple[float, numpy.typing.NDArray[numpy.floating]]
Raises:: ValueError – If the temperature is below -30 °C, outside the range covered by the empirical fits.

Note

For T >= 0 °C, a cubic polynomial is used for liquid water. For -30 <= T < 0 °C, a quadratic fit is used to represent the supercooled regime. The returned stiffness matrix C places K on the normal entries (0–2, 0–2) and zeros elsewhere.

seidart.routines.materials.compute_volume_fractions(porosity, lwc)¶

Computes volume fractions of ice, air, and water based on porosity and liquid water content.

Parameters:

porosity – Porosity fraction (0 to 1)
lwc – Liquid water content fraction (0 to 1)

Returns:

(V_ice, V_air, V_water)

seidart.routines.materials.fabric_tensor(euler_angles: ndarray[tuple[Any, ...], dtype[floating]]) → ndarray[tuple[Any, ...], dtype[floating]]¶

Compute a second-order fabric tensor from a set of orientations.

For each orientation given by Euler angles, the associated rotation matrix is formed and the (global) c-axis direction is extracted. The fabric tensor is then computed as the average of the dyadic products \(n \otimes n\) over all orientations.

Parameters:: euler_angles (numpy.typing.NDArray[numpy.floating]) – Array of Euler angle triplets with shape (M, 3), where each row contains (phi, theta, psi) in radians.
Returns:: Symmetric 3x3 fabric tensor.
Return type:: numpy.typing.NDArray[numpy.floating]
Raises:: ValueError – If euler_angles does not have shape (M, 3).

Note

The rotation matrix for each orientation is obtained via rotator_euler(). The unit vector corresponding to the local c-axis is taken as the third column of this matrix, and its outer product with itself is accumulated and averaged.

The resulting tensor summarizes the preferred orientation of the c-axis distribution and is often used as a measure of material anisotropy.

# Example: fabric tensor from random orientations
eulers = np.random.rand(100, 3)  # radians
F = fabric_tensor(eulers)

seidart.routines.materials.fujita_complex_permittivity(temperature: float, frequency: float) → float¶

Calculates the complex permittivity of ice using Fujita’s method, based on the provided temperature and frequency.

Parameters:

temperature (float) – The temperature of ice in degrees Celsius.
frequency (float) – The frequency at which to calculate permittivity, in Hz.

Returns:

The complex permittivity value for pure ice.

Return type:

float

seidart.routines.materials.get_christoffel_matrix(coefficients: Series, f0: float, direction: str | ndarray[tuple[Any, ...], dtype[floating]])¶

Build the Christoffel matrix and wave velocities for a given propagation direction.

The stiffness coefficients and density from a single row (e.g., a pandas.Series) are used to construct the 4th-order stiffness tensor. For a specified propagation direction, the Christoffel matrix is formed and its eigen-decomposition yields phase velocities and polarization vectors.

Parameters:

coefficients (pandas.Series) – Material properties for a single medium, stored as a Series with keys 'rho' and 'c11', 'c12', …, 'c66'.
f0 (float) – Reference frequency in Hz (currently unused, reserved for future frequency-dependent extensions).
direction (str or numpy.typing.NDArray[numpy.floating]) – Propagation direction, either as a string key ('x', 'y', 'z', 'xy', 'xz', 'yz') or as a 3-element NumPy array specifying the direction vector.

Returns:

A tuple (Gamma, eigenvalues, eigenvectors, velocities) where: Gamma is the 3x3 Christoffel matrix, eigenvalues are its eigenvalues, eigenvectors are the corresponding eigenvectors, and velocities are the (sorted) phase velocities obtained as the square roots of the eigenvalues (with negative values clipped to zero).

Return type:

tuple

Raises:

ValueError – If direction is an invalid string key, if the propagation vector is zero, or if required stiffness/density fields are missing.
TypeError – If direction is neither a recognized string nor a 3-element NumPy array.

Note

For string-valued directions, predefined unit vectors are used (e.g. 'x' maps to \([1, 0, 0]\), 'xy' maps to \([1, 1, 0]/\sqrt{2}\)). For array input, the vector is normalized to unit length before forming the Christoffel matrix.

The stiffness matrix in Voigt notation is first assembled from the cij fields, converted to a 3x3x3x3 tensor by voigt_to_full_tensor(), and then contracted with the propagation direction to produce the Christoffel matrix.

# Example: compute Christoffel matrix and velocities along x
row = df.iloc[0]  # pandas Series with c11..c66 and rho
Gamma, eigvals, eigvecs, velocities = get_christoffel_matrix(
    row,
    f0=10.0,
    direction="x",
)

seidart.routines.materials.get_complex_refractive_index(coefficients: Series, f0: float, direction, mu: float = 1.2566370614359173e-06)¶

Compute complex refractive indices for a given propagation direction.

Relative permittivity and (absolute) conductivity tensors are taken from a coefficient row, combined into a complex permittivity, and projected onto the plane perpendicular to the propagation direction. The resulting 2x2 operator is diagonalized to obtain complex refractive indices for the two transverse modes.

Parameters:

coefficients (pandas.Series) – Material properties for a single medium stored as a Series with entries 'e11', 'e12', 'e13', 'e22', 'e23', 'e33' for relative permittivity and 's11', 's12', 's13', 's22', 's23', 's33' for conductivity.
f0 (float) – Reference frequency in Hz.
direction (str or numpy.typing.NDArray[numpy.floating]) – Propagation direction, either as a string key ('x', 'y', 'z', 'xy', 'xz', 'yz') or as a 3-element NumPy array specifying the direction vector.
mu (float) – Magnetic permeability in H/m. Default is the vacuum permeability 4e-7 * pi. (Currently not used in the computation but included for completeness.)

Returns:

A tuple (lam, vecs, complex_index) where: lam is a 1D array of eigenvalues of the transverse operator, vecs is the corresponding 2x2 matrix of eigenvectors, and complex_index is a 1D array of complex refractive indices for the two transverse modes.

Return type:

tuple

Raises:

ValueError – If an invalid direction key is given, the propagation vector is zero, or required fields are missing.
TypeError – If direction is neither a recognized string nor a 3-element NumPy array.

Note

The complex permittivity is formed as \(\varepsilon_\mathrm{eff} = \varepsilon - i \, \sigma / (\omega \varepsilon_0)\) using relative permittivity components and conductivity, where \(\omega = 2\pi f_0\). The operator is then projected onto the plane perpendicular to the propagation direction and diagonalized.

Complex refractive indices are obtained as \(n = 1 / \sqrt{\lambda}\) for each eigenvalue \(\lambda\). Indices with non-positive real part are flipped to enforce a positive real component. Zero eigenvalues are mapped to NaN + 1j*NaN.

# Example: compute complex indices along x
row = df.iloc[0]  # pandas Series with e.. and s.. fields
lam, vecs, n_complex = get_complex_refractive_index(
    row,
    f0=10.0,
    direction="x",
)

seidart.routines.materials.get_perm(self, material)¶

Computes the permittivity and conductivity tensors for materials based on attributes contained within material and model class instances.

Parameters:: material (Material) – An instance of a class containing attributes for materials, including temperature, density, porosity, water content, and anisotropic properties.

seidart.routines.materials.get_seismic(self, material, snow_method='Hill', ice_method='Petrenko')¶

Calculates seismic stiffness coefficients based on material properties, accounting for anisotropic conditions where applicable.

Parameters:: material (Material) – The class object that contains material parameters.

seidart.routines.materials.highfreq_stability_conditions(T)¶

This is a check to make sure that no numerical instabilities or unreal wave modes are generated. See Bécache (2003) and Komatitsch (2007) for more details.

Parameters:: T (np.ndarray) – This is either a 2nd order (3-by-3) or 4th order (6-by-6) tensor
Returns:: The stability condition values
Return type:: tuple

seidart.routines.materials.ice_acidity_correction(temperature: float, frequency: float, concentration: float) → ndarray¶

Compute the real-permittivity correction tensor due to acidity in ice.

The scalar acidity increment Δε’_iso(T, f, C) from Fujita et al. (2000) is distributed isotropically and, optionally, anisotropically according to a fabric/structure tensor A and a scalar weighting beta:

Parameters:

temperature (float) – Ice temperature in degrees Celsius (≈ -40 to -10 °C for the fits).
frequency (float) – Radar frequency in Hz at which to evaluate the correction.
concentration (float) – Acid concentration C in molar units consistent with Fujita’s definition (mol acid per liter of ice, or equivalent).

Returns:

scalar Δε’ to be added to the pure-ice real permittivity.

Return type:

numpy.ndarray

seidart.routines.materials.ice_density(temperature: float, pressure: float = 0.0, method: str = None)¶

Compute the density of an ice crystal given the temperature.

Parameters:

temperature (float) – The value of the temperature in celsius
method – Specify which method to be used. Default is

seidart.routines.materials.ice_permittivity(temperature: float, density: float, center_frequency: float = None, method: str = 'fujita') → ndarray¶

Computes the complex permittivity of ice given its temperature, density, and the frequency of interest. Supports different methods of calculation.

Parameters:

temperature (float) – Temperature of the ice in degrees Celsius.
density (float) – Density of the ice in kg/m^3.
center_frequency (float, optional) – Frequency at which to compute permittivity, in Hz.
method (str) – The method used for calculating permittivity. Supports “kovacs” and “fujita”.

Returns:

The complex permittivity tensor for ice.

Return type:

np.ndarray

seidart.routines.materials.ice_stiffness_gagnon(temperature: float = None, pressure: float = 0.0) → ndarray¶

Computes the stiffness tensor for ice under specified temperature and pressure conditions based on empirical relationships.

Parameters:

temperature (float, optional) – The temperature at which to compute the stiffness tensor. Units are in Celsius.
pressure (float) – The pressure at which to compute the stiffness tensor. Units are in kbar.

Returns:

The stiffness tensor for ice.

Return type:

np.ndarray

seidart.routines.materials.ice_stiffness_petrenko(temperature: float, pressure=0) → ndarray¶

Calculate the 6×6 stiffness tensor (in Voigt notation) for ice Ih based on empirical formulas similar to those given in “Physics of Ice” by Petrenko & Whitworth.

Parameters:

temperature (float) – Temperature in °C.
pressure (float) – Pressure in kbar.

Returns:

6×6 stiffness tensor in Pa.

Return type:

np.ndarray

Note

The five independent elastic constants for hexagonal ice are C11, C12, C13, C33, and C44. The sixth constant, C66, is given by (C11 - C12)/2. The formulas below use linear (or quadratic) pressure and temperature derivatives. The coefficients here are only examples—please check with the actual values from “Physics of Ice” for your application.

seidart.routines.materials.isotropic_permittivity_tensor(temperature: float, porosity: float, water_content: float, material_name: str) → Tuple[ndarray, ndarray]¶

Computes the isotropic permittivity tensor for a material based on temperature, porosity, water content, and the material’s inherent properties.

Parameters:

temperature (float) – The temperature of the material.
porosity (float) – The porosity of the material.
water_content (float) – The water content in the material.
material_name (str) – The name of the material.

Returns:

A tuple containing the permittivity and conductivity tensors.

Return type:

Tuple[np.ndarray, np.ndarray]

seidart.routines.materials.isotropic_stiffness_tensor(pressure: float, density: float, material_limits: ndarray) → ndarray¶

Computes the isotropic stiffness tensor for a given material based on pressure, density, and predefined material properties.

Parameters:

pressure (float) – The hydrostatic pressure to which the material is subjected. Units are Pa.
density (float) – The density of the material.
material_limits (np.ndarray) – An array containing the material’s velocity limits and other relevant properties.

Returns:

The isotropic stiffness tensor for the material.

Return type:

np.ndarray

seidart.routines.materials.moduli2stiffness(K: float, G: float) → ndarray[tuple[Any, ...], dtype[floating]]¶

Construct an isotropic stiffness matrix in Voigt notation from bulk and shear moduli.

Given the bulk modulus \(K\) and shear modulus \(G\), this function returns the 6x6 stiffness matrix for an isotropic material in Voigt notation.

Parameters:

K (float) – Bulk modulus in consistent units (e.g., Pa).
G (float) – Shear modulus in consistent units (e.g., Pa).

Returns:

6x6 stiffness matrix in Voigt notation.

Return type:

numpy.typing.NDArray[numpy.floating]

Note

For isotropic media, the normal components satisfy \(C_{11} = C_{22} = C_{33} = K + 4G/3\) and \(C_{12} = C_{13} = C_{23} = K - 2G/3\), while the shear components are \(C_{44} = C_{55} = C_{66} = G\). Off-diagonal shear terms are zero.

# Example: build isotropic stiffness from K and G
K, G = 20e9, 10e9
C_iso = moduli2stiffness(K, G)

seidart.routines.materials.moduli_rock(rock_type: str, T: float, P: float) → Tuple[float, float]¶

Dispatch to rock-type-specific bulk and shear modulus models.

This function selects an empirical modulus model based on the specified rock type (e.g. basalt, granite, gneiss) and returns the corresponding bulk and shear moduli as functions of temperature and confining pressure.

Parameters:

rock_type (st) – Name of the rock type. Supported values include "basalt", "granite", and "gneiss". The comparison is case-insensitive.
T (float) – Temperature in degrees Celsius.
P (float) – Confining pressure in megapascals (MPa).

Returns:

Bulk modulus K and shear modulus G in pascals.

Return type:

tuple[float, float]

Raises:

ValueError – If rock_type is not recognized.

Note

This dispatcher computes the moduli directly using internal empirical relations for each supported rock type. Additional rock types can be added by extending the conditional branches with new parameterizations.

# Example: get moduli for basalt at 60 °C and 30 MPa
K, G = moduli_rock("basalt", T=60.0, P=30.0)

seidart.routines.materials.moduli_sand_silt_clay(P, T, lwc, porosity=0.0, composition={'clay': 0.0, 'sand': 1.0, 'silt': 0.0})¶

Computes the bulk modulus (K) and shear modulus (G) for a sand-silt-clay mixture based on pressure, temperature, and water content using Gassmann’s equation.

Parameters:

P – Overburden pressure in Pascals (Pa).
T – Temperature in degrees Celsius.
lwc – Water content fraction (0 to 1).
composition – Dictionary with fractions of sand, silt, and clay (must sum to 1).

Returns:

Bulk modulus (Pa), Shear modulus (Pa)

seidart.routines.materials.porewater_correction(temperature: float, density: float, porosity: float, liquid_water_content: float) → Tuple[float, float, float]¶

Applies corrections to the bulk density of a material based on its porosity and water content, considering temperature adjustments to the densities of air and water.

Parameters:

temperature (float) – The temperature of the material.
density (float) – The initial density of the material.
porosity (float) – The porosity percentage of the material.
liquid_water_content (float) – The percentage of the pore space filled with water.

Returns:

A tuple containing the corrected density, the density contribution from air, and the density contribution from water.

Return type:

Tuple[float, float, float]

Computes the hydrostatic pressure, temperature, and density at each grid point based on material ID, accounting for porosity and water content. The output units of pressure are in Pascals. Gammon uses kbar to compute the ice stiffness tensor.

Parameters:

im (Union[list, np.ndarray]) – An m-by-n array of integer values representing material IDs.
temp (Union[list, np.ndarray]) – Temperatures at each grid point.
rho (Union[list, np.ndarray]) – Densities at each grid point.
dz (Union[list, np.ndarray]) – Vertical spacing between grid points.
porosity (Union[list, np.ndarray], optional) – Porosities at each grid point, default is a list of zeros.
lwc (Union[list, np.ndarray], optional) – Liquid water contents at each grid point, default is a list of zeros.

Returns:

A tuple containing arrays for temperature, density, and hydrostatic pressure at each grid point.

Return type:

Tuple[np.ndarray, np.ndarray, np.ndarray]

seidart.routines.materials.read_ang(filepath: str) → ndarray¶

Reads Euler angles from a space delimited .ang file, typically associated with EBSD (Electron Backscatter Diffraction) data. Synthetic data can be built with the fabricsynth module.

Parameters:: filepath (str) – The path to the .ang file.
Returns:: An array of Euler angles extracted from the file.
Return type:: np.ndarray

Note

The .ang file is expected to contain columns for Euler angles in radians, following Bunge’s notation (z-x-z rotation), among other data related to EBSD measurements.

seidart.routines.materials.rho_water_correction(temperature: float = 0.0) → float¶

Corrects the density of water based on temperature using the empirical formula derived from the Kell equation for water >= 0C and Wagner et al. (1994) for supercoold water.

Parameters:: temperature (float) – The temperature at which to compute the water density correction.
Returns:: The corrected water density.
Return type:: float

seidart.routines.materials.rotation_matrix_to_euler_angles(rotation_matrix: ndarray[tuple[Any, ...], dtype[floating]]) → Tuple[float, float, float]¶

Convert a 3x3 rotation matrix to Euler angles.

The input orthogonal rotation matrix is converted to a set of three Euler angles (phi1, Phi, phi2) in radians, following the same convention used in the associated rotation routines.

Parameters:: rotation_matrix (numpy.typing.NDArray[numpy.floating]) – 3x3 orthogonal rotation matrix with determinant 1.
Returns:: A tuple (phi1, Phi, phi2) of Euler angles in radians.
Return type:: tuple[float, float, float]
Raises:: ValueError – If the input matrix is not 3x3 or is not a valid rotation matrix (orthogonal with determinant equal to 1).

Note

The matrix is first checked for orthogonality and unit determinant within numerical tolerance to ensure it represents a valid rotation. The angle extraction formula assumes a specific Euler-angle convention consistent with the rest of your rotation utilities.

# Example: recover Euler angles from a rotation matrix
R = rotator_euler(np.array([0.2, 0.3, 0.4]), rot="zxz")
phi1, Phi, phi2 = rotation_matrix_to_euler_angles(R)

seidart.routines.materials.rotation_matrix_to_trend_plunge_orientation(R: ndarray[tuple[Any, ...], dtype[floating]]) → Tuple[float, float, float]¶

Convert a rotation matrix to trend, plunge, and orientation angles.

This is the inverse of a trend–plunge–orientation based rotation construction. From a 3x3 rotation matrix, it recovers the trend and plunge of the transformed c-axis and the rotation (orientation) about that axis, all in radians.

Parameters:: R (numpy.typing.NDArray[numpy.floating]) – 3x3 rotation matrix.
Returns:: A tuple (trend, plunge, orientation) where all angles are in radians. trend is the azimuth of the c-axis measured in the horizontal plane, plunge is the angle of the c-axis from horizontal (positive downward), and orientation is the rotation about the c-axis.
Return type:: tuple[float, float, float]
Raises:: ValueError – If R is not a 3x3 matrix or if numerical issues prevent defining a unique orientation.

Note

The c-axis direction is obtained as the image of the local \((0, 0, 1)\) vector under R. The trend and plunge are then computed from this direction. The orientation is determined by projecting the transformed local x-axis into the plane perpendicular to the c-axis and measuring its angle relative to the similarly projected global x-axis.

The helper function rotation_matrix_to_euler_angles() is called internally, though its Euler angles are not used directly in the final trend–plunge–orientation values here.

# Example: recover trend, plunge, orientation from a rotation
R = rotator_euler(np.array([0.2, 0.3, 0.4]), rot="zxz")
trend, plunge, orientation = rotation_matrix_to_trend_plunge_orientation(R)

seidart.routines.materials.rotator(angle: float, axis: str) → ndarray[tuple[Any, ...], dtype[floating]]¶

Construct a 3x3 rotation matrix for a rotation about a coordinate axis.

The rotation is performed in 3D about one of the Cartesian axes ('x', 'y', or 'z') by the specified angle in radians.

Parameters:

angle (float) – Rotation angle in radians.
axis (str) – Axis of rotation; must be one of 'x', 'y', or 'z'.

Returns:

3x3 rotation matrix corresponding to the requested axis and angle.

Return type:

numpy.typing.NDArray[numpy.floating]

Raises:

ValueError – If axis is not one of 'x', 'y', or 'z'.

Note

The rotation matrices are defined in the standard right-handed Cartesian coordinate system. Angles follow the right-hand rule: a positive angle corresponds to a counter-clockwise rotation when looking along the positive axis toward the origin.

# Example: rotate a vector by 30 degrees about the z-axis
v = np.array([1.0, 0.0, 0.0])
R = rotator(np.deg2rad(30.0), axis="z")
v_rot = R @ v

seidart.routines.materials.rotator_euler(eul: ndarray[tuple[Any, ...], dtype[floating]], rot: str = 'zxz') → ndarray[tuple[Any, ...], dtype[floating]]¶

Generate a 3x3 rotation matrix from Euler angles for a given rotation sequence.

The three Euler angles are applied in the order specified by the rotation sequence string (e.g. "zxz", "zyx"), using rotator() to build the individual axis rotations. The default sequence is "zxz".

Parameters:

eul (numpy.typing.NDArray[numpy.floating]) – Array containing the three Euler angles in radians, in the order implied by rot.
rot (str) – Rotation sequence as a three-character string specifying the axes of rotation (e.g. "zxz", "xyx", "zyx"). Default is "zxz".

Returns:

3x3 rotation matrix derived from the Euler angles and rotation sequence.

Return type:

numpy.typing.NDArray[numpy.floating]

Raises:

ValueError – If rot does not contain exactly three axis characters.

Note

The order of rotations and the active/passive convention must be consistent with how you interpret the Euler angles elsewhere in your code. The individual axis rotations are constructed using rotator(), and combined as R = R3 @ R2 @ R1 where each Ri corresponds to one Euler angle and axis from rot.

# Example: build a zxz Euler rotation from three angles
angles = np.array([0.1, 0.2, 0.3])  # radians
R = rotator_euler(angles, rot="zxz")

seidart.routines.materials.sand_silt_clay_permittivity_conductivity(T, lwc, porosity, composition={'clay': 0.0, 'sand': 1.0, 'silt': 0.0})¶

Computes the relative permittivity (ε_r) and electrical conductivity (σ) for a sand-silt-clay mixture based on pressure, temperature, and water content.

Parameters:

P – Overburden pressure in Pascals (Pa).
T – Temperature in degrees Celsius.
lwc – Water content fraction (0 to 1).
composition – Dictionary with fractions of sand, silt, and clay (must sum to 1).

Returns:

Relative permittivity (dimensionless), Electrical conductivity (S/m)

seidart.routines.materials.scalar2tensor(scalar, A: ndarray = array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]), beta: float = 0.0) → ndarray¶

Lift a scalar into an anisotropic 3x3 tensor:

T = scalar * [ I + beta (A - I/3) ].

Parameters:

scalar – Isotropic magnitude.
A – 3x3 fabric/structure tensor (trace ≈ 1).
beta – Anisotropy weighting.

Returns:

3x3 tensor.

seidart.routines.materials.snow_conductivity(lwc: float = None, permittivity: ndarray = None, frequency: float = None, density: float = None, temperature: float = None, porosity: float = None) → ndarray¶

Computes the electrical conductivity of snow given its liquid water content (LWC), permittivity, and frequency of interest.

Parameters:

lwc (float, optional) – Liquid water content of the snow, optionally used if permittivity is not provided.
permittivity (np.ndarray, optional) – The complex permittivity of snow, if available.
frequency (float, optional) – Frequency at which conductivity is to be calculated, in Hz.

Returns:

The conductivity tensor for snow.

Return type:

np.ndarray

seidart.routines.materials.snow_permittivity(density: float = 917.0, temperature: float = 0.0, lwc: float = 0.0, porosity: float = 50.0, center_frequency: float = 10000000.0, method: str = 'shivola-tiuri') → ndarray¶

Calculates the complex permittivity of snow based on its density, temperature, liquid water content (LWC), porosity, and the chosen calculation method.

Parameters:

density (float) – Density of the snow in kg/m^3.
temperature (float) – Temperature of the snow in degrees Celsius.
lwc (float) – Liquid water content of the pore space of the snow in percentage.
porosity (float) – Porosity of the snow in percentage.
method (str) – The method to be used for calculating permittivity. Defaults to “shivola-tiuri”.

Returns:

The complex permittivity tensor for snow.

Return type:

np.ndarray

seidart.routines.materials.snow_stiffness(temperature, lwc, porosity, pressure=0.0, use_fabric=False, eulerangles=None, method=None, exponent=2, epsilon=1e-06, small_denominator=1e-12, ice_method='Petrenko')¶

Computes the effective stiffness tensor for snow using Hill, SCA, or Gassmann’s equations.

Parameters:

temperature – Temperature in Celsius.
lwc – Liquid water content (% of pore space filled with water).
porosity – Total porosity of the snow (% of volume that is pore space).
pressure – Pressure in Pascals (default: 0 Pa).
use_fabric – If True, fabric-based stiffness is used; otherwise, isotropic stiffness is computed.
eulerangles – (m x 3) array of Euler angles for fabric-based calculations.
method – Override homogenization method; options: “Hill”, “SCA”, “Gassmann”, “SCA-Gassmann”, etc.
exponent – Exponent for reducing K_ice (and optionally G_ice) by porosity.
epsilon – Small regularization for fluid shear terms in Hill or SCA-Hill mixing.
small_denominator – Threshold to clamp Gassmann’s denominator.

Returns:

(K_snow, G_snow, C_snow, rho_eff) - K_snow: Bulk modulus of snow (Pa). - G_snow: Shear modulus of snow (Pa). - C_snow: 6x6 effective stiffness tensor (Pa). - rho_eff: Effective density (kg/m³).

seidart.routines.materials.tensor2velocities(T: ndarray[tuple[Any, ...], dtype[floating]] | DataFrame, rho: float = 910.0, seismic: bool = True)¶

Compute characteristic velocities from elastic or dielectric tensors.

The input can be a Voigt-form tensor for seismic stiffness (6 or 21 components), a full matrix, or a DataFrame with named components. Depending on seismic, the function returns either seismic velocities (vpv, vph, vsv, vsh) or EM phase velocities (vx, vy, vz).

Parameters:

T (numpy.typing.NDArray[numpy.floating] or pandas.DataFrame) – Tensor representation. Accepted forms include a 1D array of length 6, a 1D array of length 21, a 3x3 or 6x6 full NumPy array, or a pandas.DataFrame with appropriate column names such as 'e11', 'e22', and 'e33' for EM data or 'c11', 'c12', 'c33', 'c44', and 'rho' for seismic data.
rho (float) – Density in kg/m^3 used for seismic velocity calculation. Ignored when a DataFrame with its own 'rho' column is provided. Default is 910.0.
seismic (bool) – If True, interpret T as a stiffness tensor and return seismic velocities. If False, interpret T as a dielectric tensor and return EM phase velocities.

Returns:

If seismic is True, returns a tuple (vpv, vph, vsv, vsh) (all in m/s). If seismic is False, returns a tuple (vx, vy, vz) of EM phase velocities (m/s) along the principal axes.

Return type:

tuple[float, float, float, float] or tuple[float, float, float]

Raises:

ValueError – If T has an unsupported shape or required fields are missing from the DataFrame.

Note

For 1D NumPy arrays, upper-triangular Voigt-style representations are expanded into full symmetric tensors before velocities are computed. For seismic tensors, the standard relationships \(v = \sqrt{C / \rho}\) are used for the relevant stiffness components. For dielectric tensors, the phase velocities are computed from the principal permittivities using \(v = c_0 \sqrt{1 / \varepsilon}\).

When a pandas.DataFrame is provided, field names such as 'e11', 'e22', 'e33' (EM) or 'c11', 'c12', 'c33', 'c44', and 'rho' (seismic) are expected.

# Example (seismic): velocities from a Voigt stiffness vector
T_voigt = np.array([
    c11, c12, c13, c22, c23, c33,
    c44, c45, c46, c55, c56,
    c66, c14, c15, c16, c24, c25, c26, c34, c35, c36
])  # length 21
vpv, vph, vsv, vsh = tensor2velocities(T_voigt, rho=2000.0, seismic=True)

# Example (EM): velocities from diagonal permittivity entries
eps_vec = np.array([eps_xx, eps_yy, eps_zz, 0.0, 0.0, 0.0])
vx, vy, vz = tensor2velocities(eps_vec, seismic=False)

seidart.routines.materials.trend_plunge_to_rotation_matrix(trend, plunge, orientation)¶

Compute the rotation matrix from a trend and plunge pair.

Parameters:

trend (float) – The trend value in degrees
plunge (float) – The plunge value in degrees
orientation (float) – The orientation/rotation angle around the axis of the plunge.

Returns:

rotation_matrix

Return type:

np.ndarray

seidart.routines.materials.voigt_to_full_tensor(C_voigt: ndarray[tuple[Any, ...], dtype[floating]]) → ndarray[tuple[Any, ...], dtype[floating]]¶

Convert a 6x6 stiffness matrix in Voigt notation to a 3x3x3x3 tensor.

The input stiffness matrix in 6x6 Voigt notation is expanded into the full 4th-order stiffness tensor \(C_{ijkl}\) using the standard Voigt index mapping and appropriate normalization factors for shear components.

Parameters:: C_voigt (numpy.typing.NDArray[numpy.floating]) – Stiffness matrix in Voigt notation with shape (6, 6).
Returns:: Full stiffness tensor with shape (3, 3, 3, 3).
Return type:: numpy.typing.NDArray[numpy.floating]
Raises:: ValueError – If C_voigt is not a 6x6 array.

Note

The mapping from index pairs \((i, j)\) to Voigt indices \(I\) follows the standard convention: (0,0)->0, (1,1)->1, (2,2)->2, (1,2)->3, (0,2)->4, (0,1)->5 (with symmetric pairs treated identically). Shear components include a factor of \(1/\sqrt{2}\) to ensure consistent energy normalization.

# Example: convert an isotropic stiffness matrix to full tensor
C_voigt = np.zeros((6, 6))
C_voigt[0, 0] = C_voigt[1, 1] = C_voigt[2, 2] = 10.0
C_voigt[3, 3] = C_voigt[4, 4] = C_voigt[5, 5] = 5.0

C_full = voigt_to_full_tensor(C_voigt)

seidart.routines.materials.vrh2(permittivity: ndarray[tuple[Any, ...], dtype[floating]], conductivity: ndarray[tuple[Any, ...], dtype[floating]], eulerangles: ndarray[tuple[Any, ...], dtype[floating]]) → Tuple[ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]]]¶

Compute Voigt-Reuss-Hill averages for 2nd-order permittivity and conductivity tensors.

The function takes a single anisotropic permittivity and conductivity tensor and a set of crystal/element orientations given by Euler angles. For each orientation, tensors are rotated, Voigt and Reuss bounds are formed, and their Hill (arithmetic) average is returned for both permittivity and conductivity.

Parameters:

permittivity (numpy.typing.NDArray[numpy.floating]) – Base permittivity tensor of shape (3, 3) in the material coordinate system.
conductivity (numpy.typing.NDArray[numpy.floating]) – Base conductivity tensor of shape (3, 3) in the material coordinate system.
eulerangles (numpy.typing.NDArray[numpy.floating]) – Array of Euler angle triplets with shape (P, 3), where P is the number of orientations. Angles are interpreted according to the convention used by rotator_euler().

Returns:

A tuple (eps_vrh, sigma_vrh) where: eps_vrh is the Voigt–Reuss–Hill averaged permittivity tensor of shape (3, 3), and sigma_vrh is the corresponding averaged conductivity tensor of shape (3, 3).

Return type:

tuple

Raises:

ValueError – If the input tensors are not 3x3 or if eulerangles does not have shape (P, 3).

Note

Voigt averages are formed by rotating the tensors into the laboratory frame for each orientation and averaging directly. Reuss averages are formed by averaging the inverses (compliance tensors) in the rotated frames and then inverting the result. The Voigt–Reuss–Hill estimate is the arithmetic mean of these two bounds.

The rotation matrices are obtained from rotator_euler(), which must be consistent with the Euler angle convention used to describe the orientations.

# Example: VRH average over a set of orientations
eps = np.eye(3) * 5.0
sigma = np.eye(3) * 1e-3
eulers = np.array([
    [0.0, 0.0, 0.0],
    [30.0, 45.0, 60.0],
])  # units must match rotator_euler

eps_vrh, sigma_vrh = vrh2(eps, sigma, eulers)

seidart.routines.materials.vrh4(C: ndarray[tuple[Any, ...], dtype[floating]], eulerangles: ndarray[tuple[Any, ...], dtype[floating]]) → ndarray[tuple[Any, ...], dtype[floating]]¶

Compute Voigt–Reuss–Hill estimates for a 4th-order tensor in Voigt notation.

The input stiffness tensor in 6x6 Voigt notation is averaged over a set of orientations given by Euler angles. For each orientation, the tensor is rotated, Voigt and Reuss bounds are formed, and their Hill (arithmetic) average is returned as an effective stiffness tensor.

Parameters:

C (numpy.typing.NDArray[numpy.floating]) – Stiffness tensor in 6x6 Voigt notation, with shape (6, 6).
eulerangles (numpy.typing.NDArray[numpy.floating]) – Array of Euler angle triplets with shape (M, 3), where M is the number of orientations. Angles must follow the same convention as used in rotator_euler().

Returns:

Voigt–Reuss–Hill averaged stiffness tensor in 6x6 Voigt notation with shape (6, 6).

Return type:

numpy.typing.NDArray[numpy.floating]

Raises:

ValueError – If C is not 6x6 or eulerangles does not have shape (M, 3).

Note

Voigt averages are obtained by rotating the stiffness tensor into the laboratory frame for each orientation using the Bond transformation matrix \(M\) and averaging the results. Reuss averages are obtained by averaging the rotated compliance tensors (inverses of stiffness) using the inverse Bond transform \(N\), and then inverting the result.

The Hill estimate is defined as the arithmetic mean of the Voigt and Reuss bounds. The functions rotator_euler() and bond() must be consistent in their use of Euler angles and coordinate conventions.

# Example: VRH averaging of a stiffness tensor over orientations
C = np.eye(6) * 1e9  # simple isotropic-like placeholder
eulers = np.array([
    [0.0, 0.0, 0.0],
    [30.0, 45.0, 60.0],
])  # units must match rotator_euler

C_vrh = vrh4(C, eulers)

seidart.routines.materials.wolff_conductivity(temperature_C: float, concentrations: dict, activation_energies: dict = {'a0': 21000.0, 'aCl': 21000.0, 'aH': 21000.0, 'aNH4': 21000.0}, a_ref: dict = {'a0': 9.0, 'aCl': 0.55, 'aH': 4.0, 'aNH4': 1.0})¶

Wolff-style scalar conductivity σ(T, C) from ions H+, NH4+, Cl-.

concentrations:: dict with keys ‘H’, ‘NH4’, ‘Cl’ giving concentrations in µeq/L.
activation_energies:: dict with keys ‘a0’, ‘aH’, ‘aNH4’, ‘aCl’ giving Ea in J/mol. If None, a shared Ea = 21 kJ/mol is used for all coefficients.

Returns:: σ in S/m (SI units).

seidart.routines.materials.wolff_fujita_complex_permittivity(temperature_C, frequency, density, acidity_M: float, concentrations, A: ndarray = array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]), beta: float = 1.0, a_ref: dict = None, activation_energies: dict = None, T_ref: float = 263.15)¶

Combine Wolff conductivity and Fujita permittivity corrections for ice.

The function computes a fabric-aware real-permittivity tensor, an acidity correction tensor, and an imaginary-permittivity tensor derived from ionic conductivity.

Parameters:

temperature_C (float) – Ice temperature in degrees Celsius.
frequency (float) – Radar frequency in Hz.
density (float) – Ice density in kg/m^3.
acidity_M (float) – Acid concentration used by ice_acidity_correction().
concentrations (dict) – Ion concentrations for wolff_conductivity(), typically with keys "H", "NH4", and "Cl".
A (numpy.ndarray) – Fabric tensor used to distribute scalar corrections.
beta (float) – Fabric anisotropy weighting.
a_ref (dict, optional) – Reference molar conductivity coefficients.
activation_energies (dict, optional) – Activation energies for each conductivity coefficient.
T_ref (float) – Reference temperature in Kelvin.

Returns:

Tuple (eps_real_tensor, eps_imag_tensor).

Return type:

tuple[numpy.ndarray, numpy.ndarray]