.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/u_particles/plot_galaxy_cluster_particles.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_u_particles_plot_galaxy_cluster_particles.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_u_particles_plot_galaxy_cluster_particles.py:


==================================
Particle Datasets: Galaxy Clusters
==================================

This example demonstrates how to generate particle datasets for a
galaxy cluster model using Pisces. To do so, we'll use the built-in
hook for generating particles, which gives us access to the
:meth:`~models.galaxy_clusters.spherical.SphericalGalaxyClusterModel.generate_particles`
method.

.. GENERATED FROM PYTHON SOURCE LINES 14-22

Setup & Imports
----------------
The first thing that needs to happen is importing the
relevant building blocks. For this, we'll need the
:class:`~models.galaxy_clusters.spherical.SphericalGalaxyClusterModel` class
which will be our model and the :class:`~profiles.density.HernquistDensityProfile`, which
we'll use to represent the total and gas density profiles of the cluster.


.. GENERATED FROM PYTHON SOURCE LINES 22-31

.. code-block:: Python

    import tempfile

    import matplotlib.pyplot as plt
    import numpy as np
    import unyt

    from pisces.models.galaxy_clusters import SphericalGalaxyClusterModel
    from pisces.profiles.density import HernquistDensityProfile








.. GENERATED FROM PYTHON SOURCE LINES 32-39

Building the Model
------------------
First things first, we need to build the model. In this case, we'll be using
basic Hernquist profiles :footcite:p:`HernquistProfile` for both the gas and total density.
This is nice for simple models because we can fix the total mass easily. We'll
create a galaxy cluster with 85\% of the total mass in dark matter and 15\% in gas.
The total mass of the cluster will be :math:`10^{15} \;{\rm M_\odot}`.

.. GENERATED FROM PYTHON SOURCE LINES 39-80

.. code-block:: Python


    # Define some parameters for the total density profile.
    total_mass = unyt.unyt_quantity(1e15, "Msun")
    gas_mass, dm_mass = total_mass * 0.15, total_mass * 0.85

    gas_scale_radius = unyt.unyt_quantity(220, "kpc")
    dm_scale_radius = unyt.unyt_quantity(200, "kpc")

    gas_rho_0 = gas_mass / (2 * np.pi * gas_scale_radius**3)
    dm_rho_0 = dm_mass / (2 * np.pi * dm_scale_radius**3)

    # Create the density profiles.
    gas_density_profile = HernquistDensityProfile(
        rho_0=gas_rho_0,
        r_s=gas_scale_radius,
    )
    total_density_profile = HernquistDensityProfile(
        rho_0=dm_rho_0,
        r_s=dm_scale_radius,
    )

    # Radial range and resolution
    rmin = unyt.unyt_quantity(1.0, "kpc")
    rmax = unyt.unyt_quantity(3.0, "Mpc")

    # Create a temporary directory to store the model file and
    # give the model a filename.
    tmpdir = tempfile.TemporaryDirectory()
    filename = f"{tmpdir.name}/cluster_temp_dens_model.h5"

    # Now we can move forward with making the model!
    model = SphericalGalaxyClusterModel.from_density_and_total_density(
        gas_density_profile,
        total_density_profile,
        filename,
        min_radius=rmin,
        max_radius=rmax,
        num_points=1000,
        overwrite=True,
    )





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    Preparing metadata...:   0%|          | 0/7 [00:00<?, ?it/s]
    Generating grid...:  14%|█▍        | 1/7 [00:00<00:00, 97541.95it/s]
    Building density...:  29%|██▊       | 2/7 [00:00<00:00, 2238.16it/s]
    Computing masses...:  43%|████▎     | 3/7 [00:00<00:00, 2184.53it/s]
    Calculating potential...:  57%|█████▋    | 4/7 [00:00<00:00, 115.34it/s]
    Calculating potential...:  71%|███████▏  | 5/7 [00:00<00:00,  7.66it/s] 
    Calculating temperature...:  71%|███████▏  | 5/7 [00:00<00:00,  7.66it/s]
    Calculating temperature...:  86%|████████▌ | 6/7 [00:18<00:03,  3.98s/it]
    Performing additional computations...:  86%|████████▌ | 6/7 [00:18<00:03,  3.98s/it]
    COMPLETE:  86%|████████▌ | 6/7 [00:18<00:03,  3.98s/it]                             
                                                       



.. GENERATED FROM PYTHON SOURCE LINES 81-91

Generating Particles
--------------------
With the model built, we can now generate the particles for the cluster. To
do this we'll use the :meth:`~models.galaxy_clusters.spherical.SphericalGalaxyClusterModel.generate_particles`
method, which allows us to specify the number of particles for each component.

Internally, this will serve to sample the particles from our density profiles,
position them in 3D space, interpolate fields onto them, and then write them to
the output file. For this example, we'll generate :math:`10^5` particles for the gas component
and :math:`10^5` particles for the dark matter component.

.. GENERATED FROM PYTHON SOURCE LINES 91-97

