Build a Polytropic Star Model

This example shows how to use the 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.

import tempfile

import matplotlib.pyplot as plt

Imports

import numpy as np
from unyt import unyt_quantity

from pisces.models.stars.polytropes import PolytropicStarModel

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

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).

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"),
        )
    )

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]

Visualize the Stellar Structure

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

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()