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

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

.. _sphx_glr_auto_examples_a_models_plot_polytropic_star.py:


=============================
Build a Polytropic Star Model
=============================

This example shows how to use the :class:`~pisces.models.stars.polytropes.PolytropicStarModel`
to generate the internal structure of a star governed by a polytropic equation of state.

We demonstrate two ways to construct a model:

1. From the total **mass and radius** of the star.
2. From the **central density and central temperature** of the core.

The resulting model outputs physical profiles for density, pressure, temperature, gravitational
potential, and mass.

.. GENERATED FROM PYTHON SOURCE LINES 17-22

.. code-block:: Python


    import tempfile

    import matplotlib.pyplot as plt








.. GENERATED FROM PYTHON SOURCE LINES 23-25

Imports
-------

.. GENERATED FROM PYTHON SOURCE LINES 25-33

.. code-block:: Python

    import numpy as np
    from unyt import unyt_quantity

    from pisces.models.stars.polytropes import PolytropicStarModel

    # Create a temp directory.
    tmpdir = tempfile.TemporaryDirectory()








.. GENERATED FROM PYTHON SOURCE LINES 34-38

Generate a Polytropic Star from Mass and Radius
-----------------------------------------------
Here we create a model for a star with 1 solar mass and 1 solar radius,
assuming a polytropic index n = 1.5 (typical for fully convective stars).

.. GENERATED FROM PYTHON SOURCE LINES 38-59

.. code-block:: Python

    mass = unyt_quantity(1, "Msun")
    radius = unyt_quantity(1, "Rsun")
    n_index = [0.2, 0.3, 0.5, 1, 1.25, 1.5, 2, 2.5, 3]

    models = []

    for n in n_index:
        filename = f"{tmpdir.name}/polytrope_example_{n}.h5"

        models.append(
            PolytropicStarModel.from_mass_and_radius(
                filename=filename,
                mass=mass,
                radius=radius,
                polytropic_index=n,
                overwrite=True,
                rmax=unyt_quantity(10, "Rsun"),
                rmin=unyt_quantity(0.01, "Rsun"),
            )
        )





