Mathematical Functions

Mathematical utility functions for rotations, random number generation, statistics, and numerical solving.

double lagrange_exp(const double&, const double&, const double&)

double dlagrange_exp(const double&, const double&, const double&)

double lagrange_pow_0(const double&, const double&, const double&, const double&, const double&)

double dlagrange_pow_0(const double&, const double&, const double&, const double&, const double&)

double lagrange_pow_1(const double&, const double&, const double&, const double&, const double&)

double dlagrange_pow_1(const double&, const double&, const double&, const double&, const double&)

double d2lagrange_pow_1(const double&, const double&, const double&, const double&, const double&)

void Newton_Raphon(const arma::vec&, const arma::vec&, const arma::mat&, arma::vec&, arma::vec&, double&)

void Fischer_Burmeister(const arma::vec&, const arma::vec&, const arma::mat&, arma::vec&, arma::vec&, double&)

void Fischer_Burmeister_limits(const arma::vec&, const arma::vec&, const arma::vec&, const arma::mat&, const arma::mat&, arma::vec&, arma::vec&, double&)

void Fischer_Burmeister_m(const arma::vec&, const arma::vec&, const arma::mat&, arma::vec&, arma::vec&, double&)

void Fischer_Burmeister_m_limits(const arma::vec&, const arma::vec&, const arma::vec&, const arma::mat&, const arma::mat&, arma::vec&, arma::vec&, double&)

arma::mat denom_FB_m(const arma::vec&, const arma::mat&, const arma::vec&)

int alea(const int &a)

Generates a random integer between 0 and a-1.

Example:

int r = alea(10);  // Returns a random integer between 0 and 9

Parameters:: a – The upper bound (exclusive) (int)
Returns:: A random integer in the range [0, a-1] (int)

int aleab(const int &a, const int &b)

Generates a random integer between a and b.

Example:

int r = aleab(5, 15);  // Returns a random integer between 5 and 15

Parameters:

a – The lower bound (inclusive) (int)
b – The upper bound (inclusive) (int)

Returns:

A random integer in the range [a, b] (int)

double alead(const double &a, const double &b)

Generates a random double between a and b.

Example:

double r = alead(0.0, 1.0);  // Returns a random double between 0.0 and 1.0

Parameters:

a – The lower bound (double)
b – The upper bound (double)

Returns:

A random double in the range [a, b] (double)

arma::vec rotate_vec(const arma::vec &v, const arma::mat &R)

Rotates a 3D vector using a rotation matrix.

Example:

vec v = {1., 0., 0.};
mat R = fillR(M_PI/4., 3);
vec v_rot = rotate_vec(v, R);

Parameters:

v – The input 3D vector (arma::vec)
R – The 3x3 rotation matrix (arma::mat)

Returns:

The rotated vector (arma::vec)

arma::vec rotate_vec(const arma::vec &v, const double &angle, const int &axis)

Rotates a 3D vector by an angle around a given axis.

Example:

vec v = {1., 0., 0.};
vec v_rot = rotate_vec(v, M_PI/4., 3); // Rotate 45° around z-axis

Parameters:

v – The input 3D vector (arma::vec)
angle – The rotation angle in radians (double)
axis – The axis of rotation: 1=x, 2=y, 3=z (int)

Returns:

The rotated vector (arma::vec)

arma::mat rotate_mat(const arma::mat &m, const arma::mat &R)

Rotates a 3x3 matrix using a rotation matrix.

