import math
from datetime import datetime
import numpy as np
from scipy import linalg as spla
from scipy.spatial.transform import Rotation
from apsg.feature._geodata import Fault, Foliation, Lineation, Pair
from apsg.helpers._math import atand
from apsg.math._matrix import Matrix3
from apsg.math._vector import Vector3
[docs]
class Rotation3(DeformationGradient3):
"""
The class to represent 3D rotation matrix.
Args:
a (3x3 array_like): Input data, that can be converted to
3x3 2D array. This includes lists, tuples and ndarrays.
Examples:
>>> R = rotation.from_axisangle(lin(120, 60), 50)
"""
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
if not np.allclose(np.dot(np.transpose(self), self), np.eye(3)):
raise TypeError("Not valid arguments for Rotation3")
# fix improper rotations
if np.allclose(-1, np.linalg.det(self)):
U, S, Vt = np.linalg.svd(self)
# Ensure a proper rotation (det=1)
coefs = U @ np.diag([1, 1, np.linalg.det(U @ Vt)]) @ Vt
self._coefs = tuple(coefs[0]), tuple(coefs[1]), tuple(coefs[2])
[docs]
def axisangle(self):
"""Return rotation as (axis, angle) tuple."""
rotvec = Rotation.from_matrix(self).as_rotvec(degrees=True)
sign = 1.0
# if rotvec[2] < 0:
# sign = -1.0
return sign * Vector3(rotvec).uv(), sign * np.linalg.norm(rotvec)
[docs]
def angle(self):
"""Return rotation angle."""
rotvec = Rotation.from_matrix(self).as_rotvec(degrees=True)
sign = 1.0
if rotvec[2] < 0:
sign = -1.0
return sign * np.linalg.norm(rotvec)
[docs]
def euler(self, seq):
"""Return rotation as Euler angles specified in degrees.
Note: Each character in seq defines one axis around which angles turns.
Args:
seq (str): sequence of axes for rotations. Up to 3 characters belonging to the set
{'X', 'Y', 'Z'} for intrinsic rotations, or {'x', 'y', 'z'} for extrinsic rotations.
Returns:
ndarray: rotation as Euler angles specified in degrees.
"""
return Rotation.from_matrix(self).as_euler(seq, degrees=True)
[docs]
def quat(self, scalar_first=False):
"""Return rotation as quaternion.
Keyword Args:
scalar_first (bool): Whether the scalar component goes first or last. Default is False
Returns:
ndarray: rotation as quaternion.
"""
return Rotation.from_matrix(self).as_quat(canonical=True)
[docs]
@classmethod
def from_pair(cls, p):
"""Return ``Rotation3`` representing rotation defined by ``Pair``.
Rotation bring x-axis to lineation and z-axis to normal to plane.
Args:
p (``Pair``): Pair object
Examples:
>>> p = pair(40, 20, 75, 16)
>>> F = defgrad.from_pair(p)
Returns:
Rotation3: ``Rotation3`` representing rotation defined by ``Pair``.
"""
try:
p = Pair(p)
except Exception:
raise TypeError("Unsupported argument for from_pair. Expecting Pair")
return cls(
np.asarray(
[
np.asarray(p.lvec),
np.asarray(p.fvec.cross(p.lvec)),
np.asarray(p.fvec),
]
).T
)
[docs]
@classmethod
def from_axisangle(cls, vector, theta):
"""Return ``Rotation3`` representing rotation around axis.
Args:
vector: Rotation axis as ``Vector3`` like object
theta: Angle of rotation in degrees
Examples:
>>> F = rotation.from_axisangle(lin(120, 30), 45)
Returns:
Rotation3: ``Rotation3`` representing rotation around axis.
"""
try:
vector = Vector3(vector)
except Exception:
raise TypeError(
"Unsupported argument for from_axisangle. Expecting Vector3"
)
rotvec = theta * np.array(vector.uv())
return cls(Rotation.from_rotvec(rotvec, degrees=True).as_matrix())
[docs]
@classmethod
def from_two_vectors(cls, v1, v2):
"""Return ``Rotation3`` representing rotation around axis perpendicular
to both vectors and rotate v1 to v2.
Args:
v1: ``Vector3`` like object
v2: ``Vector3`` like object
Examples:
>>> F = rotation.from_two_vectors(lin(120, 30), lin(210, 60))
Returns:
Rotation3: ``Rotation3`` representing rotation around axis perpendicular to both vectors.
"""
try:
v1 = Vector3(v1)
except Exception:
raise TypeError(
"Unsupported first argument for from_two_vectors. Expecting Vector3"
)
try:
v2 = Vector3(v2)
except Exception:
raise TypeError(
"Unsupported second argument for from_two_vectors. Expecting Vector3"
)
return cls.from_axisangle(v1.cross(v2), v1.angle(v2))
[docs]
@classmethod
def from_vectors_axis(cls, v1, v2, a):
"""
Return ``Rotation3`` representing rotation of vector v1 to v2 around
axis a.
If v1.angle(a) is not equal to v2.angle(b), the minimum adjustment of rotation
axis is done automatically.
Args:
v1: ``Vector3`` like object
v2: ``Vector3`` like object
a: estimated rotation axis ``Vector3`` like object
Examples:
>>> v1 = lin(130, 49)
>>> v2 = lin(209, 77)
>>> a = lin(30, 30)
>>> R = rotation.from_vectors_axis(v1, v2, a)
>>> v1.transform(R) == v2
True
>>> a_fix, theta = R.axisangle()
>>> lin(a_fix)
L:31/30
Returns:
Rotation3: ``Rotation3`` representing rotation of vector v1 to v2 around axis a.
"""
try:
v1 = Vector3(v1)
except Exception:
raise TypeError(
"Unsupported first argument for from_vectors_axis. Expecting Vector3"
)
try:
v2 = Vector3(v2)
except Exception:
raise TypeError(
"Unsupported second argument for from_vectors_axis. Expecting Vector3"
)
try:
a = Vector3(a)
except Exception:
raise TypeError(
"Unsupported third argument for from_vectors_axis. Expecting Vector3"
)
n = v1.cross(v2).cross(v1.slerp(v2, 0.5))
a_fix = a.reject(n).normalized()
v1p = v1.reject(a_fix)
v2p = v2.reject(a_fix)
return cls.from_axisangle(a_fix, v1p.angle(v2p))
[docs]
@classmethod
def from_two_pairs(cls, p1, p2, symmetry=False):
"""
Return ``Rotation3`` representing rotation of coordinates from system
defined by ``Pair`` p1 to system defined by ``Pair`` p2.
Lineation in pair define x axis and normal to foliation in pair define z axis.
Args:
p1 (``Pair``): from
p2 (``Pair``): to
Keyword Args:
symmetry (bool): If True, returns minimum angle rotation of axial pairs
Examples:
>>> p1 = pair(58, 36, 81, 34)
>>> p2 = pair(217,42, 162, 27)
>>> R = rotation.from_two_pairs(p1, p2)
>>> p1.transform(R) == p2
True
Returns:
Rotation3: ``Rotation3`` representing rotation of coordinates from system defined by ``Pair``.
"""
if symmetry:
R4 = [
cls(cls.from_pair(Pair(p2.fvec, p2.lvec)) @ cls.from_pair(p1).I),
cls(cls.from_pair(Pair(-p2.fvec, p2.lvec)) @ cls.from_pair(p1).I),
cls(cls.from_pair(Pair(p2.fvec, -p2.lvec)) @ cls.from_pair(p1).I),
cls(cls.from_pair(Pair(-p2.fvec, -p2.lvec)) @ cls.from_pair(p1).I),
]
axes, angles = zip(*[R.axisangle() for R in R4])
angles = [abs(a) for a in angles]
ix = angles.index(min(angles))
return R4[ix]
else:
return cls(cls.from_pair(p2) @ cls.from_pair(p1).I)
[docs]
@classmethod
def from_declination(cls, lat, lon, year=None, alt=0):
"""
Return ``Rotation3`` representing rotation of coordinates correcting
magnetic declination at given coordinates and given time.
Args:
lat (float): latitude
lon (float): longitude
Keyword Args:
year (float): decimal year
alt (float): altitude in km
Examples:
>>> R = rotation.from_declination(48.6, 13.2, alt=0.6)
>>> f = fol(20, 48)
>>> f.transform(R)
S:25/48
Returns:
Rotation3: ``Rotation3`` representing rotation of coordinates correcting magnetic declination.
"""
from pygeomag import GeoMag
geo_mag = GeoMag(high_resolution=True)
if year is None:
year = datetime.now().year + datetime.now().month / 12
result = geo_mag.calculate(
glat=lat, glon=lon, alt=0, time=year, allow_date_outside_lifespan=True
)
return cls.from_axisangle(Lineation(0, 90), result.d)
[docs]
@classmethod
def from_quat(cls, quat, scalar_first=False):
"""
Return ``Rotation3`` representing rotation of coordinates created
from unit norm quaternion.
Args:
quat (array_like): quaternion
Keyword Args:
scalar_first (bool): Whether the scalar component goes first or last. Default is False
Examples:
>>> q = [-0.11543715, 0.19994301, 0.39988603, 0.88701083]
>>> R = rotation.from_quat(q)
>>> f = fol(20, 48)
>>> f.transform(R)
S:82/23
Returns:
Rotation3: ``Rotation3`` representing rotation of coordinates created from unit norm quaternion.
"""
return cls(Rotation.from_quat(quat, scalar_first=False).as_matrix())
[docs]
@classmethod
def from_euler(cls, seq, angles):
"""
Return ``Rotation3`` representing rotation of coordinates created
from unit norm quaternion.
Args:
seq (str): sequence of axes for rotations. Up to 3 characters belonging to the set
{'X', 'Y', 'Z'} for intrinsic rotations, or {'x', 'y', 'z'} for extrinsic rotations.
angles (array_like): Euler angles specified in degrees. Each character in seq
defines one axis around which angles turns.
Examples:
>>> R = rotation.from_euler('zxz', [30,-64, 125])
>>> f = fol(20, 48)
>>> f.transform(R)
S:74/70
Returns:
Rotation3: ``Rotation3`` representing rotation of coordinates created from Euler angles.
"""
return cls(Rotation.from_euler(seq, angles, degrees=True).as_matrix())
[docs]
class VelocityGradient3(Matrix3):
"""
The class to represent 3D velocity gradient tensor.
Args:
a (3x3 array_like): Input data, that can be converted to
3x3 2D array. This includes lists, tuples and ndarrays.
Examples:
>>> L = velgrad(np.diag([0.1, 0, -0.1]))
"""
[docs]
@classmethod
def from_comp(cls, **kwargs):
"""Return ``VelocityGradient3`` defined by individual components. Default is zero
tensor.
Keyword Args:
xx (float): tensor component L_xx
xy (float): tensor component L_xy
xz (float): tensor component L_xz
yx (float): tensor component L_yx
yy (float): tensor component L_yy
yz (float): tensor component L_yz
zx (float): tensor component L_zx
zy (float): tensor component L_zy
zz (float): tensor component L_zz
Examples:
>>> L = velgrad.from_comp(xy=1, zy=-0.5)
>>> L
[[ 0. 1. 0. ]
[ 0. 0. 0. ]
[ 0. -0.5 0. ]]
Returns:
VelocityGradient3: ``VelocityGradient3`` defined by individual components. Default is zero tensor.
"""
xx = kwargs.get("xx", 0)
xy = kwargs.get("xy", 0)
xz = kwargs.get("xz", 0)
yx = kwargs.get("yx", 0)
yy = kwargs.get("yy", 0)
yz = kwargs.get("yz", 0)
zx = kwargs.get("zx", 0)
zy = kwargs.get("zy", 0)
zz = kwargs.get("zz", 0)
return cls([[xx, xy, xz], [yx, yy, yz], [zx, zy, zz]])
[docs]
def defgrad(self, time=1, steps=1):
"""
Return ``DeformationGradient3`` tensor accumulated after given time.
Keyword Args:
time (float): time of deformation. Default 1
steps (int): when bigger than 1, will return a list
of ``DeformationGradient3`` tensors for each timestep.
Returns:
DeformationGradient3: ``DeformationGradient3`` tensor accumulated after given time.
"""
if steps > 1: # FIX once container for matrix will be implemented
return [
DeformationGradient3(spla.expm(np.asarray(self) * t))
for t in np.linspace(0, time, steps)
]
else:
return DeformationGradient3(spla.expm(np.asarray(self) * time))
[docs]
def rate(self):
"""Return rate of deformation tensor."""
return type(self)((self + self.T) / 2)
[docs]
def spin(self):
"""Return spin tensor."""
return type(self)((self - self.T) / 2)
class Tensor3(Matrix3):
@property
def _eig(self):
if "eig" not in self._cache:
evals, evecs = np.linalg.eigh(np.asarray(self))
idx = evals.argsort()[::-1]
evals = evals[idx]
evals[np.isclose(evals, np.zeros_like(evals))] = 0
evecs = evecs[:, idx]
self._cache["eig"] = evals, evecs
return self._cache["eig"]
def eigenlins(self, which=None):
"""Return eigenvectors as ``Lineation`` objects.
Args:
which: if None returns sorted tuple of eigenlins.
If int returns given eigen lineation. Default None.
Returns:
tuple of Lineation: eigenvectors as Lineation objects.
"""
if which is None:
return tuple(Lineation(v) for v in self.eigenvectors())
else:
return Lineation(self.eigenvectors(which))
def eigenfols(self, which=None):
"""Return tuple of eigenvectors as ``Foliation`` objects.
Args:
which: if None returns sorted tuple of eigenfols.
If int returns given eigen foliation. Default None.
Returns:
tuple of Foliation: eigenvectors as Foliation objects.
"""
if which is None:
return tuple(Foliation(v) for v in self.eigenvectors())
else:
return Foliation(self.eigenvectors(which))
@property
def pair(self):
"""Return ``Pair`` representing orientation of principal axes."""
ev = self.eigenvectors()
return Pair(ev[2], ev[0])
[docs]
class Stress3(Tensor3):
"""
The class to represent 3D stress tensor.
The real eigenvalues of the stress tensor are what we call
the principal stresses. There are 3 of these in 3D, available
as properties E1, E2, and E3 in descending order of magnitude
(max, intermediate, and minimum principal stresses) with orientations
available as properties V1, V2 and V3. The minimum principal stress
is simply the eigenvalue that has the lowest magnitude. Therefore,
the maximum principal stress is the most tensile (least compressive)
and the minimum principal stress is the least tensile (most compressive).
Tensile normal stresses have positive values, and compressive normal
stresses have negative values. If the maximum principal stress is <=0
and the minimum principal stress is negative then the stresses are
completely compressive.
Note: Stress tensor has a special properties sigma1, sigma2 and sigma3
to follow common geological terminology. sigma1 is most compressive
(least tensile) while sigma3 is most tensile (least compressive).
Their orientation could be accessed with properties sigma1dir,
sigma2dir and sigma3dir.
Args:
a (3x3 array_like): Input data, that can be converted to
3x3 2D array. This includes lists, tuples and ndarrays.
Examples:
>>> S = stress([[-8, 0, 0],[0, -5, 0],[0, 0, -1]])
"""
[docs]
@classmethod
def from_comp(cls, **kwargs):
"""
Return ``Stress`` tensor. Default is zero tensor.
Note that stress tensor is always symmetrical.
Keyword Args:
xx, xy|yx, xz|zx, yy, yz|zy, zz (float): tensor components
Examples:
>>> S = stress.from_comp(xx=-5, yy=-2, zz=10, xy=1)
>>> S
Stress3
[[-5. 1. 0.]
[ 1. -2. 0.]
[ 0. 0. 10.]]
Returns:
Stress3: ``Stress`` tensor. Default is zero tensor.
"""
xx = kwargs.get("xx", 0)
xy = kwargs.get("xy", kwargs.get("yx", 0))
xz = kwargs.get("xz", kwargs.get("zx", 0))
yy = kwargs.get("yy", 0)
yz = kwargs.get("yz", kwargs.get("zy", 0))
zz = kwargs.get("zz", 0)
return cls([[xx, xy, xz], [xy, yy, yz], [xz, yz, zz]])
[docs]
@classmethod
def from_ratio(cls, r=0.5, mag=1):
"""
Return ``Stress`` tensor with given shape ration.
Keyword Args:
r (float): shape ratio between 0 and 1. Default 0.5
mag (float): magnitude of differential stress. Default 1.
Examples:
>>> S = stress.from_ratio(r=0.25, mag=10)
>>> S
Stress3
[[-5. 0. 0. ]
[ 0. -2.5 0. ]
[ 0. 0. 5. ]]
Returns:
Stress3: ``Stress`` tensor with given shape ration.
"""
xx = -mag / 2
yy = xx + r * mag
zz = mag / 2
return cls([[xx, 0, 0], [0, yy, 0], [0, 0, zz]])
@property
def mean_stress(self):
"""Mean stress."""
return self.I1 / 3
@property
def hydrostatic(self):
"""Mean hydrostatic stress tensor component."""
return type(self)(np.diag(self.mean_stress * np.ones(3)))
@property
def deviatoric(self):
"""A stress deviator tensor component."""
return type(self)(self - self.hydrostatic)
[docs]
def effective(self, fp):
"""
Return effective stress tensor reduced by fluid pressure.
Args:
fp (flot): fluid pressure
Returns:
Stress3: effective stress tensor reduced by fluid pressure.
"""
return type(self)(self + fp * Stress3())
@property
def sigma1(self):
"""A maximum principal stress (max compressive)."""
return self.E3
@property
def sigma2(self):
"""A intermediate principal stress."""
return self.E2
@property
def sigma3(self):
"""A minimum principal stress (max tensile)."""
return self.E1
@property
def sigma1dir(self):
"""Return unit length vector in direction of maximum."""
return self.V3
@property
def sigma2dir(self):
"""Return unit length vector in direction of intermediate."""
return self.V2
@property
def sigma3dir(self):
"""Return unit length vector in direction of minimum."""
return self.V1
@property
def sigma1vec(self):
"""Return maximum principal stress vector (max compressive)."""
return self.E3 * self.V3
@property
def sigma2vec(self):
"""Return intermediate principal stress vector."""
return self.E2 * self.V2
@property
def sigma3vec(self):
"""Return minimum principal stress vector (max tensile)."""
return self.E1 * self.V1
@property
def I1(self):
"""First invariant."""
return float(np.trace(self))
@property
def I2(self):
"""Second invariant."""
return float((self.I1**2 - np.trace(self**2)) / 2)
@property
def I3(self):
"""Third invariant."""
return self.det
@property
def diagonalized(self):
"""Returns diagonalized Stress tensor and orthogonal matrix R, which transforms
coordinate system to the principal one."""
return (
type(self)(np.diag(self.eigenvalues())),
DeformationGradient3(self.eigenvectors()),
)
[docs]
def cauchy(self, n):
"""
Return stress vector associated with plane given by normal vector.
Args:
n: normal given as ``Vector3`` or ``Foliation`` object
Examples:
>>> S = stress.from_comp(xx=-5, yy=-2, zz=10, xy=1)
>>> S.cauchy(fol(160, 30))
Vector3(-2.52, 0.812, 8.66)
Returns:
Vector3: stress vector associated with plane given by normal vector.
"""
return Vector3(np.dot(self, n.normalized()))
[docs]
def fault(self, n):
"""
Return ``Fault`` object derived from given by normal vector.
Args:
n: normal given as ``Vector3`` or ``Foliation`` object
Examples:
>>> S = stress.from_comp(xx=-5, yy=-2, zz=10, xy=8)
>>> S.fault(fol(160, 30))
F:160/30-141/29 +
Returns:
Fault: ``Fault`` object derived from given by normal vector.
"""
sn, tau = self.stress_comp(n)
# return Fault(sn.normalized(), tau.normalized())
return Fault(n.normalized(), -tau.normalized())
[docs]
def stress_comp(self, n):
"""Return normal and shear stress ``Vector3`` components on plane given
by normal vector.
Returns:
tuple: normal and shear stress ``Vector3`` components.
"""
t = self.cauchy(n)
sn = t.proj(n)
return sn, t - sn
[docs]
def normal_stress(self, n):
"""Return normal stress magnitude on plane given by normal vector.
Returns:
float: normal stress magnitude on plane given by normal vector.
"""
return float(np.dot(n, self.cauchy(n)))
[docs]
def shear_stress(self, n):
"""Return shear stress magnitude on plane given by normal vector.
Returns:
float: shear stress magnitude on plane given by normal vector.
"""
sn, tau = self.stress_comp(n)
return abs(tau)
[docs]
def slip_tendency(self, n, fp=0, log=False):
"""
Return slip tendency calculated as the ratio of shear stress
to normal stress acting on the plane.
Note: Providing fluid pressure effective normal stress is calculated.
Keyword Args:
fp (float): fluid pressure. Default 0
log (bool): when True, returns logarithm of slip tendency
Returns:
float: slip tendency calculated as the ratio of shear stress to normal stress.
"""
Se = self.effective(fp)
sn, tau = Se.stress_comp(n)
if log:
return np.log(abs(tau) / abs(sn))
else:
return abs(tau) / abs(sn)
[docs]
def dilation_tendency(self, n, fp=0):
"""
Return dilation tendency of the plane.
Note: Providing fluid pressure effective stress is used.
Keyword Args:
fp (float): fluid pressure. Default 0
Returns:
float: dilation tendency of the plane.
"""
Se = self.effective(fp)
sn, tau = Se.stress_comp(n)
denom = Se.sigma1 - Se.sigma3
return np.where(np.isclose(denom, 0), np.nan, (Se.sigma1 - abs(sn)) / denom)
@property
def shape_ratio(self):
"""Return shape ratio R (Gephart & Forsyth 1984)."""
return float((self.sigma1 - self.sigma2) / (self.sigma1 - self.sigma3))
[docs]
class Ellipsoid(Tensor3):
"""
The class to represent 3D ellipsoid.
See following methods and properties for additional operations.
Args:
matrix (3x3 array_like): Input data, that can be converted to
3x3 2D matrix. This includes lists, tuples and ndarrays.
Examples:
>>> E = ellipsoid([[8, 0, 0], [0, 2, 0], [0, 0, 1]])
>>> E
Ellipsoid
[[8 0 0]
[0 2 0]
[0 0 1]]
(S1:2.83, S2:1.41, S3:1)
"""
def __repr__(self) -> str:
return (
f"{Matrix3.__repr__(self)}\n"
f"(S1:{self.S1:.3g}, S2:{self.S2:.3g}, S3:{self.S3:.3g})"
)
[docs]
@classmethod
def from_defgrad(cls, F, form="left", **kwargs) -> "Ellipsoid":
"""
Return deformation tensor from ``Defgrad3``.
Args:
F: DeformationGradient3 tensor
form: 'left' or 'B' for left Cauchy–Green (Finger) deformation tensor,
'right' or 'C' for right Cauchy–Green (Green's) deformation tensor.
Default 'left'.
Returns:
Ellipsoid: deformation tensor from ``Defgrad3``.
"""
if form in ("left", "B"):
return cls(np.dot(F, np.transpose(F)), **kwargs)
elif form in ("right", "C"):
return cls(np.dot(np.transpose(F), F), **kwargs)
else:
raise TypeError("Wrong form argument")
[docs]
@classmethod
def from_stretch(cls, x=1, y=1, z=1, **kwargs) -> "Ellipsoid":
"""Return diagonal tensor defined by magnitudes of principal stretches."""
return cls([[x * x, 0, 0], [0, y * y, 0], [0, 0, z * z]], **kwargs)
@property
def kind(self) -> str:
"""Return descriptive type of ellipsoid."""
nu = self.lode
if np.allclose(self.eoct, 0):
res = "O"
elif nu < -0.75:
res = "L"
elif nu > 0.75:
res = "S"
elif nu < -0.15:
res = "LLS"
elif nu > 0.15:
res = "SSL"
else:
res = "LS"
return res
@property
def strength(self) -> float:
"""Return the Woodcock strength."""
return self.e13
@property
def shape(self) -> float:
"""Return the Woodcock shape."""
return self.K
@property
def S1(self) -> float:
"""Return the maximum principal stretch."""
return math.sqrt(self.E1)
@property
def S2(self) -> float:
"""Return the middle principal stretch."""
return math.sqrt(self.E2)
@property
def S3(self) -> float:
"""Return the minimum principal stretch."""
return math.sqrt(self.E3)
@property
def e1(self) -> float:
"""Return the maximum natural principal strain."""
return math.log(self.S1)
@property
def e2(self) -> float:
"""Return the middle natural principal strain."""
return math.log(self.S2)
@property
def e3(self) -> float:
"""Return the minimum natural principal strain."""
return math.log(self.S3)
@property
def Rxy(self) -> float:
"""Return the Rxy ratio."""
return self.S1 / self.S2 if self.S2 != 0 else float("inf")
@property
def Ryz(self) -> float:
"""Return the Ryz ratio."""
return self.S2 / self.S3 if self.S3 != 0 else float("inf")
@property
def e12(self) -> float:
"""Return the e1 - e2."""
return self.e1 - self.e2
@property
def e13(self) -> float:
"""Return the e1 - e3."""
return self.e1 - self.e3
@property
def e23(self) -> float:
"""Return the e2 - e3."""
return self.e2 - self.e3
@property
def k(self) -> float:
"""Return the strain symmetry."""
return (self.Rxy - 1) / (self.Ryz - 1) if self.Ryz > 1 else float("inf")
@property
def d(self) -> float:
"""Return the strain intensity."""
return math.sqrt((self.Rxy - 1) ** 2 + (self.Ryz - 1) ** 2)
@property
def K(self) -> float:
"""Return the strain symmetry (Ramsay, 1983)."""
return self.e12 / self.e23 if self.e23 > 0 else float("inf")
@property
def D(self) -> float:
"""Return the strain intensity."""
return self.e12**2 + self.e23**2
@property
def r(self) -> float:
"""Return the strain intensity (Watterson, 1968)."""
return self.Rxy + self.Ryz - 1
@property
def goct(self) -> float:
"""Return the natural octahedral unit shear (Nadai, 1963)."""
return 2 * math.sqrt(self.e12**2 + self.e23**2 + self.e13**2) / 3
@property
def eoct(self) -> float:
"""Return the natural octahedral unit strain (Nadai, 1963)."""
return math.sqrt(3) * self.goct / 2
@property
def lode(self) -> float:
"""Return Lode parameter (Lode, 1926)."""
return (
(2 * self.e2 - self.e1 - self.e3) / (self.e1 - self.e3)
if (self.e1 - self.e3) > 0
else 0
)
@property
def P(self) -> float:
"""Point index (Vollmer, 1990)."""
return self.E1 - self.E2
@property
def G(self) -> float:
"""Girdle index (Vollmer, 1990)."""
return 2 * (self.E2 - self.E3)
@property
def R(self) -> float:
"""Random index (Vollmer, 1990)."""
return 3 * self.E3
@property
def B(self) -> float:
"""Cylindricity index (Vollmer, 1990)."""
return self.P + self.G
@property
def Intensity(self) -> float:
"""Intensity index (Lisle, 1985)."""
return 7.5 * float(np.sum((np.array(self.eigenvalues()) - 1 / 3) ** 2))
@property
def MAD_l(self) -> float:
"""Return maximum angular deviation (MAD) of linearly distributed vectors."""
return float(atand(np.sqrt((self.E2 + self.E3) / self.E1)))
@property
def MAD_p(self) -> float:
"""Return maximum angular deviation (MAD) of planarly distributed vectors."""
return float(atand(np.sqrt(self.E3 / self.E2 + self.E3 / self.E1)))
@property
def MAD(self) -> float:
"""Return maximum angular deviation (MAD)."""
if self.shape > 1:
return self.MAD_l
else:
return self.MAD_p
[docs]
class OrientationTensor3(Ellipsoid):
"""
Represents an 3D orientation tensor, which characterize data distribution
using eigenvalue method. See (Watson 1966, Scheidegger 1965).
See following methods and properties for additional operations.
Args:
matrix (3x3 array_like): Input data, that can be converted to
3x3 2D matrix. This includes lists, tuples and ndarrays.
Array could be also ``Group`` (for backward compatibility)
Examples:
>>> ot = ortensor([[8, 0, 0], [0, 2, 0], [0, 0, 1]])
>>> ot
OrientationTensor3
[[8 0 0]
[0 2 0]
[0 0 1]]
(S1:2.83, S2:1.41, S3:1)
"""
[docs]
@classmethod
def from_features(cls, g) -> "OrientationTensor3":
"""
Return ``Ortensor`` of data in Vector3Set of features.
Args:
g (Vector3Set): Set of features
Examples:
>>> g = linset.random_fisher(position=lin(120,50))
>>> ot = ortensor.from_features(g)
>>> ot
OrientationTensor3
[[ 0.126 -0.149 -0.202]
[-0.149 0.308 0.373]
[-0.202 0.373 0.566]]
(S1:0.955, S2:0.219, S3:0.2)
>>> ot.eigenlins()
(L:119/51, L:341/31, L:237/21)
Returns:
OrientationTensor3: orientation tensor of data in Vector3Set of features.
"""
axes = np.array(g)
norms = np.linalg.norm(axes, axis=1, keepdims=True)
norms[norms == 0] = 1.0
unit_axes = axes / norms
return cls(np.dot(unit_axes.T, unit_axes) / len(unit_axes))
[docs]
@classmethod
def from_pairs(cls, p, shift=True) -> "OrientationTensor3":
"""
Return Lisle (1989) ``Ortensor`` of orthogonal data in ``PairSet``.
Lisle, R. (1989). The Statistical Analysis of Orthogonal Orientation Data.
The Journal of Geology, 97(3), 360-364.
Note: Tensor is by default shifted towards positive eigenvalues, so it
could be used as Scheidegger orientation tensor for plotting. When
original Lisle tensor is needed, set shift to False.
Args:
p: ``PairSet``
Keyword Args:
shift (bool): When True the tensor is shifted. Default True
Examples:
>>> p = pairset([pair(109, 82, 21, 10),
pair(118, 76, 30, 11),
pair(97, 86, 7, 3),
pair(109, 75, 23, 14)])
>>> ot = ortensor.from_pairs(p)
>>> ot
OrientationTensor3
[[0.577 0.192 0.029]
[0.192 0.092 0.075]
[0.029 0.075 0.332]]
(S1:0.807, S2:0.579, S3:0.114)
Returns:
OrientationTensor3: Lisle (1989) orientation tensor of orthogonal data in ``PairSet``.
"""
if shift:
return cls(
(
OrientationTensor3.from_features(p.lvec)
- OrientationTensor3.from_features(p.fvec)
+ np.eye(3)
)
/ 3
)
else:
return cls(
OrientationTensor3.from_features(p.lvec)
- OrientationTensor3.from_features(p.fvec)
)