.. 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, 107546.26it/s]
    Computing model constants...:  29%|██▊       | 2/7 [00:00<00:00, 2202.31it/s]
    Solving Lane-Emden Eq...:  43%|████▎     | 3/7 [00:00<00:00, 159.30it/s]     /home/runner/work/Pisces/Pisces/pisces/math_utils/integration.py:237: RuntimeWarning: invalid value encountered in scalar power
      dpsi_dxi = -(theta**n) * xi**2

    Computing fields...:  57%|█████▋    | 4/7 [00:00<00:00, 169.57it/s]     
    Computing gravitational potential...:  71%|███████▏  | 5/7 [00:00<00:00, 100.04it/s]
    Computing gravitational potential...:  86%|████████▌ | 6/7 [00:01<00:00,  4.95it/s] 
    Finalizing model...:  86%|████████▌ | 6/7 [00:01<00:00,  4.95it/s]                 
                                                                  

    Preparing metadata...:   0%|          | 0/7 [00:00<?, ?it/s]
    Generating grid...:  14%|█▍        | 1/7 [00:00<00:00, 139810.13it/s]
    Computing model constants...:  29%|██▊       | 2/7 [00:00<00:00, 2370.33it/s]
    Solving Lane-Emden Eq...:  43%|████▎     | 3/7 [00:00<00:00, 153.54it/s]     
    Computing fields...:  57%|█████▋    | 4/7 [00:00<00:00, 168.89it/s]     
    Computing gravitational potential...:  71%|███████▏  | 5/7 [00:00<00:00, 120.64it/s]
    Computing gravitational potential...:  86%|████████▌ | 6/7 [00:01<00:00,  5.20it/s] 
    Finalizing model...:  86%|████████▌ | 6/7 [00:01<00:00,  5.20it/s]                 
                                                                  

    Preparing metadata...:   0%|          | 0/7 [00:00<?, ?it/s]
    Generating grid...:  14%|█▍        | 1/7 [00:00<00:00, 113359.57it/s]
    Computing model constants...:  29%|██▊       | 2/7 [00:00<00:00, 2390.60it/s]
    Solving Lane-Emden Eq...:  43%|████▎     | 3/7 [00:00<00:00, 242.77it/s]     
    Computing fields...:  57%|█████▋    | 4/7 [00:00<00:00, 238.95it/s]     
    Computing gravitational potential...:  71%|███████▏  | 5/7 [00:00<00:00, 197.52it/s]
    Computing gravitational potential...:  86%|████████▌ | 6/7 [00:01<00:00,  5.48it/s] 
    Finalizing model...:  86%|████████▌ | 6/7 [00:01<00:00,  5.48it/s]                 
                                                                  

    Preparing metadata...:   0%|          | 0/7 [00:00<?, ?it/s]
    Generating grid...:  14%|█▍        | 1/7 [00:00<00:00, 119837.26it/s]
    Computing model constants...:  29%|██▊       | 2/7 [00:00<00:00, 2438.55it/s]
    Solving Lane-Emden Eq...:  43%|████▎     | 3/7 [00:00<00:00, 242.07it/s]     
    Computing fields...:  57%|█████▋    | 4/7 [00:00<00:00, 234.58it/s]     
    Computing gravitational potential...:  71%|███████▏  | 5/7 [00:00<00:00, 231.70it/s]
    Computing gravitational potential...:  86%|████████▌ | 6/7 [00:00<00:00, 23.85it/s] 
    Finalizing model...:  86%|████████▌ | 6/7 [00:00<00:00, 23.85it/s]                 
                                                                  

    Preparing metadata...:   0%|          | 0/7 [00:00<?, ?it/s]
    Generating grid...:  14%|█▍        | 1/7 [00:00<00:00, 116508.44it/s]
    Computing model constants...:  29%|██▊       | 2/7 [00:00<00:00, 2453.53it/s]
    Solving Lane-Emden Eq...:  43%|████▎     | 3/7 [00:00<00:00, 164.68it/s]     
    Computing fields...:  57%|█████▋    | 4/7 [00:00<00:00, 171.19it/s]     
    Computing gravitational potential...:  71%|███████▏  | 5/7 [00:00<00:00, 139.97it/s]
    Computing gravitational potential...:  86%|████████▌ | 6/7 [00:00<00:00,  8.17it/s] 
    Finalizing model...:  86%|████████▌ | 6/7 [00:00<00:00,  8.17it/s]                 
                                                                  

    Preparing metadata...:   0%|          | 0/7 [00:00<?, ?it/s]
    Generating grid...:  14%|█▍        | 1/7 [00:00<00:00, 135300.13it/s]
    Computing model constants...:  29%|██▊       | 2/7 [00:00<00:00, 2311.55it/s]
    Solving Lane-Emden Eq...:  43%|████▎     | 3/7 [00:00<00:00, 214.43it/s]     
    Computing fields...:  57%|█████▋    | 4/7 [00:00<00:00, 205.38it/s]     
    Computing gravitational potential...:  71%|███████▏  | 5/7 [00:00<00:00, 190.09it/s]
    Computing gravitational potential...:  86%|████████▌ | 6/7 [00:00<00:00, 10.59it/s] 
    Finalizing model...:  86%|████████▌ | 6/7 [00:00<00:00, 10.59it/s]                 
                                                                  

    Preparing metadata...:   0%|          | 0/7 [00:00<?, ?it/s]
    Generating grid...:  14%|█▍        | 1/7 [00:00<00:00, 139810.13it/s]
    Computing model constants...:  29%|██▊       | 2/7 [00:00<00:00, 2431.48it/s]
    Solving Lane-Emden Eq...:  43%|████▎     | 3/7 [00:00<00:00, 173.84it/s]     
    Computing fields...:  57%|█████▋    | 4/7 [00:00<00:00, 171.72it/s]     
    Computing gravitational potential...:  71%|███████▏  | 5/7 [00:00<00:00, 149.20it/s]
    Computing gravitational potential...:  86%|████████▌ | 6/7 [00:00<00:00, 22.89it/s] 
    Finalizing model...:  86%|████████▌ | 6/7 [00:00<00:00, 22.89it/s]                 
                                                                  

    Preparing metadata...:   0%|          | 0/7 [00:00<?, ?it/s]
    Generating grid...:  14%|█▍        | 1/7 [00:00<00:00, 102300.10it/s]
    Computing model constants...:  29%|██▊       | 2/7 [00:00<00:00, 1822.42it/s]
    Solving Lane-Emden Eq...:  43%|████▎     | 3/7 [00:00<00:00, 150.28it/s]     
    Computing fields...:  57%|█████▋    | 4/7 [00:00<00:00, 146.54it/s]     
    Computing gravitational potential...:  71%|███████▏  | 5/7 [00:00<00:00, 123.48it/s]
    Computing gravitational potential...:  86%|████████▌ | 6/7 [00:00<00:00, 19.11it/s] 
    Finalizing model...:  86%|████████▌ | 6/7 [00:00<00:00, 19.11it/s]                 
                                                                  

    Preparing metadata...:   0%|          | 0/7 [00:00<?, ?it/s]
    Generating grid...:  14%|█▍        | 1/7 [00:00<00:00, 91180.52it/s]
    Computing model constants...:  29%|██▊       | 2/7 [00:00<00:00, 2418.86it/s]
    Solving Lane-Emden Eq...:  43%|████▎     | 3/7 [00:00<00:00, 250.15it/s]     
    Computing fields...:  57%|█████▋    | 4/7 [00:00<00:00, 192.43it/s]     
    Computing gravitational potential...:  71%|███████▏  | 5/7 [00:00<00:00, 197.89it/s]
    Computing gravitational potential...:  86%|████████▌ | 6/7 [00:00<00:00, 16.28it/s] 
    Finalizing model...:  86%|████████▌ | 6/7 [00:00<00:00, 16.28it/s]                 
                                                                  



