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.GreyAtmosphere(n_streams=32)[source]#

Bases: object

Exact grey atmosphere in Chandrasekhar’s \(n\)-stream discrete-ordinates approximation.

With Gauss-Legendre nodes \(\pm\mu_i\) and weights \(a_i\) (\(i = 1..n\)), the intensities are

\[I(\tau, \mu_i) = \tfrac34 F\left[\tau + \mu_i + Q + \sum_{\alpha=1}^{n-1} \frac{L_\alpha e^{-k_\alpha\tau}}{1 + k_\alpha\mu_i}\right],\]

where the \(k_\alpha\) are the positive roots of \(\sum_i a_i/(1 - k^2\mu_i^2) = 1\) and \(Q\), \(L_\alpha\) follow from no incoming radiation at \(\tau = 0\). The Hopf function is \(q(\tau) = Q + \sum_\alpha L_\alpha e^{-k_\alpha\tau}\).

Parameters:

n_streams (int) – Number \(n\) of directions per hemisphere.

Variables:

  • Q (float) – \(q(\infty)\), which tends to Hopf’s 0.710446 as \(n\) grows.

  • L (k,) – Decay constants and amplitudes of the transient terms.

Examples

>>> atm = GreyAtmosphere()
>>> round(float(atm.hopf(0.0)), 6)  # 1/sqrt(3) in every approximation
0.57735
>>> round(atm.Q, 4)
0.7104
emergent_intensity(mu)[source]#

Emergent intensity in units of \(F/\pi = \sigma T_{\text{eff}}^4/\pi\).

The Laplace transform of \(S = \tfrac34 F(\tau + q)/\pi\) gives \(I(0, \mu) = \tfrac34\left[\mu + Q + \sum_\alpha L_\alpha/(1 + k_\alpha\mu)\right]\).

Parameters:

mu (ArrayLike)

Return type:

NDArray[double]

Returns:

ndarray

hopf(tau)[source]#

Hopf function \(q(\tau)\).

Parameters:

tau (ArrayLike)

Return type:

NDArray[double]

Returns:

ndarray

limb_darkening(mu)[source]#

Normalized limb darkening \(I(0, \mu)/I(0, 1)\).

Return type:

NDArray[double]

Parameters:

mu (ArrayLike)

temperature(tau, T_eff=1.0)[source]#

\(T(\tau) = T_{\text{eff}}\,[\tfrac34(\tau + q(\tau))]^{1/4}\).

Return type:

NDArray[double]

Parameters:

  • tau (ArrayLike)

  • T_eff (float)

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.

class physicskit.astro.SPH1D(x, v, m, u, gamma=1.4, eta=1.2, alpha=1.0, beta=2.0, fixed=None)[source]#

Bases: object

One-dimensional SPH gas with the cubic-spline kernel and Monaghan artificial viscosity.

Smoothing lengths follow the density, \(h_i = \eta\, m_i/\rho_i\), and the equations are advanced with a kick-drift-kick leapfrog at a Courant-limited time step.

Parameters:

  • x (ndarray of float, shape (n,)) – Particle positions and velocities.

  • v (ndarray of float, shape (n,)) – Particle positions and velocities.

  • m (ndarray of float or float) – Particle masses.

  • u (ndarray of float, shape (n,)) – Specific internal energies.

  • gamma (float, default 1.4) – Adiabatic index.

  • eta (float, default 1.2) – Smoothing-length factor, in units of the local particle spacing.

  • alpha (float, default 1.0, 2.0) – Artificial-viscosity coefficients.

  • beta (float, default 1.0, 2.0) – Artificial-viscosity coefficients.

  • fixed (ndarray of bool, optional) – Particles held at rest (boundaries).

Examples

>>> gas = SPH1D.sod_shock_tube(n_left=200)
>>> gas.run(0.05)
>>> bool(abs(gas.total_energy() - gas.energy_history[0]) < 1e-3 * gas.energy_history[0])
True
property pressure: NDArray[float64]#

Pressure \(P = (\gamma - 1)\rho u\).

run(t_end, cfl=0.3)[source]#

Advance to self.t + t_end with kick-drift-kick leapfrog.

Parameters:

  • t_end (float) – Time to advance by.

  • cfl (float) – Courant factor in \(\Delta t = C\min_i h_i/(c_i + |v_i|)\).

Return type:

None

classmethod sod_shock_tube(n_left=400, gamma=1.4, left=(1.0, 1.0), right=(0.125, 0.1), length=1.0, **kwargs)[source]#