The rotation is performed as \( \mathbf{m}' = \mathbf{R} \cdot \mathbf{m} \cdot \mathbf{R}^T \)

mat m = eye(3,3);
mat R = fillR(M_PI/4., 3);
mat m_rot = rotate_mat(m, R);

Parameters:

m – The input 3x3 matrix (arma::mat)
R – The 3x3 rotation matrix (arma::mat)

Returns:

The rotated matrix (arma::mat)

arma::mat rotate_mat(const arma::mat &m, const double &angle, const int &axis)

Rotates a 3x3 matrix by an angle around a given axis.

Parameters:

m – The input 3x3 matrix (arma::mat)
angle – The rotation angle in radians (double)
axis – The axis of rotation: 1=x, 2=y, 3=z (int)

Returns:

The rotated matrix (arma::mat)

arma::mat fillR(const double &angle, const int &axis, const bool &active = true)

Generates a 3x3 rotation matrix for rotation around a single axis.

Example:

mat R = fillR(M_PI/4., 3); // 45° rotation around z-axis

Parameters:

angle – The rotation angle in radians (double)
axis – The axis of rotation: 1=x, 2=y, 3=z (int)
active – If true (default), active (alibi) rotation; if false, passive (alias) rotation (bool)

Returns:

The 3x3 rotation matrix (arma::mat)

arma::mat fillR(const double &psi, const double &theta, const double &phi, const bool &active = true, const std::string &conv = "zxz")

Generates a 3x3 rotation matrix using Euler angles.

Example:

mat R = fillR(M_PI/4., M_PI/6., M_PI/3., true, "zxz");

Parameters:

psi – First Euler angle in radians (double)
theta – Second Euler angle in radians (double)
phi – Third Euler angle in radians (double)
active – If true (default), active (alibi) rotation; if false, passive (alias) rotation (bool)
conv – The Euler angle convention: “zxz” (default), “zyz”, “xyz”, etc. Use “user” for custom convention (std::string)

Returns:

The 3x3 rotation matrix (arma::mat)

arma::mat fillQS(const double &angle, const int &axis, const bool &active = true)

Generates a 6x6 rotation matrix for stress tensors in Voigt notation.

Parameters:

angle – The rotation angle in radians (double)
axis – The axis of rotation: 1=x, 2=y, 3=z (int)
active – If true (default), active rotation (bool)

Returns:

The 6x6 rotation matrix for stress (arma::mat)

arma::mat fillQS(const arma::mat &R, const bool &active = true)

Generates a 6x6 rotation matrix for stress tensors from a 3x3 rotation matrix.

Parameters:

R – The 3x3 rotation matrix (arma::mat)
active – If true (default), active rotation (bool)

Returns:

The 6x6 rotation matrix for stress (arma::mat)

arma::mat fillQE(const double &angle, const int &axis, const bool &active = true)

Generates a 6x6 rotation matrix for strain tensors in Voigt notation.

Parameters:

angle – The rotation angle in radians (double)
axis – The axis of rotation: 1=x, 2=y, 3=z (int)
active – If true (default), active rotation (bool)

Returns:

The 6x6 rotation matrix for strain (arma::mat)

arma::mat fillQE(const arma::mat &R, const bool &active = true)

Generates a 6x6 rotation matrix for strain tensors from a 3x3 rotation matrix.

Parameters:

R – The 3x3 rotation matrix (arma::mat)
active – If true (default), active rotation (bool)

Returns:

The 6x6 rotation matrix for strain (arma::mat)

arma::mat rotateL(const arma::mat &L, const double &angle, const int &axis, const bool &active = true)

Rotates a 6x6 stiffness matrix.

Example:

mat L = L_iso(E, nu, "EnuG");
mat L_rot = rotateL(L, M_PI/4., 3);

Parameters:

L – The 6x6 stiffness matrix in Voigt notation (arma::mat)
angle – The rotation angle in radians (double)
axis – The axis of rotation: 1=x, 2=y, 3=z (int)
active – If true (default), active rotation (bool)

Returns:

The rotated 6x6 stiffness matrix (arma::mat)

arma::mat rotateL(const arma::mat &L, const arma::mat &R, const bool &active = true)

Rotates a 6x6 stiffness matrix using a 3x3 rotation matrix.

Parameters:

L – The 6x6 stiffness matrix in Voigt notation (arma::mat)
R – The 3x3 rotation matrix (arma::mat)
active – If true (default), active rotation (bool)

Returns:

The rotated 6x6 stiffness matrix (arma::mat)

arma::mat rotateM(const arma::mat &M, const double &angle, const int &axis, const bool &active = true)

Rotates a 6x6 compliance matrix.

Parameters:

M – The 6x6 compliance matrix in Voigt notation (arma::mat)
angle – The rotation angle in radians (double)
axis – The axis of rotation: 1=x, 2=y, 3=z (int)
active – If true (default), active rotation (bool)

Returns:

The rotated 6x6 compliance matrix (arma::mat)

arma::mat rotateM(const arma::mat &M, const arma::mat &R, const bool &active = true)

Rotates a 6x6 compliance matrix using a 3x3 rotation matrix.

Parameters:

M – The 6x6 compliance matrix in Voigt notation (arma::mat)
R – The 3x3 rotation matrix (arma::mat)
active – If true (default), active rotation (bool)

Returns:

The rotated 6x6 compliance matrix (arma::mat)

arma::mat rotateA(const arma::mat &A, const double &angle, const int &axis, const bool &active = true)

Rotates a 6x6 strain localization matrix A.

Parameters:

A – The 6x6 strain localization matrix (arma::mat)
angle – The rotation angle in radians (double)
axis – The axis of rotation: 1=x, 2=y, 3=z (int)
active – If true (default), active rotation (bool)

Returns:

The rotated 6x6 strain localization matrix (arma::mat)

arma::mat rotateA(const arma::mat &A, const arma::mat &R, const bool &active = true)

Rotates a 6x6 strain localization matrix A using a 3x3 rotation matrix.

Parameters:

A – The 6x6 strain localization matrix (arma::mat)
R – The 3x3 rotation matrix (arma::mat)
active – If true (default), active rotation (bool)

Returns:

The rotated 6x6 strain localization matrix (arma::mat)

arma::mat rotateB(const arma::mat &B, const double &angle, const int &axis, const bool &active = true)

Rotates a 6x6 stress localization matrix B.

Parameters:

B – The 6x6 stress localization matrix (arma::mat)
angle – The rotation angle in radians (double)
axis – The axis of rotation: 1=x, 2=y, 3=z (int)
active – If true (default), active rotation (bool)

Returns:

The rotated 6x6 stress localization matrix (arma::mat)

arma::mat rotateB(const arma::mat &B, const arma::mat &R, const bool &active = true)

Rotates a 6x6 stress localization matrix B using a 3x3 rotation matrix.

Parameters:

B – The 6x6 stress localization matrix (arma::mat)
R – The 3x3 rotation matrix (arma::mat)
active – If true (default), active rotation (bool)

Returns:

The rotated 6x6 stress localization matrix (arma::mat)

arma::vec rotate_stress(const arma::vec &sigma, const double &angle, const int &axis, const bool &active = true)

Rotates a stress vector in Voigt notation.

Example:

vec sigma = {100., 50., 0., 25., 0., 0.};
vec sigma_rot = rotate_stress(sigma, M_PI/4., 3);

Parameters:

sigma – The 6-component stress vector (arma::vec)
angle – The rotation angle in radians (double)
axis – The axis of rotation: 1=x, 2=y, 3=z (int)
active – If true (default), active rotation (bool)

Returns:

The rotated stress vector (arma::vec)

arma::vec rotate_stress(const arma::vec &sigma, const arma::mat &R, const bool &active = true)

Rotates a stress vector using a 3x3 rotation matrix.

Parameters:

sigma – The 6-component stress vector (arma::vec)
R – The 3x3 rotation matrix (arma::mat)
active – If true (default), active rotation (bool)

Returns:

The rotated stress vector (arma::vec)

arma::vec rotate_strain(const arma::vec &epsilon, const double &angle, const int &axis, const bool &active = true)

Rotates a strain vector in Voigt notation.

Parameters:

epsilon – The 6-component strain vector (arma::vec)
angle – The rotation angle in radians (double)
axis – The axis of rotation: 1=x, 2=y, 3=z (int)
active – If true (default), active rotation (bool)

Returns:

The rotated strain vector (arma::vec)

arma::vec rotate_strain(const arma::vec &epsilon, const arma::mat &R, const bool &active = true)

Rotates a strain vector using a 3x3 rotation matrix.

Parameters:

epsilon – The 6-component strain vector (arma::vec)
R – The 3x3 rotation matrix (arma::mat)
active – If true (default), active rotation (bool)

Returns:

The rotated strain vector (arma::vec)

arma::mat rotate_l2g_strain(const arma::vec &E, const double &psi, const double &theta, const double &phi)

Rotates a strain tensor from local to global coordinates using Euler angles.

Parameters:

E – The 6-component strain vector in local coordinates (arma::vec)
psi – First Euler angle in radians (double)
theta – Second Euler angle in radians (double)
phi – Third Euler angle in radians (double)

Returns:

The 3x3 strain tensor in global coordinates (arma::mat)

arma::mat rotate_g2l_strain(const arma::vec &E, const double &psi, const double &theta, const double &phi)

Rotates a strain tensor from global to local coordinates using Euler angles.

Parameters:

E – The 6-component strain vector in global coordinates (arma::vec)
psi – First Euler angle in radians (double)
theta – Second Euler angle in radians (double)
phi – Third Euler angle in radians (double)

Returns:

The 3x3 strain tensor in local coordinates (arma::mat)

arma::mat rotate_l2g_stress(const arma::vec &sigma, const double &psi, const double &theta, const double &phi)

Rotates a stress tensor from local to global coordinates using Euler angles.

Parameters:

sigma – The 6-component stress vector in local coordinates (arma::vec)
psi – First Euler angle in radians (double)
theta – Second Euler angle in radians (double)
phi – Third Euler angle in radians (double)

Returns:

The 3x3 stress tensor in global coordinates (arma::mat)

arma::mat rotate_g2l_stress(const arma::vec &sigma, const double &psi, const double &theta, const double &phi)

Rotates a stress tensor from global to local coordinates using Euler angles.

Parameters:

sigma – The 6-component stress vector in global coordinates (arma::vec)
psi – First Euler angle in radians (double)
theta – Second Euler angle in radians (double)
phi – Third Euler angle in radians (double)

Returns:

The 3x3 stress tensor in local coordinates (arma::mat)

arma::mat rotate_l2g_L(const arma::mat &L, const double &psi, const double &theta, const double &phi)

Rotates a stiffness matrix from local to global coordinates using Euler angles.

Parameters:

L – The 6x6 stiffness matrix in local coordinates (arma::mat)
psi – First Euler angle in radians (double)
theta – Second Euler angle in radians (double)
phi – Third Euler angle in radians (double)

Returns:

The 6x6 stiffness matrix in global coordinates (arma::mat)

arma::mat rotate_g2l_L(const arma::mat &L, const double &psi, const double &theta, const double &phi)

Rotates a stiffness matrix from global to local coordinates using Euler angles.

Parameters:

L – The 6x6 stiffness matrix in global coordinates (arma::mat)
psi – First Euler angle in radians (double)
theta – Second Euler angle in radians (double)
phi – Third Euler angle in radians (double)

Returns:

The 6x6 stiffness matrix in local coordinates (arma::mat)

arma::mat rotate_l2g_A(const arma::mat &A, const double &psi, const double &theta, const double &phi)

Rotates a strain localization matrix from local to global coordinates using Euler angles.

Parameters:

A – The 6x6 strain localization matrix in local coordinates (arma::mat)
psi – First Euler angle in radians (double)
theta – Second Euler angle in radians (double)
phi – Third Euler angle in radians (double)

Returns:

The 6x6 strain localization matrix in global coordinates (arma::mat)

arma::mat rotate_g2l_A(const arma::mat &A, const double &psi, const double &theta, const double &phi)

Rotates a strain localization matrix from global to local coordinates using Euler angles.

Parameters:

A – The 6x6 strain localization matrix in global coordinates (arma::mat)
psi – First Euler angle in radians (double)
theta – Second Euler angle in radians (double)
phi – Third Euler angle in radians (double)

Returns:

The 6x6 strain localization matrix in local coordinates (arma::mat)

arma::mat rotate_l2g_B(const arma::mat &B, const double &psi, const double &theta, const double &phi)

Rotates a stress localization matrix from local to global coordinates using Euler angles.

Parameters:

B – The 6x6 stress localization matrix in local coordinates (arma::mat)
psi – First Euler angle in radians (double)
theta – Second Euler angle in radians (double)
phi – Third Euler angle in radians (double)

Returns:

The 6x6 stress localization matrix in global coordinates (arma::mat)

arma::mat rotate_g2l_B(const arma::mat &B, const double &psi, const double &theta, const double &phi)

Rotates a stress localization matrix from global to local coordinates using Euler angles.

Parameters:

B – The 6x6 stress localization matrix in global coordinates (arma::mat)
psi – First Euler angle in radians (double)
theta – Second Euler angle in radians (double)
phi – Third Euler angle in radians (double)

Returns:

The 6x6 stress localization matrix in local coordinates (arma::mat)

arma::mat rotate_l2g_M(const arma::mat &M, const double &psi, const double &theta, const double &phi)

Rotates a compliance matrix from local to global coordinates using Euler angles.

Parameters:

M – The 6x6 compliance matrix in local coordinates (arma::mat)
psi – First Euler angle in radians (double)
theta – Second Euler angle in radians (double)
phi – Third Euler angle in radians (double)

Returns:

The 6x6 compliance matrix in global coordinates (arma::mat)

arma::mat rotate_g2l_M(const arma::mat &M, const double &psi, const double &theta, const double &phi)

Rotates a compliance matrix from global to local coordinates using Euler angles.

Parameters:

M – The 6x6 compliance matrix in global coordinates (arma::mat)
psi – First Euler angle in radians (double)
theta – Second Euler angle in radians (double)
phi – Third Euler angle in radians (double)

Returns:

The 6x6 compliance matrix in local coordinates (arma::mat)

arma::vec quadratic(const double&, const double&, const double&)

arma::cx_vec cx_quadratic(const double&, const double&, const double&)

arma::cx_vec cx_quadratic(const arma::cx_double&, const arma::cx_double&, const arma::cx_double&)

double normal_distrib(const double &x, const double &mu, const double &sigma)

Computes the value of a normal (Gaussian) distribution at a given point.

The normal distribution is given by:

\[ f(x) = \frac{1}{\sigma \sqrt{2\pi}} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right) \]

