physicskit.astro#

physicskit.astro hub · History · Examples · API reference

physicskit.astro: stellar structure, N-body dynamics, orbital mechanics, and galactic dynamics.

Uses gravitational units with \(G=1\) by default throughout (every function that needs it accepts G as a keyword argument), the same convention physicskit.relativity uses for its own geometrized units.

class physicskit.astro.NBodySystem(positions, velocities, masses, G=1.0, softening=0.0)[source]#

Bases: object

A stateful direct-summation N-body gravitational simulation.

Parameters:

  • positions (ndarray of shape (N, 3))

  • velocities (ndarray of shape (N, 3))

  • masses (ndarray of shape (N,))

  • G (float, default=1.0)

  • softening (float, default=0.0)

simulate(dt, n_steps)[source]#

Advance n_steps steps, returning the position history.

Returns:

ndarray of shape (n_steps + 1, N, 3)

step(dt)[source]#

Advance the system’s stored state by one leapfrog step, in place.

total_angular_momentum()[source]#

Total angular momentum, \(\vec L=\sum_im_i\,\vec r_i\times\vec v_i\).

total_energy()[source]#

Kinetic plus (softened) gravitational potential energy, at the current state.

class physicskit.astro.PolytropicStar(n, K, rho_c, G=1.0)[source]#

Bases: object

A physical polytropic star built from one Lane-Emden solution.

Parameters:

  • n (float) – Polytropic index, \(P=K\rho^{1+1/n}\).

  • K (float) – Polytropic constant.

  • rho_c (float) – Central density.

  • G (float, default=1.0) – Gravitational constant.

Raises:

ValueError – If the Lane-Emden solution has no surface (\(n \ge 5\), where the polytrope’s radius is infinite).

Examples

>>> star = PolytropicStar(n=0.0, K=1.0, rho_c=1.0)
>>> round(star.mass / (4 / 3 * np.pi * star.radius**3 * star.rho_c), 3)
1.0
property alpha#

Physical length scale, \(\alpha=\sqrt{(n+1)K\rho_c^{1/n-1}/(4\pi G)}\).

property mass#

Physical stellar mass, \(M=4\pi\alpha^3\rho_c\left[-\xi_1^2\theta'(\xi_1)\right]\).

property radius#

Physical stellar radius, \(R=\alpha\xi_1\).

property xi1#

The dimensionless surface radius \(\xi_1\), the first zero of the Lane-Emden solution.

physicskit.astro.chandrasekhar_mass(mu_e=2.0)[source]#

The Chandrasekhar mass limit, in solar masses.

\[M_{\rm Ch} \approx \frac{5.83}{\mu_e^2}\ M_\odot\]

(the standard numerical prefactor from the n=3 polytrope / relativistic-electron-degeneracy treatment; see e.g. Kippenhahn, Weigert & Weiss, Stellar Structure and Evolution).

Parameters:

mu_e (float, default=2.0) – Mean molecular weight per electron (2.0 for carbon/oxygen white dwarfs).

Returns:

float – Mass limit, in solar masses.

Examples

>>> M = chandrasekhar_mass(2.0)
>>> 1.2 < M < 1.5
True
physicskit.astro.circular_velocity(r, mass_enclosed_func, G=1.0)[source]#

Circular orbital speed for any enclosed-mass profile, \(v_c(r)=\sqrt{GM(<r)/r}\).

Parameters:

  • r (array_like) – Radius.

  • mass_enclosed_func (callable) – mass_enclosed_func(r) -> mass, e.g. lambda r: nfw_enclosed_mass(r, rho_s, r_s).

  • G (float, default=1.0) – Gravitational constant.

Returns:

ndarray or float

Examples

>>> M0 = 5.0
>>> round(float(circular_velocity(2.0, lambda r: M0, G=1.0)), 6)
1.581139
physicskit.astro.gravitational_acceleration(positions, masses, G=1.0, softening=0.0)[source]#

Pairwise Plummer-softened gravitational acceleration.

\[\vec a_i = G\sum_{j\neq i} m_j\, \frac{\vec r_j-\vec r_i}{\left(|\vec r_j-\vec r_i|^2+\epsilon^2\right)^{3/2}}\]
Parameters:

  • positions (ndarray of shape (N, 3))

  • masses (ndarray of shape (N,))

  • G (float, default=1.0)

  • softening (float, default=0.0) – Plummer softening length \(\epsilon\).

Returns:

ndarray of shape (N, 3)

Examples

>>> import numpy as np
>>> pos = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]])
>>> m = np.array([1.0, 1.0])
>>> acc = gravitational_acceleration(pos, m)
>>> round(float(acc[0, 0]), 6)
1.0
physicskit.astro.hohmann_transfer(r1, r2, mu)[source]#

The two burns and transfer time for a Hohmann transfer.

Parameters:

  • r1 (float) – Initial circular orbit radius.

  • r2 (float) – Final circular orbit radius.

  • mu (float) – Gravitational parameter.

Returns:

  • delta_v1 (float) – Speed change (magnitude) to leave the initial circular orbit and enter the transfer ellipse.

  • delta_v2 (float) – Speed change (magnitude) to circularize at r2.

  • transfer_time (float) – Time to complete the transfer, half the transfer ellipse’s period.

Examples

