# plasma_plots.theory.exact

*module*

Exact solutions to verify simulations: the Riemann problem of gas dynamics (with vacuum), the
dam break, diffusion, advection, the wave equation, and pressureless flow and its caustics.

Profiles take the coordinates first and the time second, ``f(x, t, ...)``, so that
``lambda x, t: f(x, t, ...)`` (or a field of the returned result) can be passed as the
``reference=`` of plasma-plots' profile plots. Everything is plain numpy and vectorized: the
coordinates and times broadcast against each other.

Where a state can be empty (a vacuum in gas dynamics, a dry bed in shallow water) its density
or depth is zero and its velocity is ``nan`` (undefined, so plots leave a gap); momenta and
energies there are zero.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L1-L1)

## plasma_plots.theory.exact.LagrangianFlow

*class* · *dataclass*

```python
class LagrangianFlow
```

A pressureless (ballistic) flow followed along its fluid elements.

**Attributes**

- `position` (`float or numpy.ndarray`) — The Eulerian position x = q + v0(q) t of each element.
- `velocity` (`float or numpy.ndarray`) — Its velocity v0(q), constant in time.
- `density` (`float or numpy.ndarray`) — The density at x, ρ0(q)/|1 + v0'(q) t|; ``inf`` at a caustic.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L733-L749)

### plasma_plots.theory.exact.LagrangianFlow.density

*attribute* · *instance attribute*

```python
density: object
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L749-L749)

### plasma_plots.theory.exact.LagrangianFlow.position

*attribute* · *instance attribute*

```python
position: object
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L747-L747)

### plasma_plots.theory.exact.LagrangianFlow.velocity

*attribute* · *instance attribute*

```python
velocity: object
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L748-L748)

## plasma_plots.theory.exact.RiemannSolution

*class* · *dataclass*

```python
class RiemannSolution
```

The primitive variables of a gas-dynamics solution.

**Attributes**

- `density` (`float or numpy.ndarray`) — The mass density ρ.
- `velocity` (`float or numpy.ndarray`) — The velocity u; ``nan`` in a vacuum.
- `pressure` (`float or numpy.ndarray`) — The pressure p.
- `internal_energy` (`float or numpy.ndarray`) — The specific internal energy e = p/((γ − 1) ρ); ``nan`` in a vacuum.
- `gamma` (`float`) — The ratio of specific heats γ.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L90-L145)

### plasma_plots.theory.exact.RiemannSolution.density

*attribute* · *instance attribute*

```python
density: object
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L108-L108)

### plasma_plots.theory.exact.RiemannSolution.gamma

*attribute* · *instance attribute*

```python
gamma: float
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L112-L112)

### plasma_plots.theory.exact.RiemannSolution.internal_energy

*attribute* · *instance attribute*

```python
internal_energy: object
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L111-L111)

### plasma_plots.theory.exact.RiemannSolution.pressure

*attribute* · *instance attribute*

```python
pressure: object
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L110-L110)

### plasma_plots.theory.exact.RiemannSolution.velocity

*attribute* · *instance attribute*

```python
velocity: object
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L109-L109)

### plasma_plots.theory.exact.RiemannSolution.energy

*property*

```python
energy
```

The total energy density ρu²/2 + p/(γ − 1) (zero in a vacuum).

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L126-L132)

### plasma_plots.theory.exact.RiemannSolution.momentum

*property*

```python
momentum
```

The momentum density ρu (zero in a vacuum).

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L115-L123)

### plasma_plots.theory.exact.RiemannSolution.sound_speed

*property*

```python
sound_speed
```

The sound speed a = √(γ p/ρ); ``nan`` in a vacuum.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L135-L145)

## plasma_plots.theory.exact.ShallowWaterSolution

*class* · *dataclass*

```python
class ShallowWaterSolution
```

The depth and velocity of a shallow-water solution.

**Attributes**

- `depth` (`float or numpy.ndarray`) — The water depth h (zero on a dry bed).
- `velocity` (`float or numpy.ndarray`) — The depth-averaged velocity u; ``nan`` on a dry bed.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L479-L498)

### plasma_plots.theory.exact.ShallowWaterSolution.depth

