.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/c_fields/plot_elem_lap_nfw.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_c_fields_plot_elem_lap_nfw.py>`
        to download the full example code.

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

.. _sphx_glr_auto_examples_c_fields_plot_elem_lap_nfw.py:


==========================
NFW Density from Potential
==========================

Compute the NFW density from a potential profile.

.. GENERATED FROM PYTHON SOURCE LINES 9-22

This example demonstrates how to compute the dark matter density profile
of a Navarro–Frenk–White (NFW) halo by applying the Laplacian operator
to its gravitational potential in spherical coordinates.

We use the Poisson equation to numerically recover the NFW density profile:

.. math::

    \nabla^2 \Phi = 4\pi G \rho

Given the analytic form of the NFW potential :math:`\Phi(r)`, we compute the Laplacian
numerically using PyMetric and compare the inferred density :math:`\rho(r)` to its
known analytic form.

.. GENERATED FROM PYTHON SOURCE LINES 22-35

.. code-block:: Python


    import matplotlib.pyplot as plt

    # Imports
    import numpy as np

    from pymetric import DenseField, GenericGrid, SphericalCoordinateSystem

    # Characteristic density and scale radius for the NFW profile
    rho0 = 1.0
    Rs = 1.0









.. GENERATED FROM PYTHON SOURCE LINES 36-41

Define the NFW gravitational potential.

.. math::
   \Phi(r) = -\frac{4\pi G \rho_0 R_s^3}{r} \log\left(1 + \frac{r}{R_s}\right)


.. GENERATED FROM PYTHON SOURCE LINES 41-45

.. code-block:: Python

    def nfw_potential(r, _rho0=rho0, _Rs=Rs):
        return -(4 * np.pi * _rho0 * _Rs**3 / r) * np.log(1 + r / _Rs)









.. GENERATED FROM PYTHON SOURCE LINES 46-51

Define the analytic NFW density profile derived from the Poisson equation.

.. math::
   \rho(r) = \frac{\rho_0}{\left(\frac{r}{R_s}\right)(1 + \frac{r}{R_s})^2}


.. GENERATED FROM PYTHON SOURCE LINES 51-56

.. code-block:: Python

    def nfw_density(r, _rho0=rho0, _Rs=Rs):
        xi = r / _Rs
        return _rho0 / (xi * (1 + xi) ** 2)









.. GENERATED FROM PYTHON SOURCE LINES 57-61

Setup the grid using a spherical coordinate system.

We use high radial resolution to accurately capture structure in :math:`r`,
and minimal angular resolution since the profile is spherically symmetric.

.. GENERATED FROM PYTHON SOURCE LINES 61-70

.. code-block:: Python

    csys = SphericalCoordinateSystem()

    r_coord = np.geomspace(1.0, 1e4, 3000)
    theta_coord = np.linspace(0, np.pi, 10)
    phi_coord = np.linspace(0, 2 * np.pi, 10)

    grid = GenericGrid(csys, [r_coord, theta_coord, phi_coord], center="vertex")
    grid.fill_values = {"r": 1, "theta": 1, "phi": 1}  # handle r=0 and boundaries safely








.. GENERATED FROM PYTHON SOURCE LINES 71-72

Evaluate the NFW potential on the grid.

.. GENERATED FROM PYTHON SOURCE LINES 72-74

.. code-block:: Python

    field: DenseField = DenseField.from_function(nfw_potential, grid, axes=["r"])








.. GENERATED FROM PYTHON SOURCE LINES 75-76

Visualize the magnitude of the potential vs radius.

.. GENERATED FROM PYTHON SOURCE LINES 76-85

.. code-block:: Python

    plt.figure(figsize=(6, 4))
    plt.loglog(r_coord, np.abs(field[...]), label=r"$|\Phi(r)|$")
    plt.xlabel("r")
    plt.ylabel("Potential Magnitude")
    plt.title("NFW Gravitational Potential")
    plt.legend()
    plt.tight_layout()
    plt.show()




.. image-sg:: /auto_examples/c_fields/images/sphx_glr_plot_elem_lap_nfw_001.png
   :alt: NFW Gravitational Potential
   :srcset: /auto_examples/c_fields/images/sphx_glr_plot_elem_lap_nfw_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 86-91

Compute the numerical Laplacian of the field and infer the density using:

.. math::
   \rho(r) = \frac{1}{4\pi G} \nabla^2 \Phi


.. GENERATED FROM PYTHON SOURCE LINES 91-94

.. code-block:: Python

    laplacian = field.element_wise_laplacian()
    numerical_density = laplacian[...] / (4 * np.pi)








.. GENERATED FROM PYTHON SOURCE LINES 95-96

Compare the numerical density to the analytic NFW density profile.

.. GENERATED FROM PYTHON SOURCE LINES 96-105

.. code-block:: Python

    plt.figure(figsize=(6, 4))
    plt.loglog(r_coord, numerical_density, label="Numerical $\\rho(r)$", lw=2)
    plt.loglog(r_coord, nfw_density(r_coord), "--", label="Analytic $\\rho(r)$", lw=2)
    plt.xlabel("r")
    plt.ylabel("Density")
    plt.title("NFW Density from Laplacian of Potential")
    plt.legend()
    plt.tight_layout()
    plt.show()



.. image-sg:: /auto_examples/c_fields/images/sphx_glr_plot_elem_lap_nfw_002.png
   :alt: NFW Density from Laplacian of Potential
   :srcset: /auto_examples/c_fields/images/sphx_glr_plot_elem_lap_nfw_002.png
   :class: sphx-glr-single-img






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

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


.. _sphx_glr_download_auto_examples_c_fields_plot_elem_lap_nfw.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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