SOC examples updated to use new slicing

Numpy-like array slicing for fluids and elements added; get_n() function added for fluids and elements
Version changed to 1.1.3; config files updated
2025-09-03 15:57:55 +02:00 · 2025-09-03 15:57:55 +02:00 · 2025-09-03 15:57:55 +02:00 · 2025-09-03 15:57:55 +02:00 · 2025-09-03 15:57:55 +02:00 · 2025-07-28 17:11:23 +02:00
43 changed files with 2582 additions and 1240 deletions
--- a/.flake8
+++ b/.flake8
@ -14,6 +14,7 @@ exclude =
    dist,
    .conda,
    tests/autogenerated_*,
+    docs/source/api
    .venv,
    venv

--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@ -27,6 +27,13 @@ jobs:
      run: |
        python -m pip install --upgrade pip
        python -m pip install -e .[dev]
+        if [ "${{ matrix.python-version }}" = "3.13" ]; then
+          python -m pip install cffconvert
+        fi
+
+    - name: Validate CITATION.cff
+      if: ${{ matrix.python-version == '3.13' }}
+      run: cffconvert --validate

    - name: Prepare data and compile thermo database
      run: |
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@ -9,7 +9,7 @@ permissions:
  contents: write

 jobs:
-  build-and-deploy:
+  build:
    runs-on: ubuntu-latest
    steps:
      - uses: actions/checkout@v4
@ -17,24 +17,48 @@ jobs:
        uses: actions/setup-python@v3
        with:
          python-version: "3.x"
-      - name: Install dependencies
-        run: pip install sphinx sphinx_rtd_theme sphinx-autodoc-typehints myst-parser
-      - name: Generate Class List
+      - name: Build database
        run: |
-          echo "Create a dummy file to ensure the class list generation works"
-          mkdir -p src/gaspype/data
-          printf 'gapy\x00\x00\x00\x00' > src/gaspype/data/therm_data.bin
-          pip install .
-          python ./docs/source/generate_class_list.py
+          pip install pyyaml
+          python thermo_data/combine_data.py thermo_data/combined_data.yaml thermo_data/nasa9*.yaml
+          python thermo_data/compile_to_bin.py thermo_data/combined_data.yaml src/gaspype/data/therm_data.bin
+        # echo "Create a dummy file to ensure gaspype does't crash"
+        # mkdir -p src/gaspype/data
+        # printf 'gapy\x00\x00\x00\x00' > src/gaspype/data/therm_data.bin
+      - name: Install gaspype and dependencies
+        run: |
+          pip install .[doc_build]
+          python -m ipykernel install --user --name temp_kernel --display-name "Python (temp_kernel)"
+      - name: Generate Docs
+        run: python ./docs/source/generate_class_list.py
+      - name: Generate Examples
+        run: python ./docs/source/render_examples.py
      - name: Build Docs
        run: |
+          cp LICENSE docs/source/LICENSE.md
          cd docs
          sphinx-apidoc -o ./source/ ../src/ -M --no-toc
          rm ./source/*.rst
          make html
          touch ./build/html/.nojekyll
-      - name: Deploy to GitHub Pages
-        uses: JamesIves/github-pages-deploy-action@v4
+          mkdir -p ./build/html/_autogenerated
+          cp ./build/html/api/* ./build/html/_autogenerated/
+      - name: Upload artifact
+        uses: actions/upload-pages-artifact@v3
        with:
-          branch: gh-pages
-          folder: docs/build/html
+          path: docs/build/html
+
+  deploy:
+    needs: build
+    runs-on: ubuntu-latest
+    permissions:
+      contents: read
+      pages: write
+      id-token: write
+    environment:
+      name: github-pages
+      url: ${{ steps.deployment.outputs.page_url }}
+    steps:
+      - name: Deploy to GitHub Pages
+        id: deployment
+        uses: actions/deploy-pages@v4
--- a/.github/workflows/publish.yml
+++ b/.github/workflows/publish.yml
@ -10,6 +10,10 @@ jobs:
    name: Build and publish
    runs-on: ubuntu-latest

+    environment:
+      name: pypi
+      url: https://pypi.org/project/${{ github.event.repository.name }}/
+
    steps:
      - uses: actions/checkout@v3

--- a/.gitignore
+++ b/.gitignore
@ -9,8 +9,9 @@ __pycache__
 .pytest_cache
 tests/autogenerated_*.py
 docs/build/
-docs/source/modules.md
+docs/source/api/
 venv/
 .venv/
 thermo_data/combined_data.yaml
 src/gaspype/data/therm_data.bin
+*/*/build/
--- a/CITATION.cff
+++ b/CITATION.cff
@ -1,4 +1,5 @@
-cff-version: 1.2.0
+cff-version: 1.1.0
+message: "If you use this software, please cite it as below."
 title: Gaspype
 abstract: Gaspype is a performant library for thermodynamic calculations with ideal gases
 authors:
@ -8,11 +9,11 @@ authors:
    affiliation: "German Aerospace Center (DLR)"
    address: "Linder Höhe"
    city: Köln
-version: 1.0.1
-date-released: "2025-05-09"
+version: v1.1.3
+date-released: "2025-06-24"
 #identifiers:
 #  - description: This is the collection of archived snapshots of all versions of Gaspype
 #    type: doi
 #    value: ""
-license: MIT License
+license: MIT
 repository-code: "https://github.com/DLR-Institute-of-Future-Fuels/gaspype"
--- a/README.md
+++ b/README.md
@ -1,20 +1,21 @@
 # Gaspype
-The python package provides an performant library for thermodynamic calculations
+The Python package provides a performant library for thermodynamic calculations
 like equilibrium reactions for several hundred gas species and their mixtures -
-written in Python/Numpy.
+written in Python/NumPy.

 Species are treated as ideal gases. Therefore the application is limited to moderate
 pressures or high temperature applications.

-Its designed with goal to be portable to Numpy-style GPU frameworks like JAX and PyTorch.
+It is designed with goal to be portable to NumPy-style GPU frameworks like JAX and PyTorch.

-## Key features
- Tensor operations to prevent bottlenecks by the python interpreter
- Immutable types
- Elegant pythonic interface
- Great readable and compact syntax when using this package
- Good usability in Jupyter Notebook
- High performance for multidimensional fluid arrays
+## Key Features
+
+- Pure Python implementation with NumPy vectorization for high performance
+- Immutable types and comprehensive type hints for reliability
+- Intuitive, Pythonic API for both rapid prototyping and complex multidimensional models
+- Ready for Jupyter Notebook and educational use
+- Designed for future GPU support (JAX, PyTorch)
+- Ships with a comprehensive NASA9-based species database

 ## Installation
 Installation with pip:
@ -52,7 +53,7 @@ mass = fl.get_mass()
 gas_volume = fl.get_v(t=800+273.15, p=1e5)
 ```

-The arguments can be provided as numpy-arrays:
+The arguments can be provided as NumPy-arrays:

 ``` python
 import numpy as np
@ -62,7 +63,8 @@ fl.get_density(t=t_range, p=1e5)
 ```
 array([0.10122906, 0.09574625, 0.09082685, 0.08638827, 0.08236328])
 ```
-A ```fluid``` object can have multiple compositions. A multidimensional ```fluid``` object can be created for example by multiplication with a numpy array:
+A ```fluid``` object can have multiple compositions. A multidimensional ```fluid``` object
+can be created for example by multiplication with a NumPy array:

 ``` python
 fl2 = gp.fluid({'H2O': 1, 'N2': 2}) + \
@ -125,7 +127,8 @@ array([[[0.        , 0.5       , 0.5       ],
 ```

 ### Elements
-In some cases not the molecular but the atomic composition is of interest. The ```elements``` class can be used for atom based balances and works similar:
+In some cases not the molecular but the atomic composition is of interest.
+The ```elements``` class can be used for atom based balances and works similar:

 ``` python
 el = gp.elements({'N': 1, 'Cl': 2})
@ -134,7 +137,9 @@ el.get_mass()
 ```
 np.float64(0.08490700000000001)
 ```
-A ```elements``` object can be as well instantiated from a ```fluid``` object. Arithmetic operations between ```elements``` and ```fluid``` result in an ```elements``` object:
+A ```elements``` object can be as well instantiated from a ```fluid``` object.
+Arithmetic operations between ```elements``` and ```fluid``` result in
+an ```elements``` object:
 ``` python
 el2 = gp.elements(fl) + el - 0.3 * fl
 el2
@ -146,7 +151,8 @@ N                1.000e+00 mol
 O                7.000e-01 mol
 ```

-Going from an atomic composition to an molecular composition is a little bit less straight forward, since there is no universal approach. One way is to calculate the thermodynamic equilibrium for a mixture:
+Going from an atomic composition to an molecular composition is possible as well.
+One way is to calculate the thermodynamic equilibrium for a mixture:

 ``` python
 fs = gp.fluid_system('CH4, H2, CO, CO2, O2')
@ -163,7 +169,13 @@ CO2              33.07 %
 O2                0.00 %
 ```

-The ```equilibrium``` function can be called with a ```fluid``` or ```elements``` object as first argument. ```fluid``` and ```elements``` referencing a ```fluid_system``` object witch can be be set as shown above during the object instantiation. If not provided, a new one will be created automatically. Providing a ```fluid_system``` gives more control over which molecular species are included in derived ```fluid``` objects. Furthermore arithmetic operations between objects with the same ```fluid_system``` are potentially faster:
+The ```equilibrium``` function can be called with a ```fluid``` or ```elements``` object
+as first argument. ```fluid``` and ```elements``` referencing a ```fluid_system``` object
+witch can be be set as shown above during the object instantiation. If not provided,
+a new one will be created automatically. Providing a ```fluid_system``` gives more
+control over which molecular species are included in derived ```fluid``` objects.
+Furthermore arithmetic operations between objects with the same ```fluid_system```
+are potentially faster:

 ``` python
 fl3 + gp.fluid({'CH4': 1}, fs)
@ -177,7 +189,9 @@ CO2              18.07 %
 O2                0.00 %
 ```

-Especially if the ```fluid_system``` of one of the operants has not a subset of molecular species of the other ```fluid_system``` a new ```fluid_system``` will be created for the operation which might degrade performance:
+Especially if the ```fluid_system``` of one of the operants has not a subset of
+molecular species of the other ```fluid_system``` a new ```fluid_system``` will
+be created for the operation which might degrade performance:

 ``` python
 fl3 + gp.fluid({'NH3': 1})
@ -207,8 +221,8 @@ cd gaspype
 It's recommended to setup an venv:

 ```bash
-python -m venv venv
-source venv/bin/activate  # On Windows use `venv\Scripts\activate`
+python -m venv .venv
+source .venv/bin/activate  # On Windows use `.venv\Scripts\activate`
 ```

 Install the package and dev-dependencies while keeping the package files