*attribute* · *instance attribute*

```python
depth: object
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L491-L491)

### plasma_plots.theory.exact.ShallowWaterSolution.velocity

*attribute* · *instance attribute*

```python
velocity: object
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L492-L492)

### plasma_plots.theory.exact.ShallowWaterSolution.discharge

*property*

```python
discharge
```

The discharge (momentum per unit density) hu, zero on a dry bed.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L495-L498)

## plasma_plots.theory.exact.StarState

*class* · *dataclass*

```python
class StarState
```

The star region of an Euler Riemann problem, between the left and the right wave.

**Attributes**

- `pressure` (`float`) — The star pressure p* (zero if a vacuum forms).
- `velocity` (`float`) — The star (contact) velocity u*; ``nan`` if a vacuum forms or a side is a vacuum.
- `density_left` (`float`) — The density ρ*L left of the contact (zero for a vacuum).
- `density_right` (`float`) — The density ρ*R right of the contact (zero for a vacuum).
- `left_wave` (`str`) — ``"shock"`` or ``"rarefaction"`` (``"none"`` if the left state is a vacuum).
- `right_wave` (`str`) — ``"shock"`` or ``"rarefaction"`` (``"none"`` if the right state is a vacuum).
- `vacuum` (`bool`) — Whether the solution contains a vacuum (given, or generated between two rarefactions).

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L59-L87)

### plasma_plots.theory.exact.StarState.density_left

*attribute* · *instance attribute*

```python
density_left: float
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L83-L83)

### plasma_plots.theory.exact.StarState.density_right

*attribute* · *instance attribute*

```python
density_right: float
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L84-L84)

### plasma_plots.theory.exact.StarState.left_wave

*attribute* · *instance attribute*

```python
left_wave: str
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L85-L85)

### plasma_plots.theory.exact.StarState.pressure

*attribute* · *instance attribute*

```python
pressure: float
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L81-L81)

### plasma_plots.theory.exact.StarState.right_wave

*attribute* · *instance attribute*

```python
right_wave: str
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L86-L86)

### plasma_plots.theory.exact.StarState.vacuum

*attribute* · *instance attribute*

```python
vacuum: bool
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L87-L87)

### plasma_plots.theory.exact.StarState.velocity

*attribute* · *instance attribute*

```python
velocity: float
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L82-L82)

## plasma_plots.theory.exact.advected

*function*

```python
def advected(profile, x, t, velocity, period=None)
```

Compute a profile carried unchanged at constant velocity, u(x, t) = u0(x − v t).

**Parameters**

- `profile` (`callable`) — The initial profile u0, a vectorized function of the position.
- `x` (`float or array_like`) — Positions.
- `t` (`float or array_like`) — Times, broadcast against ``x``.
- `velocity` (`float`) — The advection velocity v.
- `period` (`float or (float, float)`) (default: `None`) — For a periodic domain: its length L (the domain is [0, L)) or its bounds (a, b). The shifted position is wrapped into the domain before ``profile`` is called. Default: ``None`` (not periodic).

**Returns**

- (`float or numpy.ndarray`) — u0(x − v t).

**Examples**

```pycon
>>> # −0.25 wrapped into [0, 1)
>>> round(float(advected(lambda x: x, 0.1, 0.35, 1.0, period=1.0)), 6)
0.75
>>> advected(lambda x: np.exp(-(x**2)), [0.0, 1.0], 1.0, 1.0).round(4)
array([0.3679, 1.    ])
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L624-L663)

## plasma_plots.theory.exact.caustic_time

*function*

```python
def caustic_time(velocity, q)
```

Compute the time of the first caustic (shell crossing) of a pressureless flow.

Elements with v0'(q) < 0 meet at t = −1/v0'(q), so the first caustic forms at
t_c = −1/min v0'(q), at the steepest decrease of the initial velocity.

**Parameters**

- `velocity` (`callable`) — The initial velocity v0(q), a vectorized function.
- `q` (`array_like`) — The positions where v0' is sampled (the minimum is taken over these, so resolve it).

**Returns**

- (`float`) — t_c; ``inf`` if the velocity nowhere decreases.

