Skip to content

Turbulence

turbulence

KolmogorovPowerSpectrum

Bases: Module

Kolmogorov power spectrum

Source code in src/turbulence.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
class KolmogorovPowerSpectrum(hk.Module):
    """
    Kolmogorov power spectrum
    """

    def __init__(self):
        super(KolmogorovPowerSpectrum, self).__init__()

    def __call__(self, k):
        r"""Kolmogorov power spectrum

        $$\mathcal{P}_{3D}(k)= \sigma^2 \frac{e^{-\left(k/k_{\text{inj}}\right)^2} e^{-\left(k_{\text{dis}}/k\right)^2} k^{-\alpha} }{\int 4\pi k^2  \, e^{-\left(k/k_{\text{inj}}\right)^2} e^{-\left(k_{\text{dis}}/k\right)^2} k^{-\alpha} \text{d} k}$$

        """

        sigma = hk.get_parameter("sigma", [], init=Constant(250)) #km/s
        log_inj = hk.get_parameter("log_inj", [], init=Constant(jnp.log10(300))) #kpc
        log_dis = jnp.log10(10) #kpc
        alpha = hk.get_parameter("alpha", [], init=Constant(11/3))

        k_inj = 10 ** (-log_inj)
        k_dis = 10 ** (-log_dis)

        #sigma = 10**log_sigma

        k_int = jnp.geomspace(k_inj/20, k_dis*20, 1000)
        norm = integrate.trapezoid(4*jnp.pi*k_int**3*jnp.exp(-(k_inj / k_int) ** 2) * jnp.exp(-(k_int/ k_dis) ** 2) * (k_int) ** (-alpha), x=jnp.log(k_int))
        res = jnp.where(k > 0, jnp.exp(-(k_inj / k) ** 2) * jnp.exp(-(k/ k_dis) ** 2) * (k) ** (-alpha), 0.)

        return sigma**2 * res / norm

__call__(k)

Kolmogorov power spectrum

\[\mathcal{P}_{3D}(k)= \sigma^2 \frac{e^{-\left(k/k_{\text{inj}}\right)^2} e^{-\left(k_{\text{dis}}/k\right)^2} k^{-\alpha} }{\int 4\pi k^2 \, e^{-\left(k/k_{\text{inj}}\right)^2} e^{-\left(k_{\text{dis}}/k\right)^2} k^{-\alpha} \text{d} k}\]
Source code in src/turbulence.py 
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
def __call__(self, k):
    r"""Kolmogorov power spectrum

    $$\mathcal{P}_{3D}(k)= \sigma^2 \frac{e^{-\left(k/k_{\text{inj}}\right)^2} e^{-\left(k_{\text{dis}}/k\right)^2} k^{-\alpha} }{\int 4\pi k^2  \, e^{-\left(k/k_{\text{inj}}\right)^2} e^{-\left(k_{\text{dis}}/k\right)^2} k^{-\alpha} \text{d} k}$$

    """

    sigma = hk.get_parameter("sigma", [], init=Constant(250)) #km/s
    log_inj = hk.get_parameter("log_inj", [], init=Constant(jnp.log10(300))) #kpc
    log_dis = jnp.log10(10) #kpc
    alpha = hk.get_parameter("alpha", [], init=Constant(11/3))

    k_inj = 10 ** (-log_inj)
    k_dis = 10 ** (-log_dis)

    #sigma = 10**log_sigma

    k_int = jnp.geomspace(k_inj/20, k_dis*20, 1000)
    norm = integrate.trapezoid(4*jnp.pi*k_int**3*jnp.exp(-(k_inj / k_int) ** 2) * jnp.exp(-(k_int/ k_dis) ** 2) * (k_int) ** (-alpha), x=jnp.log(k_int))
    res = jnp.where(k > 0, jnp.exp(-(k_inj / k) ** 2) * jnp.exp(-(k/ k_dis) ** 2) * (k) ** (-alpha), 0.)

    return sigma**2 * res / norm