Metadata-Version: 2.4
Name: ideal-gases
Version: 0.1.5
Summary: Compute classical and quantum Euler solutions
Keywords: riemann,euler,quantum,polylogarithm,cfd
Author-Email: "Manuel A. Diaz" <manuel.ade@gmail.com>
License-Expression: MIT
License-File: LICENSE
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: 3.14
Classifier: Topic :: Scientific/Engineering :: Physics
Project-URL: Homepage, https://github.com/wme7/ideal-gases
Project-URL: Repository, https://github.com/wme7/ideal-gases
Project-URL: Issues, https://github.com/wme7/ideal-gases/issues
Requires-Python: >=3.11
Requires-Dist: numpy>=2.0
Provides-Extra: plot
Requires-Dist: matplotlib>=3.9; extra == "plot"
Provides-Extra: progress
Requires-Dist: tqdm>=4.66; extra == "progress"
Provides-Extra: tables
Requires-Dist: scipy>=1.13; extra == "tables"
Description-Content-Type: text/markdown

# Classical and Quantum Ideal Gases

[![CI](https://github.com/wme7/ideal-gases/actions/workflows/ci.yml/badge.svg)](https://github.com/wme7/ideal-gases/actions/workflows/ci.yml)
[![Publish](https://github.com/wme7/ideal-gases/actions/workflows/publish.yml/badge.svg)](https://github.com/wme7/ideal-gases/actions/workflows/publish.yml)
[![PyPI](https://img.shields.io/pypi/v/ideal-gases)](https://pypi.org/project/ideal-gases/)

Ideal-gases is a collection of numerical solvers for:
- classical and quantum Euler inviscid gases,
- classical and quantum 1-D Navier–Stokes–Fourier (NSF) viscous gases,
- classical and quantum 1-D BGK, Shakhov, and ES-BGK solvers for rarefied gases,
- C++ kernels for the polylogarithm and fugacity inversion used to resolve the quantum equation of state.

This repository ports the MATLAB implementation found in [this NTU thesis](https://doi.org/10.6342/NTU.2015.00509) to Python 3.11. The polylog function has been ported from the MATLAB implementation to C++, fugacity inversion is a C++ kernel alongside it. For Euler equations, the Toro exact Riemann solver has been extended to support Fermi–Dirac (FD), Bose–Einstein (BE), and Maxwell–Boltzmann (MB) statistics following the example found in the anexes of [this NTU thesis](https://doi.org/10.6342/NTU.2015.00509).

## Requirements

- Python 3.11+
- [pip](https://pip.pypa.io/) (or [uv](https://docs.astral.sh/uv/))

Building from source additionally requires a C++17 compiler. See [DEVELOPER_GUIDE.md](DEVELOPER_GUIDE.md).

## Installation

```bash
pip install ideal-gases
```

For plotting (`euler plot`, interactive explorers):

```bash
pip install ideal-gases[plot]
```

For QBE-matched Shakhov (`matched_qbe`):

```bash
pip install ideal-gases[tables]
```

For progress bars (`--progress` on the numerical CLIs, and the solver `progress=` flag):

```bash
pip install ideal-gases[progress]
```

After install, the `euler`, `nsf`, `bgk`, `shakhov`, and `es` command-line tools are available.

## Interactive mode

Launch matplotlib widget explorers to build custom Riemann problems with sliders, statistic toggles (quantum), and Save/Reset controls. Y-axis limits autoscale automatically on each update.

```bash
euler interactive classical
euler interactive quantum
```

Seed the initial state from CLI flags or a JSON config (same fields as `euler solve`):

```bash
euler interactive classical --gamma 1.4 --t-end 0.5 --nx 101
euler interactive quantum --rho-l 2 --t-l 1.5 --n 3 --h 0.5
euler interactive classical --config case.json
```

Optional domain flags (`--x-min`, `--x-max`, `--x0`, `--nx`) default to an interactive Sod-tube layout (`x` in `[-10, 10]`, discontinuity at `x0=0`, `nx=1024`). Use `-f path.png` to set the **Save** button target; nothing is written until you click Save.

### Example usage

```bash
euler interactive quantum
```
Outputs a Sod shock tube problem resolved with a quantum Euler solver for all statistics. We deactivate the solutions of MB and BE to focus on the FD solution. Using the slider, we can vary the left and right states and the thermal scale parameter `h` and the number of degrees of freedom `n` of the gas.

In Fig. 5 of [Hu and Jing (2010)](https://www.researchgate.net/profile/Shi-Jin-5/publication/228568274_On_kinetic_flux_vector_splitting_schemes_for_quantum_Euler_equations/links/02e7e525327fccc0cf000000/On-kinetic-flux-vector-splitting-schemes-for-quantum-Euler-equations.pdf), a fictitious 2-d fermi gas degenerate regime is used to prove the accuracy of Kinetic Flux Vector Splitting schemes for quantum Euler equations. Using the interactive mode, we set `n` : 2 and set the left and right states ($\rho,u,\theta$). Using the `h` slider, we found that the degenerate gas is resolved approximately for `h` $\approx$ 3.71.
As show in the following figure:

![Sod shock tube](https://raw.githubusercontent.com/wme7/ideal-gases/master/figures/fermi_2d_gas_yang_hsieh_shi.png)


## Command-line mode

Compute exact solution profiles, save plots to PNG, and write CSV/JSON files with the solution fields.

### Classical Sod shock tube

```bash
euler solve classical \
  --rho-l 1 --u-l 0 --p-l 1 \
  --rho-r 0.125 --u-r 0 --p-r 0.1 \
  --t-end 0.25 --gamma 1.4 \
  --nx 101 -o sod.csv
```

### Quantum Euler

```bash
euler solve quantum \
  --rho-l 1 --u-l 0 --t-l 1 \
  --rho-r 0.125 --u-r 0 --t-r 0.25 \
  --t-end 0.20 --n 2 --h 0.1 --statistic FD \
  -o euler_fd.csv
```

Write separate files for FD, MB, and BE with `--all-statistics` (e.g. `euler_case7_FD.csv`, `euler_case7_MB.csv`, `euler_case7_BE.csv`):

```bash
euler solve quantum ... --all-statistics -o euler_case7
```

### Equilibrium inversions

Compute the fugacity from density and temperature:

```bash
euler fugacity --rho 1.0 --theta 1.0 --n 3 --h 1.0 --statistic FD
```

Recover fugacity, temperature and pressure from density and internal energy:

```bash
euler moments --rho 1.0 --e 1.5 --n 3 --h 1.0 --statistic FD
```

Use `-o result.json` to write JSON output instead of printing to stdout.

### Built-in benchmarks

```bash
euler toro 1 -o toro_test1.csv
euler list --toro

euler quantum-example 7 --all-statistics -o euler_eg7
euler list --quantum
```

### JSON config files

Define a problem in JSON and run it with `euler run` or pass `--config` to `euler solve`:

```bash
euler run --config case.json
euler solve classical --config case.json -o override.csv
```

Example `case.json`:

```json
{
  "mode": "quantum",
  "left": {"rho": 1.0, "u": 0.0, "theta": 1.0},
  "right": {"rho": 0.125, "u": 0.0, "theta": 0.25},
  "t_end": 0.20,
  "n": 2.0,
  "h": 0.1,
  "statistic": "FD",
  "all_statistics": true,
  "format": "json",
  "output": "euler_case7",
  "domain": {"x_min": 0.0, "x_max": 1.0, "x0": 0.5, "nx": 101}
}
```

Use `--format json` (or a `.json` output path) for JSON instead of CSV. CLI flags override values from the config file.

### Visualization

Save a classical Sod shock tube figure:

```bash
euler plot classical \
  --rho-l 1 --u-l 0 --p-l 1 \
  --rho-r 0.125 --u-r 0 --p-r 0.1 \
  --t-end 0.2 --gamma 1.4 --nx 101 \
  -f sod.png
```

Plot a single quantum statistic or compare FD/MB/BE:

```bash
euler plot quantum \
  --rho-l 1 --u-l 0 --t-l 1 \
  --rho-r 0.125 --u-r 0 --t-r 0.25 \
  --t-end 0.20 --n 2 --h 0.1 --statistic FD \
  -f qfd.png

euler plot quantum-example 7 --all-statistics -f eg7
```

With `--all-statistics`, `-f eg7` writes `eg7_panels.png` (3×6 grid) and `eg7_comparison.png` (overlay). Use `--layout panels|comparison|both` to select one or both (default: `both`). Add `--show` for an interactive window, or `-o` to export CSV/JSON in the same run.

### Example usage

In [Filbet, Hu and Jing (2010)](https://www.cambridge.org/core/journals/esaim-mathematical-modelling-and-numerical-analysis/article/abs/numerical-scheme-for-the-quantum-boltzmann-equation-withstiff-collision-terms/BFB7B0297D8BC201F9A2C9008F4894BC), the authors use a Sod shock tube initial condition with a fictitious 2-d fermi and bose gas to prove the accuracy of their numerical scheme in classical and quantum hydronamic regimes. These are examples 7 and 8, respectively, in the CLI plot tool.

```bash
euler plot quantum-example 7 --all-statistics -f sod_2d_gas_classical --layout comparison --show
```
yields the following plot:
![Sod shock tube](https://raw.githubusercontent.com/wme7/ideal-gases/master/figures/sod_2d_gas_classical_comparison.png)

```bash
euler plot quantum-example 8 --all-statistics -f sod_2d_gas_quantum --layout comparison --show
```
yields the following plot:
![Sod shock tube](https://raw.githubusercontent.com/wme7/ideal-gases/master/figures/sod_2d_gas_quantum_comparison.png)

## NSF, BGK, Shakhov, and ES CLIs

The `nsf`, `bgk`, `shakhov`, and `es` commands mirror `euler`’s `run` / `solve` / `plot` verbs (no interactive mode, Toro presets, or `--all-statistics`). They write cell-centered profiles. `solve` and `run` require `-o`. `plot` requires `-f` or `--show`.

```bash
nsf solve classical --kn 1e-3 --t-end 0.1 --dim 3 -o nsf.csv
nsf solve quantum --kn 1e-3 --h 1 --statistic FD -o qnsf.csv
nsf solve quantum --matched-qbe --kn 1 --statistic FD -o qnsf-matched.csv
nsf run --config case.json
nsf plot classical --kn 1e-3 --t-end 0.1 -f nsf.png
nsf plot --input nsf.csv -f nsf.png

bgk solve classical --kn 1e-3 --nv 64 -o bgk.csv
shakhov solve quantum --matched-qbe --kn 1 --statistic FD -o matched.csv

es solve classical --kn 1e-3 --t-end 0.1 --b -0.5 -o es.csv
es solve quantum --kn 1e-3 --h 1 --statistic FD -o qes.csv
es solve quantum --matched-qbe --kn 1 --statistic FD -o es-matched.csv
```

Classical Sod defaults are `(ρ,u,p) = (1,0,1) / (0.125,0,0.1)`. Quantum defaults use temperatures `(1, 0.8)` (the same `p/ρ` with `R = 1`). The spatial grid is cell-centered (`nx=200` by default). Kinetic solvers (`bgk`, `shakhov`, `es`) default to `nv=128` velocity nodes; pass `--nv` to override (the examples above that pass `--nv 64` do so on purpose). CSV files start with `# key=value` metadata, then a header `x,rho,u,T,p,q` (quantum adds `z`).

`es` fixes `dim=3`. Classical Holway `b` defaults to `-0.5` and is ignored on `es solve quantum` when `--matched-qbe` is set. QBE-matched NSF and Shakhov (`--matched-qbe`) freeze `dim=3`, `idof=0` and require FD or BE. NSF calls `nsf_matched_qbe`; Shakhov calls `shakhov_matched_qbe`. Matched ES-BGK is Fermi only (`FD`) and calls `es_matched_qbe`. BGK has no matching flag.

Example numerical JSON (`kn` is required; `nsf run` will not invent it from an Euler-only file):

```json
{
  "mode": "classical",
  "kn": 0.001,
  "dim": 3,
  "idof": 0,
  "left": {"rho": 1.0, "u": 0.0, "p": 1.0},
  "right": {"rho": 0.125, "u": 0.0, "p": 0.1},
  "t_end": 0.1,
  "domain": {"x_min": 0.0, "x_max": 1.0, "x0": 0.5, "nx": 200},
  "output": "nsf.csv"
}
```

## Python module

Import `ideal_gases` to compute classical and quantum Euler, NSF, and ES-BGK solutions in your own scripts.

### Classical Euler

```python
import numpy as np
from ideal_gases import classical_euler

x = np.linspace(0.0, 1.0, 101)
result = classical_euler(
    rho_l=1.0,
    u_l=0.0,
    p_l=1.0,
    rho_r=0.125,
    u_r=0.0,
    p_r=0.1,
    t_end=0.2,
    gamma=1.4,
    x=x,
    x0=0.5,
)
```

### Quantum Euler (FD / BE / MB)

Left and right states are given in terms of density `rho`, velocity `u`, and temperature `theta` (written `t` in the API). The solver converts these to effective pressures via the quantum EOS, then applies the Toro exact Riemann solver.

```python
import numpy as np
from ideal_gases import quantum_euler

x = np.linspace(0.0, 1.0, 101)
result = quantum_euler(
    rho_l=1.0,
    u_l=0.0,
    t_l=1.0,
    rho_r=0.125,
    u_r=0.0,
    t_r=0.25,
    t_end=0.20,
    n=2.0,          # degrees of freedom; gamma = (n+2)/n
    h=0.1,          # thermal scale parameter
    statistic="FD", # "FD", "BE", or "MB"
    x=x,
    x0=0.5,
)
```

This returns a `RiemannResult` object that contains the solution fields: `x`, `rho`, `ux`, `p`, `e`, `z` (fugacity), `t` (temperature), `mach`, `entropy`.

In the classical limit, MB statistics with `h → 0` recover the ideal-gas behaviour (pressures `p = rho * theta`).

### Classical NSF

1-D Navier–Stokes–Fourier for a polytropic ideal gas. Same Sod left/right states as the classical Euler example; `dim` in `{1, 2, 3}` is the translational dimension and `idof` (default 0, monatomic) is the internal DoF, so `γ = (n+2)/n` with `n = dim + idof`. For air-like `γ = 1.4` use `dim=3, idof=2`. Classical Euler remains a γ-law Riemann solver (`gamma=`). `kn` is the Knudsen number used by the Chapman–Enskog closure `μ = kn ρ T`. Default Prandtl number is Eucken `4γ / (9γ - 5)` (`2/3` when `n = 3`); pass `pr=` to override.

```python
import numpy as np
from ideal_gases import classical_nsf

x = np.linspace(0.0, 1.0, 101)
result = classical_nsf(
    rho_l=1.0,
    u_l=0.0,
    p_l=1.0,
    rho_r=0.125,
    u_r=0.0,
    p_r=0.1,
    t_end=0.2,
    dim=3,
    kn=0.01,
    idof=2.0,  # γ = 1.4; omit or 0 for monatomic γ = 5/3
    x=x,
    x0=0.5,
)
```

This returns a `ClassicalNSFResult` object that contains the cell-centered fields: `rho`, `u`, `t` (temperature), `p`, `q` (heat flux).

### Classical Shakhov

Reduced `(g, h)` asymptotic-preserving DVM for 1-D flow of a 2-D or 3-D polytropic gas. Time is first-order IMEX. Default closures match classical NSF (`τ = kn`, Eucken `Pr = 4γ / (9γ - 5)`). Optional `idof` (default 0) sets `γ = (n+2)/n` with `n = dim + idof`; `h` holds perpendicular translational energy plus internal energy. Optional `viscosity` and `prandtl` hooks replace those closures (for example Boltzmann CE `μ ∝ T^ω`); there is no conductivity hook.

```python
import numpy as np
from ideal_gases import classical_shakhov

x = np.linspace(0.0, 1.0, 101)
result = classical_shakhov(
    rho_l=1.0,
    u_l=0.0,
    p_l=1.0,
    rho_r=0.125,
    u_r=0.0,
    p_r=0.1,
    t_end=0.2,
    dim=3,
    kn=0.01,
    idof=2.0,  # γ = 1.4; omit or 0 for monatomic γ = 5/3
    x=x,
    x0=0.5,
    nv=64,
    xi_max=8.0,
)
```

This returns the same `ClassicalNSFResult` fields as `classical_nsf`.

### Classical BGK

Reduced `(g, h)` asymptotic-preserving DVM that relaxes to the Maxwellian (`Pr = 1`). Time is first-order IMEX. Default viscosity is `μ = kn ρ T` (`τ = kn`). Optional `idof` (default 0) is the same polytropic parameter as Shakhov/NSF. Optional `viscosity` replaces that closure; there is no Prandtl hook, because BGK cannot independently match Boltzmann heat conductivity.

```python
import numpy as np
from ideal_gases import classical_bgk

x = np.linspace(0.0, 1.0, 101)
result = classical_bgk(
    rho_l=1.0,
    u_l=0.0,
    p_l=1.0,
    rho_r=0.125,
    u_r=0.0,
    p_r=0.1,
    t_end=0.2,
    dim=3,
    kn=0.01,
    idof=2.0,  # γ = 1.4; omit or 0 for monatomic γ = 5/3
    x=x,
    x0=0.5,
    nv=64,
    xi_max=8.0,
)
```

This returns the same `ClassicalNSFResult` fields as `classical_nsf`.

### Classical ES-BGK

1-D Holway ES-BGK on a full 3-D velocity lattice (`dim=3` only). Default monatomic tuning is `b = -1/2`, `μ = kn ρ T`, and `τ = (1-b) μ / p` (Eucken `Pr = 2/3`).

```python
import numpy as np
from ideal_gases import classical_es_bgk

x = np.linspace(0.0, 1.0, 101)
result = classical_es_bgk(
    rho_l=1.0,
    u_l=0.0,
    p_l=1.0,
    rho_r=0.125,
    u_r=0.0,
    p_r=0.1,
    t_end=0.2,
    kn=0.01,
    dim=3,
    b=-0.5,
    x=x,
    x0=0.5,
    nv=64,
    xi_max=8.0,
)
```

This returns a `ClassicalESResult` (`rho`, `u`, `t`, `p`, `q`).

### Quantum NSF (FD / BE / MB)

1-D Navier–Stokes–Fourier with the quantum EOS. Same left/right states as the quantum Euler example; `dim` is the translational dimension (Euler's `n` in the monatomic case), and `idof` (default 0) sets `γ = (n+2)/n` with `n = dim + idof`. `kn` sets the Chapman–Enskog viscosity `μ = kn p(z)`. Default conductivity uses `Pr_z = Pr_Eucken R(z)` with the same Eucken factor as classical NSF. Quantum Euler uses a single `n` (set `n = dim + idof` for a polytropic comparison).

```python
import numpy as np
from ideal_gases import quantum_nsf

x = np.linspace(0.0, 1.0, 101)
result = quantum_nsf(
    rho_l=1.0,
    u_l=0.0,
    t_l=1.0,
    rho_r=0.125,
    u_r=0.0,
    t_r=0.25,
    t_end=0.20,
    dim=3,
    h=0.1,
    kn=0.01,
    idof=2.0,  # γ = 1.4; omit or 0 for monatomic
    statistic="FD", # "FD", "BE", or "MB"
    x=x,
    x0=0.5,
)
```

This returns a `QuantumNSFResult` object (`NSFResult` is an alias) that contains: `rho`, `u`, `t`, `p`, `z` (fugacity), `q` (heat flux).

### Quantum Shakhov (FD / BE / MB)

Reduced `(g, h_r)` asymptotic-preserving DVM for 1-D flow of a 2-D or 3-D polytropic gas. Time is first-order IMEX. Default closures are quantum CE (`τ = kn`, `Pr_z = Pr_Eucken R(z)`; `2/3 R(z)` when `n = 3`). Optional `idof` (default 0) sets `γ = (n+2)/n` with `n = dim + idof`; `h_r` holds perpendicular translational energy plus internal energy (Diaz extra momenta). Optional `viscosity` and `prandtl` hooks replace those closures (for example QBE tables); there is no conductivity hook.

```python
import numpy as np
from ideal_gases import quantum_shakhov

x = np.linspace(0.0, 1.0, 101)
result = quantum_shakhov(
    rho_l=1.0,
    u_l=0.0,
    t_l=1.0,
    rho_r=0.125,
    u_r=0.0,
    t_r=0.25,
    t_end=0.20,
    dim=3,
    h=0.1,
    kn=0.01,
    idof=2.0,  # γ = 1.4; omit or 0 for monatomic
    statistic="FD",  # "FD", "BE", or "MB"
    x=x,
    x0=0.5,
    nv=64,
    xi_max=8.0,
)
```

This returns the same `QuantumNSFResult` fields as `quantum_nsf`.

### Quantum QBE-matched NSF, Shakhov, and ES-BGK

1-D flow of a 3-D monatomic Fermi or Bose gas with Wu CE transport tables Φ_μ(z, θ), Φ_κ(z, θ). `nsf_matched_qbe` is `quantum_nsf` and `shakhov_matched_qbe` is `quantum_shakhov`, each with `dim=3`, `idof=0`, and closures

```text
μ = Kn √T Φ_μ(z, θ),    κ = Kn (15/4) √T Φ_κ(z, θ),    Pr_z = Pr_Q / R(z)
```

`es_matched_qbe` is the same 3-D monatomic setup for Fermi gases only. Holway `b(z) = 1 - R(z)/Pr_Q(z)` and `τ = (1-b) μ_Q / p` follow Wu Section 4. Bose gases stay on `shakhov_matched_qbe`.

`θ = T/T_a` is a collision-kernel parameter (default 0 is the s-wave limit), not the hydrodynamic temperature in √T. Requires `pip install ideal-gases[tables]`. Maxwell–Boltzmann is not tabulated.

```python
import numpy as np
from ideal_gases.matched_qbe import es_matched_qbe, nsf_matched_qbe, shakhov_matched_qbe

x = np.linspace(0.0, 1.0, 101)
kin = shakhov_matched_qbe(
    rho_l=1.0,
    u_l=0.0,
    t_l=1.0,
    rho_r=0.4,
    u_r=0.0,
    t_r=0.6,
    t_end=0.10,
    h=3.0,
    kn=0.01,
    statistic="BE",  # "FD" or "BE"
    theta=0.0,       # T/T_a
    x=x,
    x0=0.5,
    nv=64,
    xi_max=8.0,
)
nsf = nsf_matched_qbe(
    rho_l=1.0,
    u_l=0.0,
    t_l=1.0,
    rho_r=0.4,
    u_r=0.0,
    t_r=0.6,
    t_end=0.10,
    h=3.0,
    kn=0.01,
    statistic="BE",
    theta=0.0,
    x=x,
    x0=0.5,
)
es = es_matched_qbe(
    rho_l=1.0,
    u_l=0.0,
    t_l=1.0,
    rho_r=0.4,
    u_r=0.0,
    t_r=0.6,
    t_end=0.10,
    h=3.0,
    kn=0.01,
    statistic="FD",  # Fermi only
    theta=0.0,
    x=x,
    x0=0.5,
    nv=64,
    xi_max=8.0,
)
```

This returns the same `QuantumNSFResult` fields as `quantum_nsf`.

### Quantum ES-BGK (FD / BE / MB)

1-D quantum ES-BGK on a full 3-D velocity lattice (`dim=3` only). Default transport is Chapman–Enskog `μ = kn p(z)` with Holway `b = -1/2`. Pass `viscosity` for a custom closure; QBE-matched tables use `es_matched_qbe` above.

```python
import numpy as np
from ideal_gases import quantum_es_bgk

x = np.linspace(0.0, 1.0, 101)
result = quantum_es_bgk(
    rho_l=1.0,
    u_l=0.0,
    t_l=1.0,
    rho_r=0.125,
    u_r=0.0,
    t_r=0.25,
    t_end=0.20,
    dim=3,
    h=0.1,
    kn=0.01,
    statistic="FD",  # "FD", "BE", or "MB"
    x=x,
    x0=0.5,
    nv=64,
    xi_max=8.0,
)
```

This returns the same `QuantumNSFResult` fields as `quantum_nsf`.

### Quantum BGK (FD / BE / MB)

Reduced `(g, h_r)` asymptotic-preserving DVM that relaxes to the FD/BE/MB equilibrium (`Pr_z = 1`). Time is first-order IMEX. Default viscosity is `μ = kn p(z)` (`τ = kn`). Optional `idof` (default 0) is the same polytropic parameter as Shakhov/NSF. Optional `viscosity` replaces that closure (`(rho, temp, z) → μ`); there is no Prandtl hook, because BGK cannot independently match heat conductivity.

```python
import numpy as np
from ideal_gases import quantum_bgk

x = np.linspace(0.0, 1.0, 101)
result = quantum_bgk(
    rho_l=1.0,
    u_l=0.0,
    t_l=1.0,
    rho_r=0.125,
    u_r=0.0,
    t_r=0.25,
    t_end=0.20,
    dim=3,
    h=0.1,
    kn=0.01,
    idof=2.0,  # γ = 1.4; omit or 0 for monatomic
    statistic="FD",  # "FD", "BE", or "MB"
    x=x,
    x0=0.5,
    nv=64,
    xi_max=8.0,
)
```

This returns the same `QuantumNSFResult` fields as `quantum_nsf`.

### Equilibrium inversions

Given density and temperature, recover the fugacity:

```python
from ideal_gases import find_fugacity

z = find_fugacity(rho=1.0, T=1.0, dim=3, h=1.0, eta=-1)
```

Given density and internal energy, recover fugacity, temperature and pressure:

```python
from ideal_gases import find_moments

z, T, p = find_moments(rho=1.0, e=2.5, dim=3, h=1.0, eta=-1, idof=2.0)
```

The `eta` parameter selects the statistic: `-1` Fermi, `0` classical (Maxwell-Boltzmann), `+1` Bose. Optional `idof` (default 0) is the internal DoF; thermodynamic closures use `n = dim + idof`.

Given density and the mass-weighted second-moment tensor `W`, recover the ES fugacity and the tensor `λ`. `dim` is the velocity-space dimension (`2` or `3`). This path uses the C++ `find_fugacity_es` kernel and does not take `idof`.

```python
import numpy as np
from ideal_gases import find_fugacity_anisotropic

w = np.eye(3)
z, lam = find_fugacity_anisotropic(rho=1.0, w=w, dim=3, h=1.0, eta=-1)
```

### Polylogarithm module

Building from source compiles three C++ extensions: `_polylog`, `_find_fugacity`, and `_find_fugacity_es`. `polylog(n, z)` uses the Fukushima minimax Fermi–Dirac / Bose–Einstein integrals for supported half-integer orders on `z < 0` and `0 < z < 1`, with [Bhagat](https://doi.org/10.1016/S0010-4655(03)00294-7) / integer analytic branches as fallback. `find_fugacity` and `find_fugacity_anisotropic` call the fugacity kernels.

We can use the polylogarithm module on our scripts as follows:

```python
import numpy as np
from ideal_gases import polylog

polylog(2, 0.5)                         # scalar
polylog(1.5, np.linspace(0.2, 0.9, 50)) # array
```

We can plot the polylogarithm function to verify the accuracy of the implementation for integer and half-integer orders as follows:

```bash
uv run python scripts/plot_polylogarithms.py --output figures/polylogarithms.png
```

Omitting `--output` writes `PolyLogPlot.png` in the repo root. The command above yields the following plot:
![Polylogarithm](https://raw.githubusercontent.com/wme7/ideal-gases/master/figures/polylogarithms.png)

### Public API

```python
from ideal_gases import (
    G,
    ClassicalESResult,
    ClassicalNSFResult,
    QuantumNSFResult,
    RiemannResult,
    DEFAULT_HOLWAY_B_3D,
    adiabatic_index,
    classical_euler,
    classical_es_bgk,
    classical_nsf,
    classical_bgk,
    classical_shakhov,
    equilibrium_moments,
    eucken_prandtl,
    find_fugacity,
    find_fugacity_anisotropic,
    find_moments,
    polylog,
    quantum_euler,
    quantum_es_bgk,
    quantum_nsf,
    quantum_bgk,
    es_matched_qbe,
    matched_qbe,
    quantum_shakhov,
    w_from_pressure,
)
```

| Symbol | Role |
|--------|------|
| `polylog(n, z)` | Fast C++ polylogarithm (Fukushima + Bhagat/integer fallback) |
| `adiabatic_index(n)` | Returns γ = (n + 2) / n |
| `eucken_prandtl(n)` | Eucken Pr = 4γ / (9γ - 5); 2/3 at n = 3 |
| `classical_euler(...)` | Classical ideal-gas exact Euler Riemann solver |
| `quantum_euler(...)` | Quantum EOS + Toro exact Euler Riemann solver |
| `classical_nsf(...)` | 1-D classical Navier–Stokes–Fourier solver (`idof` for polytropic γ) |
| `classical_bgk(...)` | 1-D classical BGK AP DVM (`dim` in `{2, 3}`, `idof` for polytropic γ) |
| `classical_shakhov(...)` | 1-D classical Shakhov AP DVM (`dim` in `{2, 3}`, `idof` for polytropic γ) |
| `classical_es_bgk(...)` | 1-D classical Holway ES-BGK (`dim=3`, full 3-D velocity lattice) |
| `quantum_nsf(...)` | 1-D quantum Navier–Stokes–Fourier solver (`idof` for polytropic γ) |
| `quantum_bgk(...)` | 1-D quantum BGK AP DVM (`dim` in `{2, 3}`, `idof` for polytropic γ) |
| `quantum_shakhov(...)` | 1-D quantum Shakhov AP DVM (`dim` in `{2, 3}`, `idof` for polytropic γ) |
| `quantum_es_bgk(...)` | 1-D quantum ES-BGK (`dim=3`, full 3-D velocity lattice) |
| `matched_qbe` | QBE-matched module (`from ideal_gases import matched_qbe`) |
| `nsf_matched_qbe(...)` | QBE-matched 3-D monatomic NSF (`FD`/`BE`, Wu Φ_μ, Φ_κ tables); import from `ideal_gases.matched_qbe` |
| `shakhov_matched_qbe(...)` | QBE-matched 3-D monatomic Shakhov (`FD`/`BE`, Wu Φ_μ, Φ_κ tables); import from `ideal_gases.matched_qbe` |
| `es_matched_qbe(...)` | QBE-matched 3-D monatomic ES-BGK (`FD` only, Wu §4 `b(z)`); also `from ideal_gases import es_matched_qbe` or `ideal_gases.matched_qbe` |
| `RiemannResult` | Euler solution profiles on the spatial grid |
| `ClassicalNSFResult` | Classical NSF fields (`rho`, `u`, `t`, `p`, `q`) |
| `ClassicalESResult` | Classical ES-BGK fields (`rho`, `u`, `t`, `p`, `q`) |
| `QuantumNSFResult` | Quantum NSF fields (`rho`, `u`, `t`, `p`, `z`, `q`) |
| `G(n, z, eta)` | Bose / Fermi / classical partition function |
| `equilibrium_moments(z, T, ...)` | Forward map (z, T) → (ρ, e) |
| `find_fugacity(rho, T, ...)` | Invert (ρ, T) → z (C++ `_find_fugacity`) |
| `find_fugacity_anisotropic(rho, w, ...)` | Invert (ρ, W) → (z, λ) for the ES reference state (`dim` in `{2, 3}`, C++ `_find_fugacity_es`) |
| `find_moments(rho, e, ...)` | Invert (ρ, e) → (z, T, p) |
| `w_from_pressure(P, b, dim)` | Holway map `W = (1-b) p I + b P` with `p = tr(P)/D` |
| `DEFAULT_HOLWAY_B_3D` | Default Holway parameter `b = -1/2` for 3-D monatomic ES-BGK |

## License

MIT License. See [LICENSE](LICENSE) for the full text.

Copyright (c) 2026 Manuel A. Diaz

For building from source, tests, linting, CI, and releases, see [DEVELOPER_GUIDE.md](DEVELOPER_GUIDE.md).
