.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/astro/galactic_dynamics/plot_04_rubin_flat_rotation_curves.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

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

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

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

.. _sphx_glr_api_gallery_astro_galactic_dynamics_plot_04_rubin_flat_rotation_curves.py:


Rubin and Ford: flat galaxy rotation curves
===============================================

Rubin and Ford (1970) measured the Doppler shifts of emission regions
across the Andromeda galaxy, and through the 1970s they, Roberts,
Whitehurst and Bosma extended such curves far past the visible disk. If
the visible stars and gas were all the mass, the circular speed
:math:`v_c(r)=\sqrt{GM(<r)/r}` would fall as :math:`r^{-1/2}` once the
light runs out. Instead it stays flat.

This example builds a mock observed rotation curve, compares it with the
visible disk alone using
:func:`~physicskit.astro.galactic_dynamics.circular_velocity`, and shows
that the outer points require an extended halo, here an NFW halo from
:func:`~physicskit.astro.galactic_dynamics.nfw_enclosed_mass`. Units are
:math:`G=1`.

.. GENERATED FROM PYTHON SOURCE LINES 21-27

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from physicskit.astro.galactic_dynamics import circular_velocity, nfw_enclosed_mass
    from physicskit.astro.visualizers import plot_rotation_curve








.. GENERATED FROM PYTHON SOURCE LINES 28-34

Visible matter: an exponential disk
---------------------------------------
Surface brightness falls off as :math:`e^{-r/R_d}`. In the spherical
approximation its enclosed mass is
:math:`M_d[1-(1+r/R_d)e^{-r/R_d}]`, which saturates at :math:`M_d`
beyond a few scale lengths.

.. GENERATED FROM PYTHON SOURCE LINES 34-49

.. code-block:: Python

    M_disk, R_d = 60.0, 3.0


    def disk_mass(r):
        x = np.asarray(r, dtype=float) / R_d
        return M_disk * (1.0 - (1.0 + x) * np.exp(-x))


    rho_s, r_s = 0.03, 12.0


    def total_mass(r):
        return disk_mass(r) + nfw_enclosed_mass(r, rho_s, r_s)









.. GENERATED FROM PYTHON SOURCE LINES 50-55

Mock observations
---------------------
Twenty-five emission regions and 21-cm points, out to ten disk scale
lengths, with a 3% velocity error each. The "true" galaxy is disk plus
halo; the observer does not know that.

.. GENERATED FROM PYTHON SOURCE LINES 55-78

.. code-block:: Python

    rng = np.random.default_rng(1970)
    r_obs = np.linspace(1.0, 30.0, 25)
    v_true = circular_velocity(r_obs, total_mass)
    v_err = 0.03 * v_true
    v_obs = v_true + rng.normal(0.0, v_err)

    r = np.linspace(0.3, 32.0, 300)
    v_disk = circular_velocity(r, disk_mass)
    v_halo = circular_velocity(r, lambda x: nfw_enclosed_mass(x, rho_s, r_s))
    v_model = circular_velocity(r, total_mass)

    fig1, ax1 = plot_rotation_curve(r, v_model, v_observed=None)
    ax1.lines[0].set_label("disk + halo")
    ax1.errorbar(r_obs, v_obs, yerr=v_err, fmt="o", color="black", ms=4, label="observed")
    ax1.plot(r, v_disk, "--", color="firebrick", label="visible disk alone")
    ax1.plot(r, v_halo, ":", color="darkorchid", label="dark halo")
    ax1.axvspan(0, 4 * R_d, color="gold", alpha=0.15, label="optical disk")
    ax1.set_xlabel("radius r")
    ax1.set_ylabel(r"$v_c$")
    ax1.set_title("Rubin-Ford: the curve stays flat past the light")
    ax1.legend(fontsize=8, loc="lower right")
    fig1.tight_layout()




.. image-sg:: /api/gallery/astro/galactic_dynamics/images/sphx_glr_plot_04_rubin_flat_rotation_curves_001.png
   :alt: Rubin-Ford: the curve stays flat past the light
   :srcset: /api/gallery/astro/galactic_dynamics/images/sphx_glr_plot_04_rubin_flat_rotation_curves_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 79-84

The mass discrepancy grows with radius
------------------------------------------
Inverting :math:`v_c^2=GM/r` turns each observed speed into a dynamical
mass. A flat curve means :math:`M(<r)\propto r`: the enclosed mass keeps
growing long after the visible mass has stopped.

.. GENERATED FROM PYTHON SOURCE LINES 84-102

.. code-block:: Python

    M_dyn = r_obs * v_obs**2
    ratio = M_dyn / disk_mass(r_obs)
    outer = r_obs > 4 * R_d
    print(f"outer-curve slope d ln v / d ln r = {np.polyfit(np.log(r_obs[outer]), np.log(v_obs[outer]), 1)[0]:+.3f}  (Keplerian: -0.5)")
    for ri, Mi, q in zip(r_obs[::6], M_dyn[::6], ratio[::6]):
        print(f"r = {ri:5.1f}:  M_dyn = {Mi:7.1f},  M_dyn / M_visible = {q:5.2f}")

    fig2, ax2 = plt.subplots(figsize=(5.5, 3.8))
    ax2.plot(r_obs, M_dyn, "o", color="black", ms=4, label=r"dynamical $r v_c^2/G$")
    ax2.plot(r, disk_mass(r), "--", color="firebrick", label="visible disk")
    ax2.set_xlabel("radius r")
    ax2.set_ylabel("enclosed mass")
    ax2.set_title("Mass keeps rising where light stops")
    ax2.legend(fontsize=8)
    fig2.tight_layout()

    plt.show()




.. image-sg:: /api/gallery/astro/galactic_dynamics/images/sphx_glr_plot_04_rubin_flat_rotation_curves_002.png
   :alt: Mass keeps rising where light stops
   :srcset: /api/gallery/astro/galactic_dynamics/images/sphx_glr_plot_04_rubin_flat_rotation_curves_002.png
   :class: sphx-glr-single-img


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

 .. code-block:: none

    outer-curve slope d ln v / d ln r = -0.014  (Keplerian: -0.5)
    r =   1.0:  M_dyn =     4.5,  M_dyn / M_visible =  1.69
    r =   8.2:  M_dyn =   114.0,  M_dyn / M_visible =  2.50
    r =  15.5:  M_dyn =   232.9,  M_dyn / M_visible =  4.02
    r =  22.8:  M_dyn =   346.1,  M_dyn / M_visible =  5.79
    r =  30.0:  M_dyn =   398.6,  M_dyn / M_visible =  6.65




.. GENERATED FROM PYTHON SOURCE LINES 103-107

Check
-----
Beyond the optical disk the curve stays flat (slope ~0, not Keplerian
-1/2), so the dynamical mass keeps growing past the visible mass.

.. GENERATED FROM PYTHON SOURCE LINES 107-110

.. code-block:: Python

    outer_slope = np.polyfit(np.log(r_obs[outer]), np.log(v_obs[outer]), 1)[0]
    assert abs(outer_slope) < 0.1
    assert ratio[-1] > 5 and ratio[-1] > 3 * ratio[0]








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

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


.. _sphx_glr_download_api_gallery_astro_galactic_dynamics_plot_04_rubin_flat_rotation_curves.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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