Source code for ggf.matlab_funcs

"""Special functions translated from Matlab to Python"""
import numpy as np
from scipy import special
from scipy import integrate


[docs]def besselh(n, z): """Hankel function with k = 1 Parameters ---------- n: int real order z: float complex argument Notes ----- https://de.mathworks.com/help/matlab/ref/besselh.html """ return special.hankel1(n, z)
[docs]def besselj(n, z): """Bessel function of first kind Parameters ---------- n: int real order z: float complex argument Notes ----- https://de.mathworks.com/help/matlab/ref/besselj.html """ return special.jv(n, z)
[docs]def gammaln(x): """Logarithm of the absolute value of the Gamma function Notes ----- https://de.mathworks.com/help/matlab/ref/gammaln.html See Also -------- scipy.special.gammaln """ return special.gammaln(x)
[docs]def legendre(n, x): """Associated Legendre functions Parameters ---------- n: int degree x: ndarray of floats argument Notes ----- x is treated always as a row vector The statement legendre(2,0:0.1:0.2) returns the matrix ======= ======= ======= ======= / x = 0 x = 0.1 x = 0.2 ======= ======= ======= ======= m = 0 -0.5000 -0.4850 -0.4400 m = 1 0 -0.2985 -0.5879 m = 2 3.0000 2.9700 2.8800 ======= ======= ======= ======= Notes ----- https://de.mathworks.com/help/matlab/ref/legendre.html """ x = np.real_if_close(x).flatten() # flattened x = np.array(x, dtype=float) # warning if not real result = np.zeros((x.shape[0], n + 1), dtype=complex) for ii in range(x.shape[0]): # Gives us row vector a = special.lpmn(n, n, x[ii])[0].transpose()[-1] result[ii] = a return result.transpose()
[docs]def lscov(A, B, w=None): """Least-squares solution in presence of known covariance :math:`A \\cdot x = B`, that is, :math:`x` minimizes :math:`(B - A \\cdot x)^T \\cdot \\text{diag}(w) \\cdot (B - A \\cdot x)`. The matrix :math:`w` typically contains either counts or inverse variances. Parameters ---------- A: matrix or 2d ndarray input matrix B: vector or 1d ndarray input vector Notes -------- https://de.mathworks.com/help/matlab/ref/lscov.html """ # https://stackoverflow.com/questions/27128688/how-to-use-least-squares-with-weight-matrix-in-python # https://de.mathworks.com/help/matlab/ref/lscov.html if w is None: Aw = A.copy() Bw = B.T.copy() else: W = np.sqrt(np.diag(np.array(w).flatten())) Aw = np.dot(W, A) Bw = np.dot(B.T, W) # set rcond=1e-10 to prevent diverging odd indices in x # (problem specific to ggf/stress computation) x, residuals, rank, s = np.linalg.lstsq(Aw, Bw.T, rcond=1e-10) return np.array(x).flatten()
[docs]def quadl(func, a, b): """Numerically evaluate integral with scipy QUADPACK quadrature Parameters ---------- func: callable function to integrate a: float interval start b: float interval end Notes ----- https://de.mathworks.com/help/matlab/ref/quadl.html """ def wrapreal(x): return func(x).real def wrapimag(x): return func(x).imag rp = integrate.quad(wrapreal, np.real_if_close(a), np.real_if_close(b))[0] ri = integrate.quad(wrapimag, np.real_if_close(a), np.real_if_close(b))[0] return rp + 1j * ri