double p = normal_distrib(0.0, 0.0, 1.0);  // Standard normal at x=0

Parameters:

x – The point at which to evaluate the distribution (double)
mu – The mean of the distribution (double)
sigma – The standard deviation of the distribution (double)

Returns:

The probability density value at x (double)

double proba_distrib_weibull(const double &x, const double &alpha, const double &beta)

Computes the probability density of a Weibull distribution.

The Weibull probability density function is:

\[ f(x) = \frac{\alpha}{\beta} \left(\frac{x}{\beta}\right)^{\alpha-1} \exp\left(-\left(\frac{x}{\beta}\right)^\alpha\right) \]

double p = proba_distrib_weibull(1.5, 2.0, 1.0);

Parameters:

x – The point at which to evaluate the distribution (double)
alpha – The shape parameter \( \alpha > 0 \) (double)
beta – The scale parameter \( \beta > 0 \) (double)

Returns:

The probability density value at x (double)

double cumul_distrib_weibull(const double &x, const double &alpha, const double &beta)

Computes the cumulative distribution function of a Weibull distribution.

The Weibull cumulative distribution function is:

\[ F(x) = 1 - \exp\left(-\left(\frac{x}{\beta}\right)^\alpha\right) \]

double F = cumul_distrib_weibull(1.5, 2.0, 1.0);

