This module calculates analytical solutions for a Sod shock tube problem.

A shock tube consists of a long tube filled with the same gas in two different physical states. The tube is divided into two parts, separated by a diaphragm. The initial state is defined by the values for density, pressure and velocity.

When the membrane bursts, discontinuity between the two initial states breaks into leftward and rightward moving waves, separated by a contact surface.

- Region 5 is the low pressure gas which is not disturbed by the shock wave.
- Regions 3 and 4 (divided by the contact surface) contain the gas immediately behind the shock traveling at a constant speed.
- The contact surface across which the density and the temperature are discontinuous lies within this zone.
- The zone between the head and the tail of the expansion fan is noted as Region 2. In this zone, the flow properties gradually change since the expansion process is isentropic.
- Region 1 denotes the undisturbed high pressure gas.

Sod solver returns after time of evolution the following variables: 1. Positions of head and foot of rarefation wave, contact discontinuity and shock 2. Constant pressure, density and velocity for all regions except rarefaction region 3. Pressure, density and velocity sampled across the domain of interest

The usage should be straightforward:

```python from shocktubecalc import sod

positions, regions, values = sod.solve(left*state=(1, 1, 0), right*state=(0.1, 0.125, 0.),
geometry=(0., 1., 0.5), t=0.2, gamma=1.4, npts=500)
```

`solve`

functionThe `positions`

is a dictionary of region names and their positions in the pipe after t=0.2 seconds.

Printing positions for each and every region:
```
python
print('Positions:')
for desc, vals in positions.items():
print('{0:10} : {1}'.format(desc, vals))
```

The `regions`

is a dictionary of constant pressure, density and velocity for
all regions except for the rarefaction wave.
```
python
print('States:')
for desc, vals in regions.items():
print('{0:10} : {1}'.format(desc, vals))
```

The `values`

is a dictionary of *pressure* , *density* , *velocity* and *x*
arrays. If one would want to plot said values as a function of x :
```python
import matplotlib.pyplot as pl
plt.figure(1)
plt.plot(values['x'], values['p'], linewidth=1.5, color='b')
plt.ylabel('pressure')
plt.xlabel('x')
plt.axis([0, 1, 0, 1.1])

plt.figure(2) plt.plot(values['x'], values['rho'], linewidth=1.5, color='r') plt.ylabel('density') plt.xlabel('x') plt.axis([0, 1, 0, 1.1])

plt.figure(3) plt.plot(values['x'], values['u'], linewidth=1.5, color='g') plt.ylabel('velocity') plt.xlabel('x')

plt.show() ```

One can easily calculate energy, temperature, speed of sound and a Mach number:

- Energy e = p/(gamma-1)
- Temperature T = p/rho
- Speed of sound c = (gamma p /rho)^0.5
- Mach number M = u/c

```
python
e = values['p']/(gamma -1)
T = values['p']/values['rho']
c = np.sqrt(gamma *values['p']/values['rho'])
M = values['p']/c
```