Skip to content

physics

xsb_fluc.physics.density

Python module containing models for the density profile of the ICM.

BetaModel

Bases: Module

Density model which use a beta-model formula

\[n_{e}^2(r) = n_{e,0}^2 \left( 1 + \left( \frac{r}{r_{c}} \right)^{2} \right)^{-3\beta}\]
Source code in src/xsb_fluc/physics/density.py 
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
class BetaModel(hk.Module):
    r"""Density model which use a beta-model formula

    $$n_{e}^2(r) = n_{e,0}^2 \left( 1 + \left( \frac{r}{r_{c}} \right)^{2} \right)^{-3\beta}$$
    """
    def __init__(self):
        super(BetaModel, self).__init__()

    def __call__(self, r):
        """Compute the density function for a given radius.

        Parameters:
            r (jnp.array): Radius to compute the density function in R500 units

        Returns:
            (jnp.array): Density function evaluated at the given radius in cm$^{-6}$
        """
        log_ne2 = hk.get_parameter("log_ne2", [], init=Constant(-3.))
        log_r_c = hk.get_parameter("log_r_c", [], init=Constant(-1))
        beta = hk.get_parameter("beta", [], init=Constant(2/3))

        ne2 = 10 ** log_ne2
        r_c = 10 ** log_r_c

        term1 = (1. + (r / r_c) ** 2) ** (-3 * beta)

        return ne2 * term1
__call__(r)

Compute the density function for a given radius.

Parameters:

Name Type Description Default
r array

Radius to compute the density function in R500 units

 required

Returns:

Type Description
array

Density function evaluated at the given radius in cm\(^{-6}\)
Source code in src/xsb_fluc/physics/density.py 
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
def __call__(self, r):
    """Compute the density function for a given radius.

    Parameters:
        r (jnp.array): Radius to compute the density function in R500 units

    Returns:
        (jnp.array): Density function evaluated at the given radius in cm$^{-6}$
    """
    log_ne2 = hk.get_parameter("log_ne2", [], init=Constant(-3.))
    log_r_c = hk.get_parameter("log_r_c", [], init=Constant(-1))
    beta = hk.get_parameter("beta", [], init=Constant(2/3))

    ne2 = 10 ** log_ne2
    r_c = 10 ** log_r_c

    term1 = (1. + (r / r_c) ** 2) ** (-3 * beta)

    return ne2 * term1

CleanVikhlininModel

Bases: Module

Density model which use a modified Vikhlinin functional form, with alpha fixed to 0 and gamma to 3

\[n_{e}^2(r) = n_{e,0}^2 \left( 1 + \left( \frac{r}{r_{c}} \right)^{2} \right)^{-3\beta} \left(1+\left( \frac{r}{r_{c}} \right)^{\gamma} \right)^{-\frac{\epsilon}{\gamma}}\]
Source code in src/xsb_fluc/physics/density.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
38
39
40
41
42
class CleanVikhlininModel(hk.Module):
    r"""
    Density model which use a modified Vikhlinin functional form, with alpha fixed to 0 and gamma to 3

    $$n_{e}^2(r) = n_{e,0}^2 \left( 1 + \left( \frac{r}{r_{c}} \right)^{2} \right)^{-3\beta} \left(1+\left( \frac{r}{r_{c}} \right)^{\gamma} \right)^{-\frac{\epsilon}{\gamma}}$$
    """

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

    def __call__(self, r: jnp.array) -> jnp.array:
        """Compute the density function for a given radius.

        Parameters:
            r (jnp.array): Radius to compute the density function in R500 units

        Returns:
            (jnp.array): Density function evaluated at the given radius in cm$^{-6}$
        """
        log_ne2 = hk.get_parameter("log_ne2", [], init=Constant(-3.))
        log_r_c = hk.get_parameter("log_r_c", [], init=Constant(-1.))
        log_r_s = hk.get_parameter("log_r_s", [], init=Constant(-0.1))
        beta = hk.get_parameter("beta", [], init=Constant(0.6))
        epsilon = hk.get_parameter("epsilon", [], init=Constant(3.))

        gamma = 3.

        ne2 = 10 ** log_ne2
        r_c = 10 ** log_r_c
        r_s = 10 ** log_r_s

        term1 = (1. + (r / r_c) ** 2) ** (-3 * beta)
        term2 = (1. + (r / r_s) ** gamma) ** (-epsilon / gamma)

        return ne2 * term1 * term2
__call__(r)

Compute the density function for a given radius.