.. code-block:: Python

    particles = model.generate_particles(
        filename=f"{tmpdir.name}/cluster_particles.h5",
        num_particles={"gas": 100_000, "dark_matter": 100_000},
        overwrite=True,
    )





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    Generating particles:   0%|          | 0/2 [00:00<?, ?species/s]
                                                                
    Generating 100000 gas particles...

    Generating particles:   0%|          | 0/2 [00:00<?, ?species/s]

    Interpolating gas fields:   0%|          | 0/8 [00:00<?, ?fields/s]

                                                                       
    Generating particles:  50%|█████     | 1/2 [00:00<00:00,  8.70species/s]
                                                                        
    Generating 100000 dark_matter particles...

    Generating particles:  50%|█████     | 1/2 [00:00<00:00,  8.70species/s]

    Interpolating dark_matter fields:   0%|          | 0/3 [00:00<?, ?fields/s]

                                                                               
                                                                        
    DF for 'dark_matter' not found, attempting to compute it...

    Generating particles:  50%|█████     | 1/2 [00:00<00:00,  8.70species/s]

    Generating particle velocities...:   0%|          | 0/100000 [00:00<?, ?it/s]

    Generating particle velocities...:   5%|▌         | 5362/100000 [00:00<00:01, 53619.80it/s]

    Generating particle velocities...:  11%|█         | 10779/100000 [00:00<00:01, 53943.32it/s]

    Generating particle velocities...:  16%|█▌        | 16225/100000 [00:00<00:01, 54178.80it/s]

    Generating particle velocities...:  22%|██▏       | 21643/100000 [00:00<00:01, 53781.66it/s]

    Generating particle velocities...:  27%|██▋       | 27105/100000 [00:00<00:01, 54080.16it/s]

    Generating particle velocities...:  33%|███▎      | 32565/100000 [00:00<00:01, 54254.35it/s]

    Generating particle velocities...:  38%|███▊      | 38027/100000 [00:00<00:01, 54371.92it/s]

    Generating particle velocities...:  43%|████▎     | 43484/100000 [00:00<00:01, 54433.54it/s]

    Generating particle velocities...:  49%|████▉     | 48958/100000 [00:00<00:00, 54528.85it/s]

    Generating particle velocities...:  54%|█████▍    | 54420/100000 [00:01<00:00, 54554.94it/s]

    Generating particle velocities...:  60%|█████▉    | 59876/100000 [00:01<00:00, 54437.56it/s]

    Generating particle velocities...:  65%|██████▌   | 65320/100000 [00:01<00:00, 54356.15it/s]

    Generating particle velocities...:  71%|███████   | 70756/100000 [00:01<00:00, 54284.70it/s]

    Generating particle velocities...:  76%|███████▌  | 76185/100000 [00:01<00:00, 54026.85it/s]

    Generating particle velocities...:  82%|████████▏ | 81588/100000 [00:01<00:00, 53956.13it/s]

    Generating particle velocities...:  87%|████████▋ | 87004/100000 [00:01<00:00, 54014.50it/s]

    Generating particle velocities...:  92%|█████████▏| 92408/100000 [00:01<00:00, 54022.08it/s]

    Generating particle velocities...:  98%|█████████▊| 97814/100000 [00:01<00:00, 54030.64it/s]

                                                                                                
    Generating particles: 100%|██████████| 2/2 [00:02<00:00,  1.53s/species]
    Generating particles: 100%|██████████| 2/2 [00:02<00:00,  1.32s/species]




.. GENERATED FROM PYTHON SOURCE LINES 98-107

Depending on your machine, this may take a few seconds to run as the velocity sampling
is a somewhat expensive operation. Once the particles are generated, we can do any number
of things with them, such as visualizing the density field, plotting the temperature profile,
or even running a simulation with them.

For this example, we'll simply visualize the density-weighted gas temperature profile,
which is similar to what one might see in an X-ray observation of a galaxy cluster. To do this,
we'll need to extract the gas temperatures, particle positions, and weights from the dataset. We
can then use histograms to visualize the temperature distribution.

.. GENERATED FROM PYTHON SOURCE LINES 107-135

.. code-block:: Python


    # Extract the gas temperature, positions, and weights.
    temperature = particles["gas.kT"].to_value("keV")
    positions = particles["gas.particle_position"].to_value("kpc")
    density = particles["gas.density"].to_value("Msun/kpc**3")

    # Create the bins.
    bins = np.linspace(-1000, 1000, 601)

    # Create the histogram weighted by density * temperature.
    # Then create the density only histogram. We can then divide
    # these to get the mass-weighted temperature.
    hist_temp_dens, _, _ = np.histogram2d(positions[:, 0], positions[:, 1], bins=bins, weights=temperature * density)
    hist_dens, _, _ = np.histogram2d(positions[:, 0], positions[:, 1], bins=bins, weights=density)

    # Calculate the mass-weighted temperature.
    mass_weighted_temp = np.zeros_like(hist_temp_dens)
    mass_weighted_temp[hist_dens > 0] = hist_temp_dens[hist_dens > 0] / hist_dens[hist_dens > 0]

    # Create the plot.
    fig, ax = plt.subplots(figsize=(8, 6))
    c = ax.pcolormesh(bins, bins, mass_weighted_temp.T, shading="auto", cmap="inferno")
    ax.set_xlabel("X Position (kpc)")
    ax.set_ylabel("Y Position (kpc)")
    ax.set_title("Mass-Weighted Gas Temperature Distribution")
    fig.colorbar(c, ax=ax, label="Temperature (keV)")
    plt.tight_layout()
    plt.show()



.. image-sg:: /auto_examples/u_particles/images/sphx_glr_plot_galaxy_cluster_particles_001.png
   :alt: Mass-Weighted Gas Temperature Distribution
   :srcset: /auto_examples/u_particles/images/sphx_glr_plot_galaxy_cluster_particles_001.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 21.165 seconds)


.. _sphx_glr_download_auto_examples_u_particles_plot_galaxy_cluster_particles.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_galaxy_cluster_particles.ipynb <plot_galaxy_cluster_particles.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_galaxy_cluster_particles.py <plot_galaxy_cluster_particles.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_galaxy_cluster_particles.zip <plot_galaxy_cluster_particles.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_