Parameters:

x – The point at which to evaluate the CDF (double)
alpha – The shape parameter \( \alpha > 0 \) (double)
beta – The scale parameter \( \beta > 0 \) (double)

Returns:

The cumulative probability at x (double)

int tri_sum(const int &a, const int &b)

Computes the triangular sum of two integers.

Returns \( \frac{(a+b)(a+b+1)}{2} + b \), which provides a unique mapping from pairs of integers to a single integer.

Parameters:

a – First integer (int)
b – Second integer (int)

Returns:

The triangular sum (int)

double ODF_sd(const double &Theta, const double &Phi, const arma::vec &params)

Computes a classic Orientation Distribution Function (ODF).

The ODF is given by:

\[ \textrm{ODF}(\Theta) = a_1 \cos^{2p_1}(\Theta) + a_2 \cos^{2p_2+1}(\Theta) \sin^{2p_2}(\Theta) \]

Parameters:

Theta – The angle \( \Theta \) in radians (double)
Phi – The angle \( \Phi \) in radians (double)
params – Vector of parameters: \( [a_1, p_1, a_2, p_2] \) (arma::vec)

Returns:

The ODF value at the given angles (double)

double ODF_hard(const double &Theta, const double &Phi, const double &alpha, const double &beta)

Computes a hardening-type Orientation Distribution Function.