> **See Also**
>
> [`pressureless()`][plasma_plots.theory.exact.pressureless] : The flow up to (and past) the caustic.

> **References**
>
> Ya. B. Zel'dovich, "Gravitational instability: an approximate theory for large density
> perturbations", Astron. Astrophys. 5, 84 (1970).

**Examples**

```pycon
>>> q = np.linspace(-np.pi, np.pi, 1001)
>>> round(caustic_time(lambda q: -0.5 * np.sin(q), q), 6)
2.0
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L897-L931)

## plasma_plots.theory.exact.dalembert

*function*

```python
def dalembert(x, t, initial, speed, initial_rate=None)
```

Compute d'Alembert's solution of the 1-D wave equation ∂²u/∂t² = c² ∂²u/∂x².

u = [u0(x − ct) + u0(x + ct)]/2 + (1/(2c)) ∫_{x−ct}^{x+ct} v0(s) ds, for the initial profile
u0 and the initial rate v0 = ∂u/∂t at t = 0, on the whole line. The integral is done by
64-point Gauss–Legendre quadrature, exact to rounding for smooth v0 over a few wavelengths.
A pulse at rest splits into two halves moving at ±c.

**Parameters**

- `x` (`float or array_like`) — Positions.
- `t` (`float or array_like`) — Times, broadcast against ``x``.
- `initial` (`callable`) — The initial profile u0, a vectorized function of the position.
- `speed` (`float`) — The wave speed c > 0.
- `initial_rate` (`callable`) (default: `None`) — The initial rate v0, a vectorized function of the position. Default: ``None`` (at rest).

**Returns**

- (`float or numpy.ndarray`) — u(x, t).

**Raises**

- `ValueError` — If ``speed`` is not positive.

> **References**
>
> J. le Rond d'Alembert, "Recherches sur la courbe que forme une corde tenduë mise en
> vibration", Hist. Acad. R. Sci. Berlin 3, 214 (1747); L. C. Evans, Partial Differential
> Equations, 2nd ed. (AMS, 2010), section 2.4.

**Examples**

```pycon
>>> pulse = lambda x: np.exp(-(x**2) / 0.01)
>>> # the pulse split in two halves
>>> dalembert([0.0, 1.0], 1.0, pulse, 1.0).round(4)
array([0. , 0.5])
>>> # sin(ct)/c
>>> round(
...     float(dalembert(0.0, 0.5, np.zeros_like, 1.0, initial_rate=np.cos)), 6
... )
0.479426
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L669-L727)

## plasma_plots.theory.exact.dam_break

*function*

```python
def dam_break(x, t, depth=1.0, gravity=1.0, x0=0.0)
```

Compute Ritter's exact dam-break solution of the 1-D shallow-water equations on a dry bed.

Water at rest of depth h0 for x < x0 is released at t = 0 onto a dry bed. A rarefaction runs
between x0 − c0 t and the front at x0 + 2 c0 t, where c0 = √(g h0), with h = (2c0 − ξ)²/(9g)
and u = 2(c0 + ξ)/3 for ξ = (x − x0)/t.

**Parameters**

- `x` (`float or array_like`) — Positions.
- `t` (`float or array_like`) — Times ≥ 0, broadcast against ``x``.
- `depth` (`float`) (default: `1.0`) — The initial depth h0 behind the dam (x < x0). Default: ``1.0``.
- `gravity` (`float`) (default: `1.0`) — The gravitational acceleration g. Default: ``1.0``.
- `x0` (`float`) (default: `0.0`) — The position of the dam. Default: ``0.0``.

**Returns**

- (`ShallowWaterSolution`) — Depth and velocity (``nan`` on the dry bed), with the broadcast shape of ``x`` and ``t``; also the discharge hu as a property.

**Raises**

- `ValueError` — If a time is negative, or the depth not positive.

> **References**
>
> A. Ritter, "Die Fortpflanzung der Wasserwellen", Z. Ver. Dtsch. Ing. 36, 947 (1892).

**Examples**