Parameters:

Name Type Description Default
r array

Radius to compute the density function in R500 units

 required

Returns:

Type Description
array

Density function evaluated at the given radius in cm\(^{-6}\)
Source code in src/xsb_fluc/physics/density.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
def __call__(self, r: jnp.array) -> jnp.array:
    """Compute the density function for a given radius.

    Parameters:
        r (jnp.array): Radius to compute the density function in R500 units

    Returns:
        (jnp.array): Density function evaluated at the given radius in cm$^{-6}$
    """
    log_ne2 = hk.get_parameter("log_ne2", [], init=Constant(-3.))
    log_r_c = hk.get_parameter("log_r_c", [], init=Constant(-1.))
    log_r_s = hk.get_parameter("log_r_s", [], init=Constant(-0.1))
    beta = hk.get_parameter("beta", [], init=Constant(0.6))
    epsilon = hk.get_parameter("epsilon", [], init=Constant(3.))

    gamma = 3.

    ne2 = 10 ** log_ne2
    r_c = 10 ** log_r_c
    r_s = 10 ** log_r_s

    term1 = (1. + (r / r_c) ** 2) ** (-3 * beta)
    term2 = (1. + (r / r_s) ** gamma) ** (-epsilon / gamma)

    return ne2 * term1 * term2

xsb_fluc.physics.temperature

Python module containing models for the temperature profile of the ICM.

GhirardiniModel

Bases: Module

Universal temperature profile as defined in Ghirardini 2018+ in the X-COP cluster sample

Source code in src/xsb_fluc/physics/temperature.py 
 6
 7
 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 GhirardiniModel(hk.Module):
    """
    Universal temperature profile as defined in Ghirardini 2018+ in the X-COP cluster sample
    """

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

    def __call__(self, r):
        """
        Compute the temperature function for a given radius.

        Parameters:
            r (jnp.array): Radius to compute the temperature in $R_{500}$ units

        Returns:
            (jnp.array): Temperature function evaluated at the given radius in keV
        """

        T0 = 1.21
        rcool = jnp.exp(-2.78)
        rt = 0.34
        TmT0 = 0.5
        acool = 1.03
        c2 = 0.27

        T500 = 7.  # keV

        term1 = (TmT0 + (r / rcool) ** acool)
        term2 = (1 + (r / rcool) ** acool) * (1 + (r / rt) ** 2) ** c2

        return T500 * T0 * term1 / term2

__call__(r)

Compute the temperature function for a given radius.

Parameters:

Name Type Description Default
r array

Radius to compute the temperature in \(R_{500}\) units

 required

Returns:

Type Description
array

Temperature function evaluated at the given radius in keV
Source code in src/xsb_fluc/physics/temperature.py 
14
15
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, r):
    """
    Compute the temperature function for a given radius.

    Parameters:
        r (jnp.array): Radius to compute the temperature in $R_{500}$ units

    Returns:
        (jnp.array): Temperature function evaluated at the given radius in keV
    """

    T0 = 1.21
    rcool = jnp.exp(-2.78)
    rt = 0.34
    TmT0 = 0.5
    acool = 1.03
    c2 = 0.27

    T500 = 7.  # keV

    term1 = (TmT0 + (r / rcool) ** acool)
    term2 = (1 + (r / rcool) ** acool) * (1 + (r / rt) ** 2) ** c2

    return T500 * T0 * term1 / term2

xsb_fluc.physics.cooling

Python module containing structures to compute the cooling function approximation for the ICM.

CoolingFunctionGrid

