Christoph Goering committed Jan 13, 2021 1 ``````from typing import Callable, Optional, Tuple, Iterable `````` Jonas Fankhauser committed Jan 07, 2021 2 3 `````` import numpy as np `````` Jonas Fankhauser committed Feb 15, 2021 4 ``````from scipy.integrate import quad `````` Jonas Fankhauser committed Feb 02, 2021 5 6 ``````from scipy.special import hankel1, lpmv, jv, factorial `````` Jonas Fankhauser committed Jan 07, 2021 7 `````` `````` Jonas Fankhauser committed Jan 06, 2021 8 ``````from gorkov import log `````` Christoph Goering committed Dec 13, 2020 9 10 11 12 13 14 15 16 17 18 19 `````` # 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 `````` Christoph Goering committed Jan 13, 2021 20 ``````def full_range(*args: int) -> Iterable: `````` Christoph Goering committed Dec 13, 2020 21 `````` """Returns the range including the end point. `````` Christoph Goering committed Jan 13, 2021 22 `````` `````` Christoph Goering committed Dec 13, 2020 23 `````` :param args: end point, start point (optional) `````` Christoph Goering committed Jan 13, 2021 24 `````` :type args: int `````` Christoph Goering committed Dec 13, 2020 25 `````` :return: range including the end point `````` Christoph Goering committed Jan 14, 2021 26 `````` :raises: ValueError if more than 2 arguments are passed `````` Christoph Goering committed Jan 14, 2021 27 `````` :rtype: Iterable `````` Christoph Goering committed Dec 13, 2020 28 29 30 31 32 33 34 35 36 `````` """ 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.') `````` Christoph Goering committed Jan 13, 2021 37 38 ``````def cartesian_2_polar_coordinates(x: float, z: float) -> Tuple[float, float]: `````` Christoph Goering committed Dec 13, 2020 39 40 `````` """Computes from inputs the polar coordinates and returns them. `````` Christoph Goering committed Jan 13, 2021 41 `````` :param x: x coordinate `````` Christoph Goering committed Dec 13, 2020 42 `````` :type x: float `````` Christoph Goering committed Jan 13, 2021 43 `````` :param z: z coordinate `````` Christoph Goering committed Dec 13, 2020 44 45 46 47 48 49 50 51 52 53 54 `````` :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]: `````` Christoph Goering committed Jan 13, 2021 55 56 `````` """Computes from inputs velocities :attr:`v_r`, :attr:`v_theta` and the angle :attr:`theta` in the cartesian velocity components `````` Christoph Goering committed Dec 13, 2020 57 58 59 60 61 62 63 64 65 66 67 68 69 70 `````` :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 `````` Jonas Fankhauser committed Feb 02, 2021 71 72 73 74 ``````def Clebsch_Gordan_coefficient(j1: int, m1: int, j2: int, m2: int, j: int, m: int) -> float: """ Clebsch-Gordan coefficient for integer valued arguments. `````` Christoph Goering committed Dec 13, 2020 75 `````` :param j: `````` Jonas Fankhauser committed Feb 02, 2021 76 `````` :type j: int `````` Christoph Goering committed Dec 13, 2020 77 `````` :param m: `````` Jonas Fankhauser committed Feb 02, 2021 78 `````` :type m: int `````` Christoph Goering committed Dec 13, 2020 79 `````` :param j1: `````` Jonas Fankhauser committed Feb 02, 2021 80 `````` :type j1: int `````` Christoph Goering committed Dec 13, 2020 81 `````` :param m1: `````` Jonas Fankhauser committed Feb 02, 2021 82 `````` :type m1: int `````` Christoph Goering committed Dec 13, 2020 83 `````` :param j2: `````` Jonas Fankhauser committed Feb 02, 2021 84 `````` :type j2: int `````` Christoph Goering committed Dec 13, 2020 85 `````` :param m2: `````` Jonas Fankhauser committed Feb 02, 2021 86 87 `````` :type m2: int :rtype: int `````` Christoph Goering committed Dec 13, 2020 88 `````` """ `````` Jonas Fankhauser committed Feb 02, 2021 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 `````` 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 `````` Jonas Fankhauser committed Jan 06, 2021 116 117 `````` `````` Christoph Goering committed Jan 13, 2021 118 ``````def integrate(func: Callable[[float], float], lower: float, upper: float, `````` Jonas Fankhauser committed Feb 17, 2021 119 `````` rel_eps: Optional[float] = 1e-3) -> Tuple[float, float]: `````` Jonas Fankhauser committed Feb 15, 2021 120 121 `````` """ Integrates :attr:`func` between :attr:`lower` and :attr:`upper`. This methods wraps Scipy's QUADPACK implementation `quad() `````` Christoph Goering committed Jan 13, 2021 122 123 `````` `__ `````` Christoph Goering committed Jan 13, 2021 124 `````` `````` Jonas Fankhauser committed Jan 06, 2021 125 126 `````` :param func: function to be integrated. If no rel_eps is passed the function needs to be vectorized. `````` Christoph Goering committed Jan 13, 2021 127 `````` :type func: Callable[[float], float] `````` Jonas Fankhauser committed Jan 06, 2021 128 `````` :param lower: lower integration limit `````` Christoph Goering committed Jan 13, 2021 129 `````` :type lower: float `````` Jonas Fankhauser committed Jan 06, 2021 130 `````` :param upper: upper integration limit `````` Christoph Goering committed Jan 13, 2021 131 `````` :type upper: float `````` Jonas Fankhauser committed Jan 06, 2021 132 `````` :param rel_eps: estimated relative error of integration `````` Christoph Goering committed Jan 13, 2021 133 `````` :type rel_eps: Optional[float] `````` Jonas Fankhauser committed Jan 06, 2021 134 135 `````` :return: returns the definite integral of func for the boundaries lower and upper. `````` Christoph Goering committed Jan 13, 2021 136 `````` :rtype: float `````` Jonas Fankhauser committed Jan 06, 2021 137 138 `````` """ `````` Jonas Fankhauser committed Feb 15, 2021 139 140 141 142 143 144 145 146 147 148 149 `````` # 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 `````` Jonas Fankhauser committed Feb 17, 2021 150 `````` ) -> Tuple[float, float]: `````` Jonas Fankhauser committed Feb 15, 2021 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 `````` """ 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 `````` Christoph Goering committed Dec 13, 2020 201 202 `````` `````` Christoph Goering committed Jan 13, 2021 203 ``````class SpecialFunctions(object): `````` Jonas Fankhauser committed Jan 21, 2021 204 205 `````` """ This class gathers all reoccurring spherical hankel and bessel functions necessary for the gorkov package. """ `````` Christoph Goering committed Jan 13, 2021 206 `````` `````` Jonas Fankhauser committed Jan 27, 2021 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 `````` @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) `````` Christoph Goering committed Jan 13, 2021 243 244 `````` @staticmethod def jn(n: int, z: complex) -> complex: `````` Christoph Goering committed Jan 13, 2021 245 `````` """Order :attr:`n` Spherical Hankel function of the first kind `````` Christoph Goering committed Jan 13, 2021 246 247 248 `````` :param n: Order :type n: int `````` Christoph Goering committed Jan 13, 2021 249 `````` :param z: argument `````` Christoph Goering committed Jan 13, 2021 250 251 252 253 254 255 256 `````` :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: `````` Christoph Goering committed Jan 13, 2021 257 258 `````` """First derivative of order :attr:`n` Spherical Bessel function of the first `````` Christoph Goering committed Jan 13, 2021 259 260 261 262 263 264 265 266 `````` kind :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ `````` Jonas Fankhauser committed Jan 27, 2021 267 `````` return -cls.jn(n + 1, z) + (n / z) * cls.jn(n, z) `````` Christoph Goering committed Jan 13, 2021 268 269 270 `````` @classmethod def d2_jn(cls, n: int, z: complex) -> complex: `````` Jonas Fankhauser committed Jan 21, 2021 271 272 `````` """Second derivative of order :attr:`n` Spherical Bessel function of the first kind `````` Christoph Goering committed Jan 13, 2021 273 274 275 276 277 278 279 `````` :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ `````` Jonas Fankhauser committed Jan 27, 2021 280 281 `````` return (cls._a0_d2(n, z) * cls.jn(n, z) + cls._a1_d2(n, z) * cls.jn(n + 1, z)) `````` Christoph Goering committed Jan 13, 2021 282 283 284 `````` @classmethod def d3_jn(cls, n: int, z: complex) -> complex: `````` Christoph Goering committed Jan 13, 2021 285 286 `````` """Third derivative of order :attr:`n` Spherical Bessel function of the first `````` Christoph Goering committed Jan 13, 2021 287 288 289 290 291 292 293 294 `````` kind :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ `````` Jonas Fankhauser committed Jan 27, 2021 295 296 `````` return (cls._a0_d3(n, z) * cls.jn(n, z) + cls._a1_d3(n, z) * cls.jn(n + 1, z)) `````` Christoph Goering committed Jan 13, 2021 297 298 299 `````` @classmethod def d4_jn(cls, n: int, z: complex) -> complex: `````` Jonas Fankhauser committed Jan 21, 2021 300 301 `````` """Fourth derivative of order :attr:`n` Spherical Bessel function of the first kind `````` Christoph Goering committed Jan 13, 2021 302 303 304 305 306 307 308 `````` :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ `````` Jonas Fankhauser committed Jan 27, 2021 309 310 `````` return (cls._a0_d4(n, z) * cls.jn(n, z) + cls._a1_d4(n, z) * cls.jn(n + 1, z)) `````` Christoph Goering committed Jan 13, 2021 311 312 313 `````` @staticmethod def hn(n: int, z: complex) -> complex: `````` Christoph Goering committed Jan 13, 2021 314 `````` """Order :attr:`n` Spherical Hankel function of the first kind `````` Christoph Goering committed Jan 13, 2021 315 316 317 318 319 320 321 322 323 324 325 `````` :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: `````` Christoph Goering committed Jan 13, 2021 326 327 `````` """First derivative of order :attr:`n` Spherical Hankel function of the first `````` Christoph Goering committed Jan 13, 2021 328 329 330 331 332 333 334 335 `````` kind :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ `````` Jonas Fankhauser committed Jan 27, 2021 336 `````` return -cls.hn(n + 1, z) + (n / z) * cls.hn(n, z) `````` Christoph Goering committed Jan 13, 2021 337 338 339 `````` @classmethod def d2_hn(cls, n: int, z: complex) -> complex: `````` Jonas Fankhauser committed Jan 21, 2021 340 341 `````` """Second derivative of order :attr:`n` Spherical Hankel function of the first kind `````` Christoph Goering committed Jan 13, 2021 342 343 344 345 346 347 348 `````` :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ `````` Jonas Fankhauser committed Jan 27, 2021 349 350 `````` return (cls._a0_d2(n, z) * cls.hn(n, z) + cls._a1_d2(n, z) * cls.hn(n + 1, z)) `````` Christoph Goering committed Jan 13, 2021 351 352 353 `````` @classmethod def d3_hn(cls, n: int, z: complex) -> complex: `````` Jonas Fankhauser committed Jan 21, 2021 354 355 `````` """Third derivative of order :attr:`n` Spherical Hankel function of the first kind `````` Christoph Goering committed Jan 13, 2021 356 357 358 359 360 361 362 `````` :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ `````` Jonas Fankhauser committed Jan 27, 2021 363 364 `````` return (cls._a0_d3(n, z) * cls.hn(n, z) + cls._a1_d3(n, z) * cls.hn(n + 1, z)) `````` Christoph Goering committed Jan 13, 2021 365 366 367 `````` @classmethod def d4_hn(cls, n: int, z: complex) -> complex: `````` Jonas Fankhauser committed Jan 21, 2021 368 369 `````` """Fourth derivative of order :attr:`n` Spherical Hankel function of the first kind `````` Christoph Goering committed Jan 13, 2021 370 371 372 373 374 375 376 `````` :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ `````` Jonas Fankhauser committed Jan 27, 2021 377 378 `````` return (cls._a0_d4(n, z) * cls.hn(n, z) + cls._a1_d4(n, z) * cls.hn(n + 1, z)) `````` Christoph Goering committed Jan 13, 2021 379 380 381 `````` @classmethod def d5_hn(cls, n: int, z: complex) -> complex: `````` Jonas Fankhauser committed Jan 21, 2021 382 383 `````` """Fourth derivative of order :attr:`n` Spherical Hankel function of the first kind `````` Christoph Goering committed Jan 13, 2021 384 385 386 387 388 389 390 `````` :param n: order :type n: int :param z: argument :type z: complex :rtype: complex """ `````` Jonas Fankhauser committed Jan 27, 2021 391 392 `````` return (cls._a0_d5(n, z) * cls.hn(n, z) + cls._a1_d5(n, z) * cls.hn(n + 1, z)) `````` Christoph Goering committed Jan 13, 2021 393 394 395 `````` @classmethod def adaptive_derivative_jn(cls, n: int, z: complex, i: int) -> complex: `````` Christoph Goering committed Jan 13, 2021 396 397 398 399 `````` """ 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`. `````` Christoph Goering committed Jan 13, 2021 400 401 402 403 404 405 406 407 408 409 410 411 412 `````` :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: `````` Christoph Goering committed Jan 13, 2021 413 414 415 416 `````` """ 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`. `````` Christoph Goering committed Jan 13, 2021 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 `````` :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 `````` Jonas Fankhauser committed Jan 21, 2021 459 `````` def cos_poly(theta: float, coefficients: np.array) -> np.ndarray: `````` Christoph Goering committed Jan 13, 2021 460 461 462 463 464 465 466 467 468 469 470 471 472 473 `````` """ 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 `````` Jonas Fankhauser committed Jan 21, 2021 474 `````` def cos_monomial(degree: int, theta: float, coefficient: complex) -> float: `````` Christoph Goering committed Jan 13, 2021 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 `````` """ 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 `````` Jonas Fankhauser committed Jan 21, 2021 490 `````` def first_cos_poly(theta: float, coefficients: np.array) -> np.ndarray: `````` Christoph Goering committed Jan 13, 2021 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 `````` """ 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 `````` Jonas Fankhauser committed Jan 21, 2021 507 508 `````` def first_cos_monomial(degree: int, theta: float, coefficient: complex) -> float: `````` Christoph Goering committed Jan 13, 2021 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 `````` """ 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)) `````` Christoph Goering committed Dec 13, 2020 525 526 ``````if __name__ == '__main__': pass``````