Equal-mass particles for Sod’s shock tube, with the diaphragm at \(x = 0\).

The left half \(-L < x < 0\) has density and pressure left, the right half right; both are at rest. Particle spacing scales inversely with density, and the outermost few particles on each side are held fixed as walls.

Parameters:

  • n_left (int, default 400) – Particles in the left half.

  • gamma (float, default 1.4)

  • left (tuple of float) – (rho, P) of each side.

  • right (tuple of float) – (rho, P) of each side.

  • length (float, default 1.0) – Half-length \(L\) of the tube.

  • **kwargs – Passed to the constructor.

Returns:

SPH1D

total_energy()[source]#

Total kinetic plus internal energy \(\sum_i m_i(v_i^2/2 + u_i)\).

Return type:

float

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.cubic_spline_kernel(r, h)[source]#

The 1D cubic-spline (M4) kernel \(W(r, h)\), normalized to \(\int W\,dr = 1\).

\[\begin{split}W = \frac{2}{3h}\begin{cases} 1 - \tfrac32 q^2 + \tfrac34 q^3, & 0 \le q < 1, \\ \tfrac14(2 - q)^3, & 1 \le q < 2, \\ 0, & q \ge 2, \end{cases} \qquad q = |r|/h.\end{split}\]
Parameters:
Return type:

NDArray[double]

Returns:

ndarray of float

Examples

>>> r = np.linspace(-3, 3, 60001)
>>> round(float(np.trapezoid(cubic_spline_kernel(r, 1.0), r)), 6)
1.0
physicskit.astro.eddington_limb_darkening(mu)[source]#

Eddington limb darkening \(I(0,\mu)/I(0,1) = (2 + 3\mu)/5\).

Parameters:

mu (ArrayLike) – Cosine of the angle from the disc normal (1 at disc centre, 0 at the limb).

Return type:

NDArray[double]

Returns:

ndarray

physicskit.astro.eddington_temperature(tau, T_eff=1.0)[source]#

Eddington’s grey temperature law \(T^4 = \tfrac34 T_{\text{eff}}^4(\tau + 2/3)\).

Parameters:
Return type:

NDArray[double]

Returns:

ndarray

Examples

>>> float(eddington_temperature(2 / 3, T_eff=5772.0))  # T = T_eff at tau = 2/3
5772.0
physicskit.astro.emergent_intensity(source, mu, n=64)[source]#

Formal solution \(I(0, \mu) = \int_0^\infty S(\tau)e^{-\tau/\mu}\,d\tau/\mu\).

Evaluated by \(n\)-point Gauss-Laguerre quadrature in \(t = \tau/\mu\).

Parameters:

  • source (callable) – Vectorized source function \(S(\tau)\).

  • mu (ArrayLike) – Direction cosines in \((0, 1]\).

  • n (int) – Quadrature points.

Return type:

NDArray[double]

Returns:

ndarray

Examples

A linear source function \(S = a + b\tau\) emerges as \(a + b\mu\) (the Eddington-Barbier relation):

>>> round(float(emergent_intensity(lambda t: 1 + 2 * t, 0.5)), 10)
2.0
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

Smoothed-particle hydrodynamics (SPH) in one dimension.

SPH represents a gas by particles of fixed mass \(m\) that move with the flow (Lucy, Astron. J. 82, 1013 (1977); Gingold and Monaghan, MNRAS 181, 375 (1977)). Every field is a kernel-weighted sum over neighbours; the density is

\[\rho_i = \sum_j m_j W(x_i - x_j, h_i),\]

and the Euler equations become ordinary differential equations for the particles (Monaghan, Annu. Rev. Astron. Astrophys. 30, 543 (1992)),

\[\frac{dv_i}{dt} = -\sum_j m_j\left(\frac{P_i}{\rho_i^2} + \frac{P_j}{\rho_j^2} + \Pi_{ij}\right)\partial_i W_{ij}, \qquad \frac{du_i}{dt} = \frac12\sum_j m_j\left(\frac{P_i}{\rho_i^2} + \frac{P_j}{\rho_j^2} + \Pi_{ij}\right) v_{ij}\,\partial_i W_{ij},\]

with \(P = (\gamma - 1)\rho u\), the cubic-spline kernel, and Monaghan’s artificial viscosity \(\Pi_{ij}\) to capture shocks. Because the particles follow the mass, resolution concentrates where the gas is dense, which is why SPH became the standard method for star formation, galaxy mergers and cosmological gas. SPH1D evolves these equations, and SPH1D.sod_shock_tube() sets up Sod’s test problem.

class physicskit.astro.sph.SPH1D(x, v, m, u, gamma=1.4, eta=1.2, alpha=1.0, beta=2.0, fixed=None)[source]#

Bases: object

One-dimensional SPH gas with the cubic-spline kernel and Monaghan artificial viscosity.

Smoothing lengths follow the density, \(h_i = \eta\, m_i/\rho_i\), and the equations are advanced with a kick-drift-kick leapfrog at a Courant-limited time step.

Parameters:

  • x (ndarray of float, shape (n,)) – Particle positions and velocities.

  • v (ndarray of float, shape (n,)) – Particle positions and velocities.

  • m (ndarray of float or float) – Particle masses.

  • u (ndarray of float, shape (n,)) – Specific internal energies.

  • gamma (float, default 1.4) – Adiabatic index.

  • eta (float, default 1.2) – Smoothing-length factor, in units of the local particle spacing.

  • alpha (float, default 1.0, 2.0) – Artificial-viscosity coefficients.

  • beta (float, default 1.0, 2.0) – Artificial-viscosity coefficients.

  • fixed (ndarray of bool, optional) – Particles held at rest (boundaries).

Examples

>>> gas = SPH1D.sod_shock_tube(n_left=200)
>>> gas.run(0.05)
>>> bool(abs(gas.total_energy() - gas.energy_history[0]) < 1e-3 * gas.energy_history[0])
True
property pressure: NDArray[float64]#

Pressure \(P = (\gamma - 1)\rho u\).

run(t_end, cfl=0.3)[source]#

Advance to self.t + t_end with kick-drift-kick leapfrog.

Parameters:

  • t_end (float) – Time to advance by.

  • cfl (float) – Courant factor in \(\Delta t = C\min_i h_i/(c_i + |v_i|)\).

Return type:

None

classmethod sod_shock_tube(n_left=400, gamma=1.4, left=(1.0, 1.0), right=(0.125, 0.1), length=1.0, **kwargs)[source]#

Equal-mass particles for Sod’s shock tube, with the diaphragm at \(x = 0\).

The left half \(-L < x < 0\) has density and pressure left, the right half right; both are at rest. Particle spacing scales inversely with density, and the outermost few particles on each side are held fixed as walls.

Parameters:

  • n_left (int, default 400) – Particles in the left half.

  • gamma (float, default 1.4)

  • left (tuple of float) – (rho, P) of each side.

  • right (tuple of float) – (rho, P) of each side.

  • length (float, default 1.0) – Half-length \(L\) of the tube.

  • **kwargs – Passed to the constructor.

Returns:

SPH1D

total_energy()[source]#

Total kinetic plus internal energy \(\sum_i m_i(v_i^2/2 + u_i)\).

Return type:

float

physicskit.astro.sph.cubic_spline_kernel(r, h)[source]#

The 1D cubic-spline (M4) kernel \(W(r, h)\), normalized to \(\int W\,dr = 1\).

\[\begin{split}W = \frac{2}{3h}\begin{cases} 1 - \tfrac32 q^2 + \tfrac34 q^3, & 0 \le q < 1, \\ \tfrac14(2 - q)^3, & 1 \le q < 2, \\ 0, & q \ge 2, \end{cases} \qquad q = |r|/h.\end{split}\]
Parameters:
Return type:

NDArray[double]

Returns:

ndarray of float

Examples

>>> r = np.linspace(-3, 3, 60001)
>>> round(float(np.trapezoid(cubic_spline_kernel(r, 1.0), r)), 6)
1.0

Grey radiative transfer in a plane-parallel stellar atmosphere.

Radiation crossing a stellar atmosphere obeys, at optical depth \(\tau\) and direction cosine \(\mu\),

\[\mu\frac{dI}{d\tau} = I - S.\]

For a grey opacity (independent of frequency) in radiative equilibrium the source function equals the mean intensity, \(S = J\), and the flux \(F = \sigma T_{\text{eff}}^4\) is constant with depth. The solution is (Milne 1921; Hopf 1930)

\[T^4(\tau) = \tfrac34 T_{\text{eff}}^4\,[\tau + q(\tau)],\]

where the Hopf function \(q\) rises from \(1/\sqrt3\) at the surface to \(0.7104\) deep inside. Eddington’s (1926) approximation, closing the moment equations with \(K = J/3\), gives the constant \(q = 2/3\), and with it the limb darkening \(I(0, \mu)/I(0, 1) = (2 + 3\mu)/5\).