Source code in src/xsb_fluc/physics/cooling.py 
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
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
class CoolingFunctionGrid:

    def __init__(self, redshift=0.01, n_points=20, kT_span=(1., 10.), nH_span=(9e19, 2e21)):
        """
        Awfully coded class to compute the cooling function grid for a given redshift.
        The grid is saved in a .npy file in the results folder. XSPEC is not needed if the
        grid is already computed. Else it will be computed using XSPEC.

        """

        self.redshift = round(redshift, 5)
        self.n_points = n_points
        self.kT_span = kT_span
        self.nH_span = nH_span

        pars = (self.redshift, n_points, *kT_span, *nH_span)

        self.kT, self.nH = np.meshgrid(np.linspace(*kT_span, n_points),
                                       np.geomspace(*nH_span, n_points),
                                       indexing='ij')

        cooling_path = os.path.join(config.RESULTS_PATH, 'cooling_database', f'{pars}.npy')

        if not os.path.exists(cooling_path):

            print('Cooling grid not found! Computing with XSPEC')

            from .countrate import PHabsAPEC

            phabs = PHabsAPEC()

            abundance = 0.3
            norm = 1.
            countrate = phabs.countrate(self.nH / 1e22, self.kT, abundance, float(self.redshift), norm)

            # XSPEC norm factor for an apec model
            dc = FlatLambdaCDM(70, 0.3).comoving_distance(redshift).to(u.Mpc)
            factor = 1e-14 / (4 * np.pi * dc ** 2) / 1.17

            # Units of the countrate over units of the norm of APEC model
            cooling = factor * countrate * (u.count * u.s ** (-1) * u.cm ** 5)

            # Conversion to the best unit for Lambda * ne^2
            self.cooling = (cooling.to(u.count * u.cm ** 6 / u.kpc ** 3 / u.s)).value
            np.save(cooling_path, self.cooling)

        else:

            self.cooling = np.load(cooling_path)

__init__(redshift=0.01, n_points=20, kT_span=(1.0, 10.0), nH_span=(9e+19, 2e+21))

Awfully coded class to compute the cooling function grid for a given redshift. The grid is saved in a .npy file in the results folder. XSPEC is not needed if the grid is already computed. Else it will be computed using XSPEC.

Source code in src/xsb_fluc/physics/cooling.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def __init__(self, redshift=0.01, n_points=20, kT_span=(1., 10.), nH_span=(9e19, 2e21)):
    """
    Awfully coded class to compute the cooling function grid for a given redshift.
    The grid is saved in a .npy file in the results folder. XSPEC is not needed if the
    grid is already computed. Else it will be computed using XSPEC.

    """

    self.redshift = round(redshift, 5)
    self.n_points = n_points
    self.kT_span = kT_span
    self.nH_span = nH_span

    pars = (self.redshift, n_points, *kT_span, *nH_span)

    self.kT, self.nH = np.meshgrid(np.linspace(*kT_span, n_points),
                                   np.geomspace(*nH_span, n_points),
                                   indexing='ij')

    cooling_path = os.path.join(config.RESULTS_PATH, 'cooling_database', f'{pars}.npy')

    if not os.path.exists(cooling_path):

        print('Cooling grid not found! Computing with XSPEC')

        from .countrate import PHabsAPEC

        phabs = PHabsAPEC()

        abundance = 0.3
        norm = 1.
        countrate = phabs.countrate(self.nH / 1e22, self.kT, abundance, float(self.redshift), norm)

        # XSPEC norm factor for an apec model
        dc = FlatLambdaCDM(70, 0.3).comoving_distance(redshift).to(u.Mpc)
        factor = 1e-14 / (4 * np.pi * dc ** 2) / 1.17

        # Units of the countrate over units of the norm of APEC model
        cooling = factor * countrate * (u.count * u.s ** (-1) * u.cm ** 5)

        # Conversion to the best unit for Lambda * ne^2
        self.cooling = (cooling.to(u.count * u.cm ** 6 / u.kpc ** 3 / u.s)).value
        np.save(cooling_path, self.cooling)

    else:

        self.cooling = np.load(cooling_path)

CoolingFunctionModel

Bases: Module

Cooling function model using a grid of precomputed cooling function. It is fitted to the grid using a least square method on the fly. The model is given by the following formula:

\[\Psi(N_H, T) \simeq \Lambda_0 e^{- N_H \sigma} \left( \frac{T}{T_{\text{break}}}\right)^{-\alpha_{1}}\left(\frac{1}{2} + \frac{1}{2}\left(\frac{T}{T_{\text{break}}}\right)^{1/\Delta}\right)^{(\alpha_1 - \alpha_2)\Delta}\]

Note

\(\Psi(N_H, T)\) here is different from \(\Lambda(T)\) as it includes the instrumental convolution

