from typing import Callable, Optional, Tuple, Iterable import numpy as np from scipy.integrate import quad from scipy.special import hankel1, lpmv, jv, factorial from gorkov import log # Often used NumPy and SciPy Functions and Classes conj = np.conj exp = np.exp real = np.real sqrt = np.sqrt sin = np.sin cos = np.cos pi = np.pi def full_range(*args: int) -> Iterable: """Returns the range including the end point. :param args: end point, start point (optional) :type args: int :return: range including the end point :raises: ValueError if more than 2 arguments are passed :rtype: Iterable """ if len(args) == 1: return range(args[0] + 1) elif len(args) == 2: return range(args[0], args[1] + 1) else: raise ValueError('full_range takes one or two arguments.') def cartesian_2_polar_coordinates(x: float, z: float) -> Tuple[float, float]: """Computes from inputs the polar coordinates and returns them. :param x: x coordinate :type x: float :param z: z coordinate :type z: float :rtype: Tuple[float, float] """ r = sqrt(x**2 + z**2) theta = np.arctan2(x, z) return r, theta def polar_2_cartesian_velocities(v_r: float, v_theta: float, theta: float) -> Tuple[float, float]: """Computes from inputs velocities :attr:`v_r`, :attr:`v_theta` and the angle :attr:`theta` in the cartesian velocity components :param v_r: velocity in radial direction :type v_r: float :param v_theta: velocity in angle direction :type v_theta: float :param theta: angle :type theta: float :rtype: Tuple[float, float] velocity component in `x` and `z` direction """ v_x = sin(theta) * v_r + cos(theta) * v_theta v_z = cos(theta) * v_r - sin(theta) * v_theta return v_x, v_z def Clebsch_Gordan_coefficient(j1: int, m1: int, j2: int, m2: int, j: int, m: int) -> float: """ Clebsch-Gordan coefficient for integer valued arguments. :param j: :type j: int :param m: :type m: int :param j1: :type j1: int :param m1: :type m1: int :param j2: :type j2: int :param m2: :type m2: int :rtype: int """ for arg in [j1, m1, j2, m2, j, m]: if not isinstance(arg, int): raise ValueError('Arguments need to be integer valued.') if (m1 + m2) != m: return 0 elif j1 + j2 < j or abs(j1 - j2) > j: return 0 else: c = 0.0 z = 0 a = (factorial(j1 + j2 - j) * factorial(j1 - j2 + j) * factorial(-j1 + j2 + j) / factorial(j1 + j2 + j + 1.0))**0.5 b = (factorial(j1 + m1) * factorial(j1 - m1) * factorial(j2 + m2) * factorial(j2 - m2) * factorial(j + m) * factorial(j - m))**0.5 while z < (j1 - m1 + 3): denominator = (factorial(z) * factorial(j1 + j2 - j - z) * factorial(j1 - m1 - z) * factorial(j2 + m2 - z) * factorial(j - j2 + m1 + z) * factorial(j - j1 - m2 + z)) numerator = (-1)**float(z + j1 - j2 + m) if denominator != 0: c += numerator / denominator z += 1 return (-1.0)**(j1 - j2 + m) * (2.0 * j + 1.0)**0.5 * a * b * c def integrate(func: Callable[[float], float], lower: float, upper: float, rel_eps: Optional[float] = 1e-3) -> Tuple[float, float]: """ Integrates :attr:`func` between :attr:`lower` and :attr:`upper`. This methods wraps Scipy's QUADPACK implementation `quad() `__ :param func: function to be integrated. If no rel_eps is passed the function needs to be vectorized. :type func: Callable[[float], float] :param lower: lower integration limit :type lower: float :param upper: upper integration limit :type upper: float :param rel_eps: estimated relative error of integration :type rel_eps: Optional[float] :return: returns the definite integral of func for the boundaries lower and upper. :rtype: float """ # Integration integral, abs_err = quad(func, lower, upper, epsabs=0, epsrel=rel_eps, limit=100) log.debug(f'Integral = {integral},' f'absolute error estimation = {abs_err}\n') return integral, abs_err def integrate_osc(func: Callable[[float], float], lower: float, upper: float, viscous_wavelength, boundary_layer, resolution=4, roi_factor=5, rel_eps: Optional[float] = 1e-6 ) -> Tuple[float, float]: """ Special purpose integrator for oscillating integrals that integrates :attr:`func` between :attr:`lower` and :attr:`upper`. This methods wraps Scipy's QUADPACK implementation `quad() `__ Used primarily for boundary layer integration in highly elastic fluids. The methods creates a list of uniformly distributes points where local difficulties are expected. The density of points per viscous wavelength is given by :attr:`resolution`, the range is given by [:attr:`lower`, :attr:`lower` + :attr:`roi_factor`* :attr:`boundary_layer`] :param func: function to be integrated. If no rel_eps is passed the function needs to be vectorized. :type func: Callable[[float], float] :param lower: lower integration limit :type lower: float :param upper: upper integration limit :type upper: float :param viscous_wavelength: viscous wavelength in the fluid [m] :type viscous_wavelength: float :param boundary_layer: boundary layer thickness in the fluid [m] :type boundary_layer: float :param resolution: number of integration points per wavelength passed :type resolution: float :param roi_factor: multiplier for the boundary layer (see above) :type roi_factor: float :param rel_eps: estimated relative error of integration :type rel_eps: Optional[float] :return: returns the definite integral of func for the boundaries lower and upper. :rtype: float """ # Generate points dx = viscous_wavelength / resolution region_of_interest = roi_factor * boundary_layer n = int(region_of_interest / dx) points = np.linspace(lower, lower + region_of_interest, n) # Test if subdivision limit is large enough limit = 10 * n if 10 * n > 50 else 50 # Integrate integral, abs_err = quad(func, lower, upper, epsabs=0, epsrel=rel_eps, points=points, limit=limit) log.debug(f'Integral = {integral},' f'absolute error estimation = {abs_err}\n') return integral, abs_err class SpecialFunctions(object): """ This class gathers all reoccurring spherical hankel and bessel functions necessary for the gorkov package. """ @staticmethod def _a0_d2(n: int, z: complex) -> complex: return (n**2 - n - z**2) / z**2 @staticmethod def _a1_d2(n: int, z: complex) -> complex: return 2 / z @staticmethod def _a0_d3(n: int, z: complex) -> complex: return (n**3 - 3 * n**2 + 2 * n) / z**3 + (2 - n) / z @staticmethod def _a1_d3(n: int, z: complex) -> complex: return (-n**2 - n - 6) / z**2 + 1 @staticmethod def _a0_d4(n: int, z: complex) -> complex: return (n * (n - 1) * (6 + n * (n - 5)) / z ** 4 + (-2 * n * (n - 1) - 8) / z ** 2 + 1) @staticmethod def _a1_d4(n: int, z: complex) -> complex: return 8 * (3 + n + n**2) / z**3 - 4 / z @staticmethod def _a0_d5(n: int, z: complex) -> complex: return ((-4 + n) * (-3 + n) * (-2 + n) * (-1 + n) * n - 2 * (-20 + n * (2 + (-7 + n) * n)) * z**2 + (-4 + n) * z**4) / z**5 @staticmethod def _a1_d5(n: int, z: complex) -> complex: return -((n * (1 + n) * (58 + n + n**2) - 2 * n * (1 + n) * z**2 + z**4 - 20 * (-6 + z**2)) / z**4) @staticmethod def jn(n: int, z: complex) -> complex: """Order :attr:`n` Spherical Hankel function of the first kind :param n: Order :type n: int :param z: argument :type z: complex :rtype: complex """ return sqrt(pi / (2 * z)) * jv(n + 0.5, z) @classmethod def d_jn(cls, n: int, z: complex) -> complex: """First derivative of order :attr:`n` Spherical Bessel function of the first kind :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ return -cls.jn(n + 1, z) + (n / z) * cls.jn(n, z) @classmethod def d2_jn(cls, n: int, z: complex) -> complex: """Second derivative of order :attr:`n` Spherical Bessel function of the first kind :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ return (cls._a0_d2(n, z) * cls.jn(n, z) + cls._a1_d2(n, z) * cls.jn(n + 1, z)) @classmethod def d3_jn(cls, n: int, z: complex) -> complex: """Third derivative of order :attr:`n` Spherical Bessel function of the first kind :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ return (cls._a0_d3(n, z) * cls.jn(n, z) + cls._a1_d3(n, z) * cls.jn(n + 1, z)) @classmethod def d4_jn(cls, n: int, z: complex) -> complex: """Fourth derivative of order :attr:`n` Spherical Bessel function of the first kind :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ return (cls._a0_d4(n, z) * cls.jn(n, z) + cls._a1_d4(n, z) * cls.jn(n + 1, z)) @staticmethod def hn(n: int, z: complex) -> complex: """Order :attr:`n` Spherical Hankel function of the first kind :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ return sqrt(pi / (2 * z)) * hankel1(n + 0.5, z) @classmethod def d_hn(cls, n: int, z: complex) -> complex: """First derivative of order :attr:`n` Spherical Hankel function of the first kind :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ return -cls.hn(n + 1, z) + (n / z) * cls.hn(n, z) @classmethod def d2_hn(cls, n: int, z: complex) -> complex: """Second derivative of order :attr:`n` Spherical Hankel function of the first kind :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ return (cls._a0_d2(n, z) * cls.hn(n, z) + cls._a1_d2(n, z) * cls.hn(n + 1, z)) @classmethod def d3_hn(cls, n: int, z: complex) -> complex: """Third derivative of order :attr:`n` Spherical Hankel function of the first kind :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ return (cls._a0_d3(n, z) * cls.hn(n, z) + cls._a1_d3(n, z) * cls.hn(n + 1, z)) @classmethod def d4_hn(cls, n: int, z: complex) -> complex: """Fourth derivative of order :attr:`n` Spherical Hankel function of the first kind :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ return (cls._a0_d4(n, z) * cls.hn(n, z) + cls._a1_d4(n, z) * cls.hn(n + 1, z)) @classmethod def d5_hn(cls, n: int, z: complex) -> complex: """Fourth derivative of order :attr:`n` Spherical Hankel function of the first kind :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ return (cls._a0_d5(n, z) * cls.hn(n, z) + cls._a1_d5(n, z) * cls.hn(n + 1, z)) @classmethod def adaptive_derivative_jn(cls, n: int, z: complex, i: int) -> complex: """ Returns the value of the :attr:`i`-th derivative of the spherical Bessel function of order :attr:`n`. The coefficients for the derivative are calculated at runtime. There is no restriction to the magnitude of :attr:`i` other than it must be greater than `1`. :param n: order of spherical Bessel function :type n: int :param z: argument of the function :type z: complex :param i: i-th derivative is computed :type i: int :rtype: complex """ return cls.__adaptive_derivative(cls.jn, i, n, z) @classmethod def adaptive_derivative_hn(cls, n: int, z: complex, i: int) -> complex: """ Returns the value of the :attr:`i`-th derivative of the spherical Hankel function of order :attr:`n`. The coefficients for the derivative are calculated at runtime. There is no restriction to the magnitude of :attr:`i` other than it must be greater than `1`. :param n: order of spherical Hankel function :type n: int :param z: argument of the function :type z: complex :param i: i-th derivative is computed :type i: int :rtype: complex """ return cls.__adaptive_derivative(cls.hn, i, n, z) @classmethod def __adaptive_derivative(cls, func: Callable[[int, complex], complex], i: int, n: int, z: complex) -> complex: out = 0 for j in np.arange(i + 1): out += cls.__adaptive_coefficient( i, j, n, z) * func(n + j, z) return out @classmethod def __adaptive_coefficient(cls, i: int, j: int, n: int, z: complex) -> complex: if j > i or j < 0: out = 0 elif j == i: out = (-1)**i elif j == 0 and i == 1: out = n / z else: out = (n + 2 * j + 1 - i) / z out *= cls.__adaptive_coefficient(i - 1, j, n, z) out -= cls.__adaptive_coefficient(i - 1, j - 1, n, z) return out class LegendreFunctions(object): """ This class gathers all reoccurring Legendre Functions necessary for the gorkov package. """ @staticmethod def cos_poly(theta: float, coefficients: np.array) -> np.ndarray: """ Sum of Legendre polynomial with a cosine function as an argument :param theta: angle :type theta: float :param coefficients: coefficients of the polynomial :type coefficients: np.array :return: Sum of Legendre polynomial with a cosine function as an argument :rtype: np.ndarray """ return np.polynomial.legendre.legval(cos(theta), coefficients, tensor=False) @staticmethod def cos_monomial(degree: int, theta: float, coefficient: complex) -> float: """ Returns a single term of Legendre Cosine series :param degree: degree of the polynomial to be returned :type degree: int :param theta: angle :type theta: float :param coefficient: coefficients of the polynomial :type coefficient: int :return: term of Legendre polynomial with a cosine function as an argument :rtype: float """ return coefficient * lpmv(0, degree, cos(theta)) @staticmethod def first_cos_poly(theta: float, coefficients: np.array) -> np.ndarray: """ Sum of associate Legendre Polynomial of order one with a cosine as input variable :param theta: angle :type theta: float :param coefficients: coefficients of the polynomial :type coefficients: np.array :return: Associate Legendre Polynomial of order one with a cosine as input variable of order n :rtype: np.ndarray """ degrees = len(coefficients) return np.sum([coefficients[n] * lpmv(1, n, cos(theta)) for n in range(degrees)]) @staticmethod def first_cos_monomial(degree: int, theta: float, coefficient: complex) -> float: """ Returns a single term of associate Legendre Cosine series of the first order :param degree: degree of the term :type degree: int :param theta: angle :type theta: float :param coefficient: coefficient of the polynomial :type coefficient: int :return: Sum of Legendre polynomial with a cosine function as an argument :rtype: float """ return coefficient * lpmv(1, degree, cos(theta)) if __name__ == '__main__': pass