Parameters:

Theta – The angle \( \Theta \) in radians (double)
Phi – The angle \( \Phi \) in radians (double)
alpha – Parameter \( \alpha \) (double)
beta – Parameter \( \beta \) (double)

Returns:

The ODF value at the given angles (double)

double Gaussian(const double &x, const double &mu, const double &sigma, const double &ampl = 1.)

Computes a Gaussian peak function.

The Gaussian is given by:

\[ G(x) = A \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right) \]

double g = Gaussian(0.5, 0.0, 1.0, 2.0);

Parameters:

x – The position at which to evaluate (double)
mu – The center of the peak (double)
sigma – The width parameter (standard deviation) (double)
ampl – The amplitude (default = 1.0) (double)

Returns:

The Gaussian value at x (double)

double Lorentzian(const double &x, const double &x0, const double &gamma, const double &ampl = 1.)

Computes a Lorentzian (Cauchy) peak function.

The Lorentzian is given by:

\[ L(x) = \frac{A}{\pi} \frac{\gamma}{(x-x_0)^2 + \gamma^2} \]

double l = Lorentzian(0.5, 0.0, 1.0, 2.0);

Parameters:

x – The position at which to evaluate (double)
x0 – The center of the peak (double)
gamma – The half-width at half-maximum (double)
ampl – The amplitude (default = 1.0) (double)

Returns:

The Lorentzian value at x (double)

double PseudoVoigt(const double &x, const double &x0, const double &sigma, const double &gamma, const double &eta, const arma::vec &params = arma::ones(1))

Computes a Pseudo-Voigt peak function (weighted sum of Gaussian and Lorentzian).

The Pseudo-Voigt is given by:

\[ V(x) = \eta \cdot L(x) + (1-\eta) \cdot G(x) \]

Parameters:

x – The position at which to evaluate (double)
x0 – The center of the peak (double)
sigma – The Gaussian width parameter (double)
gamma – The Lorentzian width parameter (double)
eta – The mixing parameter (0 = pure Gaussian, 1 = pure Lorentzian) (double)
params – Optional additional parameters (arma::vec)

Returns:

The Pseudo-Voigt value at x (double)

double Pearson7(const double &x, const double &x0, const double &w, arma::vec &params)

Computes a Pearson VII peak function.

The Pearson VII function provides flexible peak shapes that can approximate both Gaussian and Lorentzian profiles.

Parameters:

x – The position at which to evaluate (double)
x0 – The center of the peak (double)
w – The width parameter (double)
params – Vector containing additional shape parameters (arma::vec)

Returns:

The Pearson VII value at x (double)