```pycon
>>> w = dam_break([-2.0, 0.0, 1.0, 2.1], 1.0)
>>> w.depth.round(6), w.velocity.round(6)
(array([1.      , 0.444444, 0.111111, 0.      ]), array([0.      , 0.666667, 1.333333,      nan]))
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L501-L550)

## plasma_plots.theory.exact.heat_kernel

*function*

```python
def heat_kernel(x, t, diffusivity, width=0.0, center=0.0, mass=1.0)
```

Compute a spreading Gaussian, the solution of the diffusion equation ∂u/∂t = D ∇²u.

u = M (2π σ²)^(−d/2) exp(−|x − x_c|²/(2σ²)) with σ² = σ0² + 2 D t: the fundamental solution
for σ0 = 0, and a Gaussian of initial standard deviation σ0 otherwise. The integral is M at
all times; each coordinate's variance grows as σ0² + 2Dt.

**Parameters**

- `x` (`float, array_like or tuple of array_like`) — Positions: an array in 1-D, or a tuple ``(x, y)`` or ``(x, y, z)`` of arrays that broadcast against each other in d dimensions.
- `t` (`float or array_like`) — Times ≥ 0, broadcast against the positions.
- `diffusivity` (`float`) — The diffusion coefficient D.
- `width` (`float`) (default: `0.0`) — The initial standard deviation σ0. Default: ``0.0`` (a point source at t = 0).
- `center` (`float or tuple of float`) (default: `0.0`) — The center x_c, one value per dimension (a scalar is used for all). Default: ``0.0``.
- `mass` (`float`) (default: `1.0`) — The integral M of u. Default: ``1.0``.

**Returns**

- (`float or numpy.ndarray`) — u; at t = 0 with σ0 = 0, ``inf`` at the center and 0 elsewhere.

**Raises**

- `ValueError` — If a time or ``diffusivity`` or ``width`` is negative.

> **References**
>
> L. C. Evans, Partial Differential Equations, 2nd ed. (AMS, 2010), section 2.3.

**Examples**

```pycon
>>> round(float(heat_kernel(0.0, 0.5, 1.0)), 6)  # 1/√(4π D t)
0.398942
>>> round(float(heat_kernel((1.0, 0.0), 0.25, 1.0, width=1.0)), 6)
0.076026
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L565-L621)

## plasma_plots.theory.exact.pressureless

*function*

```python
def pressureless(q, t, velocity, density=None, velocity_derivative=None)
```

