To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

functions.py 17 KB
Newer Older
Christoph Goering's avatar
Christoph Goering committed
1
from typing import Callable, Optional, Tuple, Iterable
2 3

import numpy as np
4
from scipy.integrate import quad
Jonas Fankhauser's avatar
Jonas Fankhauser committed
5 6
from scipy.special import hankel1, lpmv, jv, factorial

7

8
from gorkov import log
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's avatar
Christoph Goering committed
20
def full_range(*args: int) -> Iterable:
21
    """Returns the range including the end point.
Christoph Goering's avatar
Christoph Goering committed
22

23
    :param args: end point, start point (optional)
Christoph Goering's avatar
Christoph Goering committed
24
    :type args: int
25
    :return: range including the end point
Christoph Goering's avatar
Christoph Goering committed
26
    :raises: ValueError if more than 2 arguments are passed
Christoph Goering's avatar
Christoph Goering committed
27
    :rtype: Iterable
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's avatar
Christoph Goering committed
37 38
def cartesian_2_polar_coordinates(x: float, z: float) -> Tuple[float,
                                                               float]:
39 40
    """Computes from inputs the polar coordinates and returns them.

41
    :param x: x coordinate
42
    :type x: float
43
    :param z: z coordinate
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]:
55 56
    """Computes from inputs velocities :attr:`v_r`, :attr:`v_theta` and the
    angle :attr:`theta` in the cartesian velocity components
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's avatar
Jonas Fankhauser committed
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.
75
    :param j:
Jonas Fankhauser's avatar
Jonas Fankhauser committed
76
    :type j: int
77
    :param m:
Jonas Fankhauser's avatar
Jonas Fankhauser committed
78
    :type m: int
79
    :param j1:
Jonas Fankhauser's avatar
Jonas Fankhauser committed
80
    :type j1: int
81
    :param m1:
Jonas Fankhauser's avatar
Jonas Fankhauser committed
82
    :type m1: int
83
    :param j2:
Jonas Fankhauser's avatar
Jonas Fankhauser committed
84
    :type j2: int
85
    :param m2:
Jonas Fankhauser's avatar
Jonas Fankhauser committed
86 87
    :type m2: int
    :rtype: int
88
    """
Jonas Fankhauser's avatar
Jonas Fankhauser committed
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
116 117


Christoph Goering's avatar
Christoph Goering committed
118
def integrate(func: Callable[[float], float], lower: float, upper: float,
119
              rel_eps: Optional[float] = 1e-3) -> Tuple[float, float]:
120 121
    """ Integrates :attr:`func` between :attr:`lower` and :attr:`upper`. This
    methods wraps Scipy's QUADPACK implementation `quad()
122 123
    <https://docs.scipy.org/doc/scipy/reference/generated/
    scipy.integrate.quad.html>`__
Christoph Goering's avatar
Christoph Goering committed
124

125 126
    :param func: function to be integrated. If no rel_eps is passed the
        function needs to be vectorized.
Christoph Goering's avatar
Christoph Goering committed
127
    :type func: Callable[[float], float]
128
    :param lower: lower integration limit
Christoph Goering's avatar
Christoph Goering committed
129
    :type lower: float
130
    :param upper: upper integration limit
Christoph Goering's avatar
Christoph Goering committed
131
    :type upper: float
132
    :param rel_eps: estimated relative error of integration
Christoph Goering's avatar
Christoph Goering committed
133
    :type rel_eps: Optional[float]
134 135
    :return: returns the definite integral of func for the boundaries
        lower and upper.
Christoph Goering's avatar
Christoph Goering committed
136
    :rtype: float
137 138
    """

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
150
                  ) -> Tuple[float, float]:
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()
    <https://docs.scipy.org/doc/scipy/reference/generated/
    scipy.integrate.quad.html>`__ 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
201 202


Christoph Goering's avatar
Christoph Goering committed
203
class SpecialFunctions(object):
Jonas Fankhauser's avatar
Jonas Fankhauser committed
204 205
    """ This class gathers all reoccurring spherical hankel and bessel
    functions necessary for the gorkov package. """
Christoph Goering's avatar
Christoph Goering committed
206

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's avatar
Christoph Goering committed
243 244
    @staticmethod
    def jn(n: int, z: complex) -> complex:
245
        """Order :attr:`n` Spherical Hankel function of the first kind
Christoph Goering's avatar
Christoph Goering committed
246 247 248

        :param n: Order
        :type n: int
249
        :param z: argument
Christoph Goering's avatar
Christoph Goering committed
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:
257 258
        """First derivative of order :attr:`n` Spherical Bessel function of the
        first
Christoph Goering's avatar
Christoph Goering committed
259 260 261 262 263 264 265 266
        kind

        :param n: order
        :type n: int
        :param z: argument
        :type z: complex
        :rtype: complex
        """
267
        return -cls.jn(n + 1, z) + (n / z) * cls.jn(n, z)
Christoph Goering's avatar
Christoph Goering committed
268 269 270

    @classmethod
    def d2_jn(cls, n: int, z: complex) -> complex:
Jonas Fankhauser's avatar
Jonas Fankhauser committed
271 272
        """Second derivative of order :attr:`n` Spherical Bessel function of
        the first kind
Christoph Goering's avatar
Christoph Goering committed
273 274 275 276 277 278 279

        :param n: order
        :type n: int
        :param z: argument
        :type z: complex
        :rtype: complex
        """
280 281
        return (cls._a0_d2(n, z) * cls.jn(n, z)
                + cls._a1_d2(n, z) * cls.jn(n + 1, z))
Christoph Goering's avatar
Christoph Goering committed
282 283 284

    @classmethod
    def d3_jn(cls, n: int, z: complex) -> complex:
285 286
        """Third derivative of order :attr:`n` Spherical Bessel function of the
        first
Christoph Goering's avatar
Christoph Goering committed
287 288 289 290 291 292 293 294
        kind

        :param n: order
        :type n: int
        :param z: argument
        :type z: complex
        :rtype: complex
        """
295 296
        return (cls._a0_d3(n, z) * cls.jn(n, z)
                + cls._a1_d3(n, z) * cls.jn(n + 1, z))
Christoph Goering's avatar
Christoph Goering committed
297 298 299

    @classmethod
    def d4_jn(cls, n: int, z: complex) -> complex:
Jonas Fankhauser's avatar
Jonas Fankhauser committed
300 301
        """Fourth derivative of order :attr:`n` Spherical Bessel function of
        the first kind
Christoph Goering's avatar
Christoph Goering committed
302 303 304 305 306 307 308

        :param n: order
        :type n: int
        :param z: argument
        :type z: complex
        :rtype: complex
        """
309 310
        return (cls._a0_d4(n, z) * cls.jn(n, z)
                + cls._a1_d4(n, z) * cls.jn(n + 1, z))
Christoph Goering's avatar
Christoph Goering committed
311 312 313

    @staticmethod
    def hn(n: int, z: complex) -> complex:
314
        """Order :attr:`n` Spherical Hankel function of the first kind
Christoph Goering's avatar
Christoph Goering committed
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:
326 327
        """First derivative of order :attr:`n` Spherical Hankel function of the
        first
Christoph Goering's avatar
Christoph Goering committed
328 329 330 331 332 333 334 335
        kind

        :param n: order
        :type n: int
        :param z: argument
        :type z: complex
        :rtype: complex
        """
336
        return -cls.hn(n + 1, z) + (n / z) * cls.hn(n, z)
Christoph Goering's avatar
Christoph Goering committed
337 338 339

    @classmethod
    def d2_hn(cls, n: int, z: complex) -> complex:
Jonas Fankhauser's avatar
Jonas Fankhauser committed
340 341
        """Second derivative of order :attr:`n` Spherical Hankel function of
        the first kind
Christoph Goering's avatar
Christoph Goering committed
342 343 344 345 346 347 348

        :param n: order
        :type n: int
        :param z: argument
        :type z: complex
        :rtype: complex
        """
349 350
        return (cls._a0_d2(n, z) * cls.hn(n, z)
                + cls._a1_d2(n, z) * cls.hn(n + 1, z))
Christoph Goering's avatar
Christoph Goering committed
351 352 353

    @classmethod
    def d3_hn(cls, n: int, z: complex) -> complex:
Jonas Fankhauser's avatar
Jonas Fankhauser committed
354 355
        """Third derivative of order :attr:`n` Spherical Hankel function of
        the first kind
Christoph Goering's avatar
Christoph Goering committed
356 357 358 359 360 361 362

        :param n: order
        :type n: int
        :param z: argument
        :type z: complex
        :rtype: complex
        """
363 364
        return (cls._a0_d3(n, z) * cls.hn(n, z)
                + cls._a1_d3(n, z) * cls.hn(n + 1, z))
Christoph Goering's avatar
Christoph Goering committed
365 366 367

    @classmethod
    def d4_hn(cls, n: int, z: complex) -> complex:
Jonas Fankhauser's avatar
Jonas Fankhauser committed
368 369
        """Fourth derivative of order :attr:`n` Spherical Hankel function of
        the first kind
Christoph Goering's avatar
Christoph Goering committed
370 371 372 373 374 375 376

        :param n: order
        :type n: int
        :param z: argument
        :type z: complex
        :rtype: complex
        """
377 378
        return (cls._a0_d4(n, z) * cls.hn(n, z)
                + cls._a1_d4(n, z) * cls.hn(n + 1, z))
Christoph Goering's avatar
Christoph Goering committed
379 380 381

    @classmethod
    def d5_hn(cls, n: int, z: complex) -> complex:
Jonas Fankhauser's avatar
Jonas Fankhauser committed
382 383
        """Fourth derivative of order :attr:`n` Spherical Hankel function of
        the first kind
Christoph Goering's avatar
Christoph Goering committed
384 385 386 387 388 389 390

        :param n: order
        :type n: int
        :param z: argument
        :type z: complex
        :rtype: complex
        """
391 392
        return (cls._a0_d5(n, z) * cls.hn(n, z)
                + cls._a1_d5(n, z) * cls.hn(n + 1, z))
Christoph Goering's avatar
Christoph Goering committed
393 394 395

    @classmethod
    def adaptive_derivative_jn(cls, n: int, z: complex, i: int) -> complex:
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's avatar
Christoph Goering committed
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:
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's avatar
Christoph Goering committed
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's avatar
Jonas Fankhauser committed
459
    def cos_poly(theta: float, coefficients: np.array) -> np.ndarray:
Christoph Goering's avatar
Christoph Goering committed
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's avatar
Jonas Fankhauser committed
474
    def cos_monomial(degree: int, theta: float, coefficient: complex) -> float:
Christoph Goering's avatar
Christoph Goering committed
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's avatar
Jonas Fankhauser committed
490
    def first_cos_poly(theta: float, coefficients: np.array) -> np.ndarray:
Christoph Goering's avatar
Christoph Goering committed
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's avatar
Jonas Fankhauser committed
507 508
    def first_cos_monomial(degree: int, theta: float,
                           coefficient: complex) -> float:
Christoph Goering's avatar
Christoph Goering committed
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))


525 526
if __name__ == '__main__':
    pass