Source code in src/xsb_fluc/physics/cooling.py 
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
class CoolingFunctionModel(hk.Module):
    r"""
    Cooling function model using a grid of precomputed cooling function. It is fitted to the grid using
    a least square method on the fly. The model is given by the following formula:

    $$\Psi(N_H, T) \simeq \Lambda_0 e^{- N_H \sigma} \left( \frac{T}{T_{\text{break}}}\right)^{-\alpha_{1}}\left(\frac{1}{2} + \frac{1}{2}\left(\frac{T}{T_{\text{break}}}\right)^{1/\Delta}\right)^{(\alpha_1 - \alpha_2)\Delta}$$

    !!! note
        $\Psi(N_H, T)$ here is different from $\Lambda(T)$ as it includes the instrumental convolution
    """

    def __init__(self, coolingFunctionGrid):
        super(CoolingFunctionModel, self).__init__()
        self.grid = coolingFunctionGrid
        # Norm, sigma, kTbreak, alpha1, alpha2, delta
        X = jnp.array([0.01, 3., 1.65, 1.7, 0.14, 0.13])

        self.res = minimize(self.fitness, X, method="BFGS", tol=1e-15)
        self.pars = self.res.x

    def kT_dependency(self, kT, kT_break, alpha1, alpha2, delta):
        t = (kT / kT_break)

        return t ** (-alpha1) * (1 / 2 * (1 + t ** (1 / delta))) ** ((alpha1 - alpha2) * delta)

    def nH_dependency(self, nH, sigma):
        return jnp.exp(-nH / 1e22 * sigma)

    def model(self, nH, kT, pars):
        return pars[0] * self.nH_dependency(nH, pars[1]) * self.kT_dependency(kT, *pars[2:])

    def __call__(self, nH, kT):
        return self.model(nH, kT, self.pars)

    def fitness(self, pars):
        lsq = (self.grid.cooling - self.model(self.grid.nH, self.grid.kT, pars)) ** 2

        return jnp.sum(lsq / self.grid.cooling ** 2) ** 1 / 2

    @property
    def residual(self):
        lsq = (self.grid.cooling - self.model(self.grid.nH, self.grid.kT, self.pars)) ** 2

        return lsq ** (1 / 2) / self.grid.cooling

    def plot_residual(self):
        plt.figure()
        plt.contourf(self.grid.kT, self.grid.nH, self.residual)
        plt.colorbar()

xsb_fluc.physics.emissivity

XrayEmissivity

Bases: Module

3D Xray emissivity build with temperature, cooling function, density model. It depends on the redshift of the cluster, since the cooling function is precomputed using XSPEC. The default models are the ones used in the papers i.e. Vikhlinin for density, Ghirardini for temperature and the interpolated cooling function.

Source code in src/xsb_fluc/physics/emissivity.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
class XrayEmissivity(hk.Module):
    """
    3D Xray emissivity build with temperature, cooling function, density model.
    It depends on the redshift of the cluster, since the cooling function is precomputed using XSPEC.
    The default models are the ones used in the papers i.e. Vikhlinin for density, Ghirardini for temperature
    and the interpolated cooling function.
    """
    def __init__(self, redshift):
        super(XrayEmissivity, self).__init__()
        self.squared_density = CleanVikhlininModel()
        self.temperature = GhirardiniModel()
        self.cooling_function = CoolingFunctionModel(CoolingFunctionGrid(n_points=10, redshift=redshift))

    @classmethod
    def from_data(cls, data: "Cluster"):
        """
        Create an emissivity model from a `Cluster` object
        """
        return cls(data.z)

    def __call__(self, r, nH):
        """
        Compute the emissivity at a given radius, including $N_H$ absorption.

        Parameters:
            r (jnp.array): radius in units of $R_{500}$
            nH (jnp.array): Hydrogen density in cm$^{-2}$ in the line of sight corresponding to the radius r
        """
        return self.cooling_function(nH, self.temperature(r)) * self.squared_density(r)

__call__(r, nH)

Compute the emissivity at a given radius, including \(N_H\) absorption.

Parameters:

Name Type Description Default
r array

radius in units of \(R_{500}\)

 required
nH array

Hydrogen density in cm\(^{-2}\) in the line of sight corresponding to the radius r

 required
Source code in src/xsb_fluc/physics/emissivity.py 
28
29
30
31
32
33
34
35
36
def __call__(self, r, nH):
    """
    Compute the emissivity at a given radius, including $N_H$ absorption.

    Parameters:
        r (jnp.array): radius in units of $R_{500}$
        nH (jnp.array): Hydrogen density in cm$^{-2}$ in the line of sight corresponding to the radius r
    """
    return self.cooling_function(nH, self.temperature(r)) * self.squared_density(r)

from_data(data) classmethod

Create an emissivity model from a Cluster object

Source code in src/xsb_fluc/physics/emissivity.py 
21
22
23
24
25
26
@classmethod
def from_data(cls, data: "Cluster"):
    """
    Create an emissivity model from a `Cluster` object
    """
    return cls(data.z)

xsb_fluc.physics.projection