.. GENERATED FROM PYTHON SOURCE LINES 60-64

Visualize the Stellar Structure
-------------------------------
We now plot four key physical fields as a function of radius:
density, temperature, pressure, and gravitational field.

.. GENERATED FROM PYTHON SOURCE LINES 64-104

.. code-block:: Python

    fig, axs = plt.subplots(2, 2, figsize=(12, 8), sharex=True)
    cmap = plt.cm.cool
    norm = plt.Normalize(vmin=min(n_index), vmax=max(n_index))

    for model, n in zip(models, n_index):
        # Scale n down to 0-1 for color.
        color = cmap(norm(n))

        # Add the plots.
        axs[0, 0].plot(model.grid["r"].to("Rsun").value, model["density"].to("g/cm**3").value, color=color)
        axs[0, 1].plot(model.grid["r"].to("Rsun").value, model["temperature"].to("K").value, color=color)
        axs[1, 0].plot(model.grid["r"].to("Rsun").value, model["mass"].to("Msun").value, color=color)
        axs[1, 1].plot(model.grid["r"].to("Rsun").value, model["gravitational_field"].to("cm/s**2").value, color=color)

    # Scales
    for ax in axs.ravel():
        ax.set_xscale("log")
        ax.set_yscale("log")


    # Axis labels
    axs[0, 0].set_ylabel(r"Density [${\rm g\;cm^{-3}}$]")
    axs[0, 1].set_ylabel(r"Temperature [${\rm K}$]")
    axs[1, 0].set_ylabel(r"Enclosed Mass [${\rm M_\odot}$]")
    axs[1, 1].set_ylabel(r"Gravitational Field [${\rm cm\;s^{-2}}$]")
    for ax in axs[1]:
        ax.set_xlabel(r"Radius [${\rm R_\odot}$]")

    # Titles
    axs[0, 0].set_title("Density Profile")
    axs[0, 1].set_title("Temperature Profile")
    axs[1, 0].set_title("Mass Profile")
    axs[1, 1].set_title("Gravitational Field")

    # Add colorbar for polytropic index
    sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
    cbar = fig.colorbar(sm, ax=axs, orientation="vertical", fraction=0.025, pad=0.02)
    cbar.set_label("Polytropic Index (n)")

    plt.show()



.. image-sg:: /auto_examples/a_models/images/sphx_glr_plot_polytropic_star_001.png
   :alt: Density Profile, Temperature Profile, Mass Profile, Gravitational Field
   :srcset: /auto_examples/a_models/images/sphx_glr_plot_polytropic_star_001.png
   :class: sphx-glr-single-img






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

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


.. _sphx_glr_download_auto_examples_a_models_plot_polytropic_star.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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