>>> dv1, dv2, t = hohmann_transfer(1.0, 4.0, mu=1.0)
>>> dv1 > 0 and dv2 > 0 and t > 0
True
physicskit.astro.lane_emden(n, xi_max=1000.0, n_points=2000)[source]#

Integrate the Lane-Emden equation for a polytrope of index n.

\[\frac{1}{\xi^2}\frac{d}{d\xi}\left(\xi^2\frac{d\theta}{d\xi}\right) + \theta^n = 0, \qquad \theta(0)=1,\ \theta'(0)=0\]

the dimensionless hydrostatic-equilibrium profile of a self-gravitating sphere with equation of state \(P=K\rho^{1+1/n}\).

Parameters:

  • n (float) – Polytropic index.

  • xi_max (float, default=1000.0) – Upper integration bound. Integration normally stops earlier, at the surface; the default is large enough to reach it for every \(n \le 4.9\) (\(\xi_1 \approx 171\) there, diverging as \(n \to 5\)). For \(n \ge 5\) there is no surface and the profile is returned out to xi_max.

  • n_points (int, default=2000) – Number of points to sample the returned arrays at.

Returns:

  • xi (ndarray) – Dimensionless radius, from just above 0 to the surface (the first zero of \(\theta\)) or, if none was found, to xi_max.

  • theta (ndarray) – The Lane-Emden function \(\theta(\xi)\). theta[-1] == 0 exactly when the surface was reached.

Examples

>>> xi, theta = lane_emden(0.0)
>>> round(float(xi[-1]), 3)  # exact analytic surface is sqrt(6)
2.449
>>> xi, theta = lane_emden(4.0)
>>> round(float(xi[-1]), 3)  # tabulated xi_1 for n=4
14.972
physicskit.astro.leapfrog_step(positions, velocities, masses, dt, G=1.0, softening=0.0)[source]#

One kick-drift-kick leapfrog (symplectic) integration step.

Delegates to the shared physicskit.integrators.leapfrog_step().

Parameters:

  • positions (ndarray of shape (N, 3))

  • velocities (ndarray of shape (N, 3))

  • masses (ndarray of shape (N,))

  • dt (float) – Timestep.

  • G (float, default=1.0)

  • softening (float, default=0.0)

Returns:

new_positions, new_velocities (ndarray of shape (N, 3))

physicskit.astro.main_sequence_luminosity(mass_solar)[source]#

Empirical main-sequence mass-luminosity relation, \(L/L_\odot\approx(M/M_\odot)^{3.5}\).

Valid roughly for main-sequence stars in the 0.5-10 solar mass range; this is a pure empirical power law, not derived from first-principles stellar structure.

Parameters:

mass_solar (float) – Mass, in solar masses.

Returns:

float – Luminosity, in solar luminosities.

Examples

>>> round(main_sequence_luminosity(1.0), 6)
1.0
physicskit.astro.nfw_density(r, rho_s, r_s)[source]#

The Navarro-Frenk-White density profile.

\[\rho_{\rm NFW}(r) = \frac{\rho_s}{(r/r_s)(1+r/r_s)^2}\]
Parameters:

  • r (array_like) – Radius.

  • rho_s (float) – Characteristic density.

  • r_s (float) – Scale radius.

Returns:

ndarray or float

Examples

>>> round(float(nfw_density(1.0, rho_s=1.0, r_s=1.0)), 6)
0.25
physicskit.astro.nfw_enclosed_mass(r, rho_s, r_s, G=1.0)[source]#

Mass enclosed within radius r for an NFW halo.

\[M(<r) = 4\pi\rho_s r_s^3\left[\ln\!\left(1+\frac{r}{r_s}\right) - \frac{r/r_s}{1+r/r_s}\right]\]

G is accepted for interface consistency with the rest of this module but is not used by this formula.

Parameters:

  • r (array_like) – Radius.

  • rho_s (float) – Characteristic density.

  • r_s (float) – Scale radius.

  • G (float, default=1.0) – Unused; accepted for a consistent function signature.

Returns:

ndarray or float

Examples

>>> round(float(nfw_enclosed_mass(1.0, rho_s=1.0, r_s=1.0)), 6)
2.427159
physicskit.astro.nfw_potential(r, rho_s, r_s, G=1.0)[source]#

The NFW gravitational potential.

\[\Phi(r) = -\frac{4\pi G\rho_sr_s^3}{r}\ln\!\left(1+\frac{r}{r_s}\right)\]
Parameters:

  • r (array_like) – Radius.

  • rho_s (float) – Characteristic density.

  • r_s (float) – Scale radius.

  • G (float, default=1.0) – Gravitational constant.

Returns:

ndarray or float

physicskit.astro.orbital_elements_from_state(r_vec, v_vec, mu)[source]#

Classical (Keplerian) orbital elements from a Cartesian state vector.

Returns semi-major axis, eccentricity, inclination, right ascension of the ascending node, argument of periapsis, and true anomaly (angles in radians). For a near-equatorial orbit (\(i\approx0\) or \(\pi\)), the node vector is near zero and raan is undefined; this implementation returns raan=0.0 in that case rather than raising – a documented simplification, not a fully general equatorial-orbit treatment.

Parameters:

  • r_vec (ndarray of shape (3,)) – Position and velocity.

  • v_vec (ndarray of shape (3,)) – Position and velocity.

  • mu (float) – Gravitational parameter.

Returns:

a, e, i, raan, argp, nu (float)

Examples

>>> r_vec, v_vec = state_from_orbital_elements(1.0, 0.3, 0.5, 1.0, 0.7, 1.2, mu=1.0)
>>> a, e, i, raan, argp, nu = orbital_elements_from_state(r_vec, v_vec, mu=1.0)
>>> [round(x, 6) for x in (a, e, i, raan, argp, nu)]
[1.0, 0.3, 0.5, 1.0, 0.7, 1.2]
physicskit.astro.orbital_period(a, mu)[source]#

Kepler’s third law, \(T=2\pi\sqrt{a^3/\mu}\).

Parameters:

  • a (float) – Semi-major axis (bound orbit, a > 0).

  • mu (float) – Gravitational parameter, \(\mu=GM\).

Returns:

float

Examples

>>> round(orbital_period(1.0, 1.0) / (2 * np.pi), 6)
1.0
physicskit.astro.plot_lane_emden(xi, theta, ax=None)[source]#

Plot the Lane-Emden function \(\theta(\xi)\).

Parameters:
Returns:

  • fig (matplotlib.figure.Figure)

  • ax (matplotlib.axes.Axes)

Examples

>>> from physicskit.astro.stellar_structure import lane_emden
>>> xi, theta = lane_emden(1.5)
>>> fig, ax = plot_lane_emden(xi, theta)
>>> isinstance(fig, plt.Figure)
True
physicskit.astro.plot_nbody_trajectories(history, ax=None)[source]#

Plot the (x, y) trajectories of every body in an N-body simulation.

Parameters:
Returns:

  • fig (matplotlib.figure.Figure)

  • ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> from physicskit.astro.nbody import NBodySystem
>>> pos = np.array([[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0]])
>>> vel = np.array([[0.0, 0.5, 0.0], [0.0, -0.5, 0.0]])
>>> system = NBodySystem(pos, vel, np.array([1.0, 1.0]))
>>> history = system.simulate(dt=0.01, n_steps=50)
>>> fig, ax = plot_nbody_trajectories(history)
>>> isinstance(fig, plt.Figure)
True
physicskit.astro.plot_rotation_curve(r, v_model, ax=None, v_observed=None)[source]#

Plot a galactic rotation curve, optionally overlaid with observed data.

Parameters:

  • r (ndarray) – Radii.

  • v_model (ndarray) – Model circular velocity, e.g. from physicskit.astro.galactic_dynamics.circular_velocity().

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

  • v_observed (ndarray, optional) – Observed velocities to overlay as points.

Returns:

  • fig (matplotlib.figure.Figure)

  • ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> r = np.linspace(1, 20, 30)
>>> v = np.sqrt(1.0 / r)
>>> fig, ax = plot_rotation_curve(r, v)
>>> isinstance(fig, plt.Figure)
True
physicskit.astro.state_from_orbital_elements(a, e, i, raan, argp, nu, mu)[source]#

The inverse of orbital_elements_from_state().

Parameters:

  • a (float) – Classical orbital elements (angles in radians).

  • e (float) – Classical orbital elements (angles in radians).

  • i (float) – Classical orbital elements (angles in radians).

  • raan (float) – Classical orbital elements (angles in radians).

  • argp (float) – Classical orbital elements (angles in radians).

  • nu (float) – Classical orbital elements (angles in radians).

  • mu (float) – Gravitational parameter.

Returns:

r_vec, v_vec (ndarray of shape (3,))

physicskit.astro.vis_viva_speed(r, a, mu)[source]#

The vis-viva equation, \(v=\sqrt{\mu\left(2/r-1/a\right)}\).

Works for a > 0 (ellipse), a < 0 (hyperbola), and a = inf (parabolic escape trajectory, where the formula reduces to \(v=\sqrt{2\mu/r}\)).

Parameters:

  • r (float) – Current orbital radius.

  • a (float) – Semi-major axis.

  • mu (float) – Gravitational parameter.

Returns:

float

Examples

>>> round(vis_viva_speed(1.0, np.inf, 1.0), 6)
1.414214

Polytropic stellar models via the Lane-Emden equation.

Uses gravitational units with \(G=1\) by default – every function that needs it accepts G as a keyword argument, the same convention physicskit.relativity uses for its own geometrized units.

class physicskit.astro.stellar_structure.PolytropicStar(n, K, rho_c, G=1.0)[source]#

Bases: object

A physical polytropic star built from one Lane-Emden solution.

Parameters:

  • n (float) – Polytropic index, \(P=K\rho^{1+1/n}\).

  • K (float) – Polytropic constant.

  • rho_c (float) – Central density.

  • G (float, default=1.0) – Gravitational constant.

Raises:

ValueError – If the Lane-Emden solution has no surface (\(n \ge 5\), where the polytrope’s radius is infinite).

Examples

>>> star = PolytropicStar(n=0.0, K=1.0, rho_c=1.0)
>>> round(star.mass / (4 / 3 * np.pi * star.radius**3 * star.rho_c), 3)
1.0
property alpha#

Physical length scale, \(\alpha=\sqrt{(n+1)K\rho_c^{1/n-1}/(4\pi G)}\).

property mass#

Physical stellar mass, \(M=4\pi\alpha^3\rho_c\left[-\xi_1^2\theta'(\xi_1)\right]\).

property radius#

Physical stellar radius, \(R=\alpha\xi_1\).

property xi1#

The dimensionless surface radius \(\xi_1\), the first zero of the Lane-Emden solution.

physicskit.astro.stellar_structure.chandrasekhar_mass(mu_e=2.0)[source]#

The Chandrasekhar mass limit, in solar masses.

\[M_{\rm Ch} \approx \frac{5.83}{\mu_e^2}\ M_\odot\]

(the standard numerical prefactor from the n=3 polytrope / relativistic-electron-degeneracy treatment; see e.g. Kippenhahn, Weigert & Weiss, Stellar Structure and Evolution).

Parameters:

mu_e (float, default=2.0) – Mean molecular weight per electron (2.0 for carbon/oxygen white dwarfs).

Returns:

float – Mass limit, in solar masses.

Examples

>>> M = chandrasekhar_mass(2.0)
>>> 1.2 < M < 1.5
True
physicskit.astro.stellar_structure.lane_emden(n, xi_max=1000.0, n_points=2000)[source]#

Integrate the Lane-Emden equation for a polytrope of index n.

\[\frac{1}{\xi^2}\frac{d}{d\xi}\left(\xi^2\frac{d\theta}{d\xi}\right) + \theta^n = 0, \qquad \theta(0)=1,\ \theta'(0)=0\]

the dimensionless hydrostatic-equilibrium profile of a self-gravitating sphere with equation of state \(P=K\rho^{1+1/n}\).

Parameters:

  • n (float) – Polytropic index.

  • xi_max (float, default=1000.0) – Upper integration bound. Integration normally stops earlier, at the surface; the default is large enough to reach it for every \(n \le 4.9\) (\(\xi_1 \approx 171\) there, diverging as \(n \to 5\)). For \(n \ge 5\) there is no surface and the profile is returned out to xi_max.

  • n_points (int, default=2000) – Number of points to sample the returned arrays at.

Returns:

  • xi (ndarray) – Dimensionless radius, from just above 0 to the surface (the first zero of \(\theta\)) or, if none was found, to xi_max.

  • theta (ndarray) – The Lane-Emden function \(\theta(\xi)\). theta[-1] == 0 exactly when the surface was reached.

Examples

>>> xi, theta = lane_emden(0.0)
>>> round(float(xi[-1]), 3)  # exact analytic surface is sqrt(6)
2.449
>>> xi, theta = lane_emden(4.0)
>>> round(float(xi[-1]), 3)  # tabulated xi_1 for n=4
14.972
physicskit.astro.stellar_structure.main_sequence_luminosity(mass_solar)[source]#

Empirical main-sequence mass-luminosity relation, \(L/L_\odot\approx(M/M_\odot)^{3.5}\).

Valid roughly for main-sequence stars in the 0.5-10 solar mass range; this is a pure empirical power law, not derived from first-principles stellar structure.

Parameters:

mass_solar (float) – Mass, in solar masses.

Returns:

float – Luminosity, in solar luminosities.

Examples

>>> round(main_sequence_luminosity(1.0), 6)
1.0

Direct-summation N-body gravitational dynamics.

Uses gravitational units with \(G=1\) by default – every function accepts G as a keyword argument, the same convention physicskit.relativity uses for its own geometrized units.

class physicskit.astro.nbody.NBodySystem(positions, velocities, masses, G=1.0, softening=0.0)[source]#

Bases: object

A stateful direct-summation N-body gravitational simulation.

Parameters:

  • positions (ndarray of shape (N, 3))

  • velocities (ndarray of shape (N, 3))

  • masses (ndarray of shape (N,))

  • G (float, default=1.0)

  • softening (float, default=0.0)

simulate(dt, n_steps)[source]#

Advance n_steps steps, returning the position history.

Returns:

ndarray of shape (n_steps + 1, N, 3)

step(dt)[source]#

Advance the system’s stored state by one leapfrog step, in place.

total_angular_momentum()[source]#

Total angular momentum, \(\vec L=\sum_im_i\,\vec r_i\times\vec v_i\).

total_energy()[source]#

Kinetic plus (softened) gravitational potential energy, at the current state.

physicskit.astro.nbody.figure_eight_initial_conditions()[source]#

Initial conditions for the equal-mass planar figure-eight three-body choreography.

Three equal masses, in \(G=1\) units, chase each other forever around a single figure-eight-shaped curve – a periodic solution discovered numerically by Moore (1993) and proven to exist by Chenciner & Montgomery (2000, Annals of Mathematics 152). Unlike typical three-body configurations it is planar, collision-free, and exactly periodic, making it a clean, visually striking demo for NBodySystem.

Returns:

  • positions (ndarray of shape (3, 3))

  • velocities (ndarray of shape (3, 3))

  • masses (ndarray of shape (3,))

Examples

>>> positions, velocities, masses = figure_eight_initial_conditions()
>>> system = NBodySystem(positions, velocities, masses)
>>> history = system.simulate(dt=0.001, n_steps=100)
>>> history.shape
(101, 3, 3)
physicskit.astro.nbody.gravitational_acceleration(positions, masses, G=1.0, softening=0.0)[source]#

Pairwise Plummer-softened gravitational acceleration.

\[\vec a_i = G\sum_{j\neq i} m_j\, \frac{\vec r_j-\vec r_i}{\left(|\vec r_j-\vec r_i|^2+\epsilon^2\right)^{3/2}}\]
Parameters:

  • positions (ndarray of shape (N, 3))

  • masses (ndarray of shape (N,))

  • G (float, default=1.0)

  • softening (float, default=0.0) – Plummer softening length \(\epsilon\).

Returns:

ndarray of shape (N, 3)

Examples

>>> import numpy as np
>>> pos = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]])
>>> m = np.array([1.0, 1.0])
>>> acc = gravitational_acceleration(pos, m)
>>> round(float(acc[0, 0]), 6)
1.0
physicskit.astro.nbody.leapfrog_step(positions, velocities, masses, dt, G=1.0, softening=0.0)[source]#

One kick-drift-kick leapfrog (symplectic) integration step.

Delegates to the shared physicskit.integrators.leapfrog_step().

Parameters:

  • positions (ndarray of shape (N, 3))

  • velocities (ndarray of shape (N, 3))

  • masses (ndarray of shape (N,))

  • dt (float) – Timestep.

  • G (float, default=1.0)

  • softening (float, default=0.0)

Returns:

new_positions, new_velocities (ndarray of shape (N, 3))

Two-body Keplerian orbital mechanics: elements, state vectors, and transfers.

Uses gravitational units with \(G=1\) by default; the standard gravitational parameter \(\mu=GM\) appears directly as a function argument, the normal convention in orbital mechanics regardless of unit system (the same style physicskit.relativity uses for its own geometrized units).

physicskit.astro.orbital_mechanics.hohmann_transfer(r1, r2, mu)[source]#

The two burns and transfer time for a Hohmann transfer.

Parameters:

  • r1 (float) – Initial circular orbit radius.

  • r2 (float) – Final circular orbit radius.

  • mu (float) – Gravitational parameter.

Returns:

  • delta_v1 (float) – Speed change (magnitude) to leave the initial circular orbit and enter the transfer ellipse.

  • delta_v2 (float) – Speed change (magnitude) to circularize at r2.

  • transfer_time (float) – Time to complete the transfer, half the transfer ellipse’s period.

Examples

>>> dv1, dv2, t = hohmann_transfer(1.0, 4.0, mu=1.0)
>>> dv1 > 0 and dv2 > 0 and t > 0
True
physicskit.astro.orbital_mechanics.orbital_elements_from_state(r_vec, v_vec, mu)[source]#

Classical (Keplerian) orbital elements from a Cartesian state vector.

Returns semi-major axis, eccentricity, inclination, right ascension of the ascending node, argument of periapsis, and true anomaly (angles in radians). For a near-equatorial orbit (\(i\approx0\) or \(\pi\)), the node vector is near zero and raan is undefined; this implementation returns raan=0.0 in that case rather than raising – a documented simplification, not a fully general equatorial-orbit treatment.

Parameters:

  • r_vec (ndarray of shape (3,)) – Position and velocity.

  • v_vec (ndarray of shape (3,)) – Position and velocity.

  • mu (float) – Gravitational parameter.

Returns:

a, e, i, raan, argp, nu (float)

Examples

>>> r_vec, v_vec = state_from_orbital_elements(1.0, 0.3, 0.5, 1.0, 0.7, 1.2, mu=1.0)
>>> a, e, i, raan, argp, nu = orbital_elements_from_state(r_vec, v_vec, mu=1.0)
>>> [round(x, 6) for x in (a, e, i, raan, argp, nu)]
[1.0, 0.3, 0.5, 1.0, 0.7, 1.2]
physicskit.astro.orbital_mechanics.orbital_period(a, mu)[source]#

Kepler’s third law, \(T=2\pi\sqrt{a^3/\mu}\).

Parameters:

  • a (float) – Semi-major axis (bound orbit, a > 0).

  • mu (float) – Gravitational parameter, \(\mu=GM\).

Returns:

float

Examples

>>> round(orbital_period(1.0, 1.0) / (2 * np.pi), 6)
1.0
physicskit.astro.orbital_mechanics.state_from_orbital_elements(a, e, i, raan, argp, nu, mu)[source]#

The inverse of orbital_elements_from_state().

Parameters:

  • a (float) – Classical orbital elements (angles in radians).

  • e (float) – Classical orbital elements (angles in radians).

  • i (float) – Classical orbital elements (angles in radians).

  • raan (float) – Classical orbital elements (angles in radians).

  • argp (float) – Classical orbital elements (angles in radians).

  • nu (float) – Classical orbital elements (angles in radians).

  • mu (float) – Gravitational parameter.

Returns:

r_vec, v_vec (ndarray of shape (3,))

physicskit.astro.orbital_mechanics.vis_viva_speed(r, a, mu)[source]#

The vis-viva equation, \(v=\sqrt{\mu\left(2/r-1/a\right)}\).

Works for a > 0 (ellipse), a < 0 (hyperbola), and a = inf (parabolic escape trajectory, where the formula reduces to \(v=\sqrt{2\mu/r}\)).

Parameters:

  • r (float) – Current orbital radius.

  • a (float) – Semi-major axis.

  • mu (float) – Gravitational parameter.

Returns:

float

Examples

>>> round(vis_viva_speed(1.0, np.inf, 1.0), 6)
1.414214

The Navarro-Frenk-White dark matter halo profile and rotation curves.

Uses gravitational units with \(G=1\) by default – every function accepts G as a keyword argument, the same convention physicskit.relativity uses for its own geometrized units.

physicskit.astro.galactic_dynamics.circular_velocity(r, mass_enclosed_func, G=1.0)[source]#

Circular orbital speed for any enclosed-mass profile, \(v_c(r)=\sqrt{GM(<r)/r}\).

Parameters:

  • r (array_like) – Radius.

  • mass_enclosed_func (callable) – mass_enclosed_func(r) -> mass, e.g. lambda r: nfw_enclosed_mass(r, rho_s, r_s).

  • G (float, default=1.0) – Gravitational constant.

Returns:

ndarray or float

Examples

>>> M0 = 5.0
>>> round(float(circular_velocity(2.0, lambda r: M0, G=1.0)), 6)
1.581139
physicskit.astro.galactic_dynamics.nfw_density(r, rho_s, r_s)[source]#

The Navarro-Frenk-White density profile.

\[\rho_{\rm NFW}(r) = \frac{\rho_s}{(r/r_s)(1+r/r_s)^2}\]
Parameters:

  • r (array_like) – Radius.

  • rho_s (float) – Characteristic density.

  • r_s (float) – Scale radius.

Returns:

ndarray or float

Examples

>>> round(float(nfw_density(1.0, rho_s=1.0, r_s=1.0)), 6)
0.25
physicskit.astro.galactic_dynamics.nfw_enclosed_mass(r, rho_s, r_s, G=1.0)[source]#

Mass enclosed within radius r for an NFW halo.

\[M(<r) = 4\pi\rho_s r_s^3\left[\ln\!\left(1+\frac{r}{r_s}\right) - \frac{r/r_s}{1+r/r_s}\right]\]

G is accepted for interface consistency with the rest of this module but is not used by this formula.

Parameters:

  • r (array_like) – Radius.

  • rho_s (float) – Characteristic density.

  • r_s (float) – Scale radius.

  • G (float, default=1.0) – Unused; accepted for a consistent function signature.

Returns:

ndarray or float

Examples

>>> round(float(nfw_enclosed_mass(1.0, rho_s=1.0, r_s=1.0)), 6)
2.427159
physicskit.astro.galactic_dynamics.nfw_potential(r, rho_s, r_s, G=1.0)[source]#

The NFW gravitational potential.

\[\Phi(r) = -\frac{4\pi G\rho_sr_s^3}{r}\ln\!\left(1+\frac{r}{r_s}\right)\]
Parameters:

  • r (array_like) – Radius.

  • rho_s (float) – Characteristic density.

  • r_s (float) – Scale radius.

  • G (float, default=1.0) – Gravitational constant.

Returns:

ndarray or float

Plotting helpers: N-body trajectories, Lane-Emden profiles, and rotation curves.

Matplotlib is used throughout. Every function returns its figure object rather than calling show(), so it composes cleanly into larger figures or headless pipelines.

physicskit.astro.visualizers.animate_dynamo_wave(times, A_snapshots, B_snapshots, x, interval=80)[source]#

Animate the alpha-omega dynamo’s poloidal and toroidal fields as traveling, growing waves.

Draws \(A(x)\) and \(B(x)\) as two line plots sharing an axis, redrawn frame by frame from physicskit.astro.stellar_dynamo.simulate_alpha_omega_dynamo()’s output. In the unstable/growing regime, both the growing amplitude and the lateral migration of the wave pattern across x – the dynamo wave whose space-time trace is the solar butterfly diagram – are visible as the animation runs; the y-axis is rescaled each frame to the current amplitude so the growing wave doesn’t run off scale.

Parameters:

  • times (ndarray, shape (n_saved,))

  • A_snapshots (ndarray, shape (n_saved, nx)) – Poloidal-flux-proxy and toroidal-field snapshots.

  • B_snapshots (ndarray, shape (n_saved, nx)) – Poloidal-flux-proxy and toroidal-field snapshots.

  • x (ndarray, shape (nx,)) – Spatial grid the snapshots live on.

  • interval (int, default=80) – Delay between animation frames, in milliseconds.

Returns:

matplotlib.animation.FuncAnimation – Assign it to a variable to keep it alive, and display it with plt.show() or save it with anim.save(...).

See also

plot_dynamo_butterfly_diagram

The static space-time summary of the same wave.

Examples

>>> import numpy as np
>>> from physicskit.astro.stellar_dynamo import simulate_alpha_omega_dynamo
>>> Lx = 2 * np.pi
>>> x = np.linspace(0, Lx, 32, endpoint=False)
>>> A0 = 1e-3 * np.cos(x)
>>> B0 = np.zeros_like(A0)
>>> times, A_snaps, B_snaps = simulate_alpha_omega_dynamo(A0, B0, alpha=1.0, shear=5.0, eta=0.05, Lx=Lx, dt=0.01, n_steps=20, save_every=5)
>>> anim = animate_dynamo_wave(times, A_snaps, B_snaps, x)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.astro.visualizers.animate_nbody_trajectories(history, interval=30, trail=200, skip=1, colors=None, title=None)[source]#

Animate an N-body simulation’s real-space trajectories being traced out.

Draws each body’s (x, y) path (the z component, if any, is ignored) with a fading trail – older trail segments are more transparent – and a bright marker at the body’s current position, in a distinct color per body.

Parameters:

  • history (ndarray of shape (n_steps + 1, N, 3)) – Position history, e.g. from physicskit.astro.nbody.NBodySystem.simulate().

  • interval (int) – Delay between animation frames, in milliseconds.

  • trail (int) – Number of most-recent steps kept visible as each body’s trail.

  • skip (int) – Number of simulation steps advanced per animation frame; increase for long simulations so the animation covers them in a reasonable number of frames.

  • colors (list, optional) – Per-body colors; a tab10 sweep is used if omitted.

  • title (str | None) – Axes title.

Returns:

matplotlib.animation.FuncAnimation – Assign it to a variable to keep it alive, and display it with plt.show() or save it with anim.save(...).

Examples

>>> from physicskit.astro.nbody import NBodySystem, figure_eight_initial_conditions
>>> pos, vel, masses = figure_eight_initial_conditions()
>>> system = NBodySystem(pos, vel, masses)
>>> history = system.simulate(dt=0.002, n_steps=200)
>>> anim = animate_nbody_trajectories(history)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.astro.visualizers.animate_stellar_convection(times, omega_snapshots, T_snapshots, interval=80)[source]#

Animate 2D convective rolls: vorticity and temperature side by side.

Two panels sharing a time axis – vorticity on the left, temperature perturbation on the right – each redrawn frame by frame from physicskit.astro.stellar_dynamo.simulate_stellar_convection()’s output, showing the small initial noise roll up into churning convective plumes.

Parameters:

  • times (ndarray, shape (n_saved,))

  • omega_snapshots (ndarray, shape (n_saved, ny, nx)) – Vorticity field snapshots.

  • T_snapshots (ndarray, shape (n_saved, ny, nx)) – Temperature-perturbation field snapshots.

  • interval (int, default=80) – Delay between animation frames, in milliseconds.

Returns:

matplotlib.animation.FuncAnimation – Assign it to a variable to keep it alive, and display it with plt.show() or save it with anim.save(...).

Examples

>>> import numpy as np
>>> from physicskit.astro.stellar_dynamo import simulate_stellar_convection
>>> times, omega_snaps, T_snaps = simulate_stellar_convection(
...     24, 24, 2 * np.pi, 2 * np.pi, n_steps=20, save_every=10, seed=0)
>>> anim = animate_stellar_convection(times, omega_snaps, T_snaps)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.astro.visualizers.animate_zeldovich_collapse(q, D_values, k_vectors, amplitudes, phases, interval=80, s=2)[source]#

Animate the Zel’dovich approximation’s collapse of a Lagrangian grid into the cosmic web.

Shows a scatter plot of physicskit.astro.cosmic_web.zeldovich_position() evaluated at each growth factor in D_values, frame by frame – the initially regular grid of tracer particles progressively streaming together into filaments, sheets (“pancakes”), and dense nodes at filament intersections. Points are colored by a local point-density proxy (see _local_density_on_grid()): a fixed 2D histogram of the particle positions, shared across all frames, so density colors are directly comparable as structure grows in from D_values[0] to D_values[-1].

Parameters:

  • q (ndarray, shape (n_particles, 2)) – Lagrangian tracer positions, e.g. from physicskit.astro.cosmic_web.lagrangian_grid().

  • D_values (array_like) – Growth-factor values to animate through, one frame each.

  • k_vectors (ndarray, shape (n_modes, 2))

  • amplitudes (ndarray, shape (n_modes,))

  • phases (ndarray, shape (n_modes,)) – Displacement-potential parameters, e.g. from physicskit.astro.cosmic_web.random_displacement_potential().

  • interval (int, default=80) – Delay between animation frames, in milliseconds.

  • s (float, default=2) – Marker size.

Returns:

matplotlib.animation.FuncAnimation – Assign it to a variable to keep it alive, and display it with plt.show() or save it with anim.save(...).

See also

plot_zeldovich_snapshot

The single-frame, non-animated companion.

Examples

>>> import numpy as np
>>> from physicskit.astro.cosmic_web import lagrangian_grid, random_displacement_potential, first_caustic_time
>>> q = lagrangian_grid(20, 1.0)
>>> k_vectors, amplitudes, phases = random_displacement_potential(6, k_min=2 * np.pi, k_max=6 * np.pi, amplitude_scale=0.02, seed=0)
>>> D_collapse = first_caustic_time(q, k_vectors, amplitudes, phases)
>>> D_values = np.linspace(0.0, 1.2 * D_collapse, 10)
>>> anim = animate_zeldovich_collapse(q, D_values, k_vectors, amplitudes, phases)
>>> anim.__class__.__name__
'FuncAnimation'
physicskit.astro.visualizers.plot_dynamo_butterfly_diagram(times, B_snapshots, x, ax=None)[source]#

Static space-time (Hovmoller) plot of the alpha-omega dynamo wave: the “butterfly diagram”.

Displays \(B(x,t)\) as an image with time along the x-axis and the spatial coordinate x (standing in for stellar latitude) along the y-axis. In the unstable regime, the migrating dynamo wave shows up as diagonal stripes of alternating sign drifting across the plot – the same characteristic pattern, and the same underlying alpha-omega mechanism, behind the real Sun’s sunspot-latitude butterfly diagram.

Parameters:

  • times (ndarray, shape (n_saved,))

  • B_snapshots (ndarray, shape (n_saved, nx)) – Toroidal-field snapshots, e.g. from physicskit.astro.stellar_dynamo.simulate_alpha_omega_dynamo().

  • x (ndarray, shape (nx,)) – Spatial grid the snapshots live on.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

  • fig (matplotlib.figure.Figure)

  • ax (matplotlib.axes.Axes)

See also

animate_dynamo_wave

The animated, time-domain view of the same field.

Examples

>>> import numpy as np
>>> from physicskit.astro.stellar_dynamo import simulate_alpha_omega_dynamo
>>> Lx = 2 * np.pi
>>> x = np.linspace(0, Lx, 32, endpoint=False)
>>> A0 = 1e-3 * np.cos(x)
>>> B0 = np.zeros_like(A0)
>>> times, A_snaps, B_snaps = simulate_alpha_omega_dynamo(A0, B0, alpha=1.0, shear=5.0, eta=0.05, Lx=Lx, dt=0.01, n_steps=200, save_every=5)
>>> fig, ax = plot_dynamo_butterfly_diagram(times, B_snaps, x)
>>> isinstance(fig, plt.Figure)
True
physicskit.astro.visualizers.plot_lane_emden(xi, theta, ax=None)[source]#

Plot the Lane-Emden function \(\theta(\xi)\).

Parameters:
Returns:

  • fig (matplotlib.figure.Figure)

  • ax (matplotlib.axes.Axes)

Examples

>>> from physicskit.astro.stellar_structure import lane_emden
>>> xi, theta = lane_emden(1.5)
>>> fig, ax = plot_lane_emden(xi, theta)
>>> isinstance(fig, plt.Figure)
True
physicskit.astro.visualizers.plot_nbody_trajectories(history, ax=None)[source]#

Plot the (x, y) trajectories of every body in an N-body simulation.

Parameters:
Returns:

  • fig (matplotlib.figure.Figure)

  • ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> from physicskit.astro.nbody import NBodySystem
>>> pos = np.array([[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0]])
>>> vel = np.array([[0.0, 0.5, 0.0], [0.0, -0.5, 0.0]])
>>> system = NBodySystem(pos, vel, np.array([1.0, 1.0]))
>>> history = system.simulate(dt=0.01, n_steps=50)
>>> fig, ax = plot_nbody_trajectories(history)
>>> isinstance(fig, plt.Figure)
True
physicskit.astro.visualizers.plot_rotation_curve(r, v_model, ax=None, v_observed=None)[source]#

Plot a galactic rotation curve, optionally overlaid with observed data.

Parameters:

  • r (ndarray) – Radii.

  • v_model (ndarray) – Model circular velocity, e.g. from physicskit.astro.galactic_dynamics.circular_velocity().

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

  • v_observed (ndarray, optional) – Observed velocities to overlay as points.

Returns:

  • fig (matplotlib.figure.Figure)

  • ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> r = np.linspace(1, 20, 30)
>>> v = np.sqrt(1.0 / r)
>>> fig, ax = plot_rotation_curve(r, v)
>>> isinstance(fig, plt.Figure)
True
physicskit.astro.visualizers.plot_zeldovich_snapshot(q, D, k_vectors, amplitudes, phases, ax=None, s=2)[source]#

A single static snapshot of the Zel’dovich-evolved particle field at one growth factor.

The non-animated companion to animate_zeldovich_collapse(), using the same 2D-histogram local-density coloring.

Parameters:

  • q (ndarray, shape (n_particles, 2)) – Lagrangian tracer positions, e.g. from physicskit.astro.cosmic_web.lagrangian_grid().

  • D (float) – Growth factor at which to evaluate the Zel’dovich mapping.

  • k_vectors (ndarray, shape (n_modes, 2))

  • amplitudes (ndarray, shape (n_modes,))

  • phases (ndarray, shape (n_modes,)) – Displacement-potential parameters, e.g. from physicskit.astro.cosmic_web.random_displacement_potential().

  • ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

  • s (float, default=2) – Marker size.

Returns:

  • fig (matplotlib.figure.Figure)

  • ax (matplotlib.axes.Axes)

See also

animate_zeldovich_collapse

The animated version of this snapshot.

Examples

>>> import numpy as np
>>> from physicskit.astro.cosmic_web import lagrangian_grid, random_displacement_potential, first_caustic_time
>>> q = lagrangian_grid(20, 1.0)
>>> k_vectors, amplitudes, phases = random_displacement_potential(6, k_min=2 * np.pi, k_max=6 * np.pi, amplitude_scale=0.02, seed=0)
>>> D_collapse = first_caustic_time(q, k_vectors, amplitudes, phases)
>>> fig, ax = plot_zeldovich_snapshot(q, D_collapse, k_vectors, amplitudes, phases)
>>> isinstance(fig, plt.Figure)
True