AbelTransform

Bases: Module

Projection tool for 3D models. Relies on double exponential quadrature.

\[\text{AbelTransform}\left\{ f\right\}(x, y) = \int \text{d}z \, f(x, y, z)\]

It uses double exponential quadrature to perform the integration.

References :

Source code in src/xsb_fluc/physics/projection.py 
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
38
39
40
41
42
43
44
45
46
47
48
49
class AbelTransform(hk.Module):
    r"""
    Projection tool for 3D models. Relies on double exponential quadrature.

    $$\text{AbelTransform}\left\{ f\right\}(x, y) =  \int \text{d}z \, f(x, y, z)$$

    It uses double exponential quadrature to perform the integration.

    References :

    * [Takahasi and Mori (1974)](https://ems.press/journals/prims/articles/2686)
    * [Mori and Sugihara (2001)](https://doi.org/10.1016/S0377-0427(00)00501-X)
    * [Tanh-sinh quadrature](https://en.wikipedia.org/wiki/Tanh-sinh_quadrature) from wikipedia

    """

    def __init__(self, model:"hk.Module", n_points=71):
        """
        Take a model and project it along the line of sight.

        Parameters:
            model (hk.Module): model to project
            n_points (int): number of points to use for the quadrature
        """
        super(AbelTransform, self).__init__()
        self.model = model
        self.n = n_points
        self.t = jnp.linspace(-3, 3, self.n)

    def __call__(self, r, *args):

        r = jnp.asarray(r)[..., None]
        t = self.t[None, :]
        x = phi(t)
        dx = dphi(t)

        return 2 * jnp.trapz(self.model(jnp.sqrt(r ** 2 + x ** 2), *args) * dx, x=t, axis=-1)

__init__(model, n_points=71)

Take a model and project it along the line of sight.

Parameters:

Name Type Description Default
model Module

model to project

 required
n_points int

number of points to use for the quadrature

 71
Source code in src/xsb_fluc/physics/projection.py 
29
30
31
32
33
34
35
36
37
38
39
40
def __init__(self, model:"hk.Module", n_points=71):
    """
    Take a model and project it along the line of sight.

    Parameters:
        model (hk.Module): model to project
        n_points (int): number of points to use for the quadrature
    """
    super(AbelTransform, self).__init__()
    self.model = model
    self.n = n_points
    self.t = jnp.linspace(-3, 3, self.n)

xsb_fluc.physics.surface_brightness

XraySurfaceBrightness

Bases: Module

Xray surface brightness model from 3D emissivity

Source code in src/xsb_fluc/physics/surface_brightness.py 
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
class XraySurfaceBrightness(hk.Module):
    """
    Xray surface brightness model from 3D emissivity
    """

    def __init__(self, redshift, r_500):
        super(XraySurfaceBrightness, self).__init__()
        self.emissivity = XrayEmissivity(redshift)
        self.surface_brightness = AbelTransform(self.emissivity)
        self.r_500 = r_500.to(u.kpc).value

    @classmethod
    def from_data(cls, data: "Cluster"):
        """
        Create a surface brightness model from a `Cluster` object
        """

        return cls(data.z, data.r_500)

    def __call__(self, r, nh):
        log_bkg = hk.get_parameter("log_bkg", [], init=Constant(-5.))

        # As we integrate toward the l.o.s in r_500 units,
        # the surface brightness must be rescaled to be in good units
        return self.surface_brightness(r, nh[..., None]) * self.r_500 + 10 ** log_bkg

from_data(data) classmethod

Create a surface brightness model from a Cluster object

Source code in src/xsb_fluc/physics/surface_brightness.py 
22
23
24
25
26
27
28
@classmethod
def from_data(cls, data: "Cluster"):
    """
    Create a surface brightness model from a `Cluster` object
    """

    return cls(data.z, data.r_500)

XraySurfaceBrightnessBetaModel

Bases: Module

Xray surface brightness from a 2D surface brightness Beta-model

Source code in src/xsb_fluc/physics/surface_brightness.py 
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
class XraySurfaceBrightnessBetaModel(hk.Module):
    """
    Xray surface brightness from a 2D surface brightness Beta-model
    """

    def __call__(self, r):
        log_bkg = hk.get_parameter("log_bkg", [], init=Constant(-5.))
        log_e_0 = hk.get_parameter("log_e_0", [], init=Constant(-4))
        log_r_c = hk.get_parameter("log_r_c", [], init=Constant(-1.))
        beta = hk.get_parameter("beta", [], init=Constant(2 / 3))

        e_0 = 10**log_e_0
        r_c = 10**log_r_c

        u_beta = jnp.sqrt(jnp.pi)*jnp.exp(gammaln(3*beta - 1/2)-gammaln(3*beta))

        return u_beta*e_0*(1+(r/r_c)**2)**(1/2-3*beta) + 10 ** log_bkg