--- a/docs/media/soc_inverted.svg
+++ b/docs/media/soc_inverted.svg
@ -0,0 +1,193 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<!-- Created with Inkscape (http://www.inkscape.org/) -->
+
+<svg
+   width="164.70512mm"
+   height="86.731506mm"
+   viewBox="0 0 164.70512 86.731506"
+   version="1.1"
+   id="svg1"
+   inkscape:version="1.3.2 (091e20e, 2023-11-25, custom)"
+   sodipodi:docname="soc_inverted.svg"
+   inkscape:export-filename="soc_export.svg"
+   inkscape:export-xdpi="96"
+   inkscape:export-ydpi="96"
+   xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
+   xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
+   xmlns="http://www.w3.org/2000/svg"
+   xmlns:svg="http://www.w3.org/2000/svg">
+  <sodipodi:namedview
+     id="namedview1"
+     pagecolor="#ffffff"
+     bordercolor="#000000"
+     borderopacity="0.25"
+     inkscape:showpageshadow="2"
+     inkscape:pageopacity="0.0"
+     inkscape:pagecheckerboard="0"
+     inkscape:deskcolor="#d1d1d1"
+     inkscape:document-units="mm"
+     inkscape:zoom="1.0041055"
+     inkscape:cx="425.25411"
+     inkscape:cy="332.13641"
+     inkscape:window-width="1920"
+     inkscape:window-height="1009"
+     inkscape:window-x="1912"
+     inkscape:window-y="-8"
+     inkscape:window-maximized="1"
+     inkscape:current-layer="layer1" />
+  <defs
+     id="defs1">
+    <marker
+       style="overflow:visible"
+       id="Triangle"
+       refX="0"
+       refY="0"
+       orient="auto-start-reverse"
+       inkscape:stockid="Triangle arrow"
+       markerWidth="2"
+       markerHeight="1"
+       viewBox="0 0 1 1"
+       inkscape:isstock="true"
+       inkscape:collect="always"
+       preserveAspectRatio="none">
+      <path
+         transform="scale(0.5)"
+         style="fill:context-stroke;fill-rule:evenodd;stroke:context-stroke;stroke-width:1pt"
+         d="M 5.77,0 -2.88,5 V -5 Z"
+         id="path135" />
+    </marker>
+  </defs>
+  <g
+     inkscape:label="Ebene 1"
+     inkscape:groupmode="layer"
+     id="layer1"
+     transform="translate(-18.027099,-19.242669)">
+    <path
+       style="fill:none;stroke:#dadada;stroke-width:0.465;stroke-dasharray:none;stroke-opacity:1;marker-start:url(#Triangle);marker-end:url(#Triangle)"
+       d="M 23.663111,22.545169 V 97.81995 H 179.42973"
+       id="path1" />
+    <path
+       style="fill:none;stroke:#dadada;stroke-width:0.372829;stroke-dasharray:none;stroke-opacity:1"
+       d="M 170.85884,97.81995 V 95.304579"
+       id="path2" />
+    <circle
+       style="fill:#dddddd;fill-opacity:1;stroke:#dadada;stroke-width:0.468334;stroke-dasharray:none;stroke-opacity:1"
+       id="path3"
+       cx="40.768127"
+       cy="46.63377"
+       r="1.2658331" />
+    <circle
+       style="fill:#dddddd;fill-opacity:1;stroke:#dadada;stroke-width:0.468334;stroke-dasharray:none;stroke-opacity:1"
+       id="path3-7"
+       cx="60.487274"
+       cy="54.953339"
+       r="1.2658331" />
+    <circle
+       style="fill:#dddddd;fill-opacity:1;stroke:#dadada;stroke-width:0.468334;stroke-dasharray:none;stroke-opacity:1"
+       id="path3-6"
+       cx="94.775032"
+       cy="67.060722"
+       r="1.2658331" />
+    <circle
+       style="fill:#dddddd;fill-opacity:1;stroke:#dadada;stroke-width:0.468334;stroke-dasharray:none;stroke-opacity:1"
+       id="path3-8"
+       cx="133.16524"
+       cy="76.384407"
+       r="1.2658331" />
+    <circle
+       style="fill:#dddddd;fill-opacity:1;stroke:#dadada;stroke-width:0.468334;stroke-dasharray:none;stroke-opacity:1"
+       id="path3-1"
+       cx="23.607611"
+       cy="32.419498"
+       r="1.2658331" />
+    <circle
+       style="fill:#dddddd;fill-opacity:1;stroke:#dadada;stroke-width:0.468334;stroke-dasharray:none;stroke-opacity:1"
+       id="path3-8-3"
+       cx="170.49799"
+       cy="83.681023"
+       r="1.2658331" />
+    <path
+       style="fill:none;fill-opacity:1;stroke:#dadada;stroke-width:0.465;stroke-dasharray:none;stroke-opacity:1"
+       d="m 23.61653,32.379465 17.141782,14.300346 19.796893,8.151661 34.283565,12.111043 38.38268,9.502507 37.26475,7.266625"
+       id="path4"
+       sodipodi:nodetypes="cccccc" />
+    <text
+       xml:space="preserve"
+       style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465;stroke-dasharray:none;stroke-opacity:1"
+       x="86.083488"
+       y="105.09981"
+       id="text4"><tspan
+         sodipodi:role="line"
+         id="tspan4"
+         style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465"
+         x="86.083488"
+         y="105.09981">Relative cell lenghs</tspan></text>
+    <text
+       xml:space="preserve"
+       style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465;stroke-dasharray:none;stroke-opacity:1"
+       x="169.54897"
+       y="103.00477"
+       id="text4-1"><tspan
+         sodipodi:role="line"
+         id="tspan4-1"
+         style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465"
+         x="169.54897"
+         y="103.00477">1</tspan></text>
+    <text
+       xml:space="preserve"
+       style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465;stroke-dasharray:none;stroke-opacity:1"
+       x="22.457224"
+       y="103.68204"
+       id="text4-1-1"><tspan
+         sodipodi:role="line"
+         id="tspan4-1-2"
+         style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465"
+         x="22.457224"
+         y="103.68204">0</tspan></text>
+    <path
+       style="fill:#000000;fill-opacity:1;stroke:#dadada;stroke-width:0.3;stroke-dasharray:none;stroke-opacity:1"
+       d="M 40.710985,49.035087 V 65.438055"
+       id="path5" />
+    <path
+       style="fill:#000000;fill-opacity:1;stroke:#dadada;stroke-width:0.3;stroke-dasharray:none;stroke-opacity:1"
+       d="m 60.546682,57.469645 v 7.773293"
+       id="path5-9" />
+    <path
+       style="fill:#000000;fill-opacity:1;stroke:#dadada;stroke-width:0.3;stroke-dasharray:none;stroke-opacity:1;marker-start:url(#Triangle);marker-end:url(#Triangle)"
+       d="M 58.396033,62.880678 H 42.878103"
+       id="path5-0"
+       sodipodi:nodetypes="cc" />
+    <text
+       xml:space="preserve"
+       style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465;stroke-dasharray:none;stroke-opacity:1"
+       x="48.511959"
+       y="61.409679"
+       id="text4-1-1-5"><tspan
+         sodipodi:role="line"
+         id="tspan4-1-2-2"
+         style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465"
+         x="48.511959"
+         y="61.409679">Δz</tspan></text>
+    <text
+       xml:space="preserve"
+       style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465;stroke-dasharray:none;stroke-opacity:1"
+       x="-72.630066"
+       y="21.160755"
+       id="text4-9"
+       transform="rotate(-90)"><tspan
+         sodipodi:role="line"
+         id="tspan4-5"
+         style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465"
+         x="-72.630066"
+         y="21.160755">Current</tspan></text>
+    <path
+       style="fill:#000000;fill-opacity:1;stroke:#dadada;stroke-width:0.3;stroke-dasharray:1.2, 0.3;stroke-dashoffset:0;stroke-opacity:1"
+       d="M 94.728796,67.061137 V 97.890816"
+       id="path6" />
+    <path
+       style="fill:#000000;fill-opacity:1;stroke:#dadada;stroke-width:0.3;stroke-dasharray:1.2, 0.3;stroke-dashoffset:0;stroke-opacity:1"
+       d="M 133.2181,76.566042 V 97.753468"
+       id="path6-1"
+       sodipodi:nodetypes="cc" />
+  </g>
+</svg>
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@ -10,14 +10,14 @@ import os
 import sys
 sys.path.insert(0, os.path.abspath("../src/"))

-project = 'gaspype'
+project = 'Gaspype'
 copyright = '2025, Nicolas Kruse'
 author = 'Nicolas Kruse'

 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration

-extensions = ["sphinx.ext.autodoc", "sphinx.ext.napoleon", "sphinx_autodoc_typehints", "myst_parser"]
+extensions = ["sphinx.ext.autodoc", "sphinx.ext.napoleon", "sphinx_autodoc_typehints", "myst_parser", "sphinx.ext.autosummary"]

 templates_path = ['_templates']
 exclude_patterns = []
@ -26,7 +26,8 @@ exclude_patterns = []
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output

 # html_theme = 'alabaster'
-html_theme = 'sphinx_rtd_theme'
+html_theme = 'pydata_sphinx_theme'
 html_static_path = ['_static']

 autodoc_inherit_docstrings = True
+autoclass_content = 'both'
--- a/docs/source/floatarray.md
+++ b/docs/source/floatarray.md
@ -0,0 +1,5 @@
+# gaspype.typing.FloatArray
+
+```{eval-rst}
+.. autoclass:: gaspype.typing.FloatArray
+```
--- a/docs/source/generate_class_list.py
+++ b/docs/source/generate_class_list.py
@ -1,61 +1,86 @@
+# This script generates the source md-files for all classes and functions for the docs
+
 import importlib
 import inspect
 import fnmatch
 from io import TextIOWrapper
+import os
+
+
+def write_manual(f: TextIOWrapper, doc_files: list[str], title: str) -> None:
+    write_dochtree(f, title, doc_files)


 def write_classes(f: TextIOWrapper, patterns: list[str], module_name: str, title: str, description: str = '', exclude: list[str] = []) -> None:
-
+    """Write the classes to the file."""
    module = importlib.import_module(module_name)

    classes = [
        name for name, obj in inspect.getmembers(module, inspect.isclass)
-        if (obj.__module__ == module_name and
-            any(fnmatch.fnmatch(name, pat) for pat in patterns if pat not in exclude) and
+        if (any(fnmatch.fnmatch(name, pat) for pat in patterns if pat not in exclude) and
            obj.__doc__ and '(Automatic generated stub)' not in obj.__doc__)
    ]

-    """Write the classes to the file."""
-    f.write(f'## {title}\n\n')
    if description:
        f.write(f'{description}\n\n')

+    write_dochtree(f, title, classes)
+
    for cls in classes:
-        f.write('```{eval-rst}\n')
-        f.write(f'.. autoclass:: {module_name}.{cls}\n')
-        f.write('   :members:\n')
-        f.write('   :class-doc-from: both\n')
-        f.write('   :undoc-members:\n')
-        f.write('   :show-inheritance:\n')
-        f.write('   :inherited-members:\n')
-        if title != 'Base classes':
-            f.write('   :exclude-members: select\n')
-        f.write('```\n\n')
+        with open(f'docs/source/api/{cls}.md', 'w') as f2:
+            f2.write(f'# {module_name}.{cls}\n')
+            f2.write('```{eval-rst}\n')
+            f2.write(f'.. autoclass:: {module_name}.{cls}\n')
+            f2.write('   :members:\n')
+            f2.write('   :undoc-members:\n')
+            f2.write('   :show-inheritance:\n')
+            f2.write('   :inherited-members:\n')
+            f2.write('```\n\n')


 def write_functions(f: TextIOWrapper, patterns: list[str], module_name: str, title: str, description: str = '', exclude: list[str] = []) -> None:
-
+    """Write the classes to the file."""
    module = importlib.import_module(module_name)

    functions = [
-        name for name, obj in inspect.getmembers(module, inspect.isfunction)
-        if (obj.__module__ == module_name and
-            any(fnmatch.fnmatch(name, pat) for pat in patterns if pat not in exclude))
+        name for name, _ in inspect.getmembers(module, inspect.isfunction)
+        if (any(fnmatch.fnmatch(name, pat) for pat in patterns if pat not in exclude))
    ]

-    """Write the classes to the file."""
-    f.write(f'## {title}\n\n')
    if description:
        f.write(f'{description}\n\n')

+    write_dochtree(f, title, functions)
+
    for func in functions:
        if not func.startswith('_'):
-            f.write('```{eval-rst}\n')
-            f.write(f'.. autofunction:: {module_name}.{func}\n')
+            with open(f'docs/source/api/{func}.md', 'w') as f2:
+                f2.write(f'# {module_name}.{func}\n')
+                f2.write('```{eval-rst}\n')
+                f2.write(f'.. autofunction:: {module_name}.{func}\n')
+                f2.write('```\n\n')
+
+
+def write_dochtree(f: TextIOWrapper, title: str, items: list[str]):
+    f.write('```{toctree}\n')
+    f.write(':maxdepth: 1\n')
+    f.write(f':caption: {title}:\n')
+    #f.write(':hidden:\n')
+    for text in items:
+        if not text.startswith('_'):
+            f.write(f"{text}\n")
    f.write('```\n\n')


-with open('docs/source/modules.md', 'w') as f:
-    f.write('# Functions and classes\n\n')
+if __name__ == "__main__":
+    # Ensure the output directory exists
+    os.makedirs('docs/source/api', exist_ok=True)
+
+    with open('docs/source/api/index.md', 'w') as f:
+        f.write('# Classes and functions\n\n')
+
        write_classes(f, ['*'], 'gaspype', title='Classes')
+
        write_functions(f, ['*'], 'gaspype', title='Functions')
+
+        write_manual(f, ['../ndfloat', '../floatarray'], title='Types')
--- a/docs/source/index.md
+++ b/docs/source/index.md
@ -1,9 +1,9 @@
 ```{toctree}
-:maxdepth: 2
-:caption: Contents:
-
-readme
-modules
+:maxdepth: 1
+:hidden:
+api/index
+api/examples
+repo
 ```

 ```{include} ../../README.md
--- a/docs/source/ndfloat.md
+++ b/docs/source/ndfloat.md
@ -0,0 +1,5 @@
+# gaspype.typing.NDFloat
+
+```{eval-rst}
+.. autoclass:: gaspype.typing.NDFloat
+```
--- a/docs/source/readme.md
+++ b/docs/source/readme.md
@ -1,2 +0,0 @@
-```{include} ../../README.md
-```
--- a/docs/source/render_examples.py
+++ b/docs/source/render_examples.py
@ -0,0 +1,66 @@
+# This script converts the example md-files as jupyter notebook,
+# execute the notebook and convert the notebook back to a md-file
+# with outputs included.
+
+import subprocess
+from glob import glob
+import os
+from io import TextIOWrapper
+
+
+def run_cmd(command: list[str]):
+    result = subprocess.run(command, capture_output=True, text=True)
+
+    print('> ' + ' '.join(command))
+    print(result.stdout)
+
+    assert (not result.stderr or
+            any('RuntimeWarning: ' in line for line in result.stderr.splitlines()) or
+            any('[NbConvertApp]' in line or 'Warning' in line for line in result.stderr.splitlines())), 'ERROR: ' + result.stderr
+
+
+def run_rendering(input_path: str, output_directory: str):
+    file_name = '.'.join(os.path.basename(input_path).split('.')[:-1])
+    assert file_name
+
+    print(f'- Convert {input_path} ...')
+    run_cmd(['notedown', input_path, '--to', 'notebook', '--output', f'{output_directory}/{file_name}.ipynb', '--run'])
+    run_cmd(['jupyter', 'nbconvert', '--to', 'markdown', f'{output_directory}/{file_name}.ipynb', '--output', f'{file_name}.md'])
+    run_cmd(['python', 'tests/md_to_code.py', 'script', f'{input_path}', f'{output_directory}/{file_name}.py'])
+
+
+def write_dochtree(f: TextIOWrapper, title: str, items: list[str]):
+    f.write('```{toctree}\n')
+    f.write(':maxdepth: 1\n')
+    #f.write(':hidden:\n')
+    #f.write(f':caption: {title}:\n')
+    for text in items:
+        if not text.startswith('_'):
+            f.write(f"{text}\n")
+    f.write('```\n\n')
+
+
+def render_examples(filter: str, example_file: str):
+    files = glob(filter)
+    names = ['.'.join(os.path.basename(path).split('.')[:-1]) for path in files]
+
+    with open(example_file, 'w') as f:
+        f.write('# Gaspype examples\n\n')
+        write_dochtree(f, '', [n for n in names if n.lower() != 'readme'])
+
+        f.write('## Download Jupyter Notebooks\n\n')
+        for path, name in zip(files, names):
+            if name.lower() != 'readme':
+                run_rendering(path, 'docs/source/api')
+                notebook = name + '.ipynb'
+                f.write(f'- [{notebook}]({notebook})\n\n')
+
+        f.write('## Download plain python files\n\n')
+        for path, name in zip(files, names):
+            if name.lower() != 'readme':
+                script_name = name + '.py'
+                f.write(f'- [{script_name}]({script_name})\n\n')
+
+
+if __name__ == "__main__":
+    render_examples('examples/*.md', 'docs/source/api/examples.md')
--- a/docs/source/repo.md
+++ b/docs/source/repo.md
@ -0,0 +1,3 @@
+# Code repository
+
+Code repository is on GitHub: [github.com/DLR-Institute-of-Future-Fuels/gaspype](https://github.com/DLR-Institute-of-Future-Fuels/gaspype).
--- a/examples/README.md
+++ b/examples/README.md
@ -0,0 +1,25 @@
+# Example scripts
+
+Examples can be looked-up in the
+[documentation](https://dlr-institute-of-future-fuels.github.io/gaspype/)
+rendered with results.
+
+The gaspype examples from this directory are available in the documentation as
+downloadable Jupyter Notebooks or plain python scripts with comments.
+
+The conversion is done like the following automated by the
+[docs/source/render_examples.py](../docs/source/render_examples.py) script:
+``` bash
+# Converting markdown with code sections to Jupyter Notebook and run it:
+notedown examples/soec_methane.md --to notebook --output docs/source/api/soec_methane.ipynb --run
+
+# Converting the Jupyter Notebook to Markdown and a folder with image
+# files placed in docs/source/api/:
+jupyter nbconvert --to markdown docs/source/api/soec_methane.ipynb --output soec_methane.md
+```
+
+A new example Markdown file can be created from a Jupyter Notebook running
+the following command:
+``` bash
+jupyter nbconvert --to markdown new_example.ipynb --NbConvertApp.use_output_suffix=False --ClearOutputPreprocessor.enabled=True --output-dir examples/ --output new_example.md
+```
--- a/examples/hydrogen_production_thermodynamics.md
+++ b/examples/hydrogen_production_thermodynamics.md
@ -0,0 +1,45 @@
+# Thermodynamics of Hydrogen Production
+
+The minimal energy required to produce hydrogen from liquid water is given by the
+Higher Heating Value (HHV). The HHV is the sum of the difference
+between the enthalpies of products and educts (LHV: Lower Heating Value) and
+the Heat of Evaporation for water.
+
+```python
+import gaspype as gp
+
+lhv = gp.fluid({'H2': 1, 'O2': 1/2, 'H2O': -1}).get_H(25 + 273.15)
+dh_v = 43990  # J/mol  (heat of evaporation for water @ 25 °C)
+
+hhv = lhv + dh_v
+print(f'LHV: {lhv/1e3:.1f} kJ/mol')
+print(f'HHV: {hhv/1e3:.1f} kJ/mol')
+```
+
+Thermodynamics also defines which part of the energy must be
+provided as work (e.g., electric power) and which part can be supplied
+as heat. This depends on temperature and pressure. For generating 1 bar
+of hydrogen the temperature dependency can be calculated as follows:
+
+```python
+import numpy as np
+import matplotlib.pyplot as plt
+
+t = np.linspace(0, 2000, 128)  # 0 to 2000 °C
+p = 1e5  # Pa (=1 bar)
+
+g_products = gp.fluid({'H2': 1, 'O2': 1/2, 'H2O': 0}).get_G(t + 273.15, p)
+g_educts =   gp.fluid({'H2': 0, 'O2': 0,   'H2O': 1}).get_G(t + 273.15, p)
+
+work = g_products - g_educts  # J/mol
+heat = lhv - work  # J/mol
+
+fig, ax = plt.subplots(figsize=(6, 4), dpi=120)
+ax.set_xlabel("Temperature / °C")
+ax.set_ylabel("Energy / kWh/kg")
+k = 1e-3 / 3600 / 0.002  # Conversion factor from J/mol to kWh/kg for hydrogen
+ax.stackplot(t, k * work, k * heat, k * dh_v * np.ones_like(t))
+```
+
+Green is the heat of evaporation, orange the additional heat provided at
+the given temperature and blue the work.
--- a/examples/methane_mixtures.md
+++ b/examples/methane_mixtures.md
@ -0,0 +1,60 @@
+# Methane Mixtures
+
+This example shows equilibria of methane mixed with steam and CO2
+
+
+```python
+import gaspype as gp
+import numpy as np
+import matplotlib.pyplot as plt
+```
+
+Setting temperature and pressure:
+
+
+```python
+t = 900 + 273.15
+p = 1e5
+
+fs = gp.fluid_system(['H2', 'H2O', 'CO2', 'CO', 'CH4', 'O2'])
+```
+
+Equilibrium calculation for methane steam mixtures:
+
+
+```python
+ratio = np.linspace(0.01, 1.5, num=64)
+
+fl = gp.fluid({'CH4': 1}, fs) + ratio * gp.fluid({'H2O': 1}, fs)
+equilibrium_h2o = gp.equilibrium(fl, t, p)
+```
+
+
+```python
+fig, ax = plt.subplots(figsize=(6, 4), dpi=120)
+ax.set_xlabel("H2O/CH4")
+ax.set_ylabel("molar fraction")
+ax.set_ylim(0, 1.1)
+#ax.set_xlim(0, 100)
+ax.plot(ratio, equilibrium_h2o.get_x())
+ax.legend(fs.active_species)
+```
+
+Equilibrium calculation for methane CO2 mixtures:
+
+
+```python
+fl = gp.fluid({'CH4': 1}, fs) + ratio * gp.fluid({'CO2': 1}, fs)
+equilibrium_co2 = gp.equilibrium(fl, t, p)
+```
+
+
+```python
+fig, ax = plt.subplots(figsize=(6, 4), dpi=120)
+ax.set_xlabel("CO2/CH4")
+ax.set_ylabel("molar fraction")
+ax.set_ylim(0, 1.1)
+#ax.set_xlim(0, 100)
+ax.plot(ratio, equilibrium_co2.get_x())
+ax.legend(fs.active_species)
+```
--- a/examples/soec_syngas.md
+++ b/examples/soec_syngas.md
@ -0,0 +1,124 @@
+# SOEC Co-Electrolysis
+
+This example shows a 1D isothermal SOEC (Solid oxide electrolyzer cell) model for
+converting carbon dioxide and steam into syngas.
+
+The operating parameters chosen here are not necessarily realistic. For example,
+a utilization of 0.95 causes issues with the formation of solid carbon.
+
+```python
+import gaspype as gp
+from gaspype.constants import R, F
+import numpy as np
+import matplotlib.pyplot as plt
+```
+
+Calculate equilibrium compositions for fuel and air sides in counter flow
+along the fuel flow direction:
+```python
+utilization = 0.95
+air_dilution = 0.2
+t = 700 + 273.15  # K
+p = 1e5  # Pa
+
+fs = gp.fluid_system('H2, H2O, O2, CH4, CO, CO2')
+feed_fuel = gp.fluid({'H2O': 2, 'CO2': 1}, fs)
+
+o2_full_conv = np.sum(gp.elements(feed_fuel).get_n(['C' ,'O']) * [-1/2, 1/2])
+
+feed_air = gp.fluid({'O2': 1, 'N2': 4}) * o2_full_conv * utilization * air_dilution
+
+conversion = np.linspace(0, utilization, 128)
+perm_oxygen = o2_full_conv * conversion * gp.fluid({'O2': 1})
+
+fuel_side = gp.equilibrium(feed_fuel - perm_oxygen, t, p)
+air_side  = gp.equilibrium(feed_air  + perm_oxygen, t, p)
+```
+
+Plot compositions of the fuel and air side:
+```python
+fig, ax = plt.subplots()
+ax.set_xlabel("Conversion")
+ax.set_ylabel("Molar fraction")
+ax.plot(conversion, fuel_side.get_x(), '-')
+ax.legend(fuel_side.species)
+
+fig, ax = plt.subplots()
+ax.set_xlabel("Conversion")
+ax.set_ylabel("Molar fraction")
+ax.plot(conversion, air_side.get_x(), '-')
+ax.legend(air_side.species)
+```
+
+Calculation of the oxygen partial pressures:
+```python
+o2_fuel_side = gp.oxygen_partial_pressure(fuel_side, t, p)
+o2_air_side = air_side.get_x('O2') * p
+```
+
+Plot oxygen partial pressures:
+```python
+fig, ax = plt.subplots()
+ax.set_xlabel("Conversion")
+ax.set_ylabel("Oxygen partial pressure / Pa")
+ax.set_yscale('log')
+ax.plot(conversion, np.stack([o2_fuel_side, o2_air_side], axis=1), '-')
+ax.legend(['o2_fuel_side', 'o2_air_side'])
+```
+
+The high oxygen partial pressure at the inlet is in reality lower.
+The assumption that gas inter-diffusion in the flow direction is slower
+than the gas velocity does not hold at this very high gradient. However
+often the oxygen partial pressure is still to high to prevent oxidation of the
+cell/electrode. This can be effectively prevented by recycling small amounts of
+the hydrogen riche output gas.
+
+Calculation of the local nernst potential between fuel and air side:
+```python
+z_O2 = 4
+nernst_voltage = R*t / (z_O2*F) * np.log(o2_air_side/o2_fuel_side)
+```
+
+Plot nernst potential:
+```python
+fig, ax = plt.subplots()
+ax.set_xlabel("Conversion")
+ax.set_ylabel("Voltage / V")
+ax.plot(conversion, nernst_voltage, '-')
+print(np.min(nernst_voltage))
+```
+
+The model uses between each node a constant conversion. Because
+current density depends strongly on the position along the cell
+the constant conversion does not relate to a constant distance.
+
+![Alt text](../../media/soc_inverted.svg)
+
+To calculate the local current density (**node_current**) as well
+as the total cell current (**terminal_current**) the (relative)
+physical distance between the nodes (**dz**) must be calculated:
+```python
+cell_voltage = 1.3  # V
+ASR = 0.2  # Ohm*cm²
+
+node_current = (nernst_voltage - cell_voltage) / ASR  # A/cm² (Current density at each node)
+
+current = (node_current[1:] + node_current[:-1]) / 2  # A/cm² (Average current density between the nodes)
+
+dz = 1/current / np.sum(1/current)  # Relative distance between each node
+
+terminal_current = np.sum(current * dz) # A/cm² (Total cell current per cell area)
+
+print(f'Terminal current: {terminal_current:.2f} A/cm²')
+```
+
+Plot the local current density:
+```python
+z_position = np.concatenate([[0], np.cumsum(dz)])  # Relative position of each node
+
+fig, ax = plt.subplots()
+ax.set_xlabel("Relative cell position")
+ax.set_ylabel("Current density / A/cm²")
+ax.plot(z_position, node_current, '-')
+```
+
--- a/examples/sofc_methane.md
+++ b/examples/sofc_methane.md
@ -0,0 +1,141 @@
+# SOFC with Methane
+
+This example shows a 1D isothermal SOFC (Solid oxide fuel cell) model.
+
+The operating parameters chosen here are not necessarily realistic due to
+constraints not included in this model. Under the example condition the
+formation of solid carbon is for example very likely. A pre-reforming step
+is typically applied when operating on methane.
+
+```python
+import gaspype as gp
+from gaspype.constants import R, F
+import numpy as np
+import matplotlib.pyplot as plt
+```
+
+Calculate equilibrium compositions for fuel and air sides
+in counter flow along the fuel flow direction:
+```python
+fuel_utilization = 0.90
+air_utilization = 0.5
+t = 800 + 273.15  # K
+p = 1e5  # Pa
+
+fs = gp.fluid_system('H2, H2O, O2, CH4, CO, CO2')
+feed_fuel = gp.fluid({'CH4': 1, 'H2O': 0.1}, fs)
+
+o2_full_conv = np.sum(gp.elements(feed_fuel).get_n(['H', 'C' ,'O']) * [1/4, 1, -1/2])
+
+feed_air = gp.fluid({'O2': 1, 'N2': 4}) * o2_full_conv * fuel_utilization / air_utilization
+
+conversion = np.linspace(0, fuel_utilization, 32)
+perm_oxygen = o2_full_conv * conversion * gp.fluid({'O2': 1})
+
+fuel_side = gp.equilibrium(feed_fuel + perm_oxygen, t, p)
+air_side  = gp.equilibrium(feed_air  - perm_oxygen, t, p)
+```
+
+Plot compositions of the fuel and air side:
+```python
+fig, ax = plt.subplots()
+ax.set_xlabel("Conversion")
+ax.set_ylabel("Molar fraction")
+ax.plot(conversion, fuel_side.get_x(), '-')
+ax.legend(fuel_side.species)
+
+fig, ax = plt.subplots()
+ax.set_xlabel("Conversion")
+ax.set_ylabel("Molar fraction")
+ax.plot(conversion, air_side.get_x(), '-')
+ax.legend(air_side.species)
+```
+
+Calculation of the oxygen partial pressures:
+```python
+o2_fuel_side = gp.oxygen_partial_pressure(fuel_side, t, p)
+o2_air_side = air_side.get_x('O2') * p
+```
+
+Plot oxygen partial pressures:
+```python
+fig, ax = plt.subplots()
+ax.set_xlabel("Conversion")
+ax.set_ylabel("Oxygen partial pressure / Pa")
+ax.set_yscale('log')
+ax.plot(conversion, np.stack([o2_fuel_side, o2_air_side], axis=1), '-')
+ax.legend(['o2_fuel_side', 'o2_air_side'])
+```
+
+Calculation of the local nernst potential between fuel and air side:
+```python
+z_O2 = 4  # Number of electrons transferred per O2 molecule
+nernst_voltage = R*t / (z_O2*F) * np.log(o2_air_side/o2_fuel_side)
+```
+
+Plot nernst potential:
+```python
+fig, ax = plt.subplots()
+ax.set_xlabel("Conversion")
+ax.set_ylabel("Voltage / V")
+ax.plot(conversion, nernst_voltage, '-')
+print(np.min(nernst_voltage))
+```
+
+The model uses between each node a constant conversion. Because
+current density depends strongly on the position along the cell
+the constant conversion does not relate to a constant distance.
+
+![Alt text](../../media/soc_inverted.svg)
+
+To calculate the local current density (**node_current**) as well
+as the total cell current (**terminal_current**) the (relative)
+physical distance between the nodes (**dz**) must be calculated:
+```python
+cell_voltage = 0.77  # V
+ASR = 0.2  # Ohm*cm²
+
+node_current = (nernst_voltage - cell_voltage) / ASR  # A/cm² (Current density at each node)
+
+current = (node_current[1:] + node_current[:-1]) / 2  # A/cm² (Average current density between the nodes)
+
+dz = 1/current / np.sum(1/current)  # Relative distance between each node
+
+terminal_current = np.sum(current * dz) # A/cm² (Total cell current per cell area)
+
+print(f'Terminal current: {terminal_current:.2f} A/cm²')
+```
+
+Plot the local current density:
+```python
+z_position = np.concatenate([[0], np.cumsum(dz)])  # Relative position of each node
+
+fig, ax = plt.subplots()
+ax.set_xlabel("Relative cell position")
+ax.set_ylabel("Current density / A/cm²")
+ax.plot(z_position, node_current, '-')
+```
+
+Based on the cell current and voltage the energy balance can be calculated.
+In the following the electric cell output power (often referred to as "DC power")
+and lower heating value (LHV) are calculated. The numbers here are per cell area for
+being cell and stack size independent. The quotient of both is often referred to as
+LHV based DC efficiency.
+
+```python
+dc_power = cell_voltage * terminal_current  # W/cm²
+print(f"DC power: {dc_power:.2f} W/cm²")
+
+lhv = gp.fluid({'CH4': 1, 'H2O': -2, 'CO2': -1}).get_H(25 + 273.15)  # J/mol (LHV of methane)
+
+# LHV based chemical input power:
+lhv_power = lhv * terminal_current / (2 * z_O2 * F)  # W/cm² (two O2 per CH4 for full oxidation)
+efficiency = dc_power / lhv_power
+print(f"LHV based DC efficiency: {efficiency*100:.1f} %")
+
+# Or by shortening the therms:
+lhv_voltage = lhv / (2 * z_O2 * F)  # V
+print(f"LHV voltage: {lhv_voltage:.2f} V")
+efficiency = cell_voltage / lhv_voltage  # LHV based DC efficiency
+print(f"LHV based DC efficiency: {efficiency*100:.1f} %")
+```
--- a/examples/sulfur_oxygen_equalibrium.md
+++ b/examples/sulfur_oxygen_equalibrium.md
@ -0,0 +1,58 @@
+# Sulfur Oxygen Equilibrium
+
+This example shows equilibrium calculations for sulfur/oxygen mixtures.
+
+```python
+import gaspype as gp
+import numpy as np
+import matplotlib.pyplot as plt
+```
+
+List possible sulfur/oxygen species:
+
+```python
+gp.species(element_names = 'S, O')
+```
+
+Or more specific by using regular expressions:
+
+
+```python
+gp.species('S?[2-3]?O?[2-5]?', use_regex=True)
+```
+
+Calculation of the molar equilibrium fractions for sulfur and oxygen depending on the oxygen to sulfur ratio: 
+
+
+```python
+fs = gp.fluid_system(['S2', 'S2O', 'SO2', 'SO3', 'O2'])
+
+oxygen_ratio = np.linspace(0.5, 3, num=128)
+el = gp.elements({'S': 1}, fs) + oxygen_ratio * gp.elements({'O': 1}, fs)
+
+composition = gp.equilibrium(el, 800+273.15, 1e4)
+
+fig, ax = plt.subplots()
+ax.set_xlabel("Oxygen to sulfur ratio")
+ax.set_ylabel("Molar fraction")
+ax.plot(oxygen_ratio, composition.get_x(), '-')
+ax.legend(composition.species)
+```
+
+Calculation of the molar equilibrium fractions for sulfur and oxygen depending on temperature in °C:
+
+
+```python
+fs = gp.fluid_system(['S2', 'S2O', 'SO2', 'SO3', 'O2'])
+
+el = gp.elements({'S': 1, 'O':2.5}, fs)
+
+t_range = np.linspace(500, 1300, num=32)
+composition = gp.equilibrium(el, t_range+273.15, 1e4)
+
+fig, ax = plt.subplots()
+ax.set_xlabel("Temperature / °C")
+ax.set_ylabel("Molar fraction")
+ax.plot(t_range, composition.get_x(), '-')
+ax.legend(composition.species)
+```
--- a/examples/thermodynamic_carbon_activity.md
+++ b/examples/thermodynamic_carbon_activity.md
@ -0,0 +1,90 @@
+# Carbon Activity
+
+This example shows the calculation of the carbon activity for methane mixtures
+in thermodynamic equilibrium.
+
+```python
+import gaspype as gp
+import numpy as np
+import matplotlib.pyplot as plt
+```
+
+Setting temperatures and pressure:
+
+
+```python
+t_range = np.array([600, 700, 800, 900, 1100, 1500])  # °C
+
+p = 1e5  # Pa
+
+fs = gp.fluid_system(['H2', 'H2O', 'CO2', 'CO', 'CH4'])
+```
+
+Equilibrium calculation for methane steam mixtures:
+
+
+```python
+ratio = np.linspace(0.01, 1.5, num=128)
+
+fl = gp.fluid({'CH4': 1}, fs) + ratio * gp.fluid({'H2O': 1}, fs)
+```
+
+gaspype.carbon_activity supports currently only 0D fluids therefore we build this helper function:
+
+
+```python
+def partial_c_activity(fl: gp.fluid, t: float, p: float):
+    fls = fl.array_composition.shape
+
+    eq_fl = gp.equilibrium(fl, t, p)
+
+    ret = np.zeros(fls[0])
+    for i in range(fls[0]):
+        ret[i] = gp.carbon_activity(gp.fluid(eq_fl.array_composition[i,:], fs), t, p)
+
+    return ret
+```
+
+Now we use the helper function to calculate the carbon activitx for all
+compositions in equilibrium_h2o times all temperatures in t_range:
+
+
+```python
+carbon_activity = np.vstack([partial_c_activity(fl, tc + 273.15, p) for tc in t_range])
+```
+
+Plot carbon activities, a activity of > 1 means there is thermodynamically the formation of sold carbon favored.
+
+
+```python
+fig, ax = plt.subplots(figsize=(6, 4), dpi=120)
+ax.set_xlabel("H2O/CH4")
+ax.set_ylabel("carbon activity")
+ax.set_ylim(1e-1, 1e3)
+ax.set_yscale('log')
+ax.plot(ratio, carbon_activity.T)
+ax.hlines(1, np.min(ratio), np.max(ratio), colors='k', linestyles='dashed')
+ax.legend([f'{tc} °C' for tc in t_range])
+```
+
+Let's do the equilibrium calculation for methane CO2 mixtures as well:
+
+
+```python
+fl_co2 = gp.fluid({'CH4': 1}, fs) + ratio * gp.fluid({'CO2': 1}, fs)
+carbon_activity_co2 = np.vstack([partial_c_activity(fl_co2, tc + 273.15, p) for tc in t_range])
+```
+
+And plot carbon activities over the CO2 to CH4 ratio:
+
+
+```python
+fig, ax = plt.subplots(figsize=(6, 4), dpi=120)
+ax.set_xlabel("CO2/CH4")
+ax.set_ylabel("carbon activity")
+ax.set_ylim(1e-1, 1e3)
+ax.set_yscale('log')
+ax.plot(ratio, carbon_activity_co2.T)
+ax.hlines(1, np.min(ratio), np.max(ratio), colors='k', linestyles='dashed')
+ax.legend([f'{tc} °C' for tc in t_range])
+```
--- a/pyproject.toml
+++ b/pyproject.toml
@ -1,6 +1,6 @@
 [project]
 name = "gaspype"
-version = "1.1.0"
+version = "1.1.3"
 authors = [
  { name="Nicolas Kruse", email="nicolas.kruse@dlr.de" },
 ]
@ -18,9 +18,10 @@ dependencies = [
 ]

 [project.urls]
-Homepage = "https://github.com/DLR-Institute-of-Future-Fuels/gaspype"
+Homepage = "https://dlr-institute-of-future-fuels.github.io/gaspype/"
+documentation = "https://dlr-institute-of-future-fuels.github.io/gaspype/api/"
+Repository = "https://github.com/DLR-Institute-of-Future-Fuels/gaspype.git"
 Issues = "https://github.com/DLR-Institute-of-Future-Fuels/gaspype/issues"
-documentation = "https://dlr-institute-of-future-fuels.github.io/gaspype/"

 [build-system]
 requires = ["setuptools>=61.0", "wheel"]
@ -30,7 +31,7 @@ build-backend = "setuptools.build_meta"
 where = ["src"]

 [tool.setuptools.package-data]
-gaspype = ["data/therm_data.bin"]
+gaspype = ["data/therm_data.bin", "py.typed"]

 [project.optional-dependencies]
 dev = [
@ -41,7 +42,20 @@ dev = [
    "cantera",
    "pyyaml>=6.0.1",
    "types-PyYAML",
-    "scipy-stubs"
+    "scipy-stubs",
+    "matplotlib"
+]
+doc_build = [
+  "sphinx",
+  "pydata_sphinx_theme",
+  "sphinx-autodoc-typehints",
+  "myst-parser",
+  "pandas",
+  "matplotlib",
+  "ipykernel",
+  "jupyter",
+  "nbconvert",
+  "notedown"
 ]

 [tool.mypy]
--- a/src/gaspype/init.py
+++ b/src/gaspype/init.py
--- a/src/gaspype/_main.py
+++ b/src/gaspype/_main.py
@ -0,0 +1,820 @@
+import numpy as np
+from numpy.typing import NDArray
+from typing import Sequence, Any, TypeVar, Iterator, overload, Callable
+from math import log as ln, ceil
+from scipy.linalg import null_space
+from gaspype._phys_data import atomic_weights, db_reader
+import re
+import pkgutil
+from .constants import R, epsy, p0
+from .typing import FloatArray, NDFloat, Shape, ArrayIndices
+
+T = TypeVar('T', 'fluid', 'elements')
+
+_data = pkgutil.get_data(__name__, 'data/therm_data.bin')
+assert _data is not None, 'Could not load thermodynamic data'
+_species_db = db_reader(_data)
+
+
+def species(pattern: str = '*', element_names: str | list[str] = [], use_regex: bool = False) -> list[str]:
+    """Returns a alphabetically sorted list of all available species
+    filtered by a pattern if supplied
+
+    Args:
+        pattern: Optional filter for specific molecules. Placeholder
+            characters: # A number including non written ones: 'C#H#' matches 'CH4';
+            $ Arbitrary element name; * Any sequence of characters
+        element_names: restrict results to species that contain only the specified elements.
+            The elements can be supplied as list of strings or as comma separated string.
+        use_regex: using regular expression for the pattern
+
+    Returns:
+        List of species
+    """
+    if isinstance(element_names, str):
+        elements = {s.strip() for s in element_names.split(',')}
+    else:
+        assert isinstance(element_names, list), 'type of element_names must be list or str'
+        elements = set(element_names)
+
+    for el in elements:
+        assert el in atomic_weights, f'element {el} unknown'
+
+    if not use_regex:
+        el_pattern = '|'.join([el for el in atomic_weights.keys()])
+        pattern = pattern.replace('*', '.*')
+        pattern = pattern.replace('#', '\\d*')
+        pattern = pattern.replace('$', '(' + el_pattern + ')')
+        pattern = '^' + pattern + '(,.*)?$'
+
+    if element_names == []:
+        return [sn for sn in _species_db.names if re.fullmatch(pattern, sn)]
+    else:
+        return [
+            s.name for s in _species_db
+            if re.fullmatch(pattern, s.name) and
+            (len(elements) == 0 or set(s.composition.keys()).issubset(elements))]
+
+
+class fluid_system:
+    """A class to represent a fluid_system defined by a set of selected species.
+
+    Attributes:
+        species_names (list[str]): List of selected species in the fluid_system
+        array_molar_mass (FloatArray): Array of the molar masses of the species in the fluid_system
+        array_element_composition (FloatArray): Array of the element composition of the species in the fluid_system.
+            Dimension is: (number of species, number of elements)
+        array_atomic_mass (FloatArray): Array of the atomic masses of the elements in the fluid_system
+    """
+
+    def __init__(self, species: list[str] | str, t_min: int = 250, t_max: int = 2000):
+        """Instantiates a fluid_system.
+
+        Args:
+            species: List of species names to be available in the constructed
+                fluid_system (as list of strings or a comma separated string)
+            t_min: Lower bound of the required temperature range in Kelvin
+            t_max: Upper bound of the required temperature range in Kelvin
+        """
+        if isinstance(species, str):
+            species = [s.strip() for s in species.split(',')]
+
+        temperature_base_points = range(int(t_min), ceil(t_max))
+
+        data_shape = (len(temperature_base_points), len(species))
+        self._cp_array = np.zeros(data_shape)
+        self._h_array = np.zeros(data_shape)
+        self._s_array = np.zeros(data_shape)
+        # self._g_array = np.zeros(data_shape)
+        self._g_rt_array = np.zeros(data_shape)
+
+        self._t_offset = int(t_min)
+        self.species = species
+        self.active_species = species
+        element_compositions: list[dict[str, int]] = list()
+
+        for i, s in enumerate(species):
+            species_data = _species_db.read(s)
+            if not species_data:
+                raise Exception(f'Species {s} not found')
+            element_compositions.append(species_data.composition)
+
+            assert species_data.model == 9, 'Only NASA9 polynomials are supported'
+
+            for t1, t2, a in zip(species_data.t_range[:-1], species_data.t_range[1:], species_data.data):
+
+                for j, T in enumerate(temperature_base_points):
+                    if t2 >= T >= t1:
+                        self._cp_array[j, i] = R * (a[0]*T**-2 + a[1]*T**-1 + a[2] + a[3]*T
+                                                    + a[4]*T**2 + a[5]*T**3 + a[6]*T**4)
+                        self._h_array[j, i] = R*T * (-a[0]*T**-2 + a[1]*ln(T)/T + a[2]
+                                                     + a[3]/2*T + a[4]/3*T**2 + a[5]/4*T**3
+                                                     + a[6]/5*T**4 + a[7]/T)
+                        self._s_array[j, i] = R * (-a[0]/2*T**-2 - a[1]*T**-1 + a[2]*ln(T)
+                                                   + a[3]*T + a[4]/2*T**2 + a[5]/3*T**3
+                                                   + a[6]/4*T**4 + a[8])
+                        #self._g_array[j, i] = self._h_array[j, i] - self._s_array[j, i] * T
+                        self._g_rt_array[j, i] = (self._h_array[j, i] / T - self._s_array[j, i]) / R
+
+            # TODO: Check if temperature range is not available
+            # print(f'Warning: temperature ({T}) out of range for {s}')
+
+        self.elements: list[str] = sorted(list(set(k for ac in element_compositions for k in ac.keys())))
+        self.array_species_elements: FloatArray = np.array([[ec[el] if el in ec else 0.0 for el in self.elements] for ec in element_compositions])
+
+        self.array_atomic_mass: FloatArray = np.array([atomic_weights[el] for el in self.elements]) * 1e-3  # kg/mol
+        self.array_molar_mass: FloatArray = np.sum(self.array_atomic_mass * self.array_species_elements, axis=-1)  # kg/mol
+
+        self.array_stoichiometric_coefficients: FloatArray = np.array(null_space(self.array_species_elements.T), dtype=NDFloat).T
+
+    def get_species_h(self, t: float | FloatArray) -> FloatArray:
+        """Get the molar enthalpies for all species in the fluid system
+
+        Args:
+            t: Temperature in Kelvin (can be an array)
+
+        Returns:
+            Array with the enthalpies of each specie in J/mol
+        """
+        return lookup(self._h_array, t, self._t_offset)
+
+    def get_species_s(self, t: float | FloatArray) -> FloatArray:
+        """Get the molar entropies for all species in the fluid system
+
+        Args:
+            t: Temperature in Kelvin (can be an array)
+
+        Returns:
+            Array with the entropies of each specie in J/mol/K
+        """
+        return lookup(self._s_array, t, self._t_offset)
+
+    def get_species_cp(self, t: float | FloatArray) -> FloatArray:
+        """Get the isobaric molar heat capacity for all species in the fluid system
+
+        Args:
+            t: Temperature in Kelvin (can be an array)
+
+        Returns:
+            Array with the heat capacities of each specie in J/mol/K
+        """
+        return lookup(self._cp_array, t, self._t_offset)
+
+    # def get_species_g(self, t: float | NDArray[_Float]) -> NDArray[_Float]:
+    #     return lookup(self._g_array, t, self._t_offset)
+
+    def get_species_g_rt(self, t: float | FloatArray) -> FloatArray:
+        """Get specific gibbs free energy divided by RT for all species in the
+        fluid system (g/R/T == (h/T-s)/R )
+
+        Args:
+            t: Temperature in Kelvin (can be an array)
+
+        Returns:
+            Array of gibbs free energy divided by RT (dimensionless)
+        """
+        return lookup(self._g_rt_array, t, self._t_offset)
+
+    def get_species_references(self) -> str:
+        """Get a string with the references for all fluids of the fluid system
+
+        Returns:
+            String with the references
+        """
+        return '\n'.join([f'{s:<12}: {_species_db[s].ref_string}' for s in self.species])
+
+    def __add__(self, other: 'fluid_system') -> 'fluid_system':
+        assert isinstance(other, self.__class__)
+        return self.__class__(self.species + other.species)
+
+    def __repr__(self) -> str:
+        return ('Fluid system\n    Species:  ' + ', '.join(self.species) +
+                '\n    Elements: ' + ', '.join(self.elements))
+
+
+class fluid:
+    """A class to represent a fluid defined by a composition of
+    one or more species.
+
+    Attributes:
+        species (list[str]): List of species names in the associated fluid_system
+        array_composition (FloatArray): Array of the molar amounts of the species in the fluid
+        array_element_composition (FloatArray): Array of the element composition in the fluid
+        array_fractions (FloatArray): Array of the molar fractions of the species in the fluid
+        total (FloatArray | float): Array of the sums of the molar amount of all species
+        fs (fluid_system): Reference to the fluid_system used for this fluid
+        shape (_Shape): Shape of the fluid array
+        elements (list[str]): List of elements in the fluid_system
+    """
+
+    __array_priority__ = 100
+
+    def __init__(self, composition: dict[str, float] | list[float] | FloatArray,
+                 fs: fluid_system | None = None,
+                 shape: Sequence[int] | None = None):
+        """Instantiates a fluid.
+
+        Args:
+            composition: A dict of species names with their composition, e.g.
+                {'O2':0.5,'H2O':0.5} or a list/numpy-array of compositions.
+                The array can be multidimensional, the size of the last dimension
+                must match the number of species defined for the fluid_system.
+                The indices of the last dimension correspond to the indices in
+                the active_species list of the fluid_system.
+            fs: Reference to a fluid_system. Is optional if composition is
+                defined by a dict. If not specified a new fluid_system with
+                the components from the dict is created.
+            shape: Tuple or list for the dimensions the fluid array. Can
+                only be used if composition argument is a dict. Otherwise
+                the dimensions are specified by the composition argument.
+        """
+        if fs is None:
+            assert isinstance(composition, dict), 'fluid system must be specified if composition is not a dict'
+            fs = fluid_system(list(composition.keys()))
+
+        if isinstance(composition, list):
+            composition = np.array(composition)
+
+        if isinstance(composition, dict):
+            missing_species = [s for s in composition if s not in fs.species]
+            if len(missing_species):
+                raise Exception(f'Species {missing_species[0]} is not part of the fluid system')
+
+            species_composition = [composition[s] if s in composition.keys() else 0 for s in fs.species]
+
+            comp_array = np.array(species_composition, dtype=NDFloat)
+            if shape is not None:
+                comp_array = comp_array * np.ones(list(shape) + [len(fs.species)], dtype=NDFloat)
+
+        else:
+            assert shape is None, 'specify shape by the shape of the composition array.'
+            assert composition.shape[-1] == len(fs.species), f'composition.shape[-1] ({composition.shape[-1]}) must be {len(fs.species)}'
+            comp_array = composition
+
+        self.array_composition: FloatArray = comp_array
+        self.total: FloatArray | float = np.sum(self.array_composition, axis=-1, dtype=NDFloat)
+        self.array_fractions: FloatArray = self.array_composition / (np.expand_dims(self.total, -1) + epsy)
+        self.shape: Shape = self.array_composition.shape[:-1]
+        self.fs = fs
+        self.array_elemental_composition: FloatArray = np.dot(self.array_composition, fs.array_species_elements)
+        self.species = fs.species
+        self.elements = fs.elements
+
+    def get_composition_dict(self) -> dict[str, float]:
+        """Get a dict of the molar amount of each fluid species
+
+        Returns:
+            Returns a dict of floats with the molar amount of each fluid species in mol
+        """
+        return {s: c for s, c in zip(self.fs.species, self.array_composition)}
+
+    def get_fractions_dict(self) -> dict[str, float]:
+        """Get a dict of the molar fractions of each fluid species
+
+        Returns:
+            Returns a dict of floats with the molar fractions of each fluid species
+        """
+        return {s: c for s, c in zip(self.fs.species, self.array_fractions)}
+
+    def get_h(self, t: float | FloatArray) -> FloatArray | float:
+        """Get specific enthalpy of the fluid at the given temperature
+
+        Enthalpy is referenced to 25 °C and includes enthalpy of formation.
+        Therefore the enthalpy of H2 and O2 is 0 at 25 °C, but the enthalpy
+        of water vapor at 25 °C is −241 kJ/mol (enthalpy of formation).
+
+        Args:
+            t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
+                must be broadcastable
+
+        Returns:
+            Enthalpies in J/mol
+        """
+        return np.sum(self.fs.get_species_h(t) * self.array_fractions, axis=-1, dtype=NDFloat)
+
+    def get_H(self, t: float | FloatArray) -> FloatArray | float:
+        """Get absolute enthalpy of the fluid at the given temperature
+
+        Enthalpy is referenced to 25 °C and includes enthalpy of formation.
+        Therefore the enthalpy of H2 and O2 is 0 at 25 °C, but the enthalpy
+        of water vapor at 25 °C is −241 kJ/mol (enthalpy of formation).
+
+        Args:
+            t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
+                must be broadcastable
+
+        Returns:
+            Enthalpies in J
+        """
+        return np.sum(self.fs.get_species_h(t) * self.array_composition, axis=-1, dtype=NDFloat)
+
+    def get_s(self, t: float | FloatArray, p: float | FloatArray) -> FloatArray | float:
+        """Get molar entropy of the fluid at the given temperature and pressure
+
+        Args:
+            t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
+                must be broadcastable
+
+        Returns:
+            Entropy in J/mol/K
+        """
+        x = self.array_fractions
+        s = self.fs.get_species_s(t)
+
+        return np.sum(x * (s - R * np.log(np.expand_dims(p / p0, -1) * x + epsy)), axis=-1, dtype=NDFloat)
+
+    def get_S(self, t: float | FloatArray, p: float | FloatArray) -> FloatArray | float:
+        """Get absolute entropy of the fluid at the given temperature and pressure
+
+        Args:
+            t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
+                must be broadcastable
+
+        Returns:
+            Entropy in J/K
+        """
+        x = self.array_fractions
+        n = self.array_composition
+        s = self.fs.get_species_s(t)
+
+        return np.sum(n * (s - R * np.log(np.expand_dims(p / p0, -1) * x + epsy)), axis=-1, dtype=NDFloat)
+
+    def get_cp(self, t: float | FloatArray) -> FloatArray | float:
+        """Get molar heat capacity at constant pressure
+
+        Args:
+            t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
+                must be broadcastable
+
+        Returns:
+            Heat capacity in J/mol/K
+        """
+        return np.sum(self.fs.get_species_cp(t) * self.array_fractions, axis=-1, dtype=NDFloat)
+
+    def get_g(self, t: float | FloatArray, p: float | FloatArray) -> FloatArray | float:
+        """Get molar gibbs free energy (h - Ts)
+
+        Args:
+            t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
+                must be broadcastable
+            p: Pressures(s) in Pascal. Fluid shape and shape of the temperature
+                must be broadcastable
+
+        Returns:
+            Gibbs free energy in J/mol
+        """
+        x = self.array_fractions
+        grt = self.fs.get_species_g_rt(t)
+
+        return R * t * np.sum(x * (grt + np.log(np.expand_dims(p / p0, -1) * x + epsy)), axis=-1, dtype=NDFloat)
+
+    def get_G(self, t: float | FloatArray, p: float | FloatArray) -> FloatArray | float:
+        """Get absolute gibbs free energy (H - TS)
+
+        Args:
+            t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
+                must be broadcastable
+            p: Pressures(s) in Pascal. Fluid shape and shape of the temperature
+                must be broadcastable
+
+        Returns:
+            Gibbs free energy in J
+        """
+        x = self.array_fractions
+        n = self.array_composition
+        grt = self.fs.get_species_g_rt(t)
+
+        return R * t * np.sum(n * (grt + np.log(np.expand_dims(p / p0, -1) * x + epsy)), axis=-1, dtype=NDFloat)
+
+    def get_g_rt(self, t: float | FloatArray, p: float | FloatArray) -> FloatArray | float:
+        """Get specific gibbs free energy divided by RT: g/R/T == (h/T-s)/R
+
+        Args:
+            t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
+                must be broadcastable
+            p: Pressures(s) in Pascal. Fluid shape and shape of the temperature
+                must be broadcastable
+
+        Returns:
+            Gibbs free energy divided by RT (dimensionless)
+        """
+        x = self.array_fractions
+        grt = self.fs.get_species_g_rt(t)
+
+        return np.sum(x * (grt + np.log(np.expand_dims(p / p0, -1) * x + epsy)), axis=-1, dtype=NDFloat)
+
+    def get_v(self, t: float | FloatArray, p: float | FloatArray) -> FloatArray | float:
+        """Get Absolute fluid volume
+
+        Args:
+            t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
+                must be broadcastable
+            p: Pressure in Pa. Fluid shape and shape of the pressure
+                must be broadcastable
+
+        Returns:
+            Volume of the fluid in m³
+        """
+        return R / p * t * self.total
+
+    def get_vm(self, t: float | FloatArray, p: float | FloatArray) -> FloatArray | float:
+        """Get molar fluid volume
+
+        Args:
+            t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
+                must be broadcastable
+            p: Pressure in Pa. Fluid shape and shape of the pressure
+                must be broadcastable
+
+        Returns:
+            Molar volume of the fluid in m³/mol
+        """
+        return R / p * t
+
+    def get_mass(self) -> FloatArray | float:
+        """Get Absolute fluid mass
+
+        Returns:
+            Mass of the fluid in kg
+        """
+        return np.sum(self.array_composition * self.fs.array_molar_mass, axis=-1, dtype=NDFloat)
+
+    def get_molar_mass(self) -> FloatArray | float:
+        """Get molar fluid mass
+
+        Returns:
+            Mass of the fluid in kg/mol
+        """
+        return np.sum(self.array_fractions * self.fs.array_molar_mass, axis=-1, dtype=NDFloat)
+
+    def get_density(self, t: float | FloatArray, p: float | FloatArray) -> FloatArray | float:
+        """Get mass based fluid density
+
+        Args:
+            t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
+                must be broadcastable
+            p: Pressure in Pa. Fluid shape and shape of the pressure
+                must be broadcastable
+
+        Returns:
+            Density of the fluid in kg/m³
+        """
+        return np.sum(self.array_fractions * self.fs.array_molar_mass, axis=-1, dtype=NDFloat) / (R * t) * p
+
+    def get_x(self, species: str | list[str] | None = None) -> FloatArray:
+        """Get molar fractions of fluid species
+
+        Args:
+            species: A single species name, a list of species names or None for
+                returning the molar fractions of all species
+
+        Returns:
+            Returns an array of floats with the molar fractions of the species.
+            If the a single species name is provided the return float array has
+            the same dimensions as the fluid type. If a list or None is provided
+            the return array has an additional dimension for the species.
+        """
+        if not species:
+            return self.array_fractions
+        elif isinstance(species, str):
+            assert species in self.fs.species, f'Species {species} is not part of the fluid system'
+            return self.array_fractions[..., self.fs.species.index(species)]
+        else:
+            assert set(species) <= set(self.fs.species), f'Species {", ".join([s for s in species if s not in self.fs.species])} is/are not part of the fluid system'
+            return self.array_fractions[..., [self.fs.species.index(k) for k in species]]
+
+    def get_n(self, species: str | list[str] | None = None) -> FloatArray:
+        """Get molar amount of fluid species
+
+        Args:
+            species: A single species name, a list of species names or None for
+                returning the amount of all species
+
+        Returns:
+            Returns an array of floats with the molar amount of the species.
+            If the a single species name is provided the return float array has
+            the same dimensions as the fluid type. If a list or None is provided
+            the return array has an additional dimension for the species.
+        """
+        if not species:
+            return self.array_composition
+        elif isinstance(species, str):
+            assert species in self.fs.species, f'Species {species} is not part of the fluid system'
+            return self.array_composition[..., self.fs.species.index(species)]
+        else:
+            assert set(species) <= set(self.fs.species), f'Species {", ".join([s for s in species if s not in self.fs.species])} is/are not part of the fluid system'
+            return self.array_composition[..., [self.fs.species.index(k) for k in species]]
+
+    def __add__(self, other: T) -> T:
+        return array_operation(self, other, np.add)
+
+    def __sub__(self, other: T) -> T:
+        return array_operation(self, other, np.subtract)
+
+    def __truediv__(self, other: int | float | NDArray[Any]) -> 'fluid':
+        if isinstance(other, np.ndarray):
+            k = np.expand_dims(other, -1)
+        else:
+            k = np.array(other, dtype=NDFloat)
+        return self.__class__(self.array_composition / k, self.fs)
+
+    def __mul__(self, other: int | float | NDArray[Any]) -> 'fluid':
+        k = np.expand_dims(other, -1) if isinstance(other, np.ndarray) else other
+        return self.__class__(self.array_composition * k, self.fs)
+
+    def __rmul__(self, other: int | float | NDArray[Any]) -> 'fluid':
+        k = np.expand_dims(other, -1) if isinstance(other, np.ndarray) else other
+        return self.__class__(self.array_composition * k, self.fs)
+
+    def __neg__(self) -> 'fluid':
+        return self.__class__(-self.array_composition, self.fs)
+
+    # def __array__(self) -> FloatArray:
+    #     return self.array_composition
+
+    @overload
+    def __getitem__(self, key: str) -> FloatArray:
+        pass
+
+    @overload
+    def __getitem__(self, key: ArrayIndices) -> 'fluid':
+        pass
+
+    def __getitem__(self, key: str | ArrayIndices) -> Any:
+        if isinstance(key, str):
+            assert key in self.fs.species, f'Species {key} is not part of the fluid system'
+            return self.array_composition[..., self.fs.species.index(key)]
+        else:
+            key_tuple = key if isinstance(key, tuple) else (key,)
+            return fluid(self.array_composition[(*key_tuple, slice(None))], self.fs)
+
+    def __iter__(self) -> Iterator[dict[str, float]]:
+        assert len(self.shape) < 2, 'Cannot iterate over species with more than one dimension'
+        aec = self.array_composition.reshape(-1, len(self.fs.species))
+        return iter({s: c for s, c in zip(self.fs.species, aec[i, :])} for i in range(aec.shape[0]))
+
+    def __repr__(self) -> str:
+        if len(self.array_fractions.shape) == 1:
+            lines = [f'{s:16} {c * 100:5.2f} %' for s, c in zip(self.fs.species, self.array_fractions)]
+            return f'{"Total":16} {self.total:8.3e} mol\n' + '\n'.join(lines)
+        else:
+            array_disp = self.array_fractions.__repr__()
+            padding = int(array_disp.find('\n') / (len(self.fs.species) + 1))
+            return ('Total mol:\n' + self.total.__repr__() +
+                    '\nSpecies:\n' + ' ' * int(padding / 2) +
+                    ''.join([(' ' * (padding - len(s))) + s for s in self.fs.species]) +
+                    '\nMolar fractions:\n' + self.array_fractions.__repr__())
+
+
+class elements:
+    """Represent a fluid by composition of elements.
+
+    Attributes:
+        array_element_composition (FloatArray): Array of the element composition
+    """
+
+    __array_priority__ = 100
+
+    def __init__(self, composition: fluid | dict[str, float] | list[str] | list[float] | FloatArray,
+                 fs: fluid_system | None = None, shape: Sequence[int] | None = None):
+        """Instantiates an elements object.
+
+        Args:
+            composition: A fluid object, a dict of element names with their
+                composition, e.g.: {'O':1,'H':2} or a list/numpy-array of compositions.
+                The array can be multidimensional, the size of the last dimension
+                must match the number of elements used in the fluid_system.
+                The indices of the last dimension correspond to the indices in
+                the active_species list of the fluid_system.
+            fs: Reference to a fluid_system.
+            shape: Tuple or list for the dimensions the fluid array. Can
+                only be used if composition argument is a dict. Otherwise
+                the dimensions are specified by the composition argument.
+        """
+        if isinstance(composition, list):
+            composition = np.array(composition)
+
+        if isinstance(composition, fluid):
+            new_composition: FloatArray = np.dot(composition.array_composition, composition.fs.array_species_elements)
+            if fs:
+                self.array_elemental_composition = reorder_array(new_composition, composition.fs.elements, fs.elements)
+            else:
+                self.array_elemental_composition = new_composition
+                fs = composition.fs
+        elif isinstance(composition, dict) and fs is None:
+            fs = fluid_system(species(element_names=list(composition.keys())))
+        else:
+            assert fs, 'fluid system must be specified if composition is not specified by a fluid'
+
+        if isinstance(composition, dict):
+            missing_elements = [s for s in composition if s not in fs.elements]
+            if len(missing_elements):
+                raise Exception(f'Element {missing_elements[0]} is not part of the fluid system')
+
+            self.array_elemental_composition = np.array([composition[s] if s in composition.keys() else 0 for s in fs.elements])
+
+            if shape is not None:
+                self.array_elemental_composition = self.array_elemental_composition * np.ones(list(shape) + [len(fs.species)])
+
+        elif isinstance(composition, np.ndarray):
+            assert shape is None, 'specify shape by the shape of the composition array.'
+            assert composition.shape[-1] == len(fs.elements), f'composition.shape[-1] ({composition.shape[-1]}) must be {len(fs.elements)}'
+            self.array_elemental_composition = composition
+
+        self.shape: Shape = self.array_elemental_composition.shape[:-1]
+        self.fs = fs
+        self.elements = fs.elements
+
+    def get_elemental_composition(self) -> dict[str, float]:
+        """Get a dict of the molar amount of each element
+
+        Returns:
+            Returns a dict of floats with the molar amount of each element in mol
+        """
+        return {s: c for s, c in zip(self.fs.elements, self.array_elemental_composition)}
+
+    def get_mass(self) -> FloatArray | float:
+        """Get absolute mass of elements
+
+        Returns:
+            Mass of the fluid in kg
+        """
+        return np.sum(self.array_elemental_composition * self.fs.array_atomic_mass, axis=-1, dtype=NDFloat)
+
+    def get_n(self, elemental_species: str | list[str] | None = None) -> FloatArray:
+        """Get molar amount of elements
+
+        Args:
+            elemental_species: A single element name, a list of element names or None for
+                returning the amount of all element
+
+        Returns:
+            Returns an array of floats with the molar amount of the elements.
+            If the a single element name is provided the return float array has
+            the same dimensions as the fluid type. If a list or None is provided
+            the return array has an additional dimension for the elements.
+        """
+        if not elemental_species:
+            return self.array_elemental_composition
+        elif isinstance(elemental_species, str):
+            assert elemental_species in self.fs.elements, f'Element {elemental_species} is not part of the fluid system'
+            return self.array_elemental_composition[..., self.fs.elements.index(elemental_species)]
+        else:
+            assert set(elemental_species) <= set(self.fs.elements), f'Elements {", ".join([s for s in elemental_species if s not in self.fs.elements])} is/are not part of the fluid system'
+            return self.array_elemental_composition[..., [self.fs.elements.index(k) for k in elemental_species]]
+
+    def __add__(self, other: 'fluid | elements') -> 'elements':
+        return array_operation(self, other, np.add)
+
+    def __sub__(self, other: 'fluid | elements') -> 'elements':
+        return array_operation(self, other, np.subtract)
+
+    def __truediv__(self, other: int | float | FloatArray) -> 'elements':
+        k = np.expand_dims(other, -1) if isinstance(other, np.ndarray) else other
+        ttes = self.array_elemental_composition / k
+        return self.__class__(self.array_elemental_composition / k + ttes, self.fs)
+
+    def __mul__(self, other: int | float | FloatArray) -> 'elements':
+        k = np.expand_dims(other, -1) if isinstance(other, np.ndarray) else other
+        return self.__class__(self.array_elemental_composition * k, self.fs)
+
+    def __rmul__(self, other: int | float | FloatArray) -> 'elements':
+        k = np.expand_dims(other, -1) if isinstance(other, np.ndarray) else other
+        return self.__class__(self.array_elemental_composition * k, self.fs)
+
+    def __neg__(self) -> 'elements':
+        return self.__class__(-self.array_elemental_composition, self.fs)
+
+    def __array__(self) -> FloatArray:
+        return self.array_elemental_composition
+
+    @overload
+    def __getitem__(self, key: str) -> FloatArray:
+        pass
+
+    @overload
+    def __getitem__(self, key: ArrayIndices) -> 'elements':
+        pass
+
+    def __getitem__(self, key: str | ArrayIndices) -> Any:
+        if isinstance(key, str):
+            assert key in self.fs.elements, f'Element {key} is not part of the fluid system'
+            return self.array_elemental_composition[..., self.fs.elements.index(key)]
+        else:
+            key_tuple = key if isinstance(key, tuple) else (key,)
+            return elements(self.array_elemental_composition[(*key_tuple, slice(None))], self.fs)
+
+    def __iter__(self) -> Iterator[dict[str, float]]:
+        assert len(self.shape) < 2, 'Cannot iterate over elements with more than one dimension'
+        aec = self.array_elemental_composition.reshape(-1, len(self.fs.elements))
+        return iter({s: c for s, c in zip(self.fs.elements, aec[i, :])} for i in range(aec.shape[0]))
+
+    def __repr__(self) -> str:
+        if len(self.array_elemental_composition.shape) == 1:
+            lines = [f'{s:16} {c:5.3e} mol' for s, c in zip(self.fs.elements, self.array_elemental_composition)]
+            return '\n'.join(lines)
+        else:
+            array_disp = self.array_elemental_composition.__repr__()
+            padding = int(array_disp.find('\n') / (len(self.fs.elements) + 1))
+            return ('Elements:\n' + ' ' * int(padding / 2) +
+                    ''.join([(' ' * (padding - len(s))) + s for s in self.fs.elements]) +
+                    '\nMols:\n' + self.array_elemental_composition.__repr__())
+
+
+def lookup(prop_array: FloatArray,
+           temperature: FloatArray | float,
+           t_offset: float) -> FloatArray:
+    """linear interpolates values from the given prop_array
+
+    Args:
+        prop_array: Array of the temperature depended property
+        temperature: Absolute temperature(s) in Kelvin. Must
+            be broadcastable to prop_array.
+
+    Returns:
+        Interpolates values based on given temperature
+    """
+    t = np.array(temperature) - t_offset
+    t_lim = np.minimum(np.maximum(0, t), prop_array.shape[0] - 2)
+
+    f = np.expand_dims(t - np.floor(t_lim), axis=-1)
+
+    ti1 = t_lim.astype(int)
+    return f * prop_array[ti1 + 1, :] + (1 - f) * prop_array[ti1, :]
+
+
+def reorder_array(arr: FloatArray, old_index: list[str], new_index: list[str]) -> FloatArray:
+    """Reorder the last dimension of an array according to a provided list of species
+    names in the old oder and a list in the new order.
+
+    Args:
+        arr: Array to be reordered
+        old_index: List of species names in the current order
+        new_index: List of species names in the new order
+
+    Returns:
+        Array with the last dimension reordered
+    """
+    ret_array = np.zeros([*arr.shape[:-1], len(new_index)])
+    for i, k in enumerate(old_index):
+        ret_array[..., new_index.index(k)] = arr[..., i]
+    return ret_array
+
+
+@overload
+def array_operation(self: elements, other: elements | fluid, func: Callable[[FloatArray, FloatArray], FloatArray]) -> elements:
+    pass
+
+
+@overload
+def array_operation(self: fluid, other: T, func: Callable[[FloatArray, FloatArray], FloatArray]) -> T:
+    pass
+
+
+@overload
+def array_operation(self: T, other: fluid, func: Callable[[FloatArray, FloatArray], FloatArray]) -> T:
+    pass
+
+
+def array_operation(self: elements | fluid, other: elements | fluid, func: Callable[[FloatArray, FloatArray], FloatArray]) -> elements | fluid:
+    """Perform an array operation on two fluid or elements objects.
+    The operation is provided by a Callable that takes two arguments.
+
+    Args:
+        self: First fluid or elements object
+        other: Second fluid or elements object
+        func: Callable function to perform the operation
+
+    Returns:
+        A new fluid or elements object with the result of the
+    """
+    assert isinstance(other, elements) or isinstance(other, fluid)
+    if self.fs is other.fs:
+        if isinstance(self, elements) or isinstance(other, elements):
+            return elements(func(self.array_elemental_composition, other.array_elemental_composition), self.fs)
+        else:
+            return fluid(func(self.array_composition, other.array_composition), self.fs)
+    elif set(self.fs.species) >= set(other.fs.species):
+        if isinstance(self, elements) or isinstance(other, elements):
+            el_array = reorder_array(other.array_elemental_composition, other.fs.elements, self.fs.elements)
+            return elements(func(self.array_elemental_composition, el_array), self.fs)
+        else:
+            el_array = reorder_array(other.array_composition, other.fs.species, self.fs.species)
+            return fluid(func(self.array_composition, el_array), self.fs)
+    elif set(self.fs.species) < set(other.fs.species):
+        if isinstance(self, elements) or isinstance(other, elements):
+            el_array = reorder_array(self.array_elemental_composition, self.fs.elements, other.fs.elements)
+            return elements(func(el_array, other.array_elemental_composition), other.fs)
+        else:
+            el_array = reorder_array(self.array_composition, self.fs.species, other.fs.species)
+            return fluid(func(el_array, other.array_composition), other.fs)
+    else:
+        new_fs = fluid_system(sorted(list(set(self.fs.species) | set(other.fs.species))))
+        if isinstance(self, elements) or isinstance(other, elements):
+            el_array1 = reorder_array(self.array_elemental_composition, self.fs.elements, new_fs.elements)
+            el_array2 = reorder_array(other.array_elemental_composition, other.fs.elements, new_fs.elements)
+            return elements(func(el_array1, el_array2), new_fs)
+        else:
+            el_array1 = reorder_array(self.array_composition, self.fs.species, new_fs.species)
+            el_array2 = reorder_array(other.array_composition, other.fs.species, new_fs.species)
+            return fluid(func(el_array1, el_array2), new_fs)
--- a/src/gaspype/_operations.py
+++ b/src/gaspype/_operations.py
@ -0,0 +1,150 @@
+from math import exp
+import numpy as np
+from ._main import T, elements, fluid
+from .typing import FloatArray
+from .constants import p0, R
+from ._solver import equilibrium
+
+
+def stack(arrays: list[T], axis: int = 0) -> T:
+    """Stack a list of fluid or elements objects along a new axis
+
+    Args:
+        arrays: List of arrays
+        axis: Axis to stack the fluid objects along
+
+    Returns:
+        A new array object stacked along the new axis
+    """
+    a0 = arrays[0]
+    assert all(a.fs == a0.fs for a in arrays), 'All objects must have the same fluid system'
+    assert axis <= len(a0.shape), f'Axis must be smaller or equal to len(shape) ({len(a0.shape)})'
+    return a0.__class__(np.stack(
+        [a.array_elemental_composition if isinstance(a, elements) else a.array_composition for a in arrays],
+        axis=axis), a0.fs)
+
+
+def concat(arrays: list[T], axis: int = 0) -> T:
+    """Concatenate a list of fluid or elements objects along an existing axis
+
+    Args:
+        arrays: List of arrays
+        axis: Axis to concatenate the fluid objects along
+
+    Returns:
+        A new array object stacked along the specified axis
+    """
+    a0 = arrays[0]
+    assert all(f.fs == a0.fs for f in arrays), 'All fluid objects must have the same fluid system'
+    assert axis < len(a0.shape), f'Axis must be smaller than shape len({a0.shape})'
+    return a0.__class__(np.concatenate(
+        [a.array_elemental_composition if isinstance(a, elements) else a.array_composition for a in arrays],
+        axis=axis), a0.fs)
+
+
+def carbon_activity(f: fluid | elements, t: float, p: float) -> float:
+    """Calculate the activity of carbon in a fluid at a given temperature and pressure.
+    At a value of 1 the fluid is in equilibrium with solid graphite. At a value > 1
+    additional carbon formation is thermodynamic favored. At a value < 1 a
+    depletion of solid carbon is favored.
+
+    Args:
+        f: Fluid or elements object
+        t: Temperature in Kelvin
+        p: Pressure in Pascal
+
+    Returns:
+        The activity of carbon in the fluid
+    """
+    # Values for solid state carbon (graphite) from NIST-JANAF Tables
+    # https://janaf.nist.gov/pdf/JANAF-FourthEd-1998-Carbon.pdf
+    # https://janaf.nist.gov/pdf/JANAF-FourthEd-1998-1Vol1-Intro.pdf
+    # Polynomial is valid for T from 100 to 2500 K
+    pgef = np.array([-6.76113852E-02, 2.02187857E+00, -2.38664605E+01,
+                    1.43575462E+02, -4.51375503E+02, 6.06219665E+02])
+
+    # Gibbs free energy divided by RT for carbon
+    g_rtc = -np.sum(pgef * np.log(np.expand_dims(t, -1))**np.array([5, 4, 3, 2, 1, 0])) / R
+
+    g_rt = f.fs.get_species_g_rt(t)
+
+    x = equilibrium(f, t, p).array_fractions
+
+    i_co = f.fs.species.index('CO')
+    i_co2 = f.fs.species.index('CO2')
+    i_h2 = f.fs.species.index('H2')
+    i_h2o = f.fs.species.index('H2O')
+    i_ch4 = f.fs.species.index('CH4')
+
+    if min(x[i_co], x[i_co2]) > min([x[i_ch4], x[i_h2o], x[i_h2]]) and min(x[i_co], x[i_co2]) > 0:
+        # 2 CO -> CO2 + C(s) (Boudouard reaction)
+        lnalpha = (2 * g_rt[i_co] - (g_rt[i_co2] + g_rtc)) + np.log(
+            x[i_co]**2 / x[i_co2] * (p / p0))
+    elif min([x[i_ch4], x[i_h2o], x[i_co]]) > 1E-4:
+        # CH4 + 2 CO -> 2 H2O + 3 C(s)
+        lnalpha = ((g_rt[i_ch4] + 2 * g_rt[i_co] - 3 * g_rtc - 2 * g_rt[i_h2o]) + np.log(
+            x[i_ch4] * x[i_co]**2 / x[i_h2o]**2 * (p / p0))) / 3
+    elif min(x[i_h2], x[i_ch4]) > 0:
+        # if x[i_h2] or x[i_ch4] is small compared to the precision of the
+        # component concentrations the result can be inaccurate
+        # CH4 -> 2 H2 + C(s)
+        # CH4 + CO2 -> 2 H2 + 2 CO
+        # 2 H2O - O2 -> 2 H2
+        lnalpha = (g_rt[i_ch4] - (2 * g_rt[i_h2] + g_rtc)) + np.log(
+            x[i_ch4] / x[i_h2]**2 / (p / p0))
+    elif x[i_h2] == 0:
+        # equilibrium on carbon side
+        lnalpha = 10
+    else:
+        # equilibrium on non-carbon side
+        lnalpha = -10
+
+    return exp(lnalpha)
+
+
+def oxygen_partial_pressure(f: fluid | elements, t: float, p: float) -> FloatArray | float:
+    _oxygen_data = fluid({'O2': 1})
+
+    def get_oxygen(x: FloatArray) -> float:
+        g_rt = f.fs.get_species_g_rt(t)
+        g_rt_o2 = _oxygen_data.fs.get_species_g_rt(t)[0]
+
+        i_co = f.fs.species.index('CO') if 'C' in f.fs.elements else None
+        i_co2 = f.fs.species.index('CO2') if 'C' in f.fs.elements else None
+        i_o2 = f.fs.species.index('O2') if 'O2' in f.fs.species else None
+        i_h2o = f.fs.species.index('H2O') if 'H' in f.fs.elements else None
+        i_h2 = f.fs.species.index('H2') if 'H' in f.fs.elements else None
+        i_ch4 = f.fs.species.index('CH4') if 'CH4' in f.fs.species else None
+
+        ox_ref_val = max([float(v) for v in (x[i_h2] * x[i_h2o], x[i_co2] * x[i_co], x[i_ch4]) if v.shape == tuple()] + [0])
+
+        # print([i_o2, x[i_o2], ox_ref_val])
+
+        if i_o2 is not None and x[i_o2] > ox_ref_val:
+            # print('o O2')
+            return float(x[i_o2] * p)
+
+        elif i_ch4 is not None and x[i_ch4] > x[i_co2] * 100 and x[i_ch4] > x[i_h2o] * 100:
+            # print('o ch4')
+            # 2CH4 + O2 <--> 4H2 + 2CO
+            lnpo2 = 4 * g_rt[i_h2] + 2 * g_rt[i_co] - 2 * g_rt[i_ch4] - g_rt_o2 + np.log(x[i_h2]**4 * x[i_co]**2 / x[i_ch4]**2) - 2 * np.log(p / p0)
+
+        elif (i_co is None and i_h2 is not None) or (i_h2 is not None and i_co is not None and (x[i_h2] * x[i_h2o] > x[i_co2] * x[i_co])):
+            # print('o h2o')
+            # 2H2 + O2 <--> 2H2O
+            lnpo2 = 2 * (g_rt[i_h2o] - g_rt[i_h2] + np.log(x[i_h2o] / x[i_h2])) - g_rt_o2 - np.log(p / p0)
+
+        else:
+            assert i_co is not None
+            # print('o co2')
+            # 2CO + O2 <--> 2CO2
+            lnpo2 = 2 * (g_rt[i_co2] - g_rt[i_co] + np.log(x[i_co2] / x[i_co])) - g_rt_o2 - np.log(p / p0)
+
+        return exp(lnpo2) * p
+
+    x = equilibrium(f, t, p).array_fractions
+
+    if len(x.shape):
+        return np.apply_along_axis(get_oxygen, -1, x)
+    else:
+        return get_oxygen(x)
--- a/src/gaspype/_phys_data.py
+++ b/src/gaspype/_phys_data.py
@ -1,32 +1,44 @@
 import struct
 from typing import Generator, Iterator
+from dataclasses import dataclass


-class SpeciesData():
+def split_on_space(data: bytes, offset: int, end: int) -> Iterator[str]:
+    """Splits a byte array into ASCII strings based on spaces.
+    Args:
+        data: The byte array to split.
+        offset: The starting index in the byte array.
+        end: The ending index in the byte array.
+    Yields:
+        str: ASCII strings found in the byte array.
+    """
+    start = offset
+    for i in range(offset, end):
+        if data[i] == 0x20:  # ASCII space character
+            if start < i:
+                yield data[start:i].decode('ascii')
+            start = i + 1
+    if start < end:
+        yield data[start:end].decode('ascii')
+
+
+@dataclass
+class SpeciesData:
    """Class to hold the physical data for a species.
    Attributes:
-        comp: Dictionary of species composition with element symbols as keys and their counts as values.
+        name: Name of the species.
+        composition: Dictionary of species composition with element symbols as keys and their counts as values.
        model: Number of polynomial coefficients used in the model.
        ref_string: Reference string for the data source.
        t_range: List of temperatures nodes marking intervals.
        data: List of lists containing physical data for each temperature interval.
    """
-
-    def __init__(self, name: str, comp: dict[str, int], model: int, ref: str, t_range: list[float], data: list[list[float]]):
-        self.name = name
-        self.composition: dict[str, int] = comp
-        self.model: int = model
-        self.ref_string: str = ref
-        self.t_range: list[float] = t_range
-        self.data: list[list[float]] = data
-
-    def __repr__(self) -> str:
-        return (f"Name: {self.name}\n" +
-                f"Composition: {self.composition}\n" +
-                f"Model: {self.model}\n" +
-                f"Reference: {self.ref_string}\n" +
-                f"Temperatures: {self.t_range}\n" +
-                f"Data: {self.data}".replace('),', '),\n'))
+    name: str
+    composition: dict[str, int]
+    model: int
+    ref_string: str
+    t_range: list[float]
+    data: list[list[float]]


 class db_reader():
@ -45,8 +57,8 @@ class db_reader():
        """
        assert inp_data[:4] == b'gapy', 'Unknown data format'
        self._bin_data = inp_data
-        self._name_count = struct.unpack('<I', self._bin_data[4:8])[0]
-        species_names = self._bin_data[db_reader.header_len:(db_reader.header_len + self._name_count)].decode('ASCII').split(' ')
+        self._name_lengths = struct.unpack('<I', self._bin_data[4:8])[0]
+        species_names = split_on_space(self._bin_data, db_reader.header_len, db_reader.header_len + self._name_lengths)
        self._index = {s: i for i, s in enumerate(species_names)}

    @property
@ -82,7 +94,7 @@ class db_reader():
        if name not in self._index:
            return None

-        head_offset = self._name_count + db_reader.header_len + self._index[name] * db_reader.header_len
+        head_offset = self._name_lengths + db_reader.header_len + self._index[name] * db_reader.header_len

        head = struct.unpack('<I4B', self._bin_data[head_offset:head_offset + db_reader.header_len])

--- a/src/gaspype/_solver.py
+++ b/src/gaspype/_solver.py
@ -0,0 +1,167 @@
+from typing import Literal, Any
+from scipy.optimize import minimize, root
+import numpy as np
+from ._main import elements, fluid, fluid_system
+from .typing import NDFloat, FloatArray
+from .constants import p0, epsy
+
+
+def set_solver(solver: Literal['gibs minimization', 'system of equations']) -> None:
+    """
+    Select a solver for chemical equilibrium.
+
+    Solvers:
+        - **system of equations** (default): Finds the root for a system of
+          equations covering a minimal set of equilibrium equations and elemental balance.
+          The minimal set of equilibrium equations is derived by SVD using the null_space
+          implementation of scipy.
+
+        - **gibs minimization**: Minimizes the total Gibbs Enthalpy while keeping
+          the elemental composition constant using the SLSQP implementation of scipy
+
+    Args:
+        solver: Name of the solver
+    """
+    global _equilibrium_solver
+    if solver == 'gibs minimization':
+        _equilibrium_solver = equilibrium_gmin
+    elif solver == 'system of equations':
+        _equilibrium_solver = equilibrium_eq
+    else:
+        raise ValueError('Unknown solver')
+
+
+def get_solver() -> Literal['gibs minimization', 'system of equations']:
+    """Returns the selected solver name.
+
+    Returns:
+        Solver name
+    """
+    if _equilibrium_solver == equilibrium_gmin:
+        return 'gibs minimization'
+    else:
+        assert _equilibrium_solver == equilibrium_eq
+        return 'system of equations'
+
+
+def equilibrium_gmin(fs: fluid_system, element_composition: FloatArray, t: float, p: float) -> FloatArray:
+    """Calculate the equilibrium composition of a fluid based on minimizing the Gibbs free energy"""
+    def element_balance(n: FloatArray, fs: fluid_system, ref: FloatArray) -> FloatArray:
+        return np.dot(n, fs.array_species_elements) - ref  # type: ignore
+
+    def gibbs_rt(n: FloatArray, grt: FloatArray, p_rel: float):  # type: ignore
+        # Calculate G/(R*T)
+        return np.sum(n * (grt + np.log(p_rel * n / np.sum(n) + epsy)))
+
+    cons: dict[str, Any] = {'type': 'eq', 'fun': element_balance, 'args': [fs, element_composition]}
+    bnds = [(0, None) for _ in fs.species]
+    grt = fs.get_species_g_rt(t)
+    p_rel = p / p0
+
+    start_composition_array = np.ones_like(fs.species, dtype=float)
+    sol = np.array(minimize(gibbs_rt, start_composition_array, args=(grt, p_rel), method='SLSQP',
+                   bounds=bnds, constraints=cons, options={'maxiter': 2000, 'ftol': 1e-12})['x'], dtype=NDFloat)  # type: ignore
+
+    return sol
+
+
+def equilibrium_eq(fs: fluid_system, element_composition: FloatArray, t: float, p: float) -> FloatArray:
+    """Calculate the equilibrium composition of a fluid based on equilibrium equations"""
+    el_max = np.max(element_composition)
+    element_norm = element_composition / el_max
+    element_norm_log = np.log(element_norm + epsy)
+
+    a = fs.array_stoichiometric_coefficients
+    a_sum = np.sum(a)
+    el_matrix = fs.array_species_elements.T
+
+    # Log equilibrium constants for each reaction equation
+    b = -np.sum(fs.get_species_g_rt(t) * a, axis=1)
+
+    # Pressure corrected log equilibrium constants
+    bp = b - np.sum(a * np.log(p / p0), axis=1)
+
+    # Calculating the maximum possible amount for each species based on the elements
+    species_max = np.min(element_norm / (fs.array_species_elements + epsy), axis=1)
+    logn_start = np.log(species_max + epsy)
+
+    # global count
+    # count = 0
+
+    weighting = 100
+
+    def residuals(logn: FloatArray) -> tuple[FloatArray, FloatArray]:
+        # global count
+        # count += 1
+        # print('------', count)
+        # assert count < 100
+        n = np.exp(logn)
+        n_sum = np.sum(n)
+
+        # Residuals from equilibrium equations:
+        resid_eq = np.dot(a, logn - np.log(n_sum)) - bp
+
+        # Jacobian:
+        j_eq = a - a_sum * n / n_sum
+
+        # Residuals from elemental balance:
+        el_sum = np.dot(el_matrix, n)
+        resid_ab = weighting * (np.log(el_sum) - element_norm_log)
+
+        # print(el_sum, element_norm)
+
+        # Jacobian
+        j_ab = weighting * el_matrix * n / el_sum[:, np.newaxis]
+
+        return (np.hstack([resid_eq, resid_ab]), np.concatenate([j_eq, j_ab], axis=0))
+
+    ret = root(residuals, logn_start, jac=True, tol=1e-10)
+    n = np.exp(np.array(ret['x'], dtype=NDFloat))
+    # print(ret)
+
+    return n * el_max
+
+
+def equilibrium(f: fluid | elements, t: float | FloatArray, p: float = 1e5) -> fluid:
+    """Calculate the isobaric equilibrium composition of a fluid at a given temperature and pressure"
+
+    Args:
+        f: Fluid or elements object
+        t: Temperature in Kelvin
+        p: Pressure in Pascal
+
+    Returns:
+        A new fluid object with the equilibrium composition
+    """
+    assert isinstance(f, (fluid, elements)), 'Argument f must be a fluid or elements'
+    m_shape: int = f.fs.array_stoichiometric_coefficients.shape[0]
+    if isinstance(f, fluid):
+        if not m_shape:
+            return f
+    else:
+        if not m_shape:
+            def linalg_lstsq(array_elemental_composition: FloatArray, matrix: FloatArray) -> Any:
+                # TODO: np.dot(np.linalg.pinv(a), b) is eqivalent to lstsq(a, b).
+                # the constant np.linalg.pinv(a) can be precomputed for each fs.
+                return np.dot(np.linalg.pinv(matrix), array_elemental_composition)
+
+            # print('-->', f.array_elemental_composition.shape, f.fs.array_species_elements.transpose().shape)
+            composition = np.apply_along_axis(linalg_lstsq, -1, f.array_elemental_composition, f.fs.array_species_elements.transpose())
+            return fluid(composition, f.fs)
+
+    assert np.min(f.array_elemental_composition) >= 0, 'Input element fractions must be 0 or positive'
+    if isinstance(t, np.ndarray):
+        assert f.shape == tuple(), 'Multidimensional temperature can currently only used for 0D fluids'
+        t_composition = np.zeros(t.shape + (f.fs.array_species_elements.shape[0],))
+        for t_index in np.ndindex(t.shape):
+            t_composition[t_index] = _equilibrium_solver(f.fs, f.array_elemental_composition, float(t[t_index]), p)
+        return fluid(t_composition, f.fs)
+    else:
+        composition = np.ones(f.shape + (len(f.fs.species),), dtype=float)
+        for index in np.ndindex(f.shape):
+            # print(composition.shape, index, _equilibrium(f.fs, f._element_composition[index], t, p))
+            composition[index] = _equilibrium_solver(f.fs, f.array_elemental_composition[index], t, p)
+        return fluid(composition, f.fs)
+
+
+_equilibrium_solver = equilibrium_eq
--- a/src/gaspype/constants.py
+++ b/src/gaspype/constants.py
@ -0,0 +1,24 @@
+"""
+# Constants
+This module contains physical constants used in the gas phase calculations.
+
+kB: Boltzmann constant (J/K)
+NA: Avogadro's number (1/mol)
+R: Ideal gas constant (J/mol/K)
+F: Faraday constant (C/mol)
+p0: Standard pressure 1e5 Pa
+t0: Standard temperature 298.15 K (25 °C)
+p_atm: Standard atmosphere 1 atm = 101325 Pa
+epsy: Small value for numerical stability (1e-18)
+"""
+
+
+kB = 1.380649e-23  # J/K
+NA = 6.02214076e23  # 1/mol
+R = kB * NA  # J/mol/K
+F = 96485.3321233100  # C/mol
+p0 = 1e5  # Pa
+t0 = 273.15 + 25  # K
+p_atm = 101325  # Pa
+
+epsy = 1e-30
--- a/src/gaspype/membrane.py
+++ b/src/gaspype/membrane.py
@ -1,4 +1,7 @@
-from . import R, F, oxygen_partial_pressure, fluid, elements, FloatArray
+from .constants import R, F
+from ._operations import oxygen_partial_pressure
+from ._main import fluid, elements
+from .typing import FloatArray
 import numpy as np

 # Each O2 molecule is transported as two ions and each ion has a charged of 2e
--- a/src/gaspype/py.typed
+++ b/src/gaspype/py.typed
--- a/src/gaspype/typing.py
+++ b/src/gaspype/typing.py
@ -0,0 +1,10 @@
+from numpy import float64
+from numpy.typing import NDArray
+from typing import Sequence
+from types import EllipsisType
+
+Shape = tuple[int, ...]
+NDFloat = float64
+FloatArray = NDArray[NDFloat]
+ArrayIndex = int | slice | None | EllipsisType | Sequence[int]
+ArrayIndices = ArrayIndex | tuple[ArrayIndex, ...]
--- a/tests/benchmark_cp.py
+++ b/tests/benchmark_cp.py
@ -0,0 +1,42 @@
+import cantera as ct
+import numpy as np
+import time
+import gaspype as gp
+
+gas = ct.Solution("gri30.yaml")
+composition = {"H2": 0.3, "H2O": 0.3, "N2": 0.4}
+
+n_species = gas.n_species
+n_states = 1_000_000
+
+# Random temperatures and pressures
+temperatures = np.linspace(300.0, 2500.0, n_states)
+pressures = np.full(n_states, ct.one_atm)
+
+# Create a SolutionArray with many states at once
+states = ct.SolutionArray(gas, len(temperatures))
+
+time.sleep(0.5)
+
+# Vectorized assignment
+t0 = time.perf_counter()
+states.TPX = temperatures, pressures, composition
+cp_values = states.cp_mole
+elapsed = time.perf_counter() - t0
+
+print(f"Computed {n_states} Cp values in {elapsed:.4f} seconds (vectorized cantera)")
+print("First 5 Cp values (J/mol-K):", cp_values[:5] / 1000)
+
+
+# Vectorized fluid creation
+fluid = gp.fluid(composition)
+
+time.sleep(0.5)
+
+# Benchmark: calculate Cp for all states at once
+t0 = time.perf_counter()
+cp_values = fluid.get_cp(t=temperatures)
+elapsed = time.perf_counter() - t0
+
+print(f"Computed {n_states} Cp values in {elapsed:.4f} seconds (vectorized Gaspype)")
+print("First 5 Cp values (J/mol·K):", cp_values[:5])
--- a/tests/benchmark_cp_comp.py
+++ b/tests/benchmark_cp_comp.py
@ -0,0 +1,55 @@
+import cantera as ct
+import numpy as np
+import time
+import gaspype as gp
+
+gas = ct.Solution("gri30.yaml")
+n_species = gas.n_species
+n_states = 1_000_000
+
+# Random temperatures and pressures
+temperatures = np.linspace(300.0, 2500.0, n_states)
+pressures = np.full(n_states, ct.one_atm)
+
+# Generate random compositions for H2, H2O, N2
+rng = np.random.default_rng(seed=42)
+fractions = rng.random((n_states, 3))
+fractions /= fractions.sum(axis=1)[:, None]  # normalize
+
+# Convert to full 53-species mole fraction array
+X = np.zeros((n_states, n_species))
+X[:, gas.species_index('H2')] = fractions[:, 0]
+X[:, gas.species_index('H2O')] = fractions[:, 1]
+X[:, gas.species_index('N2')] = fractions[:, 2]
+
+# Build SolutionArray
+states = ct.SolutionArray(gas, n_states)
+
+time.sleep(0.5)
+
+# Vectorized assignment
+t0 = time.perf_counter()
+states.TPX = temperatures, pressures, X
+cp_values = states.cp_mole
+elapsed = time.perf_counter() - t0
+
+print(f"Computed {n_states} Cp values in {elapsed:.4f} seconds (vectorized cantera)")
+print("First 5 Cp values (J/mol-K):", cp_values[:5] / 1000)
+
+
+# Vectorized fluid creation
+fluid = (
+    gp.fluid({'H2': 1}) * fractions[:, 0]
+    + gp.fluid({'H2O': 1}) * fractions[:, 1]
+    + gp.fluid({'N2': 1}) * fractions[:, 2]
+)
+
+time.sleep(0.5)
+
+# Benchmark: calculate Cp for all states at once
+t0 = time.perf_counter()
+cp_values = fluid.get_cp(t=temperatures)
+elapsed = time.perf_counter() - t0
+
+print(f"Computed {n_states} Cp values in {elapsed:.4f} seconds (vectorized Gaspype)")
+print("First 5 Cp values (J/mol·K):", cp_values[:5])
--- a/tests/benchmark_equalibrium.py
+++ b/tests/benchmark_equalibrium.py
@ -0,0 +1,54 @@
+import cantera as ct
+import gaspype as gp
+import numpy as np
+import time
+from gaspype import fluid_system
+
+# -----------------------
+# Settings
+# -----------------------
+n_temps = 1000
+temps_C = np.linspace(300, 1000, n_temps)     # °C
+temperatures = temps_C + 273.15               # K
+pressure = 101325                             # Pa (1 atm)
+
+composition = {"CH4": 0.8, "H2O": 0.2}
+species_to_track = ["CH4", "H2O", "CO", "CO2", "H2", "O2"]
+
+# -----------------------
+# Cantera benchmark
+# -----------------------
+gas = ct.Solution("gri30.yaml")
+gas.TPX = temperatures[0], pressure, composition
+
+eq_cantera = np.zeros((n_temps, len(species_to_track)))
+
+time.sleep(0.5)
+t0 = time.perf_counter()
+for i, T in enumerate(temperatures):
+    gas.TP = T, pressure
+    gas.equilibrate('TP')
+    eq_cantera[i, :] = [gas.X[gas.species_index(s)] for s in species_to_track]
+elapsed_cantera = time.perf_counter() - t0
+print(f"Cantera: {elapsed_cantera:.4f} s")
+
+# -----------------------
+# Gaspype benchmark
+# -----------------------
+# Construct the fluid with composition and tracked species
+fluid = gp.fluid(composition, fs=fluid_system(species_to_track))
+
+time.sleep(0.5)
+t0 = time.perf_counter()
+eq_gaspype = gp.equilibrium(fluid, t=temperatures, p=pressure)
+elapsed_gaspype = time.perf_counter() - t0
+print(f"Gaspype: {elapsed_gaspype:.4f} s")
+
+# -----------------------
+# Compare first 5 results
+# -----------------------
+print("First 5 equilibrium compositions (mole fractions):")
+for i in range(5):
+    print(f"T = {temperatures[i]:.1f} K")
+    print("  Cantera:", eq_cantera[i])
+    print("  Gaspype :", eq_gaspype.array_composition[i])
--- a/tests/md_to_code.py
+++ b/tests/md_to_code.py
@ -0,0 +1,129 @@
+import re
+from typing import Generator, Iterable
+from dataclasses import dataclass
+import sys
+
+
+@dataclass
+class markdown_segment:
+    code_block: bool
+    language: str
+    text: str
+
+
+def convert_to(target_format: str, md_filename: str, out_filename: str, language: str = 'python'):
+    with open(md_filename, "r") as f_in, open(out_filename, "w") as f_out:
+        segments = segment_markdown(f_in)
+
+        if target_format == 'test':
+            f_out.write('\n'.join(segments_to_test(segments, language)))
+        elif target_format == 'script':
+            f_out.write('\n'.join(segments_to_script(segments, language)))
+        elif target_format == 'striped_markdown':
+            f_out.write('\n'.join(segments_to_striped_markdown(segments, language)))
+        else:
+            raise ValueError('Unknown target format')
+
+
+def segment_markdown(markdown_file: Iterable[str]) -> Generator[markdown_segment, None, None]:
+    regex = re.compile(r"(?:^```\s*(?P<language>(?:\w|-)*)$)", re.MULTILINE)
+
+    block_language: str = ''
+    code_block = False
+    line_buffer: list[str] = []
+
+    for line in markdown_file:
+        match = regex.match(line)
+        if match:
+            if line_buffer:
+                yield markdown_segment(code_block, block_language, ''.join(line_buffer))
+                line_buffer.clear()
+            block_language = match.group('language')
+            code_block = not code_block
+        else:
+            line_buffer.append(line)
+
+    if line_buffer:
+        yield markdown_segment(code_block, block_language, '\n'.join(line_buffer))
+
+
+def segments_to_script(segments: Iterable[markdown_segment], test_language: str = "python") -> Generator[str, None, None]:
+    for segment in segments:
+        if segment.code_block:
+            if segment.language == test_language:
+                yield segment.text
+
+            else:
+                for line in segment.text.splitlines():
+                    yield '#   | ' + line
+            yield ''
+
+        else:
+            for line in segment.text.strip(' \n').splitlines():
+                yield '#  ' + line
+            yield ''
+
+
+def segments_to_striped_markdown(segments: Iterable[markdown_segment], test_language: str = "python") -> Generator[str, None, None]:
+    for segment in segments:
+        if segment.code_block:
+            if segment.language == test_language:
+                yield "``` " + test_language
+                yield segment.text
+                yield "```"
+
+            elif segment.language:
+                for line in segment.text.splitlines():
+                    yield '#   | ' + line
+            yield ''
+
+        else:
+            for line in segment.text.strip(' \n').splitlines():
+                yield '#  ' + line
+            yield ''
+
+
+def segments_to_test(segments: Iterable[markdown_segment], script_language: str = "python") -> Generator[str, None, None]:
+
+    ret_block_flag = False
+
+    yield 'def run_test():'
+
+    for segment in segments:
+        if segment.code_block:
+            if segment.language == script_language:
+                lines = [line for line in segment.text.splitlines() if line.strip()]
+                ret_block_flag = lines[-1] if (not re.match(r'^[^(]*=', lines[-1]) and
+                                               not lines[-1].startswith('import ') and
+                                               not lines[-1].startswith('from ') and
+                                               not lines[-1].startswith('print(') and
+                                               not lines[-1].startswith(' ')) else None
+                # print('Last line: ', ret_block_flag, '-----------', lines[-1])
+
+                yield ''
+                yield '    print("---------------------------------------------------------")'
+                yield ''
+                if ret_block_flag:
+                    yield from ['    ' + str(line) for line in segment.text.splitlines()[:-1]]
+                    yield f'    print("-- Result (({ret_block_flag})):")'
+                    yield f'    print(({ret_block_flag}).__repr__().strip())'
+                else:
+                    yield from ['    ' + str(line) for line in segment.text.splitlines()]
+
+            elif ret_block_flag:
+                yield '    ref_str = r"""'
+                yield from [str(line) for line in segment.text.splitlines()]
+                yield '"""'
+                yield f'    print("-- Reference (({ret_block_flag})):")'
+                yield '    print(ref_str.strip())'
+                yield f'    assert ({ret_block_flag}).__repr__().strip() == ref_str.strip()'
+                ret_block_flag = False
+
+    yield '\nif __name__ == "__main__":'
+    yield '    run_test()'
+
+
+if __name__ == "__main__":
+    format = sys.argv[1]
+    assert format in ['test', 'script']
+    convert_to(sys.argv[1], sys.argv[2], sys.argv[3])
--- a/tests/test_db_reader.py
+++ b/tests/test_db_reader.py
@ -9,6 +9,7 @@ def test_db_reader():
    assert 'TEST' not in db
    assert db['HCl'].name == 'HCl'
    assert db['CH4'].composition == {'C': 1, 'H': 4}
+    assert db['C2H5OH'].composition == {'C': 2, 'H': 6, 'O': 1}
    assert db['H2O'].model == 9

    for species in db:
--- a/tests/test_doc_examples.py
+++ b/tests/test_doc_examples.py
@ -1,53 +1,28 @@
-import re
-from typing import Generator
-
-
-def convert_markdown_file(md_filename: str, out_filename: str):
-    with open(md_filename, "r") as f_in:
-        with open(out_filename, "w") as f_out:
-            f_out.write('def run_test():\n')
-            for block in markdown_to_code([line for line in f_in]):
-                f_out.write(block + '\n')
-
-
-def markdown_to_code(lines: list[str], language: str = "python") -> Generator[str, None, None]:
-    regex = re.compile(
-        r"(?P<start>^```\s*(?P<block_language>(\w|-)*)\n)(?P<code>.*?\n)(?P<end>```)",
-        re.DOTALL | re.MULTILINE,
-    )
-    blocks = [
-        (match.group("block_language"), match.group("code"))
-        for match in regex.finditer("".join(lines))
-    ]
-
-    ret_block_flag = False
-
-    for block_language, block in blocks:
-        if block_language == language:
-            lines = [line for line in block.splitlines() if line.strip()]
-            ret_block_flag = lines[-1] if '=' not in lines[-1] else None
-
-            yield ''
-            yield '    print("---------------------------------------------------------")'
-            yield ''
-            if ret_block_flag:
-                yield from ['    ' + str(line) for line in block.splitlines()[:-1]]
-            else:
-                yield from ['    ' + str(line) for line in block.splitlines()]
-            yield f'    print("-- Result (({ret_block_flag})):")'
-            yield f'    print(({ret_block_flag}).__repr__().strip())'
-
-        elif ret_block_flag:
-            yield '    ref_str = r"""'
-            yield from [str(line) for line in block.splitlines()]
-            yield '"""'
-            yield f'    print("-- Reference (({ret_block_flag})):")'
-            yield '    print(ref_str.strip())'
-            yield f'    assert ({ret_block_flag}).__repr__().strip() == ref_str.strip()'
-            ret_block_flag = False
+import md_to_code
+from glob import glob
+import importlib
+import os


 def test_readme():
-    convert_markdown_file('README.md', 'tests/autogenerated_readme.py')
+    md_to_code.convert_to('test', 'README.md', 'tests/autogenerated_readme.py')
    import autogenerated_readme
    autogenerated_readme.run_test()
+
+
+def test_example_code():
+    filter = 'examples/*.md'
+
+    files = glob(filter)
+    for path in files:
+        file_name = '.'.join(os.path.basename(path).split('.')[:-1])
+        if not file_name.lower() == 'readme':
+            print(f"> Test Example {file_name} ...")
+            md_to_code.convert_to('test', path, f'tests/autogenerated_{file_name}.py')
+            mod = importlib.import_module(f'autogenerated_{file_name}')
+            mod.run_test()
+
+
+if __name__ == "__main__":
+    test_readme()
+    test_example_code()
--- a/tests/test_edge_case_results.py
+++ b/tests/test_edge_case_results.py
@ -0,0 +1,20 @@
+import pytest
+import gaspype as gp
+
+
+def test_calculation_cold():
+    # Testing a non-equilibrium case. Result is only
+    # determined by stoichiometry.
+    fs = gp.fluid_system(['CH4', 'H2O', 'H2', 'CO2', 'CO', 'O2'])
+
+    el = gp.elements(gp.fluid({'H2': 1, 'CH4': 0.08, 'O2': 0.05}), fs=fs)
+
+    composition = gp.equilibrium(el, 300, 1e5)
+
+    print(el)
+    print(composition)
+
+    ref_result = [8.00000000e-02, 1.00000000e-01, 9.00000000e-01, 3.01246781e-23,
+                  2.90783583e-27, 3.60487456e-82]
+
+    assert composition.array_composition == pytest.approx(ref_result, abs=0.001)
--- a/tests/test_results_cantera.py
+++ b/tests/test_results_cantera.py
@ -1,7 +1,8 @@
 import gaspype as gp
 import numpy as np
 # import pytest
-import cantera as ct  # type: ignore
+import cantera as ct
+from gaspype.typing import NDFloat


 def test_equilibrium_cantera():
@ -16,7 +17,7 @@ def test_equilibrium_cantera():

    composition = gp.fluid({'H2': 1}, fs) +\
        gp.fluid({'CH4': 1}, fs) * np.linspace(0, 0.05, 30) +\
-        gp.fluid({'O2': 1}, fs) * np.expand_dims(np.linspace(0, 0.5, 30), axis=1)
+        gp.fluid({'O2': 1}, fs) * np.linspace(0, 0.5, 30)[:, None]

    t = 1495 + 273.15  # K
    p = 1e5  # Pa
@ -40,10 +41,12 @@ def test_equilibrium_cantera():
        #if flow1.X[indeces][0] > 0.01:
        #    print(flow1.X[indeces])

-    ct_result_array = np.stack(ct_results)  # type: ignore
+    ct_result_array = np.stack(ct_results, dtype=NDFloat)  # type: ignore

-    max_deviation = np.max(np.abs((gp_result_array - ct_result_array)))
-    mean_deviation = np.mean(np.abs((gp_result_array - ct_result_array)))
+    deviations = np.abs(gp_result_array - ct_result_array)

-    assert max_deviation < 0.04
-    assert mean_deviation < 2e-4
+    for dev, gp_comp_result, ct_comp_result in zip(deviations, gp_result_array, ct_result_array):
+        print(f"Composition: {gp_comp_result} / {ct_comp_result}")
+        assert np.all(dev < 0.04), f"Deviateion: {dev}"
+
+    assert np.mean(deviations) < 2e-4
--- a/tests/test_slicing.py
+++ b/tests/test_slicing.py
@ -12,28 +12,28 @@ def test_str_index():
    assert el['C'].shape == (2, 3, 4)


-def test_str_list_index():
-    assert fl[['CO2', 'H2', 'CO']].shape == (2, 3, 4, 3)
-    assert el[['C', 'H', 'O']].shape == (2, 3, 4, 3)
+def test_single_axis_int_index():
+    assert fl[0].shape == (3, 4)
+    assert fl[1].shape == (3, 4)
+    assert el[1].shape == (3, 4)
+    assert el[0].shape == (3, 4)


-def test_int_list_index():
-    assert fl[[1, 2, 0, 5]].shape == (2, 3, 4, 4)
-    assert el[[1, 2, 0, 3]].shape == (2, 3, 4, 4)
+def test_single_axis_int_list():
+    assert fl[:, [0, 1]].shape == (2, 2, 4)
+    assert el[:, [0, 1]].shape == (2, 2, 4)


-def test_mixed_list_index():
-    assert el[[1, 'H', 0, 'O']].shape == (2, 3, 4, 4)
-
-
-def test_int_index():
-    assert fl[5].shape == (2, 3, 4)
-    assert el[-1].shape == (2, 3, 4)
-
-
-def test_slice_index():
-    assert fl[0:3].shape == (2, 3, 4, 3)
-    assert fl[:].shape == (2, 3, 4, 6)
-
-    assert el[0:3].shape == (2, 3, 4, 3)
-    assert el[:].shape == (2, 3, 4, 4)
+def test_multi_axis_int_index():
+    assert fl[0, 1].shape == (4,)
+    assert fl[0, 1, 2].shape == tuple()
+    assert fl[0, 2].shape == (4,)
+    assert fl[:, 2, :].shape == (2, 4)
+    assert fl[0, [1, 2]].shape == (2, 4)
+    assert fl[..., 0].shape == (2, 3)
+    assert el[0, 1].shape == (4,)
+    assert el[0, 1, 2].shape == tuple()
+    assert el[0, 2].shape == (4,)
+    assert el[:, 2, :].shape == (2, 4)
+    assert el[0, [1, 2]].shape == (2, 4)
+    assert el[..., 0].shape == (2, 3)
--- a/tests/test_to_pandas.py
+++ b/tests/test_to_pandas.py
@ -21,3 +21,14 @@ def test_elements():

    df = pd.DataFrame(list(gp.elements(fl * np.array([1, 2, 3, 4]))))
    assert df.shape == (4, 2)
+
+
+def test_iter():
+    fl = gp.fluid({'O2': 1, 'H2': 2, 'H2O': 3})
+
+    fl2 = fl * np.array([1, 2, 3, 4])
+    for i, f in enumerate(fl2):
+        if i == 1:
+            assert f == {'O2': np.float64(2.0), 'H2': np.float64(4.0), 'H2O': np.float64(6.0)}
+        if i == 3:
+            assert f == {'O2': np.float64(4.0), 'H2': np.float64(8.0), 'H2O': np.float64(12.0)}
Author	SHA1	Message	Date
Nicolas	21c6f3cb72	SOC examples updated to use new slicing	2025-09-03 15:57:55 +02:00
Nicolas	50b32bfa8f	Numpy-like array slicing for fluids and elements added; get_n() function added for fluids and elements	2025-09-03 15:57:55 +02:00
Nicolas	a04340dd7f	Version changed to 1.1.3; config files updated	2025-09-03 15:57:55 +02:00
Nicolas	fe0fd8ec21	new lines at the ends of the benchmark scripts added	2025-09-03 15:57:55 +02:00
Nicolas	b8baa14c84	benchmark scripts added	2025-09-03 15:57:55 +02:00
Nicolas	cc1057a99b	CD for pypi: environment added	2025-07-28 17:11:23 +02:00
Nicolas	e2a6541119	Docs: Link to repo added	2025-07-28 17:11:23 +02:00
Nicolas	dd86a1a83a	doc path renamed to "api"	2025-07-28 17:11:23 +02:00
Nicolas Kruse	38709c03ac	md_to_code.py updated to handle print functions in the last line	2025-07-28 15:06:32 +02:00
Nicolas Kruse	77b80a31fc	sofc example fixes of the extension	2025-07-28 15:06:32 +02:00
Nicolas	138e33e5b6	sofc example extended	2025-07-28 15:06:32 +02:00
Nicolas	14bd6605bc	switched docs deployment method	2025-07-28 15:06:32 +02:00
Nicolas	eb171fd424	fix in SOEC example	2025-07-28 15:06:32 +02:00
Nicolas Kruse	0e01d5e7be	readme updated	2025-07-21 11:39:51 +02:00
Nicolas Kruse	b277f91f5f	example text improved	2025-07-21 11:39:51 +02:00
Nicolas Kruse	c966aa74f3	unified nomenclature in methane example	2025-07-21 11:39:51 +02:00
Nicolas	f97c8a2697	py.typed file added indication package uses type annotations	2025-07-21 11:39:51 +02:00
Nicolas Kruse	10ade221e0	fix in example description	2025-07-21 11:39:51 +02:00
Nicolas Kruse	8627b59e0b	comment lines in example removed	2025-07-21 11:39:51 +02:00
Nicolas Kruse	73627613ac	Example description fixed	2025-07-03 15:11:17 +02:00
Nicolas Kruse	ac026600c5	example added	2025-07-03 15:11:17 +02:00
Nicolas Kruse	6f3ce50999	figures in sulfur example improved	2025-07-03 15:11:17 +02:00
Nicolas Kruse	e4b33b344b	Minor README fixes	2025-07-03 15:11:17 +02:00
Nicolas Kruse	7a4c705fe1	carbon example extended	2025-07-03 15:11:17 +02:00
Nicolas Kruse	e2e09a6ea4	sofc example fixed and soec sample added	2025-06-25 11:15:27 +02:00
Nicolas Kruse	df56d0ad37	citation.cff acording to cffconvert is fixed; validation to ci added CI: missing "jobs:" readded CITATION.cff fixed and CI updated for CITATION.cff validation	2025-06-24 16:37:43 +02:00
Nicolas Kruse	9a4c6351bb	citation file updated	2025-06-24 15:14:41 +02:00
Nicolas Kruse	5d5b8ec722	codestyle for test fixed	2025-06-24 15:14:41 +02:00
Nicolas Kruse	2bf7e21e73	test with cantera reference updated	2025-06-24 15:14:41 +02:00
Nicolas Kruse	762ff6929d	test added for non-equilibrium cases	2025-06-24 15:14:41 +02:00
Nicolas Kruse	0854d28e92	Time hints added in main	2025-06-24 15:14:41 +02:00
Nicolas Kruse	5172d3721e	added weighting coefficient for better stability	2025-06-24 15:14:41 +02:00
Nicolas Kruse	fe5bf70d83	version changed to v1.1.2	2025-06-24 15:14:41 +02:00
Nicolas Kruse	63b1808341	Solver updated with start value estimation and log-error for elemental species balance	2025-06-24 15:14:41 +02:00
Nicolas Kruse	a5f806aa31	For equilibrium_eq start values are estimated now by calculating the maximum possible amount for each species based on the elements	2025-06-24 15:14:41 +02:00
Nicolas Kruse	802aa68a25	equilibrium solvers moved to separate file "_solver.py"	2025-06-24 15:14:41 +02:00
Nicolas Kruse	973aedf058	cff-version in CITATION.cff fixed	2025-06-18 14:41:41 +02:00
Nicolas Kruse	d627b87df6	added mypy conformity for numpy 2.3.0	2025-06-18 14:18:46 +02:00
Nicolas Kruse	2f5cad3d05	Docs and doc generation updated for the restructuring	2025-06-18 08:48:51 +02:00
Nicolas Kruse	93f6d88cd4	docstring fixed for species function	2025-06-18 00:46:31 +02:00
Nicolas Kruse	fd94cb2f24	__iter__ function for elements and fluid updated to pass mypy type check. test added	2025-06-18 00:19:19 +02:00
Nicolas Kruse	1e54392d54	type references fixed	2025-06-17 23:31:17 +02:00
Nicolas Kruse	944aeb3546	citation.cff updated	2025-06-17 23:22:23 +02:00
Nicolas Kruse	a12983948e	moving type aliases to typing	2025-06-17 23:18:24 +02:00
Nicolas Kruse	91ec565bf8	version number increased to 1.1.1	2025-06-17 22:57:20 +02:00
Nicolas Kruse	9e7c9a47b3	db_reader memory usage reduced by using an iterator instead of split	2025-06-17 22:56:43 +02:00
Nicolas Kruse	412827226f	File structure refactored to clean up API	2025-06-17 22:55:04 +02:00
Nicolas Kruse	5c41d04a70	test extended	2025-06-17 22:53:25 +02:00
Nicolas Kruse	3016a2db4e	Example files fixed	2025-06-06 14:17:04 +02:00
Nicolas Kruse	5f6067fc63	Example file renamed to prevent name conflict	2025-06-06 14:03:13 +02:00
Nicolas Kruse	71da33d9ad	doc CD: Ignore Warnings when running notedown command	2025-06-06 13:57:46 +02:00
Nicolas Kruse	08c81444dd	doc CD: fixed package name	2025-06-06 13:52:17 +02:00
Nicolas Kruse	6659545391	doc CD: added yaml dependency	2025-06-06 13:50:52 +02:00
Nicolas Kruse	7af6176464	doc CD: fixed invalid comment	2025-06-06 13:47:51 +02:00
Nicolas Kruse	90b3c77ae1	doc CD: replaced dummy database by a partial real database for rendering examples	2025-06-06 13:44:13 +02:00
Nicolas Kruse	a65dbee17d	added "os.makedirs('docs/source/_autogenerated', exist_ok=True)" to create _autogenerated directory if not existing	2025-06-06 13:36:55 +02:00
Nicolas Kruse	ce48083e77	Merge pull request #2 from DLR-Institute-of-Future-Fuels/example_generation_feature - Example generation feature - Code example added - Unit tests and doc generation/config updated to handle examples	2025-06-06 13:29:03 +02:00
Nicolas Kruse	dafdade833	Excluded autogenerated example code files in flake8 config, code style fixed	2025-06-06 13:24:46 +02:00
Nicolas Kruse	7671a89e14	Examples added, example readme and example unit test updated	2025-06-06 13:15:34 +02:00
Nicolas Kruse	0578b552be	Example-readme updated	2025-06-06 10:48:30 +02:00
Nicolas Kruse	d18ba0f785	Docs CD updated	2025-06-06 10:48:08 +02:00
Nicolas Kruse	1c07ffea36	docs CD: LICENSE added to docs	2025-06-06 10:13:34 +02:00
Nicolas Kruse	38eda3edcf	docs for types added	2025-06-06 10:12:55 +02:00
Nicolas Kruse	9fc5ea2dc1	Docstrings updated: types for properties and solver description added	2025-06-06 10:11:54 +02:00
Nicolas Kruse	b158a86852	doc build CD updated for rendering examples	2025-06-06 09:08:13 +02:00
Nicolas Kruse	99e543e12a	script for rendering examples added for doc generation	2025-06-06 09:07:43 +02:00
Nicolas Kruse	359120bbcf	project settings updated - e.g. dependencies for building docs and testing	2025-06-06 09:06:33 +02:00
Nicolas Kruse	6c3437f509	example unit test updated to test all examples from example directory	2025-06-06 09:05:40 +02:00
Nicolas Kruse	793b2a0ab4	doc build script and settings updated	2025-06-06 09:04:42 +02:00
Nicolas Kruse	e30b4f1d47	files rearranged	2025-06-06 09:02:52 +02:00
Nicolas Kruse	60a8f99daf	drawing for example added	2025-06-04 18:00:40 +02:00
Nicolas Kruse	cf7c5ebe95	doc_build dependencies in project file added	2025-06-04 18:00:28 +02:00
Nicolas Kruse	16d1d053d7	readme for examples added	2025-06-04 18:00:01 +02:00
Nicolas Kruse	ddd543abf0	line breaks adjusted in readme	2025-06-04 17:54:09 +02:00
Nicolas Kruse	360683a633	md_to_code updated to generate tests and human readable example code from markdown scripts	2025-06-04 17:53:37 +02:00