GreyAtmosphere solves the problem exactly in Chandrasekhar’s discrete-ordinates method (Astrophys. J. 100, 76 (1944); Radiative Transfer, 1950, ch. III), which converges to the exact solution as the number of streams grows and keeps \(q(0) = 1/\sqrt3\) exactly in every order.

class physicskit.astro.radiative_transfer.GreyAtmosphere(n_streams=32)[source]#

Bases: object

Exact grey atmosphere in Chandrasekhar’s \(n\)-stream discrete-ordinates approximation.

With Gauss-Legendre nodes \(\pm\mu_i\) and weights \(a_i\) (\(i = 1..n\)), the intensities are

\[I(\tau, \mu_i) = \tfrac34 F\left[\tau + \mu_i + Q + \sum_{\alpha=1}^{n-1} \frac{L_\alpha e^{-k_\alpha\tau}}{1 + k_\alpha\mu_i}\right],\]

where the \(k_\alpha\) are the positive roots of \(\sum_i a_i/(1 - k^2\mu_i^2) = 1\) and \(Q\), \(L_\alpha\) follow from no incoming radiation at \(\tau = 0\). The Hopf function is \(q(\tau) = Q + \sum_\alpha L_\alpha e^{-k_\alpha\tau}\).

Parameters:

n_streams (int) – Number \(n\) of directions per hemisphere.

Variables:

  • Q (float) – \(q(\infty)\), which tends to Hopf’s 0.710446 as \(n\) grows.

  • L (k,) – Decay constants and amplitudes of the transient terms.

Examples

>>> atm = GreyAtmosphere()
>>> round(float(atm.hopf(0.0)), 6)  # 1/sqrt(3) in every approximation
0.57735
>>> round(atm.Q, 4)
0.7104
emergent_intensity(mu)[source]#

Emergent intensity in units of \(F/\pi = \sigma T_{\text{eff}}^4/\pi\).

The Laplace transform of \(S = \tfrac34 F(\tau + q)/\pi\) gives \(I(0, \mu) = \tfrac34\left[\mu + Q + \sum_\alpha L_\alpha/(1 + k_\alpha\mu)\right]\).

Parameters:

mu (ArrayLike)

Return type:

NDArray[double]

Returns:

ndarray

hopf(tau)[source]#

Hopf function \(q(\tau)\).

Parameters:

tau (ArrayLike)

Return type:

NDArray[double]

Returns:

ndarray

limb_darkening(mu)[source]#

Normalized limb darkening \(I(0, \mu)/I(0, 1)\).

Return type:

NDArray[double]

Parameters:

mu (ArrayLike)

temperature(tau, T_eff=1.0)[source]#

\(T(\tau) = T_{\text{eff}}\,[\tfrac34(\tau + q(\tau))]^{1/4}\).

Return type:

NDArray[double]

Parameters:

  • tau (ArrayLike)

  • T_eff (float)

physicskit.astro.radiative_transfer.eddington_limb_darkening(mu)[source]#

Eddington limb darkening \(I(0,\mu)/I(0,1) = (2 + 3\mu)/5\).

Parameters:

mu (ArrayLike) – Cosine of the angle from the disc normal (1 at disc centre, 0 at the limb).

Return type:

NDArray[double]

Returns:

ndarray

physicskit.astro.radiative_transfer.eddington_temperature(tau, T_eff=1.0)[source]#

Eddington’s grey temperature law \(T^4 = \tfrac34 T_{\text{eff}}^4(\tau + 2/3)\).

Parameters:
Return type:

NDArray[double]

Returns:

ndarray

Examples

>>> float(eddington_temperature(2 / 3, T_eff=5772.0))  # T = T_eff at tau = 2/3
5772.0
physicskit.astro.radiative_transfer.emergent_intensity(source, mu, n=64)[source]#

Formal solution \(I(0, \mu) = \int_0^\infty S(\tau)e^{-\tau/\mu}\,d\tau/\mu\).

Evaluated by \(n\)-point Gauss-Laguerre quadrature in \(t = \tau/\mu\).

Parameters:

  • source (callable) – Vectorized source function \(S(\tau)\).

  • mu (ArrayLike) – Direction cosines in \((0, 1]\).

  • n (int) – Quadrature points.

Return type:

NDArray[double]

Returns:

ndarray

Examples

A linear source function \(S = a + b\tau\) emerges as \(a + b\mu\) (the Eddington-Barbier relation):

>>> round(float(emergent_intensity(lambda t: 1 + 2 * t, 0.5)), 10)
2.0

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
On this page
Show Source