xsb_fluc.physics.ellipse

Python module containing models for the ellipticity of projected surface brightness.

EllipseRadius

Bases: Module

Include ellipticity in 2D radius maps

Source code in src/xsb_fluc/physics/ellipse.py 
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
38
39
40
41
42
43
44
45
46
47
class EllipseRadius(hk.Module):
    """
    Include ellipticity in 2D radius maps
    """

    def __init__(self, x_c, y_c, pixel_size):
        super(EllipseRadius, self).__init__()
        self.x_c_init = x_c
        self.y_c_init = y_c
        self.pixel_size = pixel_size

    @classmethod
    def from_data(cls, data: "Cluster"):
        """
        Build the ellipticity model from a `Cluster` object.
        """

        pixel_size = (data.kpc_per_pixel/data.r_500).to(1 / u.pix).value

        return cls(data.x_c, data.y_c, pixel_size)

    def __call__(self, x_ref, y_ref):
        """
        Compute the elliptical radius for a given position.

        Returns:
            (jnp.array): Elliptical radius in unit of $R_{500}$
        """
        angle = hk.get_parameter("angle", [], init=Constant(0.))
        e = hk.get_parameter("eccentricity", [], init=Constant(0.))
        x_c = self.x_c_init*(1+hk.get_parameter("x_c", [], init=Constant(0.)))
        y_c = self.y_c_init*(1+hk.get_parameter("y_c", [], init=Constant(0.)))

        x_tilde = (x_ref-x_c)*jnp.cos(angle) - (y_ref-y_c)*jnp.sin(angle)
        y_tilde = (y_ref-y_c)*jnp.cos(angle) + (x_ref-x_c)*jnp.sin(angle)
        r = jnp.sqrt(x_tilde**2*(1-e**2)**(-1/2) + y_tilde**2*(1-e**2)**(1/2))*self.pixel_size

        return r

__call__(x_ref, y_ref)

Compute the elliptical radius for a given position.

Returns:

Type Description
array

Elliptical radius in unit of \(R_{500}\)
Source code in src/xsb_fluc/physics/ellipse.py 
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
def __call__(self, x_ref, y_ref):
    """
    Compute the elliptical radius for a given position.

    Returns:
        (jnp.array): Elliptical radius in unit of $R_{500}$
    """
    angle = hk.get_parameter("angle", [], init=Constant(0.))
    e = hk.get_parameter("eccentricity", [], init=Constant(0.))
    x_c = self.x_c_init*(1+hk.get_parameter("x_c", [], init=Constant(0.)))
    y_c = self.y_c_init*(1+hk.get_parameter("y_c", [], init=Constant(0.)))

    x_tilde = (x_ref-x_c)*jnp.cos(angle) - (y_ref-y_c)*jnp.sin(angle)
    y_tilde = (y_ref-y_c)*jnp.cos(angle) + (x_ref-x_c)*jnp.sin(angle)
    r = jnp.sqrt(x_tilde**2*(1-e**2)**(-1/2) + y_tilde**2*(1-e**2)**(1/2))*self.pixel_size

    return r

from_data(data) classmethod

Build the ellipticity model from a Cluster object.

Source code in src/xsb_fluc/physics/ellipse.py 
21
22
23
24
25
26
27
28
29
@classmethod
def from_data(cls, data: "Cluster"):
    """
    Build the ellipticity model from a `Cluster` object.
    """

    pixel_size = (data.kpc_per_pixel/data.r_500).to(1 / u.pix).value

    return cls(data.x_c, data.y_c, pixel_size)

xsb_fluc.physics.turbulence

KolmogorovPowerSpectrum

Bases: Module

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/xsb_fluc/physics/turbulence.py 
 6
 7
 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
class KolmogorovPowerSpectrum(hk.Module):
    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}$$
    """

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

    def __call__(self, k):

        log_sigma = hk.get_parameter("log_sigma", [], init=Constant(-1.))
        log_inj = hk.get_parameter("log_inj", [], init=Constant(-0.3))
        log_dis = -3.
        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 = jnp.trapz(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