Follow a 1-D pressureless flow (the Zel'dovich approximation) in Lagrangian coordinates.

Without pressure every fluid element keeps its initial velocity, x(q, t) = q + v0(q) t, and
mass conservation ρ dx = ρ0 dq gives ρ = ρ0(q)/|∂x/∂q| = ρ0(q)/|1 + v0'(q) t|. Where the
velocity decreases, elements catch up with each other: the density becomes infinite (a
caustic) at [`caustic_time()`][plasma_plots.theory.exact.caustic_time], after which the map is multivalued (shell crossing) and
the density here is that of each stream.

**Parameters**

- `q` (`float or array_like`) — The initial (Lagrangian) positions of fluid elements.
- `t` (`float or array_like`) — Times, broadcast against ``q``.
- `velocity` (`callable`) — The initial velocity v0(q), a vectorized function.
- `density` (`callable`) (default: `None`) — The initial density ρ0(q), a vectorized function. Default: ``None`` (uniform, 1).
- `velocity_derivative` (`callable`) (default: `None`) — v0'(q). Default: ``None`` (by finite differences).

**Returns**

- (`LagrangianFlow`) — The positions, velocities and densities of the elements.

> **See Also**
>
> [`pressureless_eulerian()`][plasma_plots.theory.exact.pressureless_eulerian] : The same flow as a function of the position x.
> [`caustic_time()`][plasma_plots.theory.exact.caustic_time] : The time of the first shell crossing.

> **References**
>
> Ya. B. Zel'dovich, "Gravitational instability: an approximate theory for large density
> perturbations", Astron. Astrophys. 5, 84 (1970).

**Examples**

A converging sine velocity, v0 = −0.5 sin(q), with its caustic at t = 2:

```pycon
>>> flow = pressureless([0.0, np.pi / 2], 1.0, lambda q: -0.5 * np.sin(q))
>>> flow.position.round(6), flow.density.round(6)
(array([0.      , 1.070796]), array([2., 1.]))
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L758-L813)

## plasma_plots.theory.exact.pressureless_eulerian

*function*

```python
def pressureless_eulerian(x, t, velocity, density=None, velocity_derivative=None)
```

Compute a 1-D pressureless flow as a function of the position, before shell crossing.

Inverts x = q + v0(q) t for the Lagrangian position q (by Newton's method safeguarded by
bisection, to machine precision) and returns the flow there. The inverse is unique only
before the first caustic ([`caustic_time()`][plasma_plots.theory.exact.caustic_time]); later, one of the streams is returned.

**Parameters**

- `x` (`float or array_like`) — Positions.
- `t` (`float or array_like`) — Times, broadcast against ``x``.
- `velocity` (`callable`) — The initial velocity v0(q), a vectorized function.
- `density` (`callable`) (default: `None`) — The initial density ρ0(q), a vectorized function. Default: ``None`` (uniform, 1).
- `velocity_derivative` (`callable`) (default: `None`) — v0'(q). Default: ``None`` (by finite differences).

**Returns**

- (`LagrangianFlow`) — ``position`` is x itself; the velocity and density at x.

> **See Also**
>
> [`pressureless()`][plasma_plots.theory.exact.pressureless] : The flow along the fluid elements.

> **References**
>
> Ya. B. Zel'dovich, "Gravitational instability: an approximate theory for large density
> perturbations", Astron. Astrophys. 5, 84 (1970).

**Examples**

```pycon
>>> flow = pressureless_eulerian(
...     [-1.070796, 1.070796], 1.0, lambda q: -0.5 * np.sin(q)
... )
>>> flow.density.round(5), flow.velocity.round(5)
(array([1., 1.]), array([ 0.5, -0.5]))
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L816-L894)

## plasma_plots.theory.exact.riemann_euler

*function*

```python
def riemann_euler(x, t, left, right, gamma=1.4, x0=0.0)
```

Compute the exact solution of the Riemann problem of the 1-D Euler equations.

The gas is ideal, p = (γ − 1) ρ e, and initially in the state ``left`` for x < x0 and
``right`` for x > x0. The solution is self-similar in (x − x0)/t: a left wave (shock or
rarefaction), a contact discontinuity moving at u* and a right wave. All configurations are
covered, including a vacuum on either side (gas expanding into vacuum, bounded by a front
moving at u ± 2a/(γ − 1)) and the vacuum generated between two strong rarefactions.

**Parameters**

- `x` (`float or array_like`) — Positions.
- `t` (`float or array_like`) — Times ≥ 0, broadcast against ``x``. At t = 0 the initial states are returned.
- `left` (`(float, float, float)`) — The left state (density, velocity, pressure); ``(0, u, 0)`` is a vacuum.
- `right` (`(float, float, float)`) — The right state (density, velocity, pressure); ``(0, u, 0)`` is a vacuum.
- `gamma` (`float`) (default: `1.4`) — The ratio of specific heats γ > 1. Default: ``1.4``.
- `x0` (`float`) (default: `0.0`) — The position of the initial discontinuity. Default: ``0.0``.

**Returns**

- (`RiemannSolution`) — Density, velocity, pressure and specific internal energy, with the broadcast shape of ``x`` and ``t``; also the momentum and energy densities as properties.

**Raises**

- `ValueError` — If ``gamma`` ≤ 1, a time is negative, or a state is invalid (see [`star_state()`][plasma_plots.theory.exact.star_state]).

> **See Also**
>
> [`star_state()`][plasma_plots.theory.exact.star_state] : The star pressure and velocity alone.
> [`sod_shock_tube()`][plasma_plots.theory.exact.sod_shock_tube] : The classic test case.

> **References**
>
> E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed. (Springer,
> 2009), chapter 4 (4.9 for vacuum).

**Examples**

Toro's test 1 at t = 0.25, left of, inside and right of the star region:

```pycon
>>> w = riemann_euler(
...     [-0.4, 0.1, 0.5], 0.25, (1.0, 0.0, 1.0), (0.125, 0.0, 0.1)
... )
>>> w.density.round(5), w.pressure.round(5)
(array([1.     , 0.42632, 0.125  ]), array([1.     , 0.30313, 0.1    ]))
```

Gas expanding into vacuum: the front moves at 2a/(γ − 1) = 5.9161 here.

```pycon
>>> w = riemann_euler([0.0, 3.0, 6.0], 1.0, (1.0, 0.0, 1.0), (0.0, 0.0, 0.0))
>>> w.density.round(5), w.velocity.round(5)
(array([0.40188, 0.01169, 0.     ]), array([0.98601, 3.48601,     nan]))
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L338-L431)

## plasma_plots.theory.exact.sod_shock_tube

*function*

```python
def sod_shock_tube(x, t, gamma=1.4, x0=0.5)
```

Compute the exact solution of Sod's shock tube.

The Riemann problem with (ρ, u, p) = (1, 0, 1) left and (0.125, 0, 0.1) right of x0, usually
on [0, 1] up to t = 0.2 (Toro's test 1): a left rarefaction, a contact and a right shock.

**Parameters**

- `x` (`float or array_like`) — Positions.
- `t` (`float or array_like`) — Times ≥ 0, broadcast against ``x``.
- `gamma` (`float`) (default: `1.4`) — The ratio of specific heats. Default: ``1.4``.
- `x0` (`float`) (default: `0.5`) — The position of the diaphragm. Default: ``0.5``.

**Returns**

- (`RiemannSolution`) — Density, velocity, pressure and specific internal energy.

> **See Also**
>
> [`riemann_euler()`][plasma_plots.theory.exact.riemann_euler] : Any Riemann problem.

> **References**
>
> G. A. Sod, "A survey of several finite difference methods for systems of nonlinear
> hyperbolic conservation laws", J. Comput. Phys. 27, 1 (1978).

**Examples**

```pycon
>>> w = sod_shock_tube([0.1, 0.6, 0.8, 0.9], 0.2)
>>> w.density.round(5)
array([1.     , 0.42632, 0.26557, 0.125  ])
>>> float(w.velocity[2].round(5)), float(w.pressure[2].round(5))
(0.92745, 0.30313)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L434-L473)

## plasma_plots.theory.exact.star_state

*function*

```python
def star_state(left, right, gamma=1.4)
```

Solve for the star region of the Riemann problem of the 1-D Euler equations.

The star pressure p* is the root of Toro's pressure function f_L(p) + f_R(p) + u_R − u_L,
found by Newton's method safeguarded by bisection to machine precision. A wave is a shock
if p* exceeds the pressure ahead of it and a rarefaction otherwise. If a side is a vacuum
(density and pressure zero), or if 2(a_L + a_R)/(γ − 1) ≤ u_R − u_L so that the two
rarefactions leave a vacuum between them, p* = 0.

**Parameters**

- `left` (`(float, float, float)`) — The left state (density, velocity, pressure).
- `right` (`(float, float, float)`) — The right state (density, velocity, pressure).
- `gamma` (`float`) (default: `1.4`) — The ratio of specific heats γ > 1. Default: ``1.4``.

**Returns**

- (`StarState`) — p*, u*, the densities on both sides of the contact and the types of the two waves.

**Raises**

- `ValueError` — If ``gamma`` ≤ 1, or a state has a negative density or pressure, or only one of the
two zero.

> **See Also**
>
> [`riemann_euler()`][plasma_plots.theory.exact.riemann_euler] : The whole solution.

> **References**
>
> E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed. (Springer,
> 2009), chapter 4.

**Examples**

Toro's test 1 (Sod's shock tube):

```pycon
>>> s = star_state((1.0, 0.0, 1.0), (0.125, 0.0, 0.1))
>>> (
...     round(s.pressure, 5),
...     round(s.velocity, 5),
...     round(s.density_left, 5),
...     round(s.density_right, 5),
... )
(0.30313, 0.92745, 0.42632, 0.26557)
>>> s.left_wave, s.right_wave
('rarefaction', 'shock')
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/exact.py#L209-L295)