Repository: amepproject/amep (language: python)

Issue:
Website readme optimization
Website and `README.md` optimized.
DOC: Incorrect description of 'coords' and 'other_coords' in 'order.nearest_neighbors' and 'order.k_nearest_neighbors'
### Problem description:

These two descriptions should be interchanged:

**coords** (np.ndarray of shape (N,3)) – Array of coordinates/points to search for neighbors.

**other_coords** (np.ndarray of shape (M,3) or None, optional) – Array of points to search for neighbors of. The default is None (use coords).

### Suggested content improvement:

Should read:

**coords** (np.ndarray of shape (N,3)) – Array of points to search for neighbors of.

**other_coords** (np.ndarray of shape (M,3) or None, optional) –  Array of coordinates/points to search for neighbors. The default is None (use coords).


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index d5d1b6e..538b33e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -43,3 +43,11 @@ build/
 test/data
 test/.coverage
 test/coverage.xml
+examples/data/continuum/traj.h5amep
+examples/data/continuum/#traj.h5amep
+examples/data/continuum/##temp*
+examples/data/trajs/
+examples/data/lammps/traj.h5amep
+examples/data/lammps/#traj.h5amep
+examples/data/lammps/##temp*
+examples/data/invalid
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 62ec1fd..3d95355 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,26 @@
 All notable changes to **AMEP** will be documented in this file. **AMEP** 
 adheres to [semantic versioning](https://semver.org/).
 
+## AMEP 1.0.3 (22 Oct 2024)
+
+### Bug fixes:
+
+* documentation improvements (ME, KS, KD, LH)
+* bug fix related to `trajectory.add_particle_info` and `trajectory.get_particle_info` (KD, KS)
+* LammpsReader angular momentum import bug fixed (KD)
+* bug in `evaluate.Psi6dist` related to default value of number of bins fixed (KD)
+* bug in `evaluate.VelDist` related to default value of number of bins fixed (KD)
+* faster tests - now based on example data (KS, LH, KD)
+
+### Contributors:
+
+* Lukas Hecht (LH)
+* Kay-Robert Dormann (KD)
+* Kai Luca Spanheimer (KS)
+* Mahdieh Ebrahimi (ME)
+
+
+
 ## AMEP 1.0.2 (22 Mai 2024)
 
 ### Bug fixes:
diff --git a/README.md b/README.md
index 67ecfc8..1a9e189 100644
--- a/README.md
+++ b/README.md
@@ -22,68 +22,38 @@ format. To be fast and usable on modern HPC (**H**igh **P**erformance
 **C**omputing) hardware, the methods are optimized to run also in parallel.
 
 
-# Description
+# How to cite AMEP
 
-The **AMEP** Python library provides a unified framework for handling 
-both particle-based and continuum simulation data. It is made for the analysis 
-of molecular-dynamics (MD), Brownian-dynamics (BD), and continuum simulation 
-data of condensed matter systems and active matter systems in particular. 
-**AMEP** provides a huge variety of analysis methods for both data types that 
-allow to evaluate various dynamic and static observables based on the 
-trajectories of the particles or the time evolution of continuum fields. For 
-fast and efficient data handling, **AMEP** provides a unified framework for 
-loading and storing simulation data and analysis results in a compressed, 
-HDF5-based data format. **AMEP** is written purely in Python and uses powerful 
-libraries such as NumPy, SciPy, Matplotlib, and scikit-image commonly used in 
-computational physics. Therefore, understanding, modifying, and building up on 
-the provided framework is comparatively easy. All evaluation functions are 
-optimized to run efficiently on HPC hardware to provide fast computations. To 
-plot and visualize simulation data and analysis results, **AMEP** provides an 
-optimized plotting framework based on the Matplotlib Python library, which 
-allows to easily plot and animate particles, fields, and lines. Compared to 
-other analysis libraries, the huge variety of analysis methods combined with 
-the possibility to handle both most common data types used in soft-matter 
-physics and in the active matter community in particular, enables the analysis 
-of a much broader class of simulation data including not only classical 
-molecular-dynamics or Brownian-dynamics simulations but also any kind of 
-numerical solutions of partial differential equations. The following table 
-gives an overview on the observables provided by **AMEP** and on their 
-capability of running in parallel and processing particle-based and continuum 
-simulation data.
+If you use **AMEP** for a project that leads to a scientific publication, please acknowledge 
+the use of **AMEP** within the body of your publication for example by copying or adapting 
+the following formulation:
+
+*Data analysis for this publication utilized the AMEP library [1].*
 
+> [1] L. Hecht, K.-R. Dormann, K. L. Spanheimer, M. Ebrahimi, M. Cordts, S. Mandal, 
+>     A. K. Mukhopadhyay, and B. Liebchen, AMEP: The Active Matter Evaluation Package for Python, 
+>     *arXiv [Cond-Mat.Soft]* (2024). Available at: http://arxiv.org/abs/2404.16533.
 
-| Observable | Parallel | Particles | Fields |
-|:-----------|:--------:|:---------:|:------:|
-| **Spatial Correlation Functions:** ||||
-| RDF (radial pair distribution function) | ✔ | ✔ | ❌ |
-| PCF2d (2d pair correlation function) | ✔ | ✔ | ❌ |
-| PCFangle (angular pair correlation function) | ✔ | ✔ | ❌ |
-| SFiso (isotropic static structure factor) | ✔ | ✔ | ✔ |
-| SF2d (2d static structure factor) | ✔ | ✔ | ✔ |
-| SpatialVelCor (spatial velocity correlation function) | ✔ | ✔ | ❌ |
-| PosOrderCor (positional order correlation function) | ✔ | ✔ | ❌ |
-| HexOrderCor (hexagonal order correlation function) | ✔ | ✔ | ❌ |
-| **Local Order:** ||||
-| Voronoi tesselation | ❌ | ✔ | ❌ |
-| Local density | ❌ | ✔ | ✔ |
-| Local packing fraction | ❌ | ✔ | ❌ |
-| k-atic bond order parameter | ❌ | ✔ | ❌ |
-| Next/nearest neighbor search | ❌ | ✔ | ❌ |
-| **Time Correlation Functions:** ||||
-| MSD (mean square displacement) | ❌ | ✔ | ❌ |
-| VACF (velocity autocorrelation function) | ❌ | ✔ | ❌ |
-| OACF (orientation autocorrelation function) | ❌ | ✔ | ❌ |
-| **Cluster Analysis:** ||||
-| Clustersize distribution | ❌ | ✔ | ✔ |
-| Cluster growth | ❌ | ✔ | ✔ |
-| Radius of gyration | ❌ | ✔ | ✔ |
-| Linear extension | ❌ | ✔ | ✔ |
-| Center of mass | ❌ | ✔ | ✔ |
-| Gyration tensor | ❌ | ✔ | ✔ |
-| Inertia tensor | ❌ | ✔ | ✔ |
-| **Miscellaneous:** ||||
-| Translational/rotational kinetic energy | ❌ | ✔ | ❌ |
-| Kinetic temperature | ❌ | ✔ | ❌ |
+The pre-print is freely available on [arXiv](https://arxiv.org/abs/2404.16533). To cite this reference, 
+you can use the following BibTeX entry:
+
+```bibtex
+@misc{hecht2024amep,
+    title = {AMEP: The Active Matter Evaluation Package for Python}, 
+    author = {Lukas Hecht and 
+              Kay-Robert Dormann and 
+              Kai Luca Spanheimer and 
+              Mahdieh Ebrahimi and 
+              Malte Cordts and 
+              Suvendu Mandal and 
+              Aritra K. Mukhopadhyay and 
+              Benno Liebchen},
+    year = {2024},
+    eprint = {2404.16533},
+    archivePrefix = {arXiv},
+    primaryClass = {cond-mat.soft}
+}
+```
 
 
 # Installation
@@ -146,40 +116,6 @@ Anaconda installation path.
 to download FFmpeg and to get further information on how to install FFmpeg on your machine.
 
 
-# Citation
-
-If you use **AMEP** for a project that leads to a scientific publication, please acknowledge 
-the use of **AMEP** within the body of your publication for example by copying or adapting 
-the following formulation:
-
-*Data analysis for this publication utilized the AMEP library [1].*
-
-> [1] L. Hecht, K.-R. Dormann, K. L. Spanheimer, M. Ebrahimi, M. Cordts, S. Mandal, 
->     A. K. Mukhopadhyay, and B. Liebchen, AMEP: The Active Matter Evaluation Package for Python, 
->     *arXiv [Cond-Mat.Soft]* (2024). Available at: http://arxiv.org/abs/2404.16533.
-
-The pre-print is freely available on [arXiv](https://arxiv.org/abs/2404.16533). To cite this reference, 
-you can use the following BibTeX entry:
-
-```bibtex
-@misc{hecht2024amep,
-    title = {AMEP: The Active Matter Evaluation Package for Python}, 
-    author = {Lukas Hecht and 
-              Kay-Robert Dormann and 
-              Kai Luca Spanheimer and 
-              Mahdieh Ebrahimi and 
-              Malte Cordts and 
-              Suvendu Mandal and 
-              Aritra K. Mukhopadhyay and 
-              Benno Liebchen},
-    year = {2024},
-    eprint = {2404.16533},
-    archivePrefix = {arXiv},
-    primaryClass = {cond-mat.soft}
-}
-```
-
-
 # Getting started
 
 The following example briefly demonstrates the **AMEP** workflow. A typical 
@@ -218,6 +154,70 @@ fig.savefig(rdf.name + '.pdf')
 For more detailed examples, check the [examples](https://github.com/amepproject/amep/tree/main/examples) directory.
 
 
+# Project description
+
+The **AMEP** Python library provides a unified framework for handling 
+both particle-based and continuum simulation data. It is made for the analysis 
+of molecular-dynamics (MD), Brownian-dynamics (BD), and continuum simulation 
+data of condensed matter systems and active matter systems in particular. 
+**AMEP** provides a huge variety of analysis methods for both data types that 
+allow to evaluate various dynamic and static observables based on the 
+trajectories of the particles or the time evolution of continuum fields. For 
+fast and efficient data handling, **AMEP** provides a unified framework for 
+loading and storing simulation data and analysis results in a compressed, 
+HDF5-based data format. **AMEP** is written purely in Python and uses powerful 
+libraries such as NumPy, SciPy, Matplotlib, and scikit-image commonly used in 
+computational physics. Therefore, understanding, modifying, and building up on 
+the provided framework is comparatively easy. All evaluation functions are 
+optimized to run efficiently on HPC hardware to provide fast computations. To 
+plot and visualize simulation data and analysis results, **AMEP** provides an 
+optimized plotting framework based on the Matplotlib Python library, which 
+allows to easily plot and animate particles, fields, and lines. Compared to 
+other analysis libraries, the huge variety of analysis methods combined with 
+the possibility to handle both most common data types used in soft-matter 
+physics and in the active matter community in particular, enables the analysis 
+of a much broader class of simulation data including not only classical 
+molecular-dynamics or Brownian-dynamics simulations but also any kind of 
+numerical solutions of partial differential equations. The following table 
+gives an overview on the observables provided by **AMEP** and on their 
+capability of processing particle-based and continuum 
+simulation data.
+
+
+| Observable | Particles | Fields |
+|:-----------|:---------:|:------:|
+| **Spatial Correlation Functions:** |||
+| RDF (radial pair distribution function) | ✔ | ➖ |
+| PCF2d (2d pair correlation function) | ✔ | ➖ |
+| PCFangle (angular pair correlation function) | ✔ | ➖ |
+| SFiso (isotropic static structure factor) | ✔ | ✔ |
+| SF2d (2d static structure factor) | ✔ | ✔ |
+| SpatialVelCor (spatial velocity correlation function) | ✔ | ➖ |
+| PosOrderCor (positional order correlation function) | ✔ | ➖ |
+| HexOrderCor (hexagonal order correlation function) | ✔ | ➖ |
+| **Local Order:** |||
+| Voronoi tesselation | ✔ | ➖ |
+| Local density | ✔ | ✔ |
+| Local packing fraction | ✔ | ➖ |
+| k-atic bond order parameter | ✔ | ➖ |
+| Next/nearest neighbor search | ✔ | ➖ |
+| **Time Correlation Functions:** |||
+| MSD (mean square displacement) | ✔ | ➖ |
+| VACF (velocity autocorrelation function) | ✔ | ➖ |
+| OACF (orientation autocorrelation function) | ✔ | ➖ |
+| **Cluster Analysis:** |||
+| Clustersize distribution | ✔ | ✔ |
+| Cluster growth | ✔ | ✔ |
+| Radius of gyration | ✔ | ✔ |
+| Linear extension | ✔ | ✔ |
+| Center of mass | ✔ | ✔ |
+| Gyration tensor | ✔ | ✔ |
+| Inertia tensor | ✔ | ✔ |
+| **Miscellaneous:** |||
+| Translational/rotational kinetic energy | ✔ | ➖ |
+| Kinetic temperature | ✔ | ➖ |
+
+
 # Module descriptions
 
 In the following, we provide a list of all **AMEP** modules together with a 
diff --git a/amep/_version.py b/amep/_version.py
index 447f16a..8f2e8c5 100644
--- a/amep/_version.py
+++ b/amep/_version.py
@@ -20,4 +20,4 @@
 """
 AMEP version number.
 """
-__version__ = "1.0.2"
\ No newline at end of file
+__version__ = "1.0.3"
\ No newline at end of file
diff --git a/amep/base.py b/amep/base.py
index 89e7621..86dfb74 100644
--- a/amep/base.py
+++ b/amep/base.py
@@ -36,7 +36,7 @@
 import inspect
 import logging
 
-from typing import Collection
+from typing import Collection, Iterable, Sequence
 from  io import StringIO
 from contextlib import redirect_stdout
 from tqdm import TqdmExperimentalWarning
@@ -168,7 +168,8 @@ def check_path(path: str, extension: str) -> tuple[str, str]:
 LOADMODES = [
     'lammps',
     'h5amep',
-    'field'
+    'field',
+    'gsd'
 ]
 # maximum RAM usage in GB per CPU (used for parallelized methods)
 MAXMEM = 1
@@ -899,7 +900,7 @@ def data(
 
     def __read_data(
             self, key: str, ptype: int | list | None = None,
-            pid: int | list | None = None) -> np.ndarray:
+            pid: int | Sequence | None = None) -> np.ndarray:
         r'''
         Reads a dataset from the HDF5 file either for all particles
         or for all particles of a specific particle type.
@@ -929,12 +930,10 @@ def __read_data(
 
             if pid:
                 data_ids = root['frames'][str(self.__step)]['id'][:]
-                if type(pid) != list:
+                if not isinstance(pid, Sequence):
                     pid = [pid]
-                id_list = []
-                for id in pid:
-                    index_id = int(np.where(data_ids==id)[0])
-                    id_list.append(index_id)
+                id_list = [int(np.where(data_ids == part_id)[0])
+                           for part_id in pid]
                 if key in root['frames'][str(self.__step)].keys():
                     data = root['frames'][str(self.__step)][key][id_list]
                     return data
@@ -944,28 +943,28 @@ def __read_data(
             # check if a dataset with the given key exists and read the data
             if key in root['frames'][str(self.__step)].keys():
                 data = root['frames'][str(self.__step)][key][:]
-
-                # Transform ptype so it is a list of all ptypes
+                # If no ptype is provided return all data
                 if not ptype:
                     return data
-                if type(ptype) == int:
+                # Transform ptype so it is a list of all ptypes
+                if not isinstance(ptype, Sequence):
                     ptype = [ptype]
                 mask = np.zeros(data.shape[0],)
                 # check particle type
-                for single_ptype in ptype:
-                    if single_ptype in self.ptypes:
-                        mask += types==single_ptype
-                    else:
-                        warnings.warn(
-                            f"The specified particle type {single_ptype} "\
-                            "does not exist. Returning data without type "\
-                            f"{single_ptype}."
-                        )
+                mask += sum(types == single_ptype for single_ptype
+                            in ptype if single_ptype in self.ptypes)
+                if any(single_ptype not in self.ptypes for
+                       single_ptype in ptype):
+                    for single_ptype in ptype:
+                        if single_ptype not in self.ptypes:
+                            warnings.warn(
+                                f"The specified particle type {single_ptype} "
+                                "does not exist. Returning data without type "
+                                f"{single_ptype}."
+                            )
                 return data[mask.astype(bool)]
-
-            else:
-                raise KeyError(
-                    f"The key {key} does not exist in the frame. "\
+            raise KeyError(
+                    f"The key {key} does not exist in the frame. "
                     "Returning no data!"
                 )
 
@@ -1266,42 +1265,55 @@ def __init__(self, reader):
 
         '''
         self.__reader = reader
-    def __getitem__(self,item):
+
+    def __getitem__(self, item: int | slice | Iterable[int]
+                    ) -> BaseFrame | BaseField | list[BaseField | BaseFrame]:
+        """Get an individual frame or field of a simulation.
+
+        Supports slicing as well as iterables of valid integer indices.
+        The return type depends on the type of the trajectory.
+        Also it depends if only one frame or a collection of frames is requested.
+        If a collection of frames is requested a list of frames is returned.
+
+        Parameters
+        ----------
+        item : int | slice | Iterable[int]
+
+        Returns
+        -------
+        BaseFrame | BaseField | list[BaseField | BaseFrame] 
+        """
         if isinstance(item, slice):
-            sli=range(*item.indices(len(self.__reader.steps)))
-            if self.type=="field":
-                out = [BaseField(self.__reader,index) for index in sli]
-            elif self.type=="particle":
-                out = [BaseFrame(self.__reader,index) for index in sli]
-            else:
-                out = [BaseFrame(self.__reader,index) for index in sli]
-            return out
-        elif isinstance(item, list) or isinstance(item, np.ndarray):
-            if self.type=="field":
-                out = [BaseField(self.__reader,index) for index in item]
-            elif self.type=="particle":
-                out = [BaseFrame(self.__reader,index) for index in item]
-            else:
-                out = [BaseFrame(self.__reader,index) for index in item]
-            return out
-        elif isinstance(item, int) or isinstance(item, np.integer):
-            if self.type=="field":
+            sli = range(*item.indices(len(self.__reader.steps)))
+            if self.type == "field":
+                return [BaseField(self.__reader, index) for index in sli]
+            # Any of these returns defaults to particle Frames.
+            # If we get more types of trejaectories we have to add them here
+            # with an if statement as above.
+            return [BaseFrame(self.__reader, index) for index in sli]
+        if isinstance(item, Iterable):
+            if self.type == "field":
+                return [BaseField(self.__reader, index) for index in item]
+            return [BaseFrame(self.__reader, index) for index in item]
+        if isinstance(item, (int, np.integer)):
+            if self.type == "field":
                 return BaseField(self.__reader, item)
-            elif self.type=="particle":
-                return BaseFrame(self.__reader, item)
             return BaseFrame(self.__reader, item)
-        else:
-            raise KeyError(
-                '''BaseTrajectory: Invalid key. Only integer values, 1D lists
-                and arrays, and slices are allowed.'''
-            )
+        raise KeyError('''BaseTrajectory: Invalid key. Only integer values,
+        1D lists and arrays, and slices are allowed.'''
+                       )
+
     def __iter__(self):
-        for i in range(len(self.__reader.steps)):
+        """Iterate over all frames of the trajectory."""
+        for i, _ in enumerate(self.__reader.steps):
             yield self[i]
+
     def __next__(self):
         pass
+
     def __len__(self):
         return len(self.__reader.steps)
+
     def add_author_info(
             self, author: str, key: str, value: int | float | str) -> None:
         '''
@@ -1329,6 +1341,7 @@ def add_author_info(
             if author not in root['info']['authors'].keys():
                 root['info']['authors'].create_group(author)
             root['info']['authors'][author].attrs[key] = value
+
     def get_author_info(self, author: str) -> dict:
         r'''
         Returns all information for the given author.
@@ -1352,6 +1365,7 @@ def get_author_info(self, author: str) -> dict:
                 p = dict(a for a in root['info']['authors'][author].attrs.items())
                 return p
             return {}
+
     def delete_author_info(self, author: str, key: str | None = None) -> None:
         r'''
         Deletes all information (key=None) or specific information given by
@@ -1377,6 +1391,7 @@ def delete_author_info(self, author: str, key: str | None = None) -> None:
                 del root['info']['authors'][author]
             elif type(key)==str:
                 root['info']['authors'][author].attrs.__delitem__(key)
+
     @property
     def authors(self) -> list[str]:
         r'''
@@ -1396,6 +1411,7 @@ def authors(self) -> list[str]:
                 keys = list(root['info']['authors'].keys())
                 return keys
             return []
+
     def add_software_info(self, key: str, value: str | int | float) -> None:
         r'''
         Add software information to the hdf5 trajectory file.
@@ -1418,6 +1434,7 @@ def add_software_info(self, key: str, value: str | int | float) -> None:
             if 'software' not in root['info'].keys():
                 root['info'].create_group('software')
             root['info']['software'].attrs[key] = value
+
     def delete_software_info(self, key: str | None = None) -> None:
         r'''
         Deletes all software information (key=None) or specific information
@@ -1442,6 +1459,7 @@ def delete_software_info(self, key: str | None = None) -> None:
                     root['info']['software'].attrs.__delitem__(key)
             else:
                 root['info']['software'].attrs.__delitem__(key)
+
     @property
     def software(self) -> dict:
         r'''
@@ -1731,7 +1749,16 @@ def dt(self, x):
         if type(x) == float:
             self.__reader.dt = x
     @property
-    def dim(self):
+    def dim(self) -> int:
+        '''
+        Spatial dimension of the simnulation.
+
+        Returns
+        -------
+        x : int
+            Spatial dimension.
+
+        '''
         return self.__reader.d
     @property
     def savedir(self):
@@ -1783,8 +1810,12 @@ def __init__(self) -> None:
     def __getitem__(self, key: str):
         return getattr(self, key)
 
-    def keys(self):
-        return [name for (name, value) in inspect.getmembers(
+    def keys(self) -> list[str]:
+        """The keys to the evaluation object.
+        
+        Used so Evaluation-objects can be used as dictionaries.
+        """
+        return [name for (name, _) in inspect.getmembers(
             type(self), lambda x: isinstance(x, property)
         )]
 
diff --git a/amep/continuum.py b/amep/continuum.py
index fc98996..ed0397e 100644
--- a/amep/continuum.py
+++ b/amep/continuum.py
@@ -32,6 +32,7 @@
 # =============================================================================
 # IMPORT MODULES
 # =============================================================================
+from packaging.version import Version
 import numpy as np
 import scipy.fft as fft
 
@@ -782,8 +783,12 @@ def cluster_properties(
         # set all grid points outside the cluster to zero
         y[labels != i] = 0
         # integral
-        N = np.trapz(y, x=X, axis=1)
-        N = np.trapz(N, x=Y[:, 0])
+        if Version(np.__version__)<Version("2.0.0"):
+            N = np.trapz(y, x=X, axis=1)
+            N = np.trapz(N, x=Y[:, 0])
+        else:
+            N = np.trapezoid(y, x=X, axis=1)
+            N = np.trapezoid(N, x=Y[:, 0])
         sizes.append(N)
 
     # calculate the radius of gyration and linear extension
diff --git a/amep/evaluate.py b/amep/evaluate.py
index fc20c2e..8120d7c 100644
--- a/amep/evaluate.py
+++ b/amep/evaluate.py
@@ -31,6 +31,7 @@
 # =============================================================================
 # IMPORT MODULES
 # =============================================================================
+from packaging.version import Version
 import warnings
 import numpy as np
 
@@ -50,6 +51,7 @@
 from .trajectory import ParticleTrajectory, FieldTrajectory
 from . import thermo as thermo
 
+from typing import Callable
 
 # =============================================================================
 # GENERAL FUNCTION
@@ -58,7 +60,8 @@ class Function(BaseEvaluation):
     """Apply a user-defined function to a trajectory.
     """
 
-    def __init__(self, traj, func, skip=0.0, nav=10, **kwargs):
+    def __init__(self, traj: ParticleTrajectory | FieldTrajectory, 
+        func: Callable, skip: float = 0.0, nav: int = 10, **kwargs):
         r'''Calculate a given function for a trajectory.
 
         Parameters
@@ -129,7 +132,7 @@ def __init__(self, traj, func, skip=0.0, nav=10, **kwargs):
             nr=self.__nav, indices=True, **kwargs)
         
         self.__times = self.__traj.times[self.__indices]
-    def __compute(self, ind, **kwargs):
+    def __compute(self, ind: int, **kwargs):
         r'''
         Calculation for a single frame,
 
@@ -1690,8 +1693,14 @@ class PosOrderCor(BaseEvaluation):
     """
 
     def __init__(
-            self, traj, grtdata=None, sxydata=None, skip=0.0, nav=10,
-            k0=None, order='hexagonal', dk=4.0, ptype=None, other=None,
+            self, traj: ParticleTrajectory, 
+            grtdata: BaseEvaluation | None = None,
+            sxydata: BaseEvaluation | None = None, 
+            skip: float = 0.0, nav: int = 10,
+            k0: list[np.ndarray, ...] | None = None, 
+            order: str = 'hexagonal', 
+            dk: float = 4.0, ptype: int | None = None, 
+            other: int | None = None,
             **kwargs) -> None:
         r'''
         Calculate the positional order correlation function.
@@ -1731,12 +1740,12 @@ def __init__(
         ----------
         traj : Traj
             Trajectory object with simulation data.
-        grtdata : pd.DataFrame, optional
+        grtdata : BaseEvaluation, optional
             Pair correlation function g(r,theta) as obtained from
             amep.evaluate.PCFangle. If None, the angular pair
             correlation function is calculated within this method.
             The default is None.
-        sxydata : pd.DataFrame, optional
+        sxydata : BaseEvaluation, optional
             2d structure factor as obtained from amep.evaluate.SF2d.
             If None, the 2d structure factor is calculated within this
             method. The default is None.
@@ -1840,16 +1849,22 @@ def __compute(self):
                                               (self.__k0[i][0]*np.cos(self.__theta) +\
                                               self.__k0[i][1]*np.sin(self.__theta)))
 
-            ck += np.trapz(integrand, x=self.__theta, axis=1)
-        
+            if Version(np.__version__) < Version("2.0.0"):
+                ck += np.trapz(integrand, x=self.__theta, axis=1)
+            else:
+                ck += np.trapezoid(integrand, x=self.__theta, axis=1)
+
         ck = np.abs(ck/len(self.__k0))
 
         # normalization
-        norm = np.trapz(self.__grt, x=self.__theta, axis=1)
+        if Version(np.__version__) < Version("2.0.0"):
+            norm = np.trapz(self.__grt, x=self.__theta, axis=1)
+        else:
+            norm = np.trapezoid(self.__grt, x=self.__theta, axis=1)
 
         return ck/norm
-    
-    
+
+
     def __getk(self):
         r'''
         Calculates the k vectors corresponding to the first peaks of 
@@ -2518,8 +2533,10 @@ class Psi6dist(BaseEvaluation):
     """
 
     def __init__(
-            self, traj, skip=0.0, nav=10, nbins=None, ptype=None,
-            other=None, **kwargs) -> None:
+            self, traj: ParticleTrajectory, skip: float = 0.0, 
+            nav: int = 10, nbins: int = 50, 
+            ptype: int | None = None, other: int | None = None, 
+            **kwargs) -> None:
         r'''
         Calculate the distribution of the :math:`\Psi_6`.
 
@@ -2560,7 +2577,7 @@ def __init__(
 
         Parameters
         ----------
-        traj : Traj
+        traj : ParticleTrajectory
             Trajectory object.
         skip : float, optional
             Skip this fraction at the beginning of the trajectory. 
@@ -2568,7 +2585,7 @@ def __init__(
         nav : int, optional
             Number of frames to use for the average. The default is 10.
         nbins : int, optional
-            Number of bins. The default is None.
+            Number of bins. The default is 50.
         ptype : float, optional
             Particle type. If None, all particles are used.
             The default is None.
@@ -2614,9 +2631,6 @@ def __init__(
         self.__other = other
         self.__kwargs = kwargs
         
-        if self.__nbins is None:
-            self.__nbins = int(self.__traj[0].n(ptype=self.__ptype)/100)
-        
         self.__frames, res, self.__indices = average_func(
             self.__compute, np.arange(self.__traj.nframes), skip=self.__skip,
             nr=self.__nav, indices=True)
@@ -2735,8 +2749,10 @@ class VelDist(BaseEvaluation):
     """
     
     def __init__(
-            self, traj, skip=0.0, nav=10, nbins=None, ptype=None, vmin=None,
-            vmax=None, v2min=None):
+            self, traj, skip: float = 0.0, nav: int = 10,
+            nbins: int = 50, ptype: int | None = None,
+            vmin: float | None = None, vmax: float | None = None,
+            v2min: float | None = None) -> None:
         r'''
         Calculate the distribution of velocities.
 
@@ -2745,6 +2761,14 @@ def __init__(
         and its square :math:`v^2`. It also
         takes an average over several frames (= time average).
 
+        For the :math:`v^2` distribution, logarithmic
+        bins are used. Therefore `v2min`:math:`\ge 0`
+        needs to be ensured. For the maximum value,
+        `vmax` is used as :math:`3v_{max}^2`.
+        Analogously for the maximum of :math:'v',
+        where :math:`\sqrt{3v_{max}}` is used. For the minimum
+        of :math:'v', the square root of `v2min` is used.
+
         Parameters
         ----------
         traj : Traj
@@ -2755,9 +2779,22 @@ def __init__(
         nav : int, optional
             Number of frames to use for the average. The default is 10.
         nbins : int, optional
-            Number of bins. The default is None.
+            Number of bins. The default is 50.
         ptype : float, optional
             Particle type. The default is None.
+        vmin : float | None, optional
+            Minimum value for the histogram in each spatial dimension
+            :math:`v_x, v_y, v_z`.
+            If None, then the minimum value of the last frame will be used
+        vmax : float | None, optional
+            Maximum value for the histogram in each spatial dimension
+            :math:`v_x, v_y, v_z`.
+            If None, then the maximum value of the last frame will be used
+        v2min : float | None, optional
+            Minimum value for the velocity-squared histogram.
+            This value has to be :math:`\ge 0` due to the use of
+            logarithmic bins for the :math:`v^2` distribution.
+            If None, then the minimum value of the last frame will be used.
 
         Returns
         -------
@@ -2793,9 +2830,6 @@ def __init__(
         self.__vmax  = vmax
         self.__v2min = v2min
         
-        if self.__nbins is None:
-            self.__nbins = int(self.__traj[0].n(ptype=self.__ptype)/100)
-            
         if self.__vmin is None:
             self.__vmin = np.min(self.__traj[-1].velocities(ptype=self.__ptype))
             
@@ -2852,12 +2886,22 @@ def __compute(self, ind):
         
         v2 = np.sum(vel**2, axis=1)
         
+        
         v2hist, v2bins = distribution(
             v2, nbins=self.__nbins, logbins=True, xmin=self.__v2min,
-            xmax=(3*self.__vmin**2+3*self.__vmax**2)/2)
+            xmax=3*self.__vmax**2)
         vhist, vbins = distribution(
-            np.sqrt(v2), nbins=self.__nbins, xmin=0,
-            xmax=np.sqrt((3*self.__vmin**2+3*self.__vmax**2)/2))
+            np.sqrt(v2), nbins=self.__nbins, xmin=np.sqrt(self.__v2min),
+            xmax=np.sqrt(3*self.__vmax**2))
+
+
+# (KD, 2024.07.16) was:
+        # v2hist, v2bins = distribution(
+        #     v2, nbins=self.__nbins, logbins=True, xmin=self.__v2min,
+        #     xmax=(3*self.__vmin**2+3*self.__vmax**2)/2)
+        # vhist, vbins = distribution(
+        #     np.sqrt(v2), nbins=self.__nbins, xmin=0,
+        #     xmax=np.sqrt((3*self.__vmin**2+3*self.__vmax**2)/2))
         
         return xhist, xbins, yhist, ybins, zhist, zbins, vhist, vbins, v2hist,\
                v2bins 
@@ -3191,15 +3235,26 @@ def __init__(
             if self.__xmax is None:
                 if self.__use_density:
                     # integrated density
-                    val_total_x = np.trapz(
-                        self.__traj[-1].data(self.__ftype),
-                        x = self.__traj[-1].grid[0],
-                        axis=1
-                    )
-                    val_total = np.trapz(
-                        val_total_x,
-                        x = self.__traj[-1].grid[1][:, 0]
-                    )
+                    if Version(np.__version) < Version("2.0.0"):
+                        val_total_x = np.trapz(
+                            self.__traj[-1].data(self.__ftype),
+                            x=self.__traj[-1].grid[0],
+                            axis=1
+                        )
+                        val_total = np.trapz(
+                            val_total_x,
+                            x=self.__traj[-1].grid[1][:, 0]
+                        )
+                    else:
+                        val_total_x = np.trapezoid(
+                            self.__traj[-1].data(self.__ftype),
+                            x=self.__traj[-1].grid[0],
+                            axis=1
+                        )
+                        val_total = np.trapezoid(
+                            val_total_x,
+                            x=self.__traj[-1].grid[1][:, 0]
+                        )
                     self.__xmax = val_total
                 else:
                     # total number of grid points
@@ -3412,7 +3467,7 @@ def __init__(
 
         Parameters
         ----------
-        traj : Traj
+        traj : ParticleTrajectory | FieldTrajectory
             Trajectory object.
         skip : float, optional
             Skip this fraction at the beginning of the trajectory.
@@ -4528,7 +4583,7 @@ def __init__(
 
         Parameters
         ----------
-        traj : Traj
+        traj : ParticleTrajectory
             Trajectory object with simulation data.
         skip : float, default=0.0
             Skip this fraction at the beginning of
@@ -4670,7 +4725,7 @@ def __init__(
 
         Parameters
         ----------
-        traj : Traj
+        traj : ParticleTrajectory
             Trajectory object with simulation data.
         mass : float | np.ndarray
             Function to be used for analysis.
@@ -4815,7 +4870,7 @@ def __init__(self, traj: ParticleTrajectory,
 
         Parameters
         ----------
-        traj : Traj
+        traj : ParticleTrajectory
             Trajectory object with simulation data.
         mode : "total" or "particle"
             How to return the energy, c.f. thermo.py
diff --git a/amep/functions.py b/amep/functions.py
index c6ada4a..b434e83 100644
--- a/amep/functions.py
+++ b/amep/functions.py
@@ -10,8 +10,7 @@
 # This program is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
+# GNU General Public License for more details.  #
 # You should have received a copy of the GNU General Public License
 # along with this program. If not, see <https://www.gnu.org/licenses/>.
 #
@@ -23,17 +22,19 @@
 
 .. module:: amep.functions
 
-The AMEP module :mod:`amep.functions` contains regularly used functions to fit 
+The AMEP module :mod:`amep.functions` contains regularly used functions to fit
 to data as well as fitting methods.
+The naming covention is that the fit classes are named in CamelCase while the
+mathematical functions that can be fitted are name all lowercase.
 
 """
 # =============================================================================
 # IMPORT MODULES
 # =============================================================================
-import numpy as np
 import inspect
+import numpy as np
 
-from scipy.special import erf, erfc
+from scipy.special import erf, erfc, gamma
 from .base import BaseFunction
 
 # =============================================================================
@@ -43,27 +44,35 @@ def gaussian(
         x: np.ndarray, mu: float = 0.0, sig: float = 1.0, offset: float = 0.0,
         A: float = 1.0, normalized: bool = False) -> np.ndarray:
     r'''
-    1D Gaussian.
+    Gaussian Bell curve.
+
+    Equivalent to the probability density function of a normal distribution.
+    By the central limit theorem this is a good guess for most peak shapes
+    that arise from many random processes.
+    The function is given by
+
+    .. math::
+        g(x)=A\exp\left(-\frac{\left(x-\mu\right)^2}{\sigma^2}\right)+b.
 
     Parameters
     ----------
     x : np.ndarray
-        x values.
+        :math:`x` values.
     mu : float
-        Mean value.
+        Mean value :math:`\mu`.
     sig : float
-        Standard deviation.
+        Standard deviation :math:`\sigma`.
     offset : float, optional
-        Shifts the Gaussian in y direction. The default is 0.0.
+        Shifts the Gaussian by :math:`b` in y direction. The default is 0.0.
     A : float, optional
-        Amplitude. The default is 1.0.
+        Amplitude :math:`A`. The default is 1.0.
     normalized : bool, optional
         If True, the Gaussian is normalized to unit area. The default is False.
 
     Returns
     -------
     np.ndarray
-        g(x).
+        g(x)
 
     Examples
     --------
@@ -84,7 +93,7 @@ def gaussian(
     '''
     if normalized:
         A = 1.0/(np.sqrt(2*np.pi) * sig)
-    
+
     return A*np.exp(-((x-mu)**2)/2/sig**2) + offset
 
 
@@ -93,34 +102,55 @@ def gaussian2d(
         muy: float, sigx: float, sigy: float, theta: float, offset: float
         ) -> np.ndarray:
     r'''
-    2D Gaussian.
-    
+    2D Gaussian Bell curve.
+
+
+    Equivalent to the probability density function of a normal distribution
+    in two dimensions.
+    By the central limit theorem this is a good guess for most peak shapes
+    that arise from many random processes.
+    This parametrization can include correlations
+    via the angle variable :math:`\theta`.
+    The function is given by
+
+    .. math::
+        g(x)= A\exp\left(-\left(\vec{x}-\vec{\mu}\right)^T
+                R(\Theta)^{-1}\sigma^{-2}R(\Theta)
+                \left(\vec{x}-\vec{\mu}\right)
+        \right)+b
+
+    where :math:`\vec{x}` is the vector composed the x and y coordinates,
+    :math:`\vec{\mu}` is the mean vector composed
+    of :math:`\mu_x` and :math:`\mu_y` and
+    :math:`\sigma^{-2}` is the diagonal matrix with the inverse
+    variances :math:`\sigma_x^{-2}` and :math:`\sigma_y^{-2}` as entries.
+
     Parameters
     ----------
     data : tuple
         tuple (x,y) of x and y values where x and y
-        are 1D np.ndarrays
+        are 1D np.ndarrays.
     A : float
-        amplitude (float)
+        amplitude.
     mux : float
-        mean in x direction (float)
+        mean :math:`\mu_x` in x direction.
     muy : float
-        mean in y direction (float)
+        mean :math:`\mu_y` in y direction.
     sigx : float
-        sqrt(variance) in x direction (float)
+        standard deviation :math:`\sigma_x` in x direction.
     sigy : float
-        sqrt(variance) in y direction (float)
+        standard deviation :math:`\sigma_y`  in y direction.
     theta : float
-        orientation angle of the polar axis (float)
+        Orientation angle :math:`\Theta` of the polar axis.
     offset : float
-        offset (float)
-        
+        offset :math:`b`. Shifts the output value linearly.
+
     Returns
     -------
-    np.ndarray
+    g(x): np.ndarray
         1D array of floats (flattended 2D array!)
-        
-        
+
+
     Examples
     --------
     >>> import amep
@@ -142,14 +172,15 @@ def gaussian2d(
       :align: center
 
     '''
-    (x,y) = data_tuple
+    (x, y) = data_tuple
     mux = float(mux)
-    muy = float(muy)    
+    muy = float(muy)
     a = (np.cos(theta)**2)/(2*sigx**2) + (np.sin(theta)**2)/(2*sigy**2)
     b = -(np.sin(2*theta))/(4*sigx**2) + (np.sin(2*theta))/(4*sigy**2)
     c = (np.sin(theta)**2)/(2*sigx**2) + (np.cos(theta)**2)/(2*sigy**2)
-    g = offset + A*np.exp( - (a*((x-mux)**2) + 2*b*(x-mux)*(y-muy) 
-                            + c*((y-muy)**2)))
+    g = offset + A*np.exp(a*((x-mux)**2)
+                          + 2*b*(x-mux)*(y-muy)
+                          + c*((y-muy)**2))
     return g.ravel()
 
 
@@ -204,25 +235,26 @@ def __init__(self) -> None:
           :align: center
 
         '''
-        super(Gaussian, self).__init__(3)
-        
+        super().__init__(3)
+
         self.name = 'Gaussian'
         self.keys = ['mu', 'sig', 'a']
-        
+
     def f(
             self, p: list | np.ndarray,
             x: float | np.ndarray) -> float | np.ndarray:
         r'''
         Non-normalized Gaussian function in 1d of the form
-        
+
         .. math::
-            
+
             g(x) = a\exp\left\lbrace\frac{(x-\mu)^2}{2\sigma^2}\right\rbrace.
 
         Parameters
         ----------
         p : list or np.ndarray
-            Parameters $\mu$, $\sigma$, and $a$ of the Gaussian function.
+            Parameters :math:`\mu`, :math:`\sigma`,
+            and :math:`a` of the Gaussian function.
         x : float or np.ndarray
             Value(s) at which the function is evaluated.
 
@@ -236,7 +268,14 @@ def f(
 
 
 class NormalizedGaussian(BaseFunction):
-    """Normalized one-dimensional Gaussian.
+    r"""Normalized one-dimensional Gaussian.
+
+    Has the mathematical form
+
+        .. math::
+
+            g(x) = \frac{1}{\sqrt{2\pi}\sigma}
+                \exp\left\lbrace\frac{(x-\mu)^2}{2\sigma^2}\right\rbrace.
     """
     def __init__(self) -> None:
         r'''
@@ -287,25 +326,27 @@ def __init__(self) -> None:
           :align: center
 
         '''     
-        super(NormalizedGaussian, self).__init__(2)
-        
+        super().__init__(2)
+
         self.name = 'NormalizedGaussian'
-        self.keys = ['mu','sig']
-        
+        self.keys = ['mu', 'sig']
+
     def f(
           self, p: list | np.ndarray,
           x: float | np.ndarray) -> float | np.ndarray:
         r'''
         Normalized one-dimensional Gaussian function of the form
 
-        .. math..:
+        .. math::
 
-            g(x) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left\lbrace\right\frac{(x-\mu)^2}{2\sigma^2}\rbrace.
+            g(x) = \frac{1}{\sqrt{2\pi}\sigma}
+                \exp\left\lbrace\frac{(x-\mu)^2}{2\sigma^2}\right\rbrace.
 
         Parameters
         ----------
         p : list or np.ndarray
-            Parameters $\mu$ and $\sigma$ of the normalized Gaussian function.
+            Parameters :math:`\mu` and :math:`\sigma`
+            of the normalized Gaussian function.
         x : float | np.ndarray
             Value(s) at which the function is evaluated.
 
@@ -317,8 +358,8 @@ def f(
         '''
         A = 1.0 / (np.sqrt(2*np.pi) * p[1])
         return A * np.exp(-((x-p[0])**2) / 2 / p[1]**2)
-    
-    
+
+
 class ExGaussian(BaseFunction):
     """Exponentially modified Gaussian.
     """
@@ -330,22 +371,23 @@ def __init__(self) -> None:
         -----
         The exponentially modified Gaussian function has three parameters
         :math:`\lambda, \mu`, and :math:`\sigma`, and it is defined as
-        
-        .. math..:
-            
-            g_{\rm ex}(x) = \frac{\lambda}{2} e^{\frac{\lambda}{2} (2 \mu + \lambda \sigma^2 - 2 x)}
-             {\rm erfc} \left(\frac{\mu + \lambda \sigma^2 - x}{ \sqrt{2} \sigma}\right),
-            
-        where :math:`{\rm erfc}` is the complementary error function [1]_. The 
+
+        .. math::
+
+            g_{\rm ex}(x) = \frac{\lambda}{2} e^{\frac{\lambda}{2}
+             (2 \mu + \lambda \sigma^2 - 2 x)}{\rm erfc}
+             \left(\frac{\mu + \lambda \sigma^2 - x}{ \sqrt{2} \sigma}\right),
+
+        where :math:`{\rm erfc}` is the complementary error function [1]_. The
         parameters are ordered as follows:
-        
+
             p[0] : :math:`\lambda`
-            
+
             p[1] : :math:`\mu`
-            
+
             p[2] : :math:`\sigma`
-            
-        The mean value of the distribution is given by :math:`\mu+1/\lambda` 
+
+        The mean value of the distribution is given by :math:`\mu+1/\lambda`
         and the standard deviation by :math:`\sqrt{\sigma^2 + 1/\lambda^2}`.
 
         References
@@ -393,11 +435,11 @@ def __init__(self) -> None:
           :align: center
 
         '''
-        super(ExGaussian, self).__init__(3)
-        
+        super().__init__(3)
+
         self.name = 'Exponentially modified Gaussian'
         self.keys = ['lambda', 'mu', 'sigma']
-        
+
     def f(
             self, p: list | np.ndarray,
             x: float | np.ndarray) -> float | np.ndarray:
@@ -407,8 +449,11 @@ def f(
 
         .. math..:
 
-            g_{\rm ex}(x) = \frac{\lambda}{2} e^{\frac{\lambda}{2} (2 \mu + \lambda \sigma^2 - 2 x)}
-             {\rm erfc} \left(\frac{\mu + \lambda \sigma^2 - x}{ \sqrt{2} \sigma}\right)
+            g_{\rm ex}(x) = \frac{\lambda}{2} e^{\frac{\lambda}{2}
+                                (2 \mu + \lambda \sigma^2 - 2 x)}
+                                {\rm erfc} \left(
+                                \frac{\mu + \lambda \sigma^2 - x}{
+                                \sqrt{2} \sigma}\right)
 
         Parameters
         ----------
@@ -429,11 +474,11 @@ def f(
     @property
     def mean(self) -> float:
         return self.params[1] + 1.0/self.params[0]
-    
+
     @property
     def std(self) -> float:
         return np.sqrt(self.params[2]**2 + 1.0/self.params[0]**2)
-    
+
 
 class SkewGaussian(BaseFunction):
     """Skewed Gaussian distribution.
@@ -505,11 +550,11 @@ def __init__(self) -> None:
           :align: center
 
         '''
-        super(SkewGaussian, self).__init__(3)
-        
+        super().__init__(3)
+
         self.name = 'Skew Normal Distribution'
         self.keys = ['mu', 'sigma', 'alpha']
-        
+
     def f(
             self, p: list | np.ndarray,
             x: float | np.ndarray) -> float | np.ndarray:
@@ -538,25 +583,45 @@ def f(
         '''
         return 1.0/np.sqrt(2*np.pi)/p[1] * np.exp(-(x-p[0])**2/2/p[1]**2) *\
             (1 + erf(p[2]*(x-p[0])/np.sqrt(2)/p[1]))
-    
+
     @property
     def mean(self) -> float:
         delta = self.params[2]/np.sqrt(1+self.params[2]**2)
         return self.params[0] + self.params[1]*delta*np.sqrt(2/np.pi)
-    
+
     @property
     def std(self) -> float:
         delta = self.params[2]/np.sqrt(1+self.params[2]**2)
         return np.sqrt(self.params[1]**2*(1-2*delta**2/np.pi))
-    
-    
+
+
 class MaxwellBoltzmann(BaseFunction):
     """Maxwell-Boltzmann distribution.
     """
 
-    def __init__(self, d=1, m=None):
+    def __init__(self,
+                 d: int = 1,
+                 m: float = 1.):
         r'''Maxwell-Boltzmann velocity distribution.
 
+        Initializes a function object of a Maxwell-Boltzmann distribution.
+        This distribution describes velocities of free particles in
+        thermal equilibrium.
+        It is usually derived in three dimensions where it is
+
+        .. math::
+
+            f(v) = {\left[\frac{m}{2\pi{}k_bT}\right]}^\frac{3}{2}4\pi{}
+            v^2\exp\left(-\frac{mv^2}{k_bT}\right)
+
+        in arbitrary spatial dimensions :math:`d` it is given by
+
+        .. math::
+
+            f(v) =  \frac{2}{\Gamma\left(\frac{d}{2}\right)}
+                    {\left[\frac{m}{k_bT}\right]}^\frac{d}{2}
+                    v^{d-1}\exp\left(-\frac{mv^2}{k_bT}\right)
+
         Parameters
         ----------
         d : int, optional
@@ -593,36 +658,47 @@ def __init__(self, d=1, m=None):
           :align: center
 
         '''
-        super(MaxwellBoltzmann, self).__init__(1)
+        super().__init__(1)
 
         self.__d = d
         self.__m = m
 
-        if self.__m is None:
-            self.__m = 1
-
         self.name = 'Maxwell-Boltzmann Distribution'
         self.keys = ['kbT']
 
-    def f(self, p, x):
+    def f(self,
+          p: list[float] | np.ndarray,
+          x: float | np.ndarray
+          ) -> float | np.ndarray:
         r"""Calculate the value of the Maxwell-Boltzmann distribution.
 
+        This distribution describes velocities of free particles in
+        thermal equilibrium.
+        It is usually derived in three dimensions where it is
+
+        .. math::
+
+            f(v) = {\left[\frac{m}{2\pi{}k_bT}\right]}^\frac{3}{2}4\pi{}
+            v^2\exp\left(-\frac{mv^2}{k_bT}\right)
+
+        in arbitrary spatial dimensions :math:`d` it is given by
+
+        .. math::
+
+            f(v) =  \frac{2}{\Gamma\left(\frac{d}{2}\right)}
+                    {\left[\frac{m}{k_bT}\right]}^\frac{d}{2}
+                    v^{d-1}\exp\left(-\frac{mv^2}{k_bT}\right)
+
         Parameters
         ----------
         p: list | np.ndarray
             Parameters :math:`k_{\rm B}T`.
         x: float | np.ndarray
-            Value(s) at which the function is evaluated.
+            Velocities(s) :math:`v` at which the function is evaluated.
 
         """
-        if self.__d == 1:
-            return (self.__m/(2*np.pi*p[0]))**(1./2.) * np.exp(
-                    -self.__m*x**2/(2*p[0]))
-        elif self.__d == 2:
-            return (self.__m/p[0]) * x * np.exp(-self.__m*x**2/(2*p[0]))
-        elif self.__d == 3:
-            return (self.__m/(2*np.pi*p[0]))**(3./2.) * 4*np.pi*x**2 * np.exp(
-                    -self.__m*x**2/(2*p[0]))
+        return ((self.__m/(2*p[0]))**(self.__d/2)/gamma(self.__d/2)*2 *
+                np.exp(-self.__m*x**2/(2*p[0]))*x**(self.__d-1))
 
     @property
     def mean(self):
@@ -636,15 +712,36 @@ def std(self):
 class Gaussian2d(BaseFunction):
     """Two-dimensional Gaussian.
     """
-    
+
     def __init__(self):
         r'''
-        Two-dimensional Gaussian function.
+        Initializes a function object of Two-dimensional Gaussian function.
+
+        Equivalent to the probability density function of a normal distribution
+        in two dimensions.
+        By the central limit theorem this is a good guess for most peak shapes
+        that arise from many random processes.
+        This parametrization can include correlations
+        via the angle variable :math:`\theta`.
+        The function is given by
+
+        .. math::
+            g(x)= A\exp\left(-\left(\vec{x}-\vec{\mu}\right)^T
+                    R(\Theta)^{-1}\sigma^{-2}R(\Theta)
+                    \left(\vec{x}-\vec{\mu}\right)
+            \right)+b
+
+        where :math:`\vec{x}` is the vector composed the x and y coordinates,
+        :math:`\vec{\mu}` is the mean vector composed
+        of :math:`\mu_x` and :math:`\mu_y` and
+        :math:`\sigma^{-2}` is the diagonal matrix with the inverse
+        variances :math:`\sigma_x^{-2}` and :math:`\sigma_y^{-2}` as entries.
+
 
         Returns
         -------
         None.
-        
+
         Examples
         --------
         >>> import amep
@@ -666,51 +763,64 @@ def __init__(self):
         .. image:: /_static/images/functions/functions-Gaussian2d.png
           :width: 400
           :align: center
-        
+
         '''
-        super(Gaussian2d, self).__init__(6)
-    
+        super().__init__(6)
+
         self.name = 'Two-dimensional Gaussian'
         self.keys = ['a', 'mux', 'muy', 'sigx', 'sigy', 'theta']
-        
+
     def f(self, p, x):
         r'''
         2D Gaussian.
-        
+        The function is given by
+
+        .. math::
+            g(x)= A\exp\left(-\left(\vec{x}-\vec{\mu}\right)^T
+                    R(\Theta)^{-1}\sigma^{-2}R(\Theta)
+                    \left(\vec{x}-\vec{\mu}\right)
+            \right)+b
+
+        where :math:`\vec{x}` is the vector composed the x and y coordinates,
+        :math:`\vec{\mu}` is the mean vector composed
+        of :math:`\mu_x` and :math:`\mu_y` and
+        :math:`\sigma^{-2}` is the diagonal matrix with the inverse
+        variances :math:`\sigma_x^{-2}` and :math:`\sigma_y^{-2}` as entries.
+
         Parameters
         ----------
         p: list
-            Parameters.
+            Parameters :math:`(A,\mu_x,\mu_y,\sigma_x,\sigma_y,\Theta)`.
         x: np.ndarray
             x-values as 2d array of shape (N,2).
-            
+
         Returns
         -------
         np.ndarray
             1D array of shape (N,).
-            
-            
+
+
         Examples
         --------
         >>> 
         '''
-        X = x[:,0]
-        Y = x[:,1]
+        X = x[:, 0]
+        Y = x[:, 1]
         a = (np.cos(p[5])**2)/(2*p[3]**2) + (np.sin(p[5])**2)/(2*p[4]**2)
         b = -(np.sin(2*p[5]))/(4*p[3]**2) + (np.sin(2*p[5]))/(4*p[4]**2)
         c = (np.sin(p[5])**2)/(2*p[3]**2) + (np.cos(p[5])**2)/(2*p[4]**2)
-        g = p[0]*np.exp( - (a*((X-p[1])**2) + 2*b*(X-p[1])*(Y-p[2]) 
-                                + c*((Y-p[2])**2)))
+        g = p[0]*np.exp(- (a*((X-p[1])**2) + 2*b*(X-p[1])*(Y-p[2])
+                           + c*((Y-p[2])**2)))
         return g
-    
+
     @property
     def mean(self):
         pass
-    
+
     @property
     def std(self):
         pass
-    
+
 
 class Fit(BaseFunction):
     """Fitting a 1d function.
@@ -777,19 +887,23 @@ def __init__(self, g, **kwargs):
         len_diff = len(inspected_object[0]) - len(inspected_object[3])
         fct_call_params_defaults = inspected_object[3]
         fct_def_params = inspected_object[0][len_diff:]
-        # 2: Get user defined function to get all parameters and their default values, then get the parameters from
+        # 2: Get user defined function to get all parameters and
+        # their default values, then get the parameters from
         # the actual function call and categorize parameters
         # - their function call values if they are mentioned
-        # - their default values if they are mentioned as "=None" in the function call
+        # - their default values if they are mentioned as
+        # "=None" in the function call
         # - fit to them if neither of the above is true
 
         # Extract all params from def and call of fct
         self.fct_def_params = fct_def_params
-        fct_call_params = [i for i in kwargs.keys()]
-        fct_call_values = [i for i in kwargs.values()]
-        # Get the length difference between arguments and kwargs, this allows us to match name to value
+        fct_call_params = list(kwargs.keys())
+        fct_call_values = list(kwargs.values())
+        # Get the length difference between arguments and kwargs,
+        # this allows us to match name to value
 
-        # Determine which params have what values in fit and which are to be fitted
+        # Determine which params have what values
+        # in fit and which are to be fitted
         fit_params_input = []
         fit_params_to_determine = []
         for k, param in enumerate(fct_def_params):
@@ -807,7 +921,6 @@ def __init__(self, g, **kwargs):
         self.defaults = inspected_object[3]
         self.g = g
 
-
         self.fit_params_input = fit_params_input
         self.keys = fit_params_to_determine
 
diff --git a/amep/load.py b/amep/load.py
index eddeff8..d552a35 100644
--- a/amep/load.py
+++ b/amep/load.py
@@ -95,22 +95,27 @@ def traj(
     Minimal example for loading lammps trajectory from path `"data"`:
 
     >>> import amep
-    >>> traj = amep.load.traj("data")
-    >>> traj = amep.load.traj("data", reload=True) # data will be read again
+    >>> traj = amep.load.traj("/examples/data/lammps/")
+    >>> traj = amep.load.traj("/examples/data/lammps/", reload=True) # data will be read again
 
 
     Shortcuts to load an h5amep file directly:
 
-    >>> traj = amep.load.traj("data/traj.h5amep")
+    >>> traj = amep.load.traj("/examples/data/lammps.h5amep")
     >>> traj = amep.load.traj(
     ...     "data", trajfile="traj.h5amep", mode="h5amep"
     ... )
     >>> 
 
+    
+    Example for loading continuum data:
+
+    >>> traj = amep.load.traj("/examples/data/continuum/", mode="field", dumps="field_*")
+
 
     Fix for working with remote files and insufficient access rights:
-    Save the h5amep file locally while importing/loading the simulation data
-    for the first time:
+    This can be solved by saving the h5amep file locally while
+    importing/loading the simulation data from the remote directory.
 
     >>> traj = amep.load.traj(
     ...     "/path/to/sftp:host=192.168.255.255/remote/directory/",
@@ -221,6 +226,7 @@ def traj(
             **kwargs
         )
         return FieldTrajectory(reader)
+
     # here one has to check both the amep version with which the file has been
     # created (reader.version) and the data type (particles or fields) -
     # the latter is needed to decide whether a ParticleTrajectory or a
diff --git a/amep/order.py b/amep/order.py
index 1a1dd06..01a8f75 100644
--- a/amep/order.py
+++ b/amep/order.py
@@ -703,7 +703,7 @@ def voronoi_density(
         Coordinates of the particles for the voronoi tessellation.
         If None, the user has to specify the Voronoi tessellation with the
         keyword `vor` (and ids for pbc=True).
-    box_boundary : np.ndarray of shape (3,2) or None
+    box_boundary : np.ndarray of shape (3,2) | None, optional
         Boundary of the simulation box in the form of
         `np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]])`.
         If None, the user has to specify the Voronoi tessellation with the
@@ -1017,13 +1017,13 @@ def nearest_neighbors(
     Parameters
     ----------
     coords : np.ndarray of shape (N,3)
-        Array of coordinates/points to search for neighbors.
+        Array of coordinates to search for neighbors of.
     box_boundary : np.ndarray of shape (3,2)
         Boundary of the simulation box in the form of
         `np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]])`.
     other_coords : np.ndarray of shape (M,3) or None, optional
-        Array of points to search for neighbors of. The default is None (use 
-        coords).
+        Array of coordinates of the particles (neighbors).
+        The default is None (use coords).
     rmax : float, optional
         Distance cutoff. Only neighbors with distance smaller than the cutoff
         are returned. The default is 1.0.
diff --git a/amep/pbc.py b/amep/pbc.py
index f67eb9f..024b7c0 100644
--- a/amep/pbc.py
+++ b/amep/pbc.py
@@ -369,9 +369,9 @@ def mirror_points(
     ----------
         coordinates: np.ndarray
             coordinate frame (Nx3 array of floats; required)
-        box_boundary : np.ndarray
-            Boundaries of the simulation box given in the form
-            np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]]).
+	box_boundary : np.ndarray of shape (3,2)
+            Boundary of the simulation box in the form of
+            `np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]])`.
         width: float or None
             Width of the periodic images relative to the box dimensions.
             None means the original and all periodic images are returned;
@@ -550,8 +550,8 @@ def pbc_diff(
         First vector.
     v2 : np.ndarray
         Second vector. 
-    box_boundary : np.ndarray
-        Boundaries of the simulation box given in the form
+   box_boundary : np.ndarray of shape (3,2)
+        Boundary of the simulation box in the form of
         `np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]])`.
     pbc : bool, optional
         If True, periodic boundary conditions are considered.
@@ -581,8 +581,8 @@ def pbc_diff_rect(v1, v2, box_boundary):
         First vector.
     v2 : np.ndarray
         Second vector. 
-    box_boundary : np.ndarray
-        Boundaries of the simulation box given in the form
+    box_boundary : np.ndarray of shape (3,2)
+        Boundary of the simulation box in the form of
         `np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]])`.
         
     Returns
@@ -912,7 +912,7 @@ def __dis_chunk(
         Coordinates to which the distances should be calculated.
     otree : KDTree
         KDTree of coordinates from which the distances are calculated.
-    box_boundary : np.ndarray
+    box_boundary : np.ndarray of shape (3,2)
         Boundary of the simulation box in the form of
         `np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]])`.
     pbc : bool
diff --git a/amep/plot.py b/amep/plot.py
index 3e1f6cf..1e01628 100644
--- a/amep/plot.py
+++ b/amep/plot.py
@@ -539,7 +539,7 @@ def new(figsize: tuple[float, float] = None,
         **kwargs
         ) -> tuple[plt.Figure, mpl.axes.Axes | np.ndarray]:
     r'''
-    Create a new matplotlib.pyplot figure.
+    Create a new matplotlib.pyplot figure (uses matplotlib.pyplot.subplots).
 
     Parameters
     ----------
@@ -579,6 +579,14 @@ def new(figsize: tuple[float, float] = None,
     axes : AxesSubplot
         matplotlib.pyplot axes objects..
 
+    Examples
+    --------
+    >>> import amep
+    >>> import numpy as np
+    >>> fig, axs = amep.plot.new()
+    >>> x = np.linspace(0,10,1000)
+    >>> axs.plot(x, np.sin(x))
+
     '''
     gridspec = {}
     if width_ratios:
@@ -875,9 +883,9 @@ def box(axis: mpl.axes.Axes, box_boundary: np.ndarray, **kwargs) -> None:
     ----------
     axis : AxisSubplot
         Matplotlib.pyplot AxisSubplot object.
-    box_boundary : np.ndarray
-        Boundaries of the simulation box.
-        [[xmin, xmax], [ymin, ymax], [zmin, zmax]]
+    box_boundary : np.ndarray of shape (3,2)
+        Boundary of the simulation box in the form of
+        `np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]])`.
     **kwargs : 
         Forwarded to axis.plot.
 
@@ -901,7 +909,7 @@ def box(axis: mpl.axes.Axes, box_boundary: np.ndarray, **kwargs) -> None:
 
 
 def particles(
-        ax: mpl.axes, coords: np.ndarray, box_boundary: np.ndarray,
+        ax: mpl.axes.Axes, coords: np.ndarray, box_boundary: np.ndarray,
         radius: np.ndarray | float, scalefactor: float = 1.0,
         values: np.ndarray | None = None,
         cmap: list | str = 'viridis', set_ax_limits: bool = True,
@@ -927,9 +935,9 @@ def particles(
         The matplotlib axes object where the particles will be plotted.
     coords : np.ndarray
         An array of coordinates for the particles.
-    box_boundary : np.ndarray
-        The boundary of the simulation box.
-        [[xmin, xmax], [ymin, ymax], [zmin, zmax]]
+    box_boundary : np.ndarray of shape (3,2)
+        Boundary of the simulation box in the form of
+        `np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]])`.
     radius : float or np.ndarray
         The radius of the particles.
     scalefactor : float, optional
@@ -1073,7 +1081,7 @@ def particles(
 
 
 def field(
-        ax: mpl.axes, density: np.ndarray, X: np.ndarray, Y: np.ndarray,
+        ax: mpl.axes.Axes, density: np.ndarray, X: np.ndarray, Y: np.ndarray,
         cmap: list | str = 'plasma', box_boundary: np.ndarray | None = None,
         vmin: float | None = None, vmax: float | None = None,
         cscale: str = 'lin', verbose: bool = False, **kwargs
@@ -1095,10 +1103,10 @@ def field(
     cmap: list or str, optional
         A list representing the colormap to use for the particles or
         the name of a matplotlib colormap. The default is 'viridis'.
-    box_boundary : np.ndarray | None, optional
-        The boundary of the simulation box.
+    box_boundary : np.ndarray of shape (3,2) | None, optional
+        Boundary of the simulation box in the form of
+        `np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]])`.
         Needed if `set_ax_limits = True`.
-        [[xmin, xmax], [ymin, ymax], [zmin, zmax]]
     vmin : float or None, optional
         Lower limit for the colormap. The default is None.
     vmax : float or None, optional
@@ -1511,7 +1519,7 @@ def animate_trajectory(
     ...     traj, "./figures/plot/plot-animate_trajectory_1.gif",
     ...     formatter=lambda x: amep.plot.format_axis(x, direction="out"),
     ...     vmin=0.0, vmax=3.0, ftype='c', xlabel=r"$x$", ylabel=r"$y$",
-    ...     cbar_label=r"$c(x,y)$"
+    ...     cbar_label=r"$c(x,y)$", fps=5
     ... )
     >>>
     
@@ -1525,7 +1533,7 @@ def animate_trajectory(
     ...     traj, "./figures/plot/plot-animate_trajectory_2.gif",
     ...     painter=lambda x,p: (x.velocities(ptype=p)**2).sum(axis=1),
     ...     xlabel='x', ylabel='y', cbar_label=r'$|\vec{v}|^2$',
-    ...     vmin=1e0, vmax=1e6, cscale="log"
+    ...     vmin=1e0, vmax=1e6, cscale="log", fps=5
     ... )
     >>> 
 
diff --git a/amep/reader.py b/amep/reader.py
index 363115a..be74704 100644
--- a/amep/reader.py
+++ b/amep/reader.py
@@ -37,7 +37,6 @@
 import warnings
 import h5py
 import string
-
 import numpy as np
 
 from tqdm.autonotebook import tqdm
@@ -125,6 +124,7 @@ def version(self):
         '''
         return self.__version
 
+
 class LammpsReader(BaseReader):
     '''Reads LAMMPS simulation data and writes it to an hdf5 file.
     '''
@@ -502,6 +502,17 @@ def __init__(
                         else:
                             frame['torque'][:] = torque
 
+                        if 'angmom' not in frame.keys():
+                            frame.create_dataset('angmom',
+                                                 (N, 3),
+                                                 data=angmom,
+                                                 dtype=DTYPE,
+                                                 compression=COMPRESSION,
+                                                 shuffle=SHUFFLE,
+                                                 fletcher32=FLETCHER)
+                        else:
+                            frame['angmom'][:] = angmom
+
                         # spatial dimension
                         frame.attrs['d'] = d
 
@@ -562,6 +573,7 @@ def __init__(
         finally:
             pass
 
+
     def __sorter(self, item):
         r'''
         Returns the time step of a dump file that is given
@@ -580,7 +592,6 @@ def __sorter(self, item):
         return float(basename.split(key)[1].split('.')[0])
 
 
-
     def __get_timestep_from_logfile(self):
         r'''
         Reads the time step from the log.lammps file.
@@ -1083,4 +1094,4 @@ def __detect_delimiter(
         
         if delimiter == []:
             return None
-        return delimiter
+        return delimiter
\ No newline at end of file
diff --git a/amep/spatialcor.py b/amep/spatialcor.py
index 024b042..b75a19f 100644
--- a/amep/spatialcor.py
+++ b/amep/spatialcor.py
@@ -1392,9 +1392,9 @@ def sfiso(
     ----------
     coords : np.ndarray of shape (N,3)
         Coordinates.
-    np.ndarray of shape (3,2)
+    box_boundary : np.ndarray of shape (3,2)
         Boundary of the simulation box in the form of
-        `np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]])`
+        `np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]])`.
     N : int
         Total number of particles.
     qmax : float, optional
@@ -1650,9 +1650,9 @@ def sf2d(
     ----------
     coords : np.ndarray
         Coordinate frame (3D).
-    box_boundary : np.ndarray
-        Boundaries of the simulation box given in the form
-        np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]]).
+    box_boundary : np.ndarray of shape (3,2)
+        Boundary of the simulation box in the form of
+        `np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]])`.
     other_coords : np.ndarray, optional
         Coordinate frame of the other species to which the pair correlation
         is calculated. The default is None (uses coords).
diff --git a/amep/trajectory.py b/amep/trajectory.py
index 3b8379d..964be10 100644
--- a/amep/trajectory.py
+++ b/amep/trajectory.py
@@ -27,7 +27,7 @@
 are containers for a time series of simulation data (particle-based data or 
 continuum fields). The following classes are included:
 
-    ParticleTracetory :
+    ParticleTrajetory :
         Time-series of per-particle data such as coordinates, velocities,
         forces, etc.
     FieldTrajectory :
@@ -59,6 +59,8 @@ def __init__(self, reader) -> None:
         Creates a trajectory object containing data frames for multiple
         time steps.
 
+        Can be viewed as a Sequence of :class:`.BaseFrame`.
+
         Parameters
         ----------
         reader : BaseReader
@@ -180,7 +182,7 @@ def add_particle_info(
         Returns
         -------
         None.
-        
+
         Examples
         --------
         >>> import amep
@@ -190,9 +192,11 @@ def add_particle_info(
         >>> print(traj.get_particle_info())
         {1: {'name': 'active'}, 2: {'name': 'passive'}}
         >>> traj.delete_particle_info(None)
-        >>> 
+        >>>
 
         '''
+        if not isinstance(ptype, int):
+            raise TypeError("Trajectory.add_particle_info(): The particle type <ptype> must be an integer.")
         with h5py.File(os.path.join(self.reader.savedir, self.reader.filename), 'a') as root:
             if str(ptype) not in root['particles'].keys():
                 root['particles'].create_group(str(ptype))
@@ -213,7 +217,7 @@ def get_particle_info(self, ptype: int | None = None) -> dict:
         -------
         p : dict
             Parameters.
-        
+
         Examples
         --------
         >>> import amep
@@ -223,9 +227,11 @@ def get_particle_info(self, ptype: int | None = None) -> dict:
         >>> print(traj.get_particle_info())
         {1: {'name': 'active'}, 2: {'name': 'passive'}}
         >>> traj.delete_particle_info(None)
-        >>> 
+        >>>
 
         '''
+        if not isinstance(ptype, int) and ptype is not None:
+            raise TypeError("Trajectory.get_particle_info(): The particle type <ptype> must be an integer.")
         with h5py.File(os.path.join(self.reader.savedir, self.reader.filename), 'r') as root:
             if ptype is None:
                 p = {}
@@ -320,6 +326,8 @@ def __init__(self, reader):
         Creates a trajectory object containing data frames for multiple
         time steps.
 
+        Can be viewed as a Sequence of :class:`.BaseField`.
+
         Parameters
         ----------
         reader : BaseReader
diff --git a/amep/utils.py b/amep/utils.py
index 3e4e166..2e4b33b 100644
--- a/amep/utils.py
+++ b/amep/utils.py
@@ -31,6 +31,7 @@
 # =============================================================================
 import os
 import time
+from packaging.version import Version
 
 import numpy as np
 
@@ -309,9 +310,9 @@ def profile(
     ----------
     coords : np.ndarray
         Coordinate frame.
-    box_boundary : np.ndarray
-        Boundaries of the simulation box.
-        [[xmin,xmax],[ymin,ymax],[zmin,zmax]]
+    box_boundary : np.ndarray of shape (3,2)
+        Boundary of the simulation box in the form of
+        `np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]])`.
     values : np.ndarray
         Per-particle values (1D array).
     binwidth : float, optional
@@ -1188,10 +1189,14 @@ def sq_from_gr(
     '''
     if twod:
         ydata = ((gr - 1) * r).reshape(-1, 1) * special.jv(0, r.reshape(-1, 1) * q.reshape(1, -1))
-        return np.trapz(x=r, y=ydata, axis=0) * (2 * np.pi * rho) + 1
+        if Version(np.__version__) < Version("2.0.0"):
+            return np.trapz(x=r, y=ydata, axis=0) * (2 * np.pi * rho) + 1
+        return np.trapezoid(x=r, y=ydata, axis=0) * (2 * np.pi * rho) + 1
     else:
         ydata = ((gr - 1) * r).reshape(-1, 1) * np.sin(r.reshape(-1, 1) * q.reshape(1, -1))
-        return np.trapz(x=r, y=ydata, axis=0) * (4 * np.pi * rho / q) + 1    
+        if Version(np.__version__) < Version("2.0.0"):
+            return np.trapz(x=r, y=ydata, axis=0) * (4 * np.pi * rho / q) + 1
+        return np.trapezoid(x=r, y=ydata, axis=0) * (4 * np.pi * rho / q) + 1
     
 
 def sq_from_sf2d(
@@ -1273,7 +1278,10 @@ def msd_from_vacf(vacf, times):
         Time values as 1D array.
 
     '''
-    msd = 2*np.array([np.trapz((times[i]-times[:i+1])*vacf[:i+1], x=times[:i+1]) for i in range(len(times))])
+    if Version(np.__version__) < Version("2.0.0"):
+        msd = 2*np.array([np.trapz((times[i]-times[:i+1])*vacf[:i+1], x=times[:i+1]) for i in range(len(times))])
+    else:
+        msd = 2*np.array([np.trapezoid((times[i]-times[:i+1])*vacf[:i+1], x=times[:i+1]) for i in range(len(times))])
     return msd, times
 
 
@@ -1385,7 +1393,9 @@ def domain_length(
         s_fac = s_fac[q<=qmax]
         q = q[q<=qmax]
     
-    return 2*np.pi*(np.trapz(s_fac, x=q)/np.trapz(s_fac*q, x=q))
+    if Version(np.__version__) < Version("2.0.0"):
+        return 2*np.pi*(np.trapz(s_fac, x=q)/np.trapz(s_fac*q, x=q))
+    return 2*np.pi*(np.trapezoid(s_fac, x=q)/np.trapezoid(s_fac*q, x=q))
 
 
 # =============================================================================
diff --git a/doc/source/_templates/custom-class-template.rst b/doc/source/_templates/custom-class-template.rst
index d63879f..1a1ccd3 100644
--- a/doc/source/_templates/custom-class-template.rst
+++ b/doc/source/_templates/custom-class-template.rst
@@ -5,6 +5,7 @@
 .. autoclass:: {{ objname }}
    :members:                                    
    :show-inheritance:                          
+   :special-members: __getitem__
    :inherited-members:                         
 
    {% block methods %}
diff --git a/doc/source/conf.py b/doc/source/conf.py
index b2145f1..56891cf 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -9,7 +9,7 @@
 project = 'AMEP'
 copyright = '2023-2024, Lukas Hecht, Kay-Robert Dormann, Kai Luca Spanheimer'
 author = 'Lukas Hecht, Kay-Robert Dormann, Kai Luca Spanheimer'
-release = '1.0.2'
+release = '1.0.3'
 
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration
@@ -51,7 +51,7 @@
                 "url": "https://github.com/amepproject/amep",  # required
                 # Icon class (if "type": "fontawesome"), or path to local image (if "type": "local")
                 "icon": "fa-brands fa-square-github",
-                # The type of image to be used (see below for details)
+                # The type of image to be used
                 "type": "fontawesome",
             },
             {
@@ -61,8 +61,16 @@
                 "url": "https://pypi.org/project/amep/",  # required
                 # Icon class (if "type": "fontawesome"), or path to local image (if "type": "local")
                 "icon": "fa-brands fa-python",
-                # The type of image to be used (see below for details)
+                # The type of image to be used
                 "type": "fontawesome",
             }
-            ]
-        }
+        ],
+        #"switcher": {
+        #    "json_url": "https://amepproject.de/switcher.json",
+        #    "version_match": release
+        #},
+        #"check_switcher": True,
+        #"navbar_persistent": ["search-button.html", "theme-switcher.html"],
+        #"navbar_end": ["version-switcher.html", "icon-links.html"]
+        # "navbar_start": ["navbar_logo", "version-switcher"]
+    }
diff --git a/doc/source/user_guide/howto_cite.rst b/doc/source/howto_cite.rst
similarity index 95%
rename from doc/source/user_guide/howto_cite.rst
rename to doc/source/howto_cite.rst
index d423787..0dac878 100644
--- a/doc/source/user_guide/howto_cite.rst
+++ b/doc/source/howto_cite.rst
@@ -1,6 +1,6 @@
-=====================
-How to cite **AMEP**
-=====================
+============
+How to Cite
+============
 
 If you use **AMEP** for a project that leads to a scientific publication, please acknowledge 
 the use of **AMEP** within the body of your publication for example by copying or adapting 
diff --git a/doc/source/index.rst b/doc/source/index.rst
index 26b512e..8d63828 100644
--- a/doc/source/index.rst
+++ b/doc/source/index.rst
@@ -3,8 +3,8 @@
   :align: center
 
 
-AMEP documentation
-==================
+AMEP 1.0.3 documentation
+========================
 
 .. include:: ../../README.md
    :parser: myst_parser.sphinx_
@@ -17,7 +17,7 @@ Table of Contents
    :maxdepth: 2
 
    gettingstarted/index
-   datastructures/index
+   howto_cite
    user_guide/index
    api
 
diff --git a/doc/source/datastructures/field_data.rst b/doc/source/user_guide/datastructures/field_data.rst
similarity index 100%
rename from doc/source/datastructures/field_data.rst
rename to doc/source/user_guide/datastructures/field_data.rst
diff --git a/doc/source/datastructures/function_flow.rst b/doc/source/user_guide/datastructures/function_flow.rst
similarity index 100%
rename from doc/source/datastructures/function_flow.rst
rename to doc/source/user_guide/datastructures/function_flow.rst
diff --git a/doc/source/datastructures/index.rst b/doc/source/user_guide/datastructures/index.rst
similarity index 100%
rename from doc/source/datastructures/index.rst
rename to doc/source/user_guide/datastructures/index.rst
diff --git a/doc/source/user_guide/index.rst b/doc/source/user_guide/index.rst
index e2b9bfd..fa39e0f 100644
--- a/doc/source/user_guide/index.rst
+++ b/doc/source/user_guide/index.rst
@@ -15,7 +15,7 @@ your research.
 
 .. toctree::
 
-   howto_cite
+   datastructures/index
 
 .. toctree::
 
diff --git a/examples/data/continuum.h5amep b/examples/data/continuum.h5amep
index e5c4a60..9ba7b7c 100644
Binary files a/examples/data/continuum.h5amep and b/examples/data/continuum.h5amep differ
diff --git a/examples/data/lammps.h5amep b/examples/data/lammps.h5amep
index 1ea36ce..51b79f2 100644
Binary files a/examples/data/lammps.h5amep and b/examples/data/lammps.h5amep differ
diff --git a/pyproject.toml b/pyproject.toml
index 04eeaa2..b45a8e2 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [project]
 name = "amep"
-version = "1.0.2"
+version = "1.0.3"
 license = {file="LICENSE"}
 authors = [
     {name = "Lukas Hecht", email = "lukas.hecht@pkm.tu-darmstadt.de"},
@@ -33,7 +33,7 @@ dependencies = [
     "h5py >= 3.7.0",
     "matplotlib >= 3.6.2",
     "numba >= 0.56.4",
-    "numpy >= 1.21.6",
+    "numpy >= 1.26.6",
     "scipy >= 1.10.0",
     "scikit-image >= 0.20.0",
     "tqdm >= 4.65.0"


Test patch:
diff --git a/test/test_base.py b/test/test_base.py
index 82fbfc6..e4837c8 100644
--- a/test/test_base.py
+++ b/test/test_base.py
@@ -23,88 +23,40 @@
 # =============================================================================
 # IMPORT MODULES
 # =============================================================================
-import numpy as np
+from pathlib import Path
 import unittest
 import os
+import numpy as np
 import amep
 
 # =============================================================================
 # GLOBAL CONFIG.
 # =============================================================================
-DATADIR = './data'
-FIELDDIR = './data/field'
+DATADIR = Path('../examples/data/')
+FIELDDIR = DATADIR/"continuum"
 RNG = np.random.default_rng(1234)
 
 # for field data creation
-BOX = np.array([[-10,10], [-20,20], [-0.5,0.5]])
-SHAPE = np.array([40,80,2])
-X,Y,Z = np.meshgrid(
-    np.linspace(BOX[0,0], BOX[0,1], SHAPE[0]),
-    np.linspace(BOX[1,0], BOX[1,1], SHAPE[1]),
-    np.linspace(BOX[2,0], BOX[2,1], SHAPE[2]),
-    indexing = 'ij'
+BOX = np.array([[0.5, 49.5], [0.5, 49.5], [-0., 0.]])
+SHAPE = np.array([50, 50])
+X, Y = np.meshgrid(
+        np.linspace(BOX[0, 0], BOX[0, 1], SHAPE[0]),
+        np.linspace(BOX[1, 0], BOX[1, 1], SHAPE[1]),
+        indexing='xy'
 )
-COORDS = [X.flatten(), Y.flatten(), Z.flatten()]
-TIMESTEPS = [0, 1000, 2000, 3000]
+COORDS = [X.flatten(), Y.flatten()]
+TIMESTEPS = np.arange(0, 10001, 1000)
 DT = 1.0
-FIELDS = ['rho', 'c', 'alpha', 'beta']
-
-# create directories if they do not exist
-if not os.path.isdir('./data'):
-    os.mkdir('./data')
-if not os.path.isdir(FIELDDIR):
-    os.mkdir(FIELDDIR)
-
-# =============================================================================
-# TEST DATA GENERATORS
-# =============================================================================
-def create_field_data():
-    
-    # generate grid.txt file
-    with open(
-            os.path.join(FIELDDIR, "grid.txt"),
-            "w",
-            encoding="utf-8"
-    ) as wfile:
-        wfile.write(f'BOX:\n{BOX[0,0]}\t{BOX[0,1]}\n{BOX[1,0]}\t{BOX[1,1]}\n{BOX[2,0]}\t{BOX[2,1]}\n')
-        wfile.write('SHAPE:\n' + '\t'.join(str(s) for s in SHAPE) + '\n')
-        wfile.write('COORDINATES:\tX\tY\tZ\n')
-        wfile.write('\n'.join('\t'.join(
-            str(COORDS[i][j]) for i in range(3)
-        ) for j in range(len(COORDS[0]))))
-    
-    # generate dump files with random data
-    for i, step in enumerate(TIMESTEPS):
-        with open(
-                os.path.join(FIELDDIR, f'field_{step}.txt'),
-                "w",
-                encoding="utf-8"
-        ) as wfile:
-            wfile.write(f'TIMESTEP:\n{step}\nDATA:\t'+'\t'.join(FIELDS)+'\n')
-            wfile.write('\n'.join('\t'.join(
-                str(RNG.random()) for _ in FIELDS
-            ) for _ in range(len(COORDS[0]))))
-            
-    # generate a pseudo simulation script
-    with open(
-        os.path.join(FIELDDIR, "continuum.run"),
-        "w",
-        encoding="utf-8"
-    ) as wfile:
-        wfile.write('line 1\n')
-        wfile.write('line 2\n')
-        wfile.write('line 3\n')
-        wfile.write('line 4\n')
-        wfile.write('line 5')
+FIELDS = ['p', 'c']
 
 
 # =============================================================================
 # BASEFIELD TESTS
 # =============================================================================
 class TestBaseField(unittest.TestCase):
-    
+    """Testing Suite for all rudimetary functions of BaseField."""
     @classmethod
-    def setUpClass(self):
+    def setUpClass(cls):
         """
         Basic setup. Generating test data.
 
@@ -114,23 +66,23 @@ def setUpClass(self):
 
         """
         # generate some test field data
-        create_field_data()
-        
+        # create_field_data()
+
         # load the data
-        self.traj = amep.load.traj(
+        cls.traj = amep.load.traj(
             FIELDDIR,
             mode='field',
-            delimiter='\t',
+            delimiter=' ',
             dumps='field_*.txt',
             reload=True,
             timestep=DT
         )
-        
+
         # get frame
-        self.frame = self.traj[0]
-    
+        cls.frame = cls.traj[0]
+
     def test_data(self):
-        # check shape and data access
+        """Check shape and data access."""
         for key in FIELDS:
             self.assertEqual(
                 self.frame.data(key).shape,
@@ -138,16 +90,16 @@ def test_data(self):
                 f'''Data {key} has wrong shape {self.frame.data(key).shape}.
                     Require shape {SHAPE}.'''
             )
-        data = self.frame.data(FIELDS[0], FIELDS[1], FIELDS[2])
+        data = self.frame.data(FIELDS[0], FIELDS[1])
         self.assertEqual(
-            data.shape, 
-            (3,) + tuple(SHAPE),
-            f'''Loaded data with keys {[FIELDS[0], FIELDS[1], FIELDS[2]]}
+            data.shape,
+            (2,) + tuple(SHAPE),
+            f'''Loaded data with keys {[FIELDS[0], FIELDS[1]]}
                 has invalid shape {data.shape}. Require shape {(3,) + SHAPE}'''
         )
-    
+
     def test_step(self):
-        # step
+        """Check Data access for step method."""
         for i in range(self.traj.nframes):
             self.assertEqual(
                 self.traj[i].step,
@@ -155,45 +107,50 @@ def test_step(self):
                 f'''Incorrect time step: got {self.traj[i].step} instead of
                     {TIMESTEPS[i]} at index {i}.'''
             )
+
     def test_center(self):
-        # center
-        center = (BOX[:,0]+BOX[:,1]) / 2.0
+        """Check Data access for center method."""
+        center = (BOX[:, 0]+BOX[:, 1]) / 2.0
         comparison = self.frame.center == center
         self.assertTrue(
             comparison.all(),
             f'''Incorrect center: got {self.frame.center} instead of
                 {center}.'''
         )
+
     def test_dim(self):
-        # dim
+        """Check Data access for dimension method."""
         self.assertEqual(
             self.frame.dim,
-            3,
+            2,
             f'''Incorrect spatial dimension: got {self.frame.dim}
-                instead of 3.'''
+                instead of 2.'''
         )
+
     def test_box(self):
-        # box
+        """Check Data access for box method."""
         comparison = self.frame.box == BOX
         self.assertTrue(
             comparison.all(),
             f'''Incorrect simulation box: got {self.frame.box} instead
                 of {BOX}.'''
         )
+
     def test_volume(self):
-        # volume
-        volume = np.prod(np.diff(BOX))
+        """Check Data access for Volume method."""
+        volume = np.prod(np.diff(BOX[:2]))
         comparison = self.frame.volume == volume
         self.assertTrue(
             comparison.all(),
             f'''Incorrect simulation box volume: got {self.frame.volume}
                 instead of {volume} with simulation box {BOX}.'''
         )
+
     def test_grid(self):
-        # grid
+        """Check Data access for grid method."""
         self.assertEqual(
             len(self.frame.grid),
-            len(COORDS),
+            2,
             f'''Incorrect grid dimension: got {len(self.frame.grid)}
                 instead of {len(COORDS)}.'''
         )
@@ -205,39 +162,34 @@ def test_grid(self):
         )
         comparisonx = np.isclose(self.frame.grid[0], X)
         comparisony = np.isclose(self.frame.grid[1], Y)
-        comparisonz = np.isclose(self.frame.grid[2], Z)
         self.assertTrue(
             comparisonx.all(),
-            f'''Wrong x coordinates. Found {len(X[comparisonx==0])}
+            f'''Wrong x coordinates. Found {len(X[comparisonx == 0])}
                 invalid entries out of {len(X.flatten())}.'''
         )
         self.assertTrue(
             comparisony.all(),
-            f'''Wrong y coordinates. Found {len(Y[comparisony==0])}
+            f'''Wrong y coordinates. Found {len(Y[comparisony == 0])}
                 invalid entries out of {len(Y.flatten())}.'''
         )
-        self.assertTrue(
-            comparisonz.all(),
-            f'''Wrong z coordinates. Found {len(Z[comparisonz==0])}
-                invalid entries out of {len(Z.flatten())}.'''
-        )
+
     def test_keys(self):
-        # keys
+        """Check Data access for keys method."""
         keys = self.frame.keys
         for key in keys:
             self.assertTrue(
                 key in FIELDS,
                 f'''Incorrect key {key}. Allowed keys are {FIELDS}.'''
             )
-    
+
 
 # =============================================================================
 # BASETRAJECTORY TESTS
 # =============================================================================
 class TestBaseTrajectory(unittest.TestCase):
-    
+    """Testing Suite for all rudimetary functions of BaseTrajectory."""
     @classmethod
-    def setUpClass(self):
+    def setUpClass(cls):
         """
         Basic setup. Generating test data.
 
@@ -247,34 +199,40 @@ def setUpClass(self):
 
         """
         # generate some test field data
-        create_field_data()
-        
+        # create_field_data()
         # create a trajectory object
-        self.traj = amep.load.traj(
+        cls.traj = amep.load.traj(
             FIELDDIR,
             mode='field',
-            delimiter='\t',
+            delimiter=' ',
             dumps='field_*.txt',
             reload=True,
             timestep=DT
         )
-        
         # add author information
-        self.traj.add_author_info(
+        cls.traj.add_author_info(
             'Author A', 'email', 'A@beispiel.de')
-        self.traj.add_author_info(
+        cls.traj.add_author_info(
             'Author A', 'affiliation', 'Institute A')
-        self.traj.add_author_info(
+        cls.traj.add_author_info(
             'Author B', 'email', 'B@beispiel.de')
-        self.traj.add_author_info(
+        cls.traj.add_author_info(
             'Author B', 'affiliation', 'Institute B')
-        
+
         # add software information
-        self.traj.add_software_info('version', '4May2023')
-        self.traj.add_software_info('web', 'lammps.org')
-        self.traj.add_software_info('name', 'LAMMPS')
-        
+        cls.traj.add_software_info('version', '4May2023')
+        cls.traj.add_software_info('web', 'lammps.org')
+        cls.traj.add_software_info('name', 'LAMMPS')
+
     def test_get_item(self):
+        """
+        Test access to items in the trajectory.
+
+        Returns
+        -------
+        None.
+
+        """
         self.assertTrue(
             isinstance(self.traj[0], amep.base.BaseField),
             f'''traj[0] returns an incorrect item type. Got
@@ -289,20 +247,20 @@ def test_get_item(self):
             len(self.traj[2:]),
             len(TIMESTEPS)-2,
             f'''traj[2:] has incorrect length. Got {len(self.traj[2:])}
-            instead of {len(TIMESTEPS)-2}.'''    
+            instead of {len(TIMESTEPS)-2}.'''
         )
         self.assertTrue(
-            isinstance(self.traj[[0,2]], list),
-            f'''traj[[0,2]] returns an incorrect type. Got {type(self.traj[[0,2]])}
-            instead of list.'''
+            isinstance(self.traj[[0, 2]], list),
+            f'''traj[[0,2]] returns an incorrect type.
+                Got {type(self.traj[[0, 2]])} instead of list.'''
         )
         self.assertEqual(
-            len(self.traj[[0,2]]),
+            len(self.traj[[0, 2]]),
             2,
-            f'''traj[[0,2]] has incorrect length. Got {len(self.traj[[0,2]])}
-            instead of 2.'''    
+            f'''traj[[0,2]] has incorrect length. Got {len(self.traj[[0, 2]])}
+            instead of 2.'''
         )
-        for frame,i in zip(self.traj[[0,2]], [0,2]):
+        for frame, i in zip(self.traj[[0, 2]], [0, 2]):
             self.assertEqual(
                 frame.step,
                 TIMESTEPS[i],
@@ -311,9 +269,17 @@ def test_get_item(self):
             )
 
     def test_add_author_info(self):
+        """
+        Test adding, reading and removing author info to the trajectory.
+
+        Returns
+        -------
+        None.
+
+        """
         # add info
         self.traj.add_author_info('Author A', 'test', 'test')
-        
+
         # check if info has been stored
         self.assertIn(
             'test',
@@ -322,11 +288,19 @@ def test_add_author_info(self):
             The new key "test" is not in the stored information given by
             {self.traj.get_author_info('Author A').keys()}.'''
         )
-        
+
         # delete info
         self.traj.delete_author_info('Author A', 'test')
 
     def test_get_author_info(self):
+        """
+        Test reading author from to the trajectory.
+
+        Returns
+        -------
+        None.
+
+        """
         self.assertTrue(
             isinstance(self.traj.get_author_info('Author A'), dict),
             f'''Wrong type returned. Got
@@ -344,9 +318,16 @@ def test_get_author_info(self):
             f'''Missing key "affiliation". Stored keys are
             {self.traj.get_author_info('Author A').keys()}.'''
         )
-    
+
     def test_delete_author_info(self):
-        # delete author information
+        """
+        Test deleting author from to the trajectory.
+
+        Returns
+        -------
+        None.
+
+        """
         self.traj.delete_author_info('Author B', 'email')
         # check if it has been deleted correctly
         self.assertNotIn(
@@ -356,8 +337,16 @@ def test_delete_author_info(self):
         )
         # add the information again
         self.traj.add_author_info('Author B', 'email', 'B@beispiel.de')
-    
+
     def test_authors(self):
+        """
+        Test acessing author list from to the trajectory.
+
+        Returns
+        -------
+        None.
+
+        """
         self.assertTrue(
             isinstance(self.traj.authors, list),
             f'''Invalid type. Got {type(self.traj.authors)} instead of list'''
@@ -378,8 +367,16 @@ def test_authors(self):
             self.traj.authors,
             f'''Author "Author B" not in {self.traj.authors}.'''
         )
-    
+
     def test_add_software_info(self):
+        """
+        Test adding software info from to the trajectory.
+
+        Returns
+        -------
+        None.
+
+        """
         # add info
         self.traj.add_software_info('test', 'v102')
         # check if info has been added correctly
@@ -397,8 +394,16 @@ def test_add_software_info(self):
         )
         # delete software info
         self.traj.delete_software_info('test')
-    
+
     def test_delete_software_info(self):
+        """
+        Test deleting software info from to the trajectory.
+
+        Returns
+        -------
+        None.
+
+        """
         # delete software information
         self.traj.delete_software_info('version')
         # check if it has been deleted correctly
@@ -410,8 +415,16 @@ def test_delete_software_info(self):
         )
         # add deleted information again
         self.traj.add_software_info('version', '4May2023')
-    
+
     def test_software(self):
+        """
+        Test acessing software info from to the trajectory.
+
+        Returns
+        -------
+        None.
+
+        """
         self.assertTrue(
             isinstance(self.traj.software, dict),
             f'''Invalid type. Got {type(self.traj.software)} instead of dict'''
@@ -428,8 +441,16 @@ def test_software(self):
             f'''Got wrong software name {self.traj.software['name']} instead of
             "LAMMPS".'''
         )
-    
+
     def test_info(self):
+        """
+        Test acessing general info from to the trajectory.
+
+        Returns
+        -------
+        None.
+
+        """
         self.assertTrue(
             isinstance(self.traj.info, dict),
             f'''Invalid type. Got {type(self.traj.info)} instead of dict.'''
@@ -449,8 +470,16 @@ def test_info(self):
             self.traj.info,
             f'''The key "params" is missing. Got {self.traj.info}.'''
         )
-        
+
     def test_params(self):
+        """
+        Test acessing params of the trajectory.
+
+        Returns
+        -------
+        None.
+
+        """
         self.assertTrue(
             isinstance(self.traj.params, dict),
             f'''Invalid type. Got {type(self.traj.params)} instead of dict.'''
@@ -462,8 +491,16 @@ def test_params(self):
                 self.traj.params,
                 f'''Missing parameter "{p}". Got {self.traj.params}.'''
             )
-    
+
     def test_add_param(self):
+        """
+        Test adding, acessing and deleting params of the trajectory.
+
+        Returns
+        -------
+        None.
+
+        """
         # add parameter
         self.traj.add_param('A', 1)
         # check if parameter has been added correctly
@@ -480,8 +517,16 @@ def test_add_param(self):
         )
         # delete parameter
         self.traj.delete_param('A')
-    
+
     def test_delete_param(self):
+        """
+        Test deleting params of the trajectory.
+
+        Returns
+        -------
+        None.
+
+        """
         # delete parameter
         self.traj.delete_param('dt')
         # check if parameter has been deleted correctly
@@ -492,83 +537,121 @@ def test_delete_param(self):
         )
         # add pararmeter again
         self.traj.add_param('dt', 1.0)
-    
+
     def test_add_script(self):
-        # add script
-        self.traj.add_script(os.path.join(FIELDDIR, "continuum.run"))
+        """
+        Test adding a script to the trajectory.
+
+        Returns
+        -------
+        None.
+
+        """
+        self.traj.add_script(os.path.join(FIELDDIR, "solver.py"))
         # check if script has been added correctly
         self.assertIn(
-            'continuum.run',
+            'solver.py',
             self.traj.scripts,
-            f'''The script "continuum.run" has not been added to the 
+            f'''The script "solver.py" has not been added to the
             trajectory. Included files are {self.traj.scripts}.'''
         )
         # delete script again
-        self.traj.delete_script('continuum.run')
-    
+        self.traj.delete_script('solver.py')
+
     def test_get_script(self):
+        """
+        Test adding, acessing and deleting a script in the trajectory.
+
+        Returns
+        -------
+        None.
+
+        """
         # add script
-        self.traj.add_script(os.path.join(FIELDDIR, "continuum.run"))
+        self.traj.add_script(os.path.join(FIELDDIR, "solver.py"))
         # clean data directory
-        if os.path.exists(os.path.join(DATADIR, "continuum.run")):
-            os.remove(os.path.join(DATADIR, "continuum.run"))
+        (DATADIR/"solver.py").unlink(missing_ok=True)
         # test get_script method
         self.assertTrue(
-            isinstance(self.traj.get_script('continuum.run'), list),
-            f'''Invalid type. Got {type(self.traj.get_script('continuum.run'))}
+            isinstance(self.traj.get_script('solver.py'), list),
+            f'''Invalid type. Got {type(self.traj.get_script('solver.py'))}
             instead of list.'''
         )
         self.assertEqual(
-            len(self.traj.get_script('continuum.run')), 
-            5,
+            len(self.traj.get_script('solver.py')),
+            199,
             f'''Invalid length. Got
-            {len(self.traj.get_script('continuum.run'))} instead of 5.'''
+            {len(self.traj.get_script('solver.py'))} instead of 5.'''
         )
         # write script to file
         _ = self.traj.get_script(
-            'continuum.run',
-            store = True, 
-            directory = DATADIR
+            'solver.py',
+            store=True,
+            directory=DATADIR
         )
         self.assertTrue(
-            os.path.exists(os.path.join(DATADIR, "continuum.run")),
-            '''Script "continuum.run" has not been stored.'''
+            os.path.exists(os.path.join(DATADIR, "solver.py")),
+            '''Script "solver.py" has not been stored.'''
         )
         # delete script again
-        self.traj.delete_script('continuum.run')
+        self.traj.delete_script('solver.py')
         # clean data directory again
-        if os.path.exists(os.path.join(DATADIR, "continuum.run")):
-            os.remove(os.path.join(DATADIR, "continuum.run"))
+        if os.path.exists(os.path.join(DATADIR, "solver.py")):
+            os.remove(os.path.join(DATADIR, "solver.py"))
 
     def test_delete_script(self):
+        """
+        Test deleting a script in the trajectory.
+
+        Returns
+        -------
+        None.
+
+        """
         # add script
-        self.traj.add_script(os.path.join(FIELDDIR, "continuum.run"))
+        self.traj.add_script(os.path.join(FIELDDIR, "solver.py"))
         # delete script again
-        self.traj.delete_script('continuum.run')
+        self.traj.delete_script('solver.py')
         # check if script has been deleted correctly
         self.assertNotIn(
-            'continuum.run',
+            'solver.py',
             self.traj.scripts,
-            f'''Script "continuum.run" has not been deleted. Available scripts
+            f'''Script "solver.py" has not been deleted. Available scripts
             are {self.traj.scripts}.'''
         )
-    
+
     def test_scripts(self):
+        """
+        Test adding, acessing and deleting a script in the trajectory.
+
+        Returns
+        -------
+        None.
+
+        """
         # add script
-        self.traj.add_script(os.path.join(FIELDDIR, "continuum.run"))
+        self.traj.add_script(os.path.join(FIELDDIR, "solver.py"))
         self.assertEqual(
             len(self.traj.scripts),
             1,
             f'''Invalid length. Got {len(self.traj.scripts)} instead of 1.'''
         )
         self.assertIn(
-            'continuum.run',
+            'solver.py',
             self.traj.scripts,
-            f'''Missing file "continuum.run". Included files are
+            f'''Missing file "solver.py". Included files are
             {self.traj.scripts}.'''
         )
 
     def test_properties(self):
+        """
+        Test acessing properties in the trajectory.
+
+        Returns
+        -------
+        None.
+
+        """
         self.assertEqual(
             self.traj.version, amep.__version__,
             f'''Invalid version. Got {self.traj.version} instead of
@@ -581,7 +664,7 @@ def test_properties(self):
         )
         self.assertTrue(
             isinstance(self.traj.reader, amep.reader.ContinuumReader),
-            f'''Invalid type. Got {type(self.traj.reader)} instead of 
+            f'''Invalid type. Got {type(self.traj.reader)} instead of
             amep.reader.ContinuumReader.'''
         )
         self.assertEqual(
@@ -602,14 +685,12 @@ def test_properties(self):
         self.assertEqual(
             self.traj.nframes,
             len(TIMESTEPS),
-            f'''Invalid number of frames. Got {self.traj.nframes} instead of 
+            f'''Invalid number of frames. Got {self.traj.nframes} instead of
             {len(TIMESTEPS)}.'''
         )
         # check steps
-        invalid_steps = False
-        for step in self.traj.steps:
-            if step not in TIMESTEPS:
-                invalid_steps = True
+        invalid_steps = any(step not in np.array(TIMESTEPS)*DT
+                            for step in self.traj.steps)
         self.assertFalse(
             invalid_steps,
             f'''Invalid step detected. Got {self.traj.steps} instead of
@@ -618,17 +699,15 @@ def test_properties(self):
         self.assertEqual(
             len(self.traj.steps),
             len(TIMESTEPS),
-            f'''Invalid length of steps. Got {len(self.traj.steps)} instead of 
+            f'''Invalid length of steps. Got {len(self.traj.steps)} instead of
             {len(TIMESTEPS)}.'''
         )
         # check times
-        invalid_times = False
-        for time in self.traj.times:
-            if time not in np.array(TIMESTEPS)*DT:
-                invalid_times = True
+        invalid_times = any(time not in np.array(TIMESTEPS)*DT
+                            for time in self.traj.times)
         self.assertFalse(
             invalid_times,
-            f'''Invalid time detected. Got {self.traj.times} instead of 
+            f'''Invalid time detected. Got {self.traj.times} instead of
             {np.array(TIMESTEPS)*DT}.'''
         )
         self.assertEqual(
@@ -644,7 +723,7 @@ def test_properties(self):
         )
         self.assertEqual(
             self.traj.dim,
-            3,
+            2,
             f'''Invalid dimension. Got {self.traj.dim} instead of 3.'''
         )
         self.assertEqual(
@@ -653,14 +732,16 @@ def test_properties(self):
             f'''Invalid savedir. Got {self.traj.savedir} instead of
             {FIELDDIR}.'''
         )
-        
+
+
 # =============================================================================
 # BASEFUNCTION TESTS
 # =============================================================================
 class TestBaseFunction(unittest.TestCase):
-    
+    """Test generation and base utilities of the Function object."""
+
     @classmethod
-    def setUpClass(self):
+    def setUpClass(cls):
         """
         Basic setup. Generating test data.
 
@@ -670,94 +751,142 @@ def setUpClass(self):
 
         """
         # create test data
-        self.x = np.linspace(0,10,100)
-        self.y = 10*self.x + 5
-        
+        cls.x = np.linspace(0, 10, 100)
+        cls.y = 10*cls.x + 5
+
     def test_init(self):
+        """
+        Test initialisation of a BaseFunction.
+
+        Returns
+        -------
+        None.
+
+        """
         # initialize function object
-        Function = amep.base.BaseFunction(2)
+        function = amep.base.BaseFunction(2)
         # check number of parameters
         self.assertEqual(
-            Function.nparams,
+            function.nparams,
             2,
-            f"""Invalid number of parameters. Got {Function.nparams}
+            f"""Invalid number of parameters. Got {function.nparams}
             instead of 2."""
         )
         # check length of key and value lists
         self.assertEqual(
-            len(Function.keys),
+            len(function.keys),
             2,
-            f"""Invalid number of keys. Got {len(Function.keys)}
+            f"""Invalid number of keys. Got {len(function.keys)}
             instead of 2."""
         )
         self.assertEqual(
-            len(Function.errors),
+            len(function.errors),
             2,
             f"""Invalid length of fit error list. Got
-            {len(Function.errors)} instead of 2."""
+            {len(function.errors)} instead of 2."""
         )
         self.assertEqual(
-            len(Function.params),
+            len(function.params),
             2,
             f"""Invalid length of fit results list. Got
-            {len(Function.params)} instead of 2."""
+            {len(function.params)} instead of 2."""
         )
-    
+
     def test_fit(self):
+        """
+        Test fit routines of a BaseFunction.
+
+        Returns
+        -------
+        None.
+
+        """
         # initialize function object
-        Function = amep.base.BaseFunction(2)
+        function = amep.base.BaseFunction(2)
         # fit the test data
-        Function.fit(self.x, self.y)
+        function.fit(self.x, self.y)
         # check fit parameters
         self.assertTrue(
-            np.isclose(Function.params[0], 10.0),
-            f"""Invalid parameter {Function.keys[0]}. Got {Function.params[0]}
+            np.isclose(function.params[0], 10.0),
+            f"""Invalid parameter {function.keys[0]}. Got {function.params[0]}
             instead of 10.0."""
         )
         self.assertTrue(
-            np.isclose(Function.params[1], 5.0),
-            f"""Invalid parameter {Function.keys[1]}. Got {Function.params[1]}
+            np.isclose(function.params[1], 5.0),
+            f"""Invalid parameter {function.keys[1]}. Got {function.params[1]}
             instead of 5.0."""
         )
-    
+
     def test_generate(self):
+        """
+        Test generation routines of a BaseFunction.
+
+        Returns
+        -------
+        None.
+
+        """
         # initialize function object
-        Function = amep.base.BaseFunction(2)
+        function = amep.base.BaseFunction(2)
         # fit the test data
-        Function.fit(self.x, self.y)
+        function.fit(self.x, self.y)
         # generate y data
-        y = Function.generate(self.x)
+        y = function.generate(self.x)
         self.assertTrue(
             np.all(np.isclose(y, self.y)),
             """The generate method returned wrong data."""
         )
-    
+
     def test_properties(self):
+        """
+        Test property access in a BaseFunction.
+
+        Returns
+        -------
+        None.
+
+        """
         # initialize function object
-        Function = amep.base.BaseFunction(2)
+        function = amep.base.BaseFunction(2)
         # fit the test data
-        Function.fit(self.x, self.y)
+        function.fit(self.x, self.y)
         # check data type
         self.assertTrue(
-            isinstance(Function.keys, list),
-            f"""Invalid type for property keys. Got {type(Function.keys)}
+            isinstance(function.keys, list),
+            f"""Invalid type for property keys. Got {type(function.keys)}
             instead of list."""
         )
         self.assertTrue(
-            isinstance(Function.errors, np.ndarray),
-            f"""Invalid type for property errors. Got {type(Function.errors)}
+            isinstance(function.errors, np.ndarray),
+            f"""Invalid type for property errors. Got {type(function.errors)}
             instead of np.ndarray."""
         )
         self.assertTrue(
-            isinstance(Function.params, np.ndarray),
-            f"""Invalid type for property params. Got {type(Function.params)}
+            isinstance(function.params, np.ndarray),
+            f"""Invalid type for property params. Got {type(function.params)}
             instead of np.ndarray."""
         )
         self.assertTrue(
-            isinstance(Function.results, dict),
-            f"""Invalid type for property results. Got {type(Function.results)}
+            isinstance(function.results, dict),
+            f"""Invalid type for property results. Got {type(function.results)}
             instead of dict."""
         )
-    
-if __name__ == '__main__':
-    unittest.main()
+
+
+# =============================================================================
+# UTILS TESTS
+# =============================================================================
+class TestHelperFunctions(unittest.TestCase):
+    """Test generation and base utilities of the Function object."""
+    def test_check_path(self):
+        """"""
+        self.assertRaises(ValueError, amep.base.check_path,
+                          "file.txt", ".csv")
+        self.assertRaises(FileNotFoundError, amep.base.check_path,
+                          "not_there/file.txt", ".txt")
+        self.assertEqual(amep.base.check_path("../test", ""),
+                         ("..", "test"))
+        self.assertEqual(amep.base.check_path("../test", ".csv"),
+                         ("../test", ""))
+        self.assertEqual(amep.base.check_path("test.csv", ".csv"),
+                         (os.getcwd(), "test.csv"))
diff --git a/test/test_continuum.py b/test/test_continuum.py
index 680d0e4..34942c3 100644
--- a/test/test_continuum.py
+++ b/test/test_continuum.py
@@ -21,132 +21,46 @@
 
 Including it's main class, readers and methods."""
 import unittest
-from zipfile import ZipFile
 from pathlib import Path
-from random import random
-from requests import get
 from matplotlib import use
 from amep.load import traj
 from amep import plot
 from amep.continuum import identify_clusters
 use("Agg")
-DATA_DIR = Path("./data")
+DATA_DIR = Path("../examples/data/")
 N_LINES = 100
 COORDS = ("X", "Y", "Z")
 FIELDS = ("c", "rho", "alpha", "omega")
 TIMES = (0, 0.4, 0.5, 0.9, 2.3)
 
-SERVER_URL: str = "https://kuno.fkp.physik.tu-darmstadt.de/d/a3d9887b8a5747e0a56e/files/?p=/"
-URL_END: str = "&dl=1"
-FIELD_DIR: Path = DATA_DIR/"fields"
-PLOT_DIR: Path = DATA_DIR/"plots"
-FILE_NAMES: tuple[str, ...] = ("fields",)
-FILE_END: str = ".zip"
-
-
-def get_test_files():
-    """Gets and extracts needed Simulation data.
-
-    In order to make it accessible for the tests.
-    """
-    for file_name in FILE_NAMES:
-        sim_dir = DATA_DIR
-        zip_path = sim_dir/Path(f"{file_name}{FILE_END}")
-        file_url = f"{SERVER_URL}{file_name}{FILE_END}{URL_END}"
-        if not sim_dir.is_dir():
-            sim_dir.mkdir()
-        try:
-            print(f"Downloading {file_name}")
-            with open(zip_path, "wb") as ofile:
-                ofile.write(get(file_url, allow_redirects=True).content)
-            print(f"unziping {file_name}")
-            with ZipFile(zip_path, "r") as z_file:
-                z_file.extractall(sim_dir)
-            zip_path.unlink()
-        except:
-            print(f"Downloading and Unziping {file_name} failed!")
-
-
-def create_test_files():
-    """Creates grid and some field data files with random data"""
-    with open(DATA_DIR/Path("grid.txt"), "w", encoding="utf-8") as wfile:
-        wfile.write("BOX:\n10\t10\n10\t10\n0.5\t0.5\nCOORDINATES:\t" +
-                    "\t".join(COORDS)+"\n")
-        wfile.write("\n".join("\t".join(str(random())for _ in COORDS)
-                              for _ in range(N_LINES)))
-    for index, time in enumerate(TIMES):
-        with open(DATA_DIR/Path(f"field_{index}.txt"), "w",
-                  encoding="utf-8") as wfile:
-            wfile.write(f"TIMESTEP:\n{time}\nDATA:\t"+"\t".join(FIELDS)+"\n")
-            wfile.write("\n".join("\t".join(str(random()) for _ in FIELDS)
-                                  for _ in range(N_LINES)))
+FIELD_DIR: Path = DATA_DIR/"continuum"
 
 
 class TestFieldMethods(unittest.TestCase):
     """A test case for all field and continuum data methods"""
-    @classmethod
-    def setUpClass(cls):
-        DATA_DIR.mkdir(exist_ok=True)
-        PLOT_DIR.mkdir(exist_ok=True)
-        if not FIELD_DIR.is_dir():
-            get_test_files()
-        else:
-            print("Already here")
-        # if not all(path.isfile(filius) for filius in
-        #            (path.join(DATA_DIR, "grid.txt"),
-        #             *(f"field_{index}" for index, _ in enumerate(TIMES)))):
-        #     DATA_DIR.mkdir(parents=True, exist_ok=True)
-        #     create_test_files()
-        # else:
-        #     print("All needed  Field-Data is here")
-
-    def test_creation(self):
-        """Tests the creation of AMEP file from field data"""
-        # reader=ContinuumReader("../data/continuum_data")
-        # field=BaseField(reader,0)
-        trajectory = traj(str(FIELD_DIR/"cahn-hilliard-model"),
-                          mode="field", reload=True,
-                          start=0.0, stop=1.0, nth=1,
-                          delimiter=" ",
-                          dumps="dump_*.txt")
-        print(trajectory.nframes)
-        field = trajectory[1]
-        print("Box:", field.box)
-        print("Center:", field.center)
-        print("Dimension:", field.dim)
-        print("Grid:", field.grid)
-        print("Step: ", field.step)
-        print("Time: ", field.time)
-        print("Volume", field.volume)
-        print(field.data())
 
-    def test_plottability(self):
-        """Try plotting a field"""
-        trajectory = traj(str(FIELD_DIR/"cahn-hilliard-model"),
-                          mode="field")
-        field = trajectory[-1]
-        fig, axe = plot.new()
-        print(field.grid)
-        axe.pcolormesh(field.grid[0],
-                       field.grid[1],
-                       field.data("c"))
-        fig.savefig(PLOT_DIR/"Field_test.pdf")
+    def test_coords_to_density(self):
+        """Test coords to density.
+            TO BE IMPLEMENTED
+        """
+        pass
 
-    def test_cluster(self):
-        """Test the cluster detection with both methods."""
-        trajectory = traj(str(FIELD_DIR/"cahn-hilliard-model"),
-                          mode="field", reload=False)
-        fig, axe = plot.new(ncols=3)
-        plot.field(axe[0], trajectory[-1].data("c"), *trajectory[-1].grid)
-        _, labels = identify_clusters(trajectory[-1].data("c"))
-        plot.field(axe[1], labels, *trajectory[-1].grid)
-        _, labels = identify_clusters(
-            trajectory[-1].data("c"),
-            cutoff=0.7, method="threshold"
-        )
-        plot.field(axe[2], labels, *trajectory[-1].grid)
-        fig.savefig(PLOT_DIR/"Cluster_test.pdf")
+    def test_density_estimators(self):
+        """Test density estimators.
+            TO BE IMPLEMENTED
+        """
+        pass
 
+    def test_sf2d(self):
+        """Test 2d structure factor.
+            TO BE IMPLEMENTED
+        """
+        pass
 
-if __name__ == '__main__':
-    unittest.main()
\ No newline at end of file
+    def test_cluster_methods(self):
+        """Test the cluster detection with both methods.
+            TO BE IMPLEMENTED when I know how to check cluster
+            detection in a sane way.
+            Maybe combine it with coords to density.
+        """
+        pass
diff --git a/test/test_evaluate.py b/test/test_evaluate.py
index 9e11af5..edd3fb9 100644
--- a/test/test_evaluate.py
+++ b/test/test_evaluate.py
@@ -21,73 +21,105 @@
 
 Including it's main class, readers and methods."""
 import unittest
-from zipfile import ZipFile
 from pathlib import Path
-from random import random
-from requests import get
-import numpy as np
 from matplotlib import use
 from amep.load import traj
-from amep.evaluate import ClusterGrowth, ClusterSizeDist
+from amep.evaluate import ClusterGrowth, ClusterSizeDist, Function, SpatialVelCor, RDF, PCF2d, PCFangle, SF2d
 use("Agg")
-DATA_DIR = Path("./data")
+DATA_DIR = Path("../examples/data/")
 
-SERVER_URL: str = "https://kuno.fkp.physik.tu-darmstadt.de/d/a3d9887b8a5747e0a56e/files/?p=/"
-URL_END: str = "&dl=1"
-FIELD_DIR: Path = DATA_DIR/"fields"
-PARTICLE_DIR: Path = DATA_DIR/"particles"
+FIELD_DIR: Path = DATA_DIR/"continuum"
+PARTICLE_DIR: Path = DATA_DIR/"lammps"
+RESULT_DIR: Path = DATA_DIR/"results"
 PLOT_DIR: Path = DATA_DIR/"plots"
-FILE_NAMES: tuple[str, ...] = ("fields", "particles")
-FILE_END: str = ".zip"
-
-
-def get_test_files():
-    """Gets and extracts needed Simulation data.
-
-    In order to make it accessible for the tests.
-    """
-    for file_name in FILE_NAMES:
-        sim_dir = DATA_DIR
-        zip_path = sim_dir/Path(f"{file_name}{FILE_END}")
-        file_url = f"{SERVER_URL}{file_name}{FILE_END}{URL_END}"
-        sim_dir.mkdir(exist_ok=True)
-        try:
-            print(f"Downloading {file_name}")
-            with open(zip_path, "wb") as ofile:
-                ofile.write(get(file_url, allow_redirects=True).content)
-            print(f"unziping {file_name}")
-            with ZipFile(zip_path, "r") as z_file:
-                z_file.extractall(sim_dir)
-            zip_path.unlink()
-        except:
-            print(f"Downloading and Unziping {file_name} failed!")
 
 
 class TestEvaluateMethods(unittest.TestCase):
     """A test case for all field and continuum data methods"""
     @classmethod
     def setUpClass(cls):
-        DATA_DIR.mkdir(exist_ok=True)
-        PLOT_DIR.mkdir(exist_ok=True)
-        if not FIELD_DIR.is_dir() or not PARTICLE_DIR.is_dir():
-            get_test_files()
-        else:
-            print("Already here")
-        cls.field_trajs = [traj(field) for field in FIELD_DIR.iterdir()]
-        cls.particle_trajs = [traj(field) for field in PARTICLE_DIR.iterdir()]
-
-    # def test_clustersize_dist(self):
-    #     print(ClusterSizeDist(self.field_trajs[0], ftype="c").frames)
+        """Set up needed data"""
+        RESULT_DIR.mkdir(exist_ok=True)
+        cls.field_traj = traj(DATA_DIR/"continuum.h5amep")
+        cls.particle_traj = traj(DATA_DIR/"lammps.h5amep")
 
     def test_cluster_growth(self):
-        self.assertTrue((ClusterGrowth(self.field_trajs[1], scale=1.5, cutoff=0.8,
+        """Test the cluster growth methods.
+
+        Due to their weighting they have a given order they take
+        for arbitrary trajectories.
+        This order gets tested here.
+        Since ClusterGrowth calls ClusterSizeDist we don't
+        have to check it separately.
+        """
+        self.assertTrue((ClusterGrowth(self.field_traj, scale=1.5,
+                                       cutoff=0.8,
                                        ftype="c", mode="mean").frames.sum() <=
-                        ClusterGrowth(self.field_trajs[1], scale=1.5, cutoff=0.8,
-                                      ftype="c",
+                        ClusterGrowth(self.field_traj, scale=1.5,
+                                      cutoff=0.8, ftype="c",
                                       mode="weighted mean").frames.sum()))
-        self.assertTrue((ClusterGrowth(self.field_trajs[1],
+        self.assertTrue((ClusterGrowth(self.field_traj,
                                        ftype="c",
                                        mode="largest").frames >= 0).all())
 
-if __name__ == '__main__':
-    unittest.main()
\ No newline at end of file
+    def test_energy_methods(self):
+        """Test the energy methods.
+        TO BE IMPLEMENTED
+        """
+        pass
+
+    def test_function(self):
+        """Test arbitray function evaluation.
+        TO BE IMPLEMENTED
+        """
+        traj = self.particle_traj
+        ftraj = self.field_traj
+        def msd(frame, start=None):
+            vec = start.unwrapped_coords() - frame.unwrapped_coords()
+            return (vec ** 2).sum(axis=1).mean()
+        msd_eval = Function(
+            self.particle_traj, msd,
+            nav=self.particle_traj.nframes,
+            start=self.particle_traj[0]
+            )
+        msd_eval.name = "msd"
+        msd_eval.save(RESULT_DIR/"msd_eval.h5")
+        pcf2d = PCF2d(traj,
+                      nav=2, nxbins=50, nybins=50,
+                      njobs=4, skip=0.9
+                      )
+        pcf2d.save(RESULT_DIR/"pcf2d.h5")
+        pcfangle = PCFangle(
+            traj, nav=2, ndbins=50, nabins=50,
+            njobs=4, rmax=8.0, skip=0.9
+            )
+        pcfangle.save(RESULT_DIR/"pcfangle.h5")
+        psf2d = SF2d(traj, skip=0.9, nav=2)
+        psf2d.save(RESULT_DIR/"sf2d_eval.h5", database=True, name="particles")
+        fsf2d = SF2d(ftraj, skip=0.9, nav=2, ftype="c")
+        fsf2d.save(RESULT_DIR/"sf2d_eval.h5", database=True, name="field")
+
+    def test_order_evaluations(self):
+        """Test order parameter evaluation.
+        TO BE IMPLEMENTED
+        """
+        pass
+
+    def test_correlation(self):
+        """Test order parameter evaluation.
+        TO BE IMPLEMENTED
+        """
+        svc = SpatialVelCor(self.particle_traj, skip=0.9, nav=2, njobs=4)
+
+        svc.save(RESULT_DIR/"svc.h5")
+        rdfcalc = RDF(
+                self.particle_traj,
+                nav=2, nbins=1000,
+                skip=0.9, njobs=4)
+        rdfcalc.save(RESULT_DIR/'rdf.h5')
+
+    def test_transforms(self):
+        """Test order parameter evaluation.
+        TO BE IMPLEMENTED
+        """
+        pass
diff --git a/test/test_functions.py b/test/test_functions.py
index 6b1e9a9..92448ec 100644
--- a/test/test_functions.py
+++ b/test/test_functions.py
@@ -23,11 +23,11 @@
 import unittest
 from pathlib import Path
 from numpy import linspace, average
+from numpy.random import normal
 from amep.functions import Gaussian, gaussian
 from amep.plot import new
 
-DATA_DIR = Path("./data/")
-PLOT_DIR = DATA_DIR/"plots"
+DATA_DIR = Path("../examples/data/")
 
 
 class TestMathFunctions(unittest.TestCase):
@@ -45,11 +45,4 @@ def test_gaussian(self):
         amp = max(y_vec)
         gauss = Gaussian()
         gauss.fit(x_vec, y_vec, p0=[mean, sig, amp])
-        fig, axs = new(figsize=(5, 4))
-        axs.plot(x_vec, y_vec, label="data")
-        new_x = x_vec*0.3 + 2
-        axs.plot(new_x, gauss.generate(new_x), label="fit")
-        fig.savefig(PLOT_DIR/Path("test_gaussian.pdf"))
-        
-if __name__ == '__main__':
-    unittest.main()
+        # Check if fit and generated data are close.
diff --git a/test/test_load.py b/test/test_load.py
index 5d27c16..209f447 100644
--- a/test/test_load.py
+++ b/test/test_load.py
@@ -24,147 +24,92 @@
 # IMPORT MODULES
 # =============================================================================
 import unittest
-from zipfile import ZipFile
-from requests import get
 from pathlib import Path
 import amep
 
-SERVER_URL: str = "https://kuno.fkp.physik.tu-darmstadt.de/d/a3d9887b8a5747e0a56e/files/?p=/"
-URL_END: str = "&dl=1"
-DATA_DIR: Path = Path("./data")
-PARTICLE_DIR: Path = DATA_DIR/"particles"
-FIELD_DIR: Path = DATA_DIR/"fields"
-TRAJ_DIR: Path = DATA_DIR/"trajs"
+DATA_DIR: Path = Path("../examples/data/")
+PARTICLE_DIR: Path = DATA_DIR/"lammps"
+FIELD_DIR: Path = DATA_DIR/"continuum"
 FILE_NAMES: tuple[str, ...] = ("particles",
                                "fields"
                                )
-FILE_END: str = ".zip"
-
-
-def get_test_files():
-    """Gets and extracts needed Simulation data in order to make it accessible for the tests"""
-    for file_name in FILE_NAMES:
-        sim_dir = DATA_DIR
-        zip_path = sim_dir/Path(f"{file_name}{FILE_END}")
-        file_url = f"{SERVER_URL}{file_name}{FILE_END}{URL_END}"
-        if not sim_dir.is_dir():
-            sim_dir.mkdir()
-        if (sim_dir/file_name).is_dir():
-            print("Directory already there.")
-        else:
-            try:
-                print(f"Downloading {file_name}")
-                with open(zip_path, "wb") as ofile:
-                    ofile.write(get(file_url, allow_redirects=True).content)
-                print(f"unziping {file_name}")
-                with ZipFile(zip_path, "r") as z_file:
-                    z_file.extractall(sim_dir)
-                zip_path.unlink()
-            except:
-                print(f"Downloading and Unziping {file_name} failed!")
 
 
 # =============================================================================
 # MAIN READER TEST
 # =============================================================================
 class TestLoad(unittest.TestCase):
-
-    @classmethod
-    def setUpClass(cls):
-        """
-        Basic setup. Generating test data.
-
-        Returns
-        -------
-        None.
-
-        """
-        DATA_DIR.mkdir(exist_ok=True)
-        TRAJ_DIR.mkdir(exist_ok=True)
-        if not PARTICLE_DIR.is_dir() or not FIELD_DIR.is_dir():
-            get_test_files()
+    '''A test case for all load functions for AMEP save files.'''
 
     def test_traj(self):
         # test lammps
         traj = amep.load.traj(
-            str(PARTICLE_DIR/"different-sized-particles"),
+            PARTICLE_DIR,
             mode='lammps',
             dumps='dump*.txt'
         )
-        traj = amep.load.traj(
-            str(PARTICLE_DIR/"different-sized-particles"),
-            savedir=TRAJ_DIR,
-            mode='lammps',
-            dumps='dump*.txt',
-            trajfile='lammps.h5amep'
-        )
-
         # test fields
         traj = amep.load.traj(
-            str(FIELD_DIR/'cahn-hilliard-model'),
+            FIELD_DIR,
             mode='field',
-            dumps='dump_*.txt',
-            delimiter='\t',
+            dumps='field_*.txt',
+            delimiter=' ',
         )
-        traj = amep.load.traj(
-            str(FIELD_DIR/'cahn-hilliard-model'),
-            savedir='./data/trajs',
+        allframes = amep.load.traj(
+            FIELD_DIR,
             mode='field',
-            dumps='dump_*.txt',
-            trajfile='field.h5amep',
+            dumps='field_*.txt',
+            trajfile='traj.h5amep',
             delimiter=' ',
-            reload = True
-        )
-        allframes = traj.nframes
-        traj = amep.load.traj(
-            str(FIELD_DIR/'cahn-hilliard-model'),
-            savedir='./data/trajs',
+            reload=True,
+            timestep=1
+        ).nframes
+        halfframes = amep.load.traj(
+            FIELD_DIR,
             mode='field',
-            dumps='dump_*.txt',
-            trajfile='field.h5amep',
+            dumps='field_*.txt',
+            trajfile='traj.h5amep',
             delimiter=' ',
-            reload = True,
-            start = 0.5
-        )
-        halfframes = traj.nframes
-        traj = amep.load.traj(
-            str(FIELD_DIR/'cahn-hilliard-model'),
-            savedir='./data/trajs',
+            reload=True,
+            start=0.5,
+            timestep=1
+        ).nframes
+        quarterframes = amep.load.traj(
+            FIELD_DIR,
             mode='field',
-            dumps='dump_*.txt',
-            trajfile='field.h5amep',
+            dumps='field_*.txt',
+            trajfile='traj.h5amep',
             delimiter=' ',
-            reload = True,
-            nth = 4
-        )
-        quarterframes = traj.nframes
-        traj = amep.load.traj(
-            str(FIELD_DIR/'cahn-hilliard-model'),
-            savedir='./data/trajs',
+            reload=True,
+            nth=4,
+            timestep=1
+        ).nframes
+        onefifthframes = amep.load.traj(
+            FIELD_DIR,
             mode='field',
-            dumps='dump_*.txt',
-            trajfile='field.h5amep',
+            dumps='field_*.txt',
+            trajfile='traj.h5amep',
             delimiter=' ',
-            reload = True,
-            stop = 0.2
-        )
-        onefifthframes = traj.nframes
+            reload=True,
+            stop=0.2,
+            timestep=1
+        ).nframes
         self.assertEqual(
             halfframes,
             int(allframes/2)+1,
-            "Invalid number of frames with start=0.5. "\
+            "Invalid number of frames with start=0.5. "
             f"Got {halfframes} instead of {int(allframes/2)+1}."
         )
         self.assertEqual(
             quarterframes,
             int(allframes/4)+1,
-            "Invalid number of frames with nth=4. "\
+            "Invalid number of frames with nth=4. "
             f"Got {halfframes} instead of {int(allframes/4)+1}."
         )
         self.assertEqual(
             onefifthframes,
             int(allframes/5),
-            "Invalid number of frames with stop=0.2. "\
+            "Invalid number of frames with stop=0.2. "
             f"Got {halfframes} instead of {int(allframes/5)}."
         )
 
@@ -173,6 +118,3 @@ def test_evaluation(self):
 
     def test_database(self):
         pass
-
-if __name__ == '__main__':
-    unittest.main()
diff --git a/test/test_order.py b/test/test_order.py
index ad2cb7f..10a8e6e 100644
--- a/test/test_order.py
+++ b/test/test_order.py
@@ -23,17 +23,19 @@
 # =============================================================================
 # IMPORT MODULES
 # =============================================================================
-import numpy as np
 import unittest
+import numpy as np
 import amep
 
+
 # =============================================================================
 # MAIN PBC TEST
 # =============================================================================
 class TestOrder(unittest.TestCase):
+    '''A test case for the calculation of order parameters.'''
 
     @classmethod
-    def setUpClass(self):
+    def setUpClass(cls):
         """
         Basic setup. Generating test data.
 
@@ -42,9 +44,10 @@ def setUpClass(self):
         None.
 
         """
-        self.coords = np.array([[0,0,0],[1,0,0],[0,1,0],[-1,0,0],[0,-2,0]])
-        self.box = np.array([[-4,4],[-4,4],[-0.5,0.5]])
-        self.center = np.array([0,0,0])
+        cls.coords = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0],
+                               [-1, 0, 0], [0, -2, 0]])
+        cls.box = np.array([[-4, 4], [-4, 4], [-0.5, 0.5]])
+        cls.center = np.array([0, 0, 0])
 
     def test_local_number_density(self):
         """
@@ -57,7 +60,7 @@ def test_local_number_density(self):
         """
         # calculate local number density with fixed particle radius
         ld = amep.order.local_number_density(
-            self.coords, self.box, 0.5, rmax = 2.0, pbc = False, enforce_nd = 2
+            self.coords, self.box, 0.5, rmax=2.0, pbc=False, enforce_nd=2
         )
         # expected results for the given coords
         expected = np.array(
@@ -68,26 +71,26 @@ def test_local_number_density(self):
         # check result
         self.assertTrue(
             compare.all(),
-            'Incorrect local number density for calculation with fixed '\
+            'Incorrect local number density for calculation with fixed '
             f'particle radius. Got {ld} instead of {expected}.'
         )
-        
+
         # calculate local number density with different particle radii
-        radius = np.array([0.5,0.5,0.5,0.5,1.5])
+        radius = np.array([0.5, 0.5, 0.5, 0.5, 1.5])
         ld = amep.order.local_number_density(
             self.coords, self.box, radius,
-            rmax = 2.0, pbc = False, enforce_nd = 2
+            rmax=2.0, pbc=False, enforce_nd=2
         )
         # expected results for the given coords
         expected = np.array(
-            [0.35809862, 0.31204799, 0.3315728 , 0.31204799, 0.16137229]
+            [0.35809862, 0.31204799, 0.3315728, 0.31204799, 0.16137229]
         )
         # comparison
         compare = ld.round(3) == expected.round(3)
         # check result
         self.assertTrue(
             compare.all(),
-            'Incorrect local number density for calculation with different '\
+            'Incorrect local number density for calculation with different '
             f'particle radii. Got {ld} instead of {expected}.'
         )
 
@@ -103,11 +106,11 @@ def test_local_mass_density(self):
         # calculate local mass density with fixed mass
         ld1 = amep.order.local_mass_density(
             self.coords, self.box, 0.5, 1.0,
-            rmax = 2.0, pbc = False, enforce_nd = 2
+            rmax=2.0, pbc=False, enforce_nd=2
         )
         ld2 = amep.order.local_mass_density(
             self.coords, self.box, 0.5, 2.0,
-            rmax = 2.0, pbc = False, enforce_nd = 2
+            rmax=2.0, pbc=False, enforce_nd=2
         )
         # expected results for the given coords
         expected = np.array(
@@ -118,7 +121,7 @@ def test_local_mass_density(self):
         # check result
         self.assertTrue(
             compare.all(),
-            'Incorrect local mass density for calculation with fixed '\
+            'Incorrect local mass density for calculation with fixed '
             f'mass. Got {ld1} instead of {expected}.'
         )
         # check if twice the mass leads to twice the mass density
@@ -126,29 +129,29 @@ def test_local_mass_density(self):
         # check result
         self.assertTrue(
             compare.all(),
-            'Doubling the mass does not lead to twice the local mass density.'\
+            'Doubling the mass does not lead to twice the local mass density.'
             f' Got {ld2} instead of {2*ld1}.'
         )
         # calculate local number density with different masses
-        mass = np.array([5.0,1.5,0.5,0.5,1.5])
-        radius = np.array([1.0,0.5,0.5,0.5,1.5])
+        mass = np.array([5.0, 1.5, 0.5, 0.5, 1.5])
+        radius = np.array([1.0, 0.5, 0.5, 0.5, 1.5])
         ld = amep.order.local_mass_density(
             self.coords, self.box, radius, mass,
-            rmax = 2.0, pbc = False, enforce_nd = 2
+            rmax=2.0, pbc=False, enforce_nd=2
         )
         # expected results for the given coords
         expected = np.array(
-            [0.65651414, 0.62722693, 0.6167254 , 0.58743819, 0.36031597]
+            [0.65651414, 0.62722693, 0.6167254, 0.58743819, 0.36031597]
         )
         # comparison
         compare = ld.round(3) == expected.round(3)
         # check result
         self.assertTrue(
             compare.all(),
-            'Incorrect local mass density for calculation with different '\
+            'Incorrect local mass density for calculation with different '
             f'particle radii and masses. Got {ld} instead of {expected}.'
         )
-            
+
     def test_local_packing_fraction(self):
         """
         Tests the `local_packing_fraction` function of `amep.order`.
@@ -160,7 +163,7 @@ def test_local_packing_fraction(self):
         """
         # calculate local packing fraction with fixed particle radius
         ld = amep.order.local_packing_fraction(
-            self.coords, self.box, 0.5, rmax = 2.0, pbc = False, enforce_nd = 2
+            self.coords, self.box, 0.5, rmax=2.0, pbc=False, enforce_nd=2
         )
         # expected results for the given coords
         expected = np.array(
@@ -171,24 +174,24 @@ def test_local_packing_fraction(self):
         # check result
         self.assertTrue(
             compare.all(),
-            'Incorrect local packing fraction for calculation with fixed '\
+            'Incorrect local packing fraction for calculation with fixed '
             f'particle radius. Got {ld} instead of {expected}.'
         )
         # result should be pi/4 times local number density
         lnd = amep.order.local_number_density(
-            self.coords, self.box, 0.5, rmax = 2.0, pbc = False, enforce_nd = 2
+            self.coords, self.box, 0.5, rmax=2.0, pbc=False, enforce_nd=2
         )
         compare = ld.round(3) == (lnd*np.pi/4).round(3)
         self.assertTrue(
             compare.all(),
-            'Local packing fraction and pi/4 times local number density do '\
+            'Local packing fraction and pi/4 times local number density do '
             f'not match. Got phi={ld} instead of rho*pi/4={lnd*np.pi/4}.'
         )
         # calculate local packing fraction with different particle radii
-        radius = np.array([0.5,0.5,0.5,0.5,1.0])
+        radius = np.array([0.5, 0.5, 0.5, 0.5, 1.0])
         ld = amep.order.local_packing_fraction(
             self.coords, self.box, radius,
-            rmax = 2.0, pbc = False, enforce_nd = 2
+            rmax=2.0, pbc=False, enforce_nd=2
         )
         # expected results for the given coords
         expected = np.array(
@@ -199,10 +202,6 @@ def test_local_packing_fraction(self):
         # check result
         self.assertTrue(
             compare.all(),
-            'Incorrect local packing fraction for calculation with different '\
+            'Incorrect local packing fraction for calculation with different '
             f'particle radii. Got {ld} instead of {expected}.'
         )
-
-
-if __name__ == '__main__':
-    unittest.main()
diff --git a/test/test_particle_methods.py b/test/test_particle_methods.py
index 086b3a7..dfab7f5 100644
--- a/test/test_particle_methods.py
+++ b/test/test_particle_methods.py
@@ -18,38 +18,13 @@
 # Contact: Lukas Hecht (lukas.hecht@pkm.tu-darmstadt.de)
 # =============================================================================
 """Test all the read and evaluate methods for particle based simulations"""
-import unittest
 from pathlib import Path
-from zipfile import ZipFile
-from requests import get
+import unittest
+import numpy as np
 from amep.load import traj
-SERVER_URL: str = "https://kuno.fkp.physik.tu-darmstadt.de/d/a3d9887b8a5747e0a56e/files/?p=/"
 URL_END: str = "&dl=1"
-DATA_DIR: Path = Path("./data/")
-PARTICLE_DIR: Path = DATA_DIR/"particles"
-FILE_NAMES: tuple[str, ...] = ("particles",
-                               )
-FILE_END: str = ".zip"
-
-
-def get_test_files():
-    """Gets and extracts needed Simulation data in order to make it accessible for the tests"""
-    for file_name in FILE_NAMES:
-        sim_dir = DATA_DIR
-        zip_path = sim_dir/Path(f"{file_name}{FILE_END}")
-        file_url = f"{SERVER_URL}{file_name}{FILE_END}{URL_END}"
-        if not sim_dir.is_dir():
-            sim_dir.mkdir()
-        try:
-            print(f"Downloading {file_name}")
-            with open(zip_path, "wb") as ofile:
-                ofile.write(get(file_url, allow_redirects=True).content)
-            print(f"unziping {file_name}")
-            with ZipFile(zip_path, "r") as z_file:
-                z_file.extractall(sim_dir)
-            zip_path.unlink()
-        except:
-            print(f"Downloading and Unziping {file_name} failed!")
+DATA_DIR: Path = Path("../examples/data/")
+PARTICLE_DIR: Path = DATA_DIR/"lammps"
 
 
 class TestParticleMethods(unittest.TestCase):
@@ -57,78 +32,63 @@ class TestParticleMethods(unittest.TestCase):
 
     @classmethod
     def setUpClass(cls):
-        """Checks whether data needed for tests is already there.
-
-        And gets it if not.
+        """Provides test case with instance of loaded particle data.
         """
-        DATA_DIR.mkdir(exist_ok=True)
-        if not PARTICLE_DIR.is_dir():
-            get_test_files()
-        cls.trajs = [traj(str(file_name)) for file_name in PARTICLE_DIR.iterdir()]
+        cls.traj = traj(PARTICLE_DIR)
 
     def test_class_methods(self):
         """Testing all the class methods of the particle trajectory"""
-        for trace in self.trajs:
-            first_frame = trace[0]
-            first_frames = trace[:3]
-            self.assertEqual(first_frame, first_frames[0],
-                             "indexing went wrong")
-        print([trace.reader.savedir for trace in self.trajs])
-        test_traj = [trace for trace in self.trajs if trace.reader.savedir ==
-                     str(PARTICLE_DIR/"kob-andersen-liquid")][0]
-        test_traj.add_param("T", 0.65)
-        test_traj.add_author_info("KLS", "name", "Kai Luca Spanheimer")
-        test_traj.add_author_info("KLS", "email",
+        first_frame = self.traj[0]
+        first_frames = self.traj[:3]
+        self.assertEqual(first_frame, first_frames[0], "indexing went wrong")
+        self.traj.add_param("T", 0.65)
+        self.traj.add_author_info("KLS", "name", "Kai Luca Spanheimer")
+        self.traj.add_author_info("KLS", "email",
                                   "kai.luca@pkm.tu-darmstadt.de")
-        test_traj.add_script(PARTICLE_DIR/Path("kob-andersen-liquid/production.run"))
-        test_traj.delete_param("T")
-        test_traj.nojump()
+        self.traj.add_script(PARTICLE_DIR/"abps.run")
+        self.traj.delete_param("T")
+        self.traj.nojump()
 
     def test_frame_methods(self):
-        for trace in self.trajs:
-            first_frame = trace[0]
-            second_frame = trace[1]
-            third_frame = trace[2]
-            self.assertNotEqual(first_frame.step, second_frame.step)
-            self.assertNotEqual(first_frame.time, second_frame.time)
-            self.assertEqual(first_frame.dim, second_frame.dim)
-            self.assertEqual(first_frame.volume, second_frame.volume)
-            self.assertEqual(third_frame.keys, second_frame.keys)
-            # print(trace.info)
-            # self.assertEqual(first_frame.center,second_frame.center)
-            # self.assertEqual(first_frame.box-second_frame.box,np.zeros_like(second_frame.box))
-            first_frame.density()
-            first_frame.n()
-            first_frame.coords()
-            """
-            with self.assertRaises(KeyError, msg=""):
-                first_frame.unwrapped_coords()
-                first_frame.nojump_coords()
-            """
-            first_frame.velocities()
-            first_frame.forces()
-            first_frame.orientations()
-            first_frame.omegas()
-            first_frame.torque()
-            first_frame.radius()
-            first_frame.ids()
-            first_frame.data("x", "vx", "omegay", "id")
-            for ptype in first_frame.ptypes:
-                ptypes = [ptype,]
-                first_frame.density(ptype=ptypes)
-                first_frame.n(ptype=ptypes)
-                first_frame.coords(ptype=ptypes)
-                # first_frame.unwrapped_coords(ptype=ptypes)
-                # first_frame.nojump_coords(ptype=ptypes)
-                first_frame.velocities(ptype=ptypes)
-                first_frame.forces(ptype=ptypes)
-                first_frame.orientations(ptype=ptypes)
-                first_frame.omegas(ptype=ptypes)
-                first_frame.torque(ptype=ptypes)
-                first_frame.radius(ptype=ptypes)
-                first_frame.ids(ptype=ptypes)
-                first_frame.data("x", "vx",
-                                 "omegay", "id", ptype=ptypes)
-
-if __name__ == '__main__':
-    unittest.main()
\ No newline at end of file
+        """Test all the relevant frame methods."""
+        first_frame = self.traj[0]
+        second_frame = self.traj[1]
+        third_frame = self.traj[2]
+        self.assertNotEqual(first_frame.step, second_frame.step)
+        self.assertNotEqual(first_frame.time, second_frame.time)
+        self.assertEqual(first_frame.dim, second_frame.dim)
+        self.assertEqual(first_frame.volume, second_frame.volume)
+        self.assertEqual(third_frame.keys, second_frame.keys)
+        np.testing.assert_equal(first_frame.center, second_frame.center)
+        np.testing.assert_equal(first_frame.box-second_frame.box,
+                                np.zeros_like(second_frame.box))
+        first_frame.density()
+        first_frame.n()
+        first_frame.coords()
+        # with self.assertRaises(KeyError, msg=""):
+        #     first_frame.unwrapped_coords()
+        #     first_frame.nojump_coords()
+        first_frame.velocities()
+        first_frame.forces()
+        first_frame.orientations()
+        first_frame.omegas()
+        first_frame.torque()
+        first_frame.radius()
+        first_frame.ids()
+        first_frame.data("x", "vx", "omegay", "id")
+        for ptype in first_frame.ptypes:
+            ptypes = [ptype,]
+            first_frame.density(ptype=ptypes)
+            first_frame.n(ptype=ptypes)
+            first_frame.coords(ptype=ptypes)
+            # first_frame.unwrapped_coords(ptype=ptypes)
+            # first_frame.nojump_coords(ptype=ptypes)
+            first_frame.velocities(ptype=ptypes)
+            first_frame.forces(ptype=ptypes)
+            first_frame.orientations(ptype=ptypes)
+            first_frame.omegas(ptype=ptypes)
+            first_frame.torque(ptype=ptypes)
+            first_frame.radius(ptype=ptypes)
+            first_frame.ids(ptype=ptypes)
+            first_frame.data("x", "vx",
+                             "omegay", "id", ptype=ptypes)
diff --git a/test/test_pbc.py b/test/test_pbc.py
index 8fa9ba0..980f03e 100644
--- a/test/test_pbc.py
+++ b/test/test_pbc.py
@@ -23,17 +23,18 @@
 # =============================================================================
 # IMPORT MODULES
 # =============================================================================
-import numpy as np
 import unittest
+import numpy as np
 import amep
 
+
 # =============================================================================
 # MAIN PBC TEST
 # =============================================================================
 class TestPbc(unittest.TestCase):
 
     @classmethod
-    def setUpClass(self):
+    def setUpClass(cls):
         """
         Basic setup. Generating test data.
 
@@ -42,7 +43,7 @@ def setUpClass(self):
         None.
 
         """
-        self.coords = np.array([
+        cls.coords = np.array([
             [2.0, 2.0, 0.0],
             [-1.0, 1.0, 0.0],
             [2.5, 3.0, 0.0],
@@ -50,7 +51,7 @@ def setUpClass(self):
             [0.0, 6.0, 0.0],
             [0.7, 0.7, 0.0]
         ])
-        self.box = np.array([[-5, 5], [-5, 5], [-0.5, 0.5]])
+        cls.box = np.array([[-5, 5], [-5, 5], [-0.5, 0.5]])
 
     def test_fold(self):
         """
@@ -303,22 +304,16 @@ def test_mirror_points(self):
             [-8.0, -13.0, -1.25]
         ])
         # comparison
-        failed2d = []
-        for mc in amep.pbc.mirror_points(coords2d, box, enforce_nd=2):
-            found = False
-            for c in mirror_coords2d:
-                if (mc.round(3)==c.round(3)).all():
-                    found = True
-            if not found:
-                failed2d.append(mc)
-        failed3d = []
-        for mc in amep.pbc.mirror_points(coords3d, box):
-            found = False
-            for c in mirror_coords3d:
-                if (mc.round(3)==c.round(3)).all():
-                    found = True
-            if not found:
-                failed3d.append(mc)
+        failed2d = [mc for mc in amep.pbc.mirror_points(coords2d,
+                                                        box,
+                                                        enforce_nd=2)
+                    if not any((mc.round(3) == c.round(3)).all()
+                               for c in mirror_coords2d)
+                    ]
+        failed3d = [mc for mc in amep.pbc.mirror_points(coords3d, box)
+                    if not any((mc.round(3) == c.round(3)).all()
+                               for c in mirror_coords3d)
+                    ]
         # check 2d data indentification
         shape = amep.pbc.mirror_points(coords2d, box, enforce_nd=2).shape
         self.assertEqual(
@@ -397,7 +392,7 @@ def test_mirror_points(self):
             f'''thickness=5.0, width=0.1 failure. width should be preferred
                 before thickness. Expected shape (1,3). Got {shape}.'''
         )
-        
+
     def test_pbc_diff(self):
         # generate test data
         coords1 = np.array([
@@ -444,7 +439,7 @@ def test_pbc_diff(self):
                 result={pbcdiff_to_test},\n
                 box={box}.'''
         )
-    
+
     def test_pbc_diff_rect(self):
         # generate test data
         coords1 = np.array([
@@ -475,19 +470,15 @@ def test_pbc_diff_rect(self):
                 result={pbcdiff_to_test},\n
                 box={box}.'''
         )
-        
+
     def test_kdtree(self):
         pass
-    
+
     def test_find_pairs(self):
         pass
-    
+
     def test_distance_matrix(self):
         pass
-    
+
     def test_distances(self):
         pass       
-
-
-if __name__ == '__main__':
-    unittest.main()
diff --git a/test/test_plot.py b/test/test_plot.py
index a47e7fb..4944715 100644
--- a/test/test_plot.py
+++ b/test/test_plot.py
@@ -19,17 +19,16 @@
 # =============================================================================
 """Test units for all functions in the plot part of amep"""
 import unittest
-from requests import get
 from pathlib import Path
-from zipfile import ZipFile
 from random import choice
 from numpy import linspace, exp
 import numpy as np
 from matplotlib.pyplot import close
 from amep import plot, load
 from amep.trajectory import FieldTrajectory, ParticleTrajectory
+from amep.base import BaseFrame
 
-DATA_DIR = Path("./data")
+DATA_DIR = Path("../examples/data")
 FIELD_DIR = DATA_DIR/"fields"
 PLOT_DIR = DATA_DIR/"plots"
 LATEX_TESCASES = {"<": r"\textless ",
@@ -39,43 +38,15 @@
                   "_": r"\_",
                   "3^5>9^2": r"3\^{}5\textgreater 9\^{}2"
                   }
-SERVER_URL: str = "https://kuno.fkp.physik.tu-darmstadt.de"
-WEB_FOLDER: str = "d/a3d9887b8a5747e0a56e/files/?p=%2F"
-FILE_NAMES: tuple[str] = ("fields", "particles")
-FILE_END: str = ".zip"
-URL_END: str = "&dl=1"
-
-
-def get_test_files():
-    for file_name in FILE_NAMES:
-        zip_path = DATA_DIR/f"{file_name}{FILE_END}"
-        file_url = f"{SERVER_URL}/{WEB_FOLDER}{file_name}{FILE_END}{URL_END}"
-        # try:
-        print(f"Downloading {file_name}")
-        with open(zip_path, "wb") as ofile:
-            ofile.write(get(file_url, allow_redirects=True).content)
-        print(f"unziping {file_name}")
-        print(zip_path)
-        with ZipFile(zip_path, "r") as z_file:
-            z_file.extractall(DATA_DIR)
-            zip_path.unlink()
-        # except Exception as exc:
-        #     print(f"Something failed while getting {file_name}.\nNamely{exc}")
 
 
 class TestPlotMethods(unittest.TestCase):
     """A test Case for the plot functions"""
     @classmethod
     def setUpClass(cls):
-        DATA_DIR.mkdir(exist_ok=True)
         PLOT_DIR.mkdir(exist_ok=True)
-
-        if not all((DATA_DIR/file_name).is_dir() for file_name in FILE_NAMES):
-            get_test_files()
-        cls.trajs = [*(load.traj(str(directory)) for directory
-                       in (DATA_DIR/"particles").iterdir()),
-                     *(load.traj(str(directory), mode="field") for directory
-                       in (DATA_DIR/"fields").iterdir())
+        cls.trajs = [load.traj(fil) for fil in DATA_DIR.iterdir()
+                     if ".h5amep" in fil.suffixes
                      ]
 
     def test_to_latex(self):
@@ -152,24 +123,25 @@ def test_animation(self):
         """Test the animation facilities of AMEP."""
         f_trajectory = choice([traj for traj in self.trajs
                                if isinstance(traj, FieldTrajectory)])
-        OUT_FIELD = PLOT_DIR/"field_vid.gif"
+        out_field = PLOT_DIR/"field_vid.gif"
         p_trajectory = choice([traj for traj in self.trajs
                                if isinstance(traj, ParticleTrajectory)])
-        OUT_PARTICLES = PLOT_DIR/"particle_vid.gif"
-        plot.animate_trajectory(f_trajectory, OUT_FIELD)
-        plot.animate_trajectory(p_trajectory, OUT_PARTICLES)
+        out_particles = PLOT_DIR/"particle_vid.gif"
+        plot.animate_trajectory(f_trajectory, out_field, ftype="c", nth=10)
+        plot.animate_trajectory(p_trajectory, out_particles, nth=50)
 
     def test_ll_video(self):
+        """Test the low level video interface."""
         trajectory = choice([traj for traj in self.trajs
                              if isinstance(traj, FieldTrajectory)])
         fig, axie = plot.new((1, 1))
         image = axie.imshow(trajectory[0].data("c"))
         fig.colorbar(image)
 
-        def make_frame(frame: np.ndarray) -> list:
+        def make_frame(frame: BaseFrame) -> list:
             now = frame.data("c")
-            nowmin = np.min(now)
-            nowmax = np.max(now)
+            nowmin: float = np.min(now)
+            nowmax: float = np.max(now)
             image.set_array(now)
             image.set_clim(nowmin, nowmax)
             return [image,]
@@ -178,6 +150,3 @@ def make_frame(frame: np.ndarray) -> list:
                           PLOT_DIR/"ll_vid.gif",
                           output_codec="libx264", bitrate=500000)
         close(fig)
-
-if __name__ == '__main__':
-    unittest.main()
diff --git a/test/test_reader.py b/test/test_reader.py
index 97bc279..025176b 100644
--- a/test/test_reader.py
+++ b/test/test_reader.py
@@ -23,30 +23,32 @@
 # =============================================================================
 # IMPORT MODULES
 # =============================================================================
-import numpy as np
+from pathlib import Path
 import unittest
+import numpy as np
 import amep
 import os
 
 # =============================================================================
 # GLOBAL CONFIG.
 # =============================================================================
-FIELDDIR = './data/field'
-LAMMPSDIR = './data/lammps'
-SAVEDIR = './data/trajs'
-EMPTYDIR = './data/emtpy'
-INVALIDFIELDDIR = './data/field/invalid'
+DATADIR = Path('../examples/data/')
+LAMMPSDIR = DATADIR/'lammps'
+FIELDDIR = DATADIR/'continuum'
+EMPTYDIR = DATADIR/'empty'
+SAVEDIR = DATADIR/'trajs'
+INVALIDFIELDDIR = DATADIR/'invalid'
 
 RNG = np.random.default_rng(1234)
 
 # for field data creation
-BOX = np.array([[-10,10], [-20,20], [-0.5,0.5]])
-SHAPE = np.array([40,80,2])
-X,Y,Z = np.meshgrid(
-    np.linspace(BOX[0,0], BOX[0,1], SHAPE[0]),
-    np.linspace(BOX[1,0], BOX[1,1], SHAPE[1]),
-    np.linspace(BOX[2,0], BOX[2,1], SHAPE[2]),
-    indexing = 'ij'
+BOX = np.array([[-10, 10], [-20, 20], [-0.5, 0.5]])
+SHAPE = np.array([40, 80, 2])
+X, Y, Z = np.meshgrid(
+        np.linspace(BOX[0, 0], BOX[0, 1], SHAPE[0]),
+        np.linspace(BOX[1, 0], BOX[1, 1], SHAPE[1]),
+        np.linspace(BOX[2, 0], BOX[2, 1], SHAPE[2]),
+        indexing='ij'
 )
 COORDS = [X.flatten(), Y.flatten(), Z.flatten()]
 TIMESTEPS = [0, 1000, 2000, 3000]
@@ -54,57 +56,13 @@
 DATA = ['id', 'x', 'y', 'z', 'vx', 'vy', 'fx', 'fy', 'mass', 'radius']
 NATOMS = 100
 
-# for lammps data creation
-
-
-# create directories if they do not exist
-if not os.path.isdir('./data'):
-    os.mkdir('./data')
-if not os.path.isdir(FIELDDIR):
-    os.mkdir(FIELDDIR)
-if not os.path.isdir(LAMMPSDIR):
-    os.mkdir(LAMMPSDIR)
-if not os.path.isdir(SAVEDIR):
-    os.mkdir(SAVEDIR)
-if not os.path.isdir(EMPTYDIR):
-    os.mkdir(EMPTYDIR)
-if not os.path.isdir(INVALIDFIELDDIR):
-    os.mkdir(INVALIDFIELDDIR)
-
-
 
 # =============================================================================
 # TEST DATA GENERATORS
 # =============================================================================
-def create_field_data():
-
-    # generate grid.txt file
-    with open(
-            os.path.join(FIELDDIR, "grid.txt"),
-            "w",
-            encoding="utf-8"
-    ) as wfile:
-        wfile.write(f'BOX:\n{BOX[0,0]}\t{BOX[0,1]}\n{BOX[1,0]}\t{BOX[1,1]}\n{BOX[2,0]}\t{BOX[2,1]}\n')
-        wfile.write('SHAPE:\n' + '\t'.join(str(s) for s in SHAPE) + '\n')
-        wfile.write('COORDINATES:\tX\tY\tZ\n')
-        wfile.write('\n'.join('\t'.join(
-            str(COORDS[i][j]) for i in range(3)
-        ) for j in range(len(COORDS[0]))))
-
-    # generate dump files with random data
-    for i, step in enumerate(TIMESTEPS):
-        with open(
-                os.path.join(FIELDDIR, f'field_{step}.txt'),
-                "w",
-                encoding="utf-8"
-        ) as wfile:
-            wfile.write(f'TIMESTEP:\n{step}\nDATA:\t'+'\t'.join(FIELDS)+'\n')
-            wfile.write('\n'.join('\t'.join(
-                str(RNG.random()) for _ in FIELDS
-            ) for _ in range(len(COORDS[0]))))
-
 def create_invalid_field_data():
-    
+    """Generate a set of invalid field data."""
+
     # INVALID BOX
     # generate grid.txt file
     with open(
@@ -112,7 +70,7 @@ def create_invalid_field_data():
             "w",
             encoding="utf-8"
     ) as wfile:
-        wfile.write(f'BOX:\n{BOX[0,0]}\t{BOX[0,1]}\n{BOX[1,0]}\t{BOX[1,1]}\n')
+        wfile.write(f'BOX:\n{BOX[0, 0]}\t{BOX[0, 1]}\n{BOX[1, 0]}\t{BOX[1, 1]}\n')
         wfile.write('SHAPE:\n' + '\t'.join(str(s) for s in SHAPE) + '\n')
         wfile.write('COORDINATES:\tX\tY\tZ\n')
         wfile.write('\n'.join('\t'.join(
@@ -128,7 +86,7 @@ def create_invalid_field_data():
         wfile.write('\n'.join('\t'.join(
             str(RNG.random()) for _ in FIELDS
         ) for _ in range(len(COORDS[0]))))
-    
+
     # INVALID SHAPE
     # generate grid.txt file
     with open(
@@ -136,7 +94,7 @@ def create_invalid_field_data():
             "w",
             encoding="utf-8"
     ) as wfile:
-        wfile.write(f'BOX:\n{BOX[0,0]}\t{BOX[0,1]}\n{BOX[1,0]}\t{BOX[1,1]}\n{BOX[2,0]}\t{BOX[2,1]}\n')
+        wfile.write(f'BOX:\n{BOX[0, 0]}\t{BOX[0, 1]}\n{BOX[1, 0]}\t{BOX[1, 1]}\n{BOX[2, 0]}\t{BOX[2, 1]}\n')
         wfile.write('COORDINATES:\tX\tY\tZ\n')
         wfile.write('\n'.join('\t'.join(
             str(COORDS[i][j]) for i in range(3)
@@ -151,7 +109,7 @@ def create_invalid_field_data():
         wfile.write('\n'.join('\t'.join(
             str(RNG.random()) for _ in FIELDS
         ) for _ in range(len(COORDS[0]))))
-        
+
     # INVALID COORDINATES
     # generate grid.txt file
     with open(
@@ -159,7 +117,7 @@ def create_invalid_field_data():
             "w",
             encoding="utf-8"
     ) as wfile:
-        wfile.write(f'BOX:\n{BOX[0,0]}\t{BOX[0,1]}\n{BOX[1,0]}\t{BOX[1,1]}\n{BOX[2,0]}\t{BOX[2,1]}\n')
+        wfile.write(f'BOX:\n{BOX[0, 0]}\t{BOX[0, 1]}\n{BOX[1, 0]}\t{BOX[1, 1]}\n{BOX[2, 0]}\t{BOX[2, 1]}\n')
         wfile.write('SHAPE:\n' + '\t'.join(str(s) for s in SHAPE) + '\n')
         wfile.write('COORDINATES:X\tY\tZ\n')
         wfile.write('\n'.join('\t'.join(
@@ -175,7 +133,7 @@ def create_invalid_field_data():
         wfile.write('\n'.join('\t'.join(
             str(RNG.random()) for _ in FIELDS
         ) for _ in range(len(COORDS[0]))))
-        
+
     # INVALID TIMESTEP
     # generate grid.txt file
     with open(
@@ -183,7 +141,7 @@ def create_invalid_field_data():
             "w",
             encoding="utf-8"
     ) as wfile:
-        wfile.write(f'BOX:\n{BOX[0,0]}\t{BOX[0,1]}\n{BOX[1,0]}\t{BOX[1,1]}\n{BOX[2,0]}\t{BOX[2,1]}\n')
+        wfile.write(f'BOX:\n{BOX[0, 0]}\t{BOX[0, 1]}\n{BOX[1, 0]}\t{BOX[1, 1]}\n{BOX[2, 0]}\t{BOX[2, 1]}\n')
         wfile.write('SHAPE:\n' + '\t'.join(str(s) for s in SHAPE) + '\n')
         wfile.write('COORDINATES:\tX\tY\tZ\n')
         wfile.write('\n'.join('\t'.join(
@@ -199,7 +157,7 @@ def create_invalid_field_data():
         wfile.write('\n'.join('\t'.join(
             str(RNG.random()) for _ in FIELDS
         ) for _ in range(len(COORDS[0]))))
-        
+
     # INVALID DATA
     # generate grid.txt file
     with open(
@@ -207,7 +165,7 @@ def create_invalid_field_data():
             "w",
             encoding="utf-8"
     ) as wfile:
-        wfile.write(f'BOX:\n{BOX[0,0]}\t{BOX[0,1]}\n{BOX[1,0]}\t{BOX[1,1]}\n{BOX[2,0]}\t{BOX[2,1]}\n')
+        wfile.write(f'BOX:\n{BOX[0, 0]}\t{BOX[0, 1]}\n{BOX[1, 0]}\t{BOX[1, 1]}\n{BOX[2, 0]}\t{BOX[2, 1]}\n')
         wfile.write('SHAPE:\n' + '\t'.join(str(s) for s in SHAPE) + '\n')
         wfile.write('COORDINATES:\tX\tY\tZ\n')
         wfile.write('\n'.join('\t'.join(
@@ -225,32 +183,13 @@ def create_invalid_field_data():
         ) for _ in range(len(COORDS[0]))))
 
 
-
-def create_lammps_data():
-    
-    # generate dump files with random data
-    for i, step in enumerate(TIMESTEPS):
-        with open(
-                os.path.join(LAMMPSDIR, f'dump{step}.txt'),
-                "w",
-                encoding="utf-8"
-        ) as wfile:
-            wfile.write(f'ITEM: TIMESTEP\n{step}\n')
-            wfile.write(f'ITEM: NUMBER OF ATOMS\n{NATOMS}\n')
-            wfile.write('ITEM: BOX BOUNDS pp pp pp\n')
-            wfile.write(f'{BOX[0,0]} {BOX[0,1]}\n{BOX[1,0]} {BOX[1,1]}\n{BOX[2,0]} {BOX[2,1]}\n')
-            wfile.write('ITEM: ATOMS '+' '.join(DATA)+'\n')
-            wfile.write('\n'.join('\t'.join(
-                str(RNG.random()) for _ in DATA
-            ) for _ in range(NATOMS)))
-
 # =============================================================================
 # MAIN READER TESTS
 # =============================================================================
 class TestContinuumReader(unittest.TestCase):
-    
+
     @classmethod
-    def setUpClass(self):
+    def setUpClass(cls):
         """
         Basic setup. Generating test data.
 
@@ -259,154 +198,156 @@ def setUpClass(self):
         None.
 
         """
-        # create test data
-        create_field_data()
-        
-        # create corrupted field data
+        INVALIDFIELDDIR.mkdir(exist_ok=True)
+        SAVEDIR.mkdir(exist_ok=True)
+        EMPTYDIR.mkdir(exist_ok=True)
         create_invalid_field_data()
-    
 
-    
     def test_init(self):
+        """Test constructor method."""
         # read data
+        (FIELDDIR/"#temp#traj.h5amep").touch()
         _ = amep.reader.ContinuumReader(
             FIELDDIR,
             FIELDDIR,
-            trajfile = 'traj.h5amep',
-            deleteold = False,
-            dumps = 'field_*.txt',
-            gridfile = 'grid.txt',
-            delimiter = '\t',
-            timestep = 1.0
+            trajfile='#temp#traj.h5amep',
+            deleteold=False,
+            dumps='field_*.txt',
+            gridfile='grid.txt',
+            delimiter=' ',
+            timestep=1.0,
+            nth=None
         )
         # test exceptions for invalid/corrupted data
         with self.assertRaises(
                 Exception,
-                msg = 'No exception raised for invalid box.'
+                msg='No exception raised for invalid box.'
         ):
             _ = amep.reader.ContinuumReader(
                 INVALIDFIELDDIR,
                 INVALIDFIELDDIR,
-                trajfile = 'field-box.h5amep',
-                deleteold = False,
-                dumps = 'field-box_*.txt',
-                gridfile = 'grid-box.txt',
-                delimiter = '\t',
-                timestep = 0.5
+                trajfile='field-box.h5amep',
+                deleteold=False,
+                dumps='field-box_*.txt',
+                gridfile='grid-box.txt',
+                delimiter=' ',
+                timestep=0.5
             )
         with self.assertRaises(
                 Exception,
-                msg = 'No exception raised for invalid shape.'
+                msg='No exception raised for invalid shape.'
         ):
             _ = amep.reader.ContinuumReader(
                 INVALIDFIELDDIR,
                 INVALIDFIELDDIR,
-                trajfile = 'field-shape.h5amep',
-                deleteold = False,
-                dumps = 'field-shape_*.txt',
-                gridfile = 'grid-shape.txt',
-                delimiter = '\t',
-                timestep = 0.5
+                trajfile='field-shape.h5amep',
+                deleteold=False,
+                dumps='field-shape_*.txt',
+                gridfile='grid-shape.txt',
+                delimiter=' ',
+                timestep=0.5
             )
         with self.assertRaises(
                 Exception,
-                msg = 'No exception raised for invalid coordinates.'
+                msg='No exception raised for invalid coordinates.'
         ):
             _ = amep.reader.ContinuumReader(
                 INVALIDFIELDDIR,
                 INVALIDFIELDDIR,
-                trajfile = 'field-coords.h5amep',
-                deleteold = False,
-                dumps = 'field-coords_*.txt',
-                gridfile = 'grid-coords.txt',
-                delimiter = '\t',
-                timestep = 0.5
+                trajfile='field-coords.h5amep',
+                deleteold=False,
+                dumps='field-coords_*.txt',
+                gridfile='grid-coords.txt',
+                delimiter=' ',
+                timestep=0.5
             )
         with self.assertRaises(
                 Exception,
-                msg = 'No exception raised for invalid timestep.'
+                msg='No exception raised for invalid timestep.'
         ):
             _ = amep.reader.ContinuumReader(
                 INVALIDFIELDDIR,
                 INVALIDFIELDDIR,
-                trajfile = 'field-step.h5amep',
-                deleteold = False,
-                dumps = 'field-step_*.txt',
-                gridfile = 'grid-step.txt',
-                delimiter = '\t',
-                timestep = 0.5
+                trajfile='field-step.h5amep',
+                deleteold=False,
+                dumps='field-step_*.txt',
+                gridfile='grid-step.txt',
+                delimiter=' ',
+                timestep=0.5
             )
         with self.assertRaises(
                 Exception,
-                msg = 'No exception raised for invalid data.'
+                msg='No exception raised for invalid data.'
         ):
             _ = amep.reader.ContinuumReader(
                 INVALIDFIELDDIR,
                 INVALIDFIELDDIR,
-                trajfile = 'field-data.h5amep',
-                deleteold = False,
-                dumps = 'field-data_*.txt',
-                gridfile = 'grid-data.txt',
-                delimiter = '\t',
-                timestep = 0.5
-            ) 
-            
+                trajfile='field-data.h5amep',
+                deleteold=False,
+                dumps='field-data_*.txt',
+                gridfile='grid-data.txt',
+                delimiter=' ',
+                timestep=0.5
+            )
+
     def test_savedir(self):
+        """Test saving method"""
         _ = amep.reader.ContinuumReader(
             FIELDDIR,
             SAVEDIR,
-            trajfile = 'field.h5amep',
-            deleteold = False,
-            dumps = 'field_*.txt',
-            gridfile = 'grid.txt',
-            delimiter = '\t',
-            timestep = 0.5
+            trajfile='traj.h5amep',
+            deleteold=False,
+            dumps='field_*.txt',
+            gridfile='grid.txt',
+            delimiter=' ',
+            timestep=0.5
         )
         self.assertTrue(
-            os.path.exists(os.path.join(SAVEDIR, 'field.h5amep')),
-            f'''Savedir error: {os.path.join(SAVEDIR, 'field.h5amep')}
+            os.path.exists(os.path.join(SAVEDIR, 'traj.h5amep')),
+            f'''Savedir error: {os.path.join(SAVEDIR, 'traj.h5amep')}
             does not exist.'''
         )
-    
+
     def test_trajfile(self):
         _ = amep.reader.ContinuumReader(
             FIELDDIR,
             FIELDDIR,
-            trajfile = 'field.h5amep',
-            deleteold = False,
-            dumps = 'field_*.txt',
-            gridfile = 'grid.txt',
-            delimiter = '\t',
-            timestep = 0.5
+            trajfile='traj.h5amep',
+            deleteold=False,
+            dumps='field_*.txt',
+            gridfile='grid.txt',
+            delimiter=' ',
+            timestep=0.5
         )
         self.assertTrue(
-            os.path.exists(os.path.join(FIELDDIR, 'field.h5amep')),
-            f'''trajfile error: {os.path.join(FIELDDIR, 'field.h5amep')}
+            os.path.exists(os.path.join(FIELDDIR, 'traj.h5amep')),
+            f'''trajfile error: {os.path.join(FIELDDIR, 'traj.h5amep')}
             does not exist.'''
         )
-    
+
     def test_deleteold(self):
+        """Test backup during creation."""
         # create trajfile
         _ = amep.reader.ContinuumReader(
             FIELDDIR,
             SAVEDIR,
-            trajfile = 'field2.h5amep',
-            deleteold = False,
-            dumps = 'field_*.txt',
-            gridfile = 'grid.txt',
-            delimiter = '\t',
-            timestep = 0.5
+            trajfile='field2.h5amep',
+            deleteold=False,
+            dumps='field_*.txt',
+            gridfile='grid.txt',
+            delimiter=' ',
+            timestep=0.5
         )
         # recreate trajfile and create backup
         _ = amep.reader.ContinuumReader(
             FIELDDIR,
             SAVEDIR,
-            trajfile = 'field2.h5amep',
-            deleteold = False,
-            dumps = 'field_*.txt',
-            gridfile = 'grid.txt',
-            delimiter = '\t',
-            timestep = 0.5
+            trajfile='field2.h5amep',
+            deleteold=False,
+            dumps='field_*.txt',
+            gridfile='grid.txt',
+            delimiter=' ',
+            timestep=0.5
         )
         self.assertTrue(
             os.path.exists(os.path.join(SAVEDIR, '#field2.h5amep')),
@@ -417,95 +358,96 @@ def test_deleteold(self):
         _ = amep.reader.ContinuumReader(
             FIELDDIR,
             SAVEDIR,
-            trajfile = 'field2.h5amep',
-            deleteold = True,
-            dumps = 'field_*.txt',
-            gridfile = 'grid.txt',
-            delimiter = '\t',
-            timestep = 0.5
+            trajfile='field2.h5amep',
+            deleteold=True,
+            dumps='field_*.txt',
+            gridfile='grid.txt',
+            delimiter=' ',
+            timestep=0.5
         )
         self.assertFalse(
             os.path.exists(os.path.join(SAVEDIR, '#field2.h5amep')),
-            f'''Backup error: backup file 
+            f'''Backup error: backup file
             {os.path.join(SAVEDIR, 'field2.h5amep')} has not been deleted.'''
         )
-    
+
     def test_dumps(self):
+        """Test dump methods."""
         with self.assertRaises(
                 Exception,
-                msg = f'''No exception raised for empty directory {EMPTYDIR}.'''
-        ):
+                msg=f'''No exception raised for empty directory {EMPTYDIR}.'''
+                ):
             _ = amep.reader.ContinuumReader(
                 EMPTYDIR,
                 EMPTYDIR,
-                trajfile = 'field.h5amep',
-                deleteold = False,
-                dumps = 'field_*.txt',
-                gridfile = 'grid.txt',
-                delimiter = '\t',
-                timestep = 0.5
-            )
+                trajfile='traj.h5amep',
+                deleteold=False,
+                dumps='field_*.txt',
+                gridfile='grid.txt',
+                delimiter=' ',
+                timestep=0.5
+                )
         self.assertFalse(
-            os.path.exists(os.path.join(EMPTYDIR, 'field.h5amep')),
-            f'''Trajectory file {os.path.join(EMPTYDIR, 'field.h5amep')}
+            os.path.exists(os.path.join(EMPTYDIR, 'traj.h5amep')),
+            f'''Trajectory file {os.path.join(EMPTYDIR, 'traj.h5amep')}
             has been created although there is no data in {EMPTYDIR}.'''
         )
-    
+
     def test_gridfile(self):
         with self.assertRaises(
                 FileNotFoundError,
-                msg = 'No exception raised for wrong grid file name.'
+                msg='No exception raised for wrong grid file name.'
         ):
             _ = amep.reader.ContinuumReader(
                 FIELDDIR,
                 FIELDDIR,
-                trajfile = 'traj.h5amep',
-                deleteold = False,
-                dumps = 'field_*.txt',
-                gridfile = 'test.txt',
-                delimiter = '\t',
-                timestep = 1.0
+                trajfile='traj.h5amep',
+                deleteold=False,
+                dumps='field_*.txt',
+                gridfile='test.txt',
+                delimiter=' ',
+                timestep=1.0
             )
-    
+
     def test_delimiter(self):
         with self.assertRaises(
                 ValueError,
-                msg = 'No ValueError raised for wrong delimiter.'
+                msg='No ValueError raised for wrong delimiter.'
         ):
             _ = amep.reader.ContinuumReader(
                 FIELDDIR,
                 FIELDDIR,
-                trajfile = 'field.h5amep',
-                deleteold = False,
-                dumps = 'field_*.txt',
-                gridfile = 'grid.txt',
-                delimiter = ' ',
-                timestep = 0.5
+                trajfile='traj.h5amep',
+                deleteold=False,
+                dumps='field_*.txt',
+                gridfile='grid.txt',
+                delimiter=',',
+                timestep=0.5
             )
-    
+
     def test_timestep(self):
         # test warning
-        with self.assertLogs(level = "WARNING"):
+        with self.assertLogs(level="WARNING"):
             reader = amep.reader.ContinuumReader(
                 FIELDDIR,
                 FIELDDIR,
-                trajfile = 'field.h5amep',
-                deleteold = False,
-                dumps = 'field_*.txt',
-                gridfile = 'grid.txt',
-                delimiter = '\t',
-                timestep = None
+                trajfile='traj.h5amep',
+                deleteold=False,
+                dumps='field_*.txt',
+                gridfile='grid.txt',
+                delimiter=' ',
+                timestep=None
             )
         # test value
         reader = amep.reader.ContinuumReader(
             FIELDDIR,
             FIELDDIR,
-            trajfile = 'field.h5amep',
-            deleteold = False,
-            dumps = 'field_*.txt',
-            gridfile = 'grid.txt',
-            delimiter = '\t',
-            timestep = 0.5
+            trajfile='traj.h5amep',
+            deleteold=False,
+            dumps='field_*.txt',
+            gridfile='grid.txt',
+            delimiter=' ',
+            timestep=0.5
         )
         self.assertEqual(
             reader.dt,
@@ -515,81 +457,69 @@ def test_timestep(self):
 
 
 class TestLammpsReader(unittest.TestCase):
-    
-    @classmethod
-    def setUpClass(self):
-        """
-        Basic setup. Generating test data.
-
-        Returns
-        -------
-        None.
 
-        """
-        create_lammps_data()
-        
     def test_savedir(self):
         _ = amep.reader.LammpsReader(
             LAMMPSDIR,
             SAVEDIR,
-            trajfile = 'lammps.h5amep',
-            deleteold = False,
-            dumps = 'dump*.txt',
+            trajfile='traj.h5amep',
+            deleteold=False,
+            dumps='dump*.txt',
         )
         self.assertTrue(
-            os.path.exists(os.path.join(SAVEDIR, 'lammps.h5amep')),
-            f'''Savedir error: {os.path.join(SAVEDIR, 'lammps.h5amep')}
+            os.path.exists(os.path.join(SAVEDIR, 'traj.h5amep')),
+            f'''Savedir error: {os.path.join(SAVEDIR, 'traj.h5amep')}
             does not exist.'''
         )
-        
+
     def test_dumps(self):
         with self.assertRaises(
                 Exception,
-                msg = f'''No exception raised for empty directory {EMPTYDIR}.'''
+                msg=f'''No exception raised for empty directory {EMPTYDIR}.'''
         ):
             _ = amep.reader.LammpsReader(
                 EMPTYDIR,
                 EMPTYDIR,
-                trajfile = 'lammps.h5amep',
-                deleteold = False,
-                dumps = 'dump*.txt'
+                trajfile='traj.h5amep',
+                deleteold=False,
+                dumps='dump*.txt'
             )
         self.assertFalse(
-            os.path.exists(os.path.join(EMPTYDIR, 'lammps.h5amep')),
-            f'''Trajectory file {os.path.join(EMPTYDIR, 'lammps.h5amep')}
+            os.path.exists(os.path.join(EMPTYDIR, 'traj.h5amep')),
+            f'''Trajectory file {os.path.join(EMPTYDIR, 'traj.h5amep')}
             has been created although there is no data in {EMPTYDIR}.'''
         )
-        
+
     def test_trajfile(self):
         _ = amep.reader.LammpsReader(
             LAMMPSDIR,
             LAMMPSDIR,
-            trajfile = 'lammps.h5amep',
-            deleteold = False,
-            dumps = 'dump*.txt'
+            trajfile='traj.h5amep',
+            deleteold=False,
+            dumps='dump*.txt'
         )
         self.assertTrue(
-            os.path.exists(os.path.join(LAMMPSDIR, 'lammps.h5amep')),
-            f'''trajfile error: {os.path.join(LAMMPSDIR, 'lammps.h5amep')}
+            os.path.exists(os.path.join(LAMMPSDIR, 'traj.h5amep')),
+            f'''trajfile error: {os.path.join(LAMMPSDIR, 'traj.h5amep')}
             does not exist.'''
         )
-    
+
     def test_deleteold(self):
         # create trajfile
         _ = amep.reader.LammpsReader(
             LAMMPSDIR,
             SAVEDIR,
-            trajfile = 'lammps2.h5amep',
-            deleteold = False,
-            dumps = 'dump*.txt'
+            trajfile='lammps2.h5amep',
+            deleteold=False,
+            dumps='dump*.txt'
         )
         # recreate trajfile and create backup
         _ = amep.reader.LammpsReader(
             LAMMPSDIR,
             SAVEDIR,
-            trajfile = 'lammps2.h5amep',
-            deleteold = False,
-            dumps = 'dump*.txt'
+            trajfile='lammps2.h5amep',
+            deleteold=False,
+            dumps='dump*.txt'
         )
         self.assertTrue(
             os.path.exists(os.path.join(SAVEDIR, '#lammps2.h5amep')),
@@ -600,20 +530,25 @@ def test_deleteold(self):
         _ = amep.reader.LammpsReader(
             LAMMPSDIR,
             SAVEDIR,
-            trajfile = 'lammps2.h5amep',
-            deleteold = True,
-            dumps = 'dump*.txt'
+            trajfile='lammps2.h5amep',
+            deleteold=True,
+            dumps='dump*.txt'
         )
         self.assertFalse(
             os.path.exists(os.path.join(SAVEDIR, '#lammps2.h5amep')),
             f'''Backup error: backup file 
             {os.path.join(SAVEDIR, 'lammps2.h5amep')} has not been deleted.'''
         )
-    
+
+
 class TestH5amepReader(unittest.TestCase):
-    
+    """Testcase for the standard reader.
+
+    TO BE IMPLEMENTED.
+    """
+
     @classmethod
-    def setUpClass(self):
+    def setUpClass(cls):
         """
         Basic setup. Generating test data.
 
@@ -623,22 +558,28 @@ def setUpClass(self):
 
         """
         pass
-    
+
+    def test_trajfile(self):
+        pass
+
+
 class TestGromacsReader(unittest.TestCase):
-    
+    """Testcase for GROMACS data.
+
+    TO BE IMPLEMENTED.
+    """
     @classmethod
     def setUpClass(self):
         """
         Basic setup. Generating test data.
 
+        TO BE IMPLEMENTED.
         Returns
         -------
         None.
 
         """
         pass
-    
-    
-if __name__ == '__main__':
-    unittest.main()
-    
+
+    def test_trajfile(self):
+        pass
diff --git a/test/test_spatialcor.py b/test/test_spatialcor.py
index 18b22cb..2bce23a 100644
--- a/test/test_spatialcor.py
+++ b/test/test_spatialcor.py
@@ -32,9 +32,11 @@
 # MAIN TIMECOR TEST
 # =============================================================================
 class TestSpatialcor(unittest.TestCase):
+    """Test case for spatial correlation functions.
+    """
 
     @classmethod
-    def setUpClass(self):
+    def setUpClass(cls):
         """
         Basic setup. Generating test data.
 
@@ -46,47 +48,46 @@ def setUpClass(self):
         # initialize random number generator with seed 0
         rng = np.random.default_rng(0)
         # generate random coordinates
-        self.coords = np.zeros((1000,3))
-        self.coords[:,:2] = rng.uniform(
-            low = -10,
-            high = 10,
-            size = (1000,2)
+        cls.coords = np.zeros((1000, 3))
+        cls.coords[:, :2] = rng.uniform(
+            low=-10,
+            high=10,
+            size=(1000, 2)
         )
         # create box
-        self.box = np.array(
-            [[-10,10],
-             [-10,10],
-             [-0.5,0.5]]
+        cls.box = np.array(
+            [[-10, 10],
+             [-10, 10],
+             [-0.5, 0.5]]
         )
 
     def test_rdf(self):
+        """Calculate and compare radial distribution functions.
+
+        Calculate rdf by different means and compare the results."""
         # calculate rdf with mode diff
         rdf_diff, _ = amep.spatialcor.rdf(
             self.coords,
             self.box,
-            nbins = 50,
-            pbc = True,
-            rmax = 5.0,
-            mode = 'diff'
+            nbins=50,
+            pbc=True,
+            rmax=5.0,
+            mode='diff'
         )
         # calculate rdf with mode kdtree
         rdf_kdtree, _ = amep.spatialcor.rdf(
             self.coords,
             self.box,
-            nbins = 50,
-            pbc = True,
-            rmax = 5.0,
-            mode = 'kdtree'
+            nbins=50,
+            pbc=True,
+            rmax=5.0,
+            mode='kdtree'
         )
         # compare results
         compare = rdf_diff.round(3) == rdf_kdtree.round(3)
         self.assertTrue(
             compare.all(),
-            'The rdf calculation with mode diff and mode kdtree are not the '\
-            'same. Got a summed difference of '\
+            'The rdf calculation with mode diff and mode kdtree are not the '
+            'same. Got a summed difference of '
             f'{np.abs(rdf_diff-rdf_kdtree).sum()}'
         )
-        
-        
-if __name__ == '__main__':
-    unittest.main()
\ No newline at end of file
diff --git a/test/test_timecor.py b/test/test_timecor.py
index cda43bf..d6029ff 100644
--- a/test/test_timecor.py
+++ b/test/test_timecor.py
@@ -23,8 +23,8 @@
 # =============================================================================
 # IMPORT MODULES
 # =============================================================================
-import numpy as np
 import unittest
+import numpy as np
 import amep
 
 
@@ -32,9 +32,10 @@
 # MAIN TIMECOR TEST
 # =============================================================================
 class TestTimecor(unittest.TestCase):
+    """Test case for time correlation functions."""
 
     @classmethod
-    def setUpClass(self):
+    def setUpClass(cls):
         """
         Basic setup. Generating test data.
 
@@ -43,18 +44,19 @@ def setUpClass(self):
         None.
 
         """
-        self.start_coords = np.array([
+        cls.start_coords = np.array([
             [1.0, 2.0, 3.0],
             [-3.0, 0.0, 5.0],
             [0.0, -2.0, 5.0]
         ])
-        self.end_coords = np.array([
+        cls.end_coords = np.array([
             [2.0, 4.0, 3.0],
             [-1.0, 1.0, 5.0],
             [4.0, -2.0, 4.0]
         ])
 
     def test_msd(self):
+        """Test calculation of the mean square displacement."""
         msd = amep.timecor.msd(
             self.start_coords,
             self.end_coords
@@ -63,12 +65,13 @@ def test_msd(self):
             np.isclose(msd, 9.0),
             f'Incorrect MSD value. Got {msd} instead of 9.0.'
         )
-    
+
     def test_acf(self):
+        """Test calculation of the autocorellation function."""
         # 1d test
         acf = amep.timecor.acf(
-            self.start_coords[:,0],
-            self.end_coords[:,0]
+            self.start_coords[:, 0],
+            self.end_coords[:, 0]
         )
         self.assertTrue(
             np.isclose(acf, 5/3),
@@ -76,8 +79,8 @@ def test_acf(self):
         )
         # 2d test
         acf = amep.timecor.acf(
-            self.start_coords[:,0:2],
-            self.end_coords[:,0:2]
+            self.start_coords[:, 0:2],
+            self.end_coords[:, 0:2]
         )
         self.assertTrue(
             np.isclose(acf, 4+5/3),
@@ -93,8 +96,9 @@ def test_acf(self):
             f'''Incorrect ACF value for 3D data. Got {acf}
             instead of {4+5/3+54/3}.'''
         )
-    
+
     def test_isf(self):
+        """Test calculation of the isotropic structure factor."""
         isf = amep.timecor.isf(self.start_coords, self.end_coords, 2*np.pi)
         self.assertTrue(
             np.isclose(isf, 0.05625849532694646),
@@ -106,8 +110,4 @@ def test_isf(self):
             np.isclose(isf, 0.032969797157540344),
             f'''Incorrect ISF value for k={3*np.pi}.
             Got {isf} instead of 0.05625849532694646.'''
-        )        
-        
-        
-if __name__ == '__main__':
-    unittest.main()
\ No newline at end of file
+        )
diff --git a/test/test_trajectory.py b/test/test_trajectory.py
index d6fb6c8..04181ad 100644
--- a/test/test_trajectory.py
+++ b/test/test_trajectory.py
@@ -23,16 +23,18 @@
 # =============================================================================
 # IMPORT MODULES
 # =============================================================================
-import numpy as np
+from pathlib import Path
 import unittest
-import amep
 import os
+import numpy as np
+import amep
 
 # =============================================================================
 # GLOBAL CONFIG.
 # =============================================================================
-FIELDDIR = './data/field'
-LAMMPSDIR = './data/lammps'
+DATADIR = Path("../examples/data")
+FIELDDIR = DATADIR/'continuum'
+LAMMPSDIR = DATADIR/'lammps'
 
 RNG = np.random.default_rng(1234)
 
@@ -51,64 +53,6 @@
 DATA = ['id', 'x', 'y', 'z', 'vx', 'vy', 'fx', 'fy', 'mass', 'radius']
 NATOMS = 100
 
-# for lammps data creation
-
-
-# create directories if they do not exist
-if not os.path.isdir('./data'):
-    os.mkdir('./data')
-if not os.path.isdir(FIELDDIR):
-    os.mkdir(FIELDDIR)
-if not os.path.isdir(LAMMPSDIR):
-    os.mkdir(LAMMPSDIR)
-
-# =============================================================================
-# TEST DATA GENERATORS
-# =============================================================================
-def create_field_data():
-    
-    # generate grid.txt file
-    with open(
-            os.path.join(FIELDDIR, "grid.txt"),
-            "w",
-            encoding="utf-8"
-    ) as wfile:
-        wfile.write(f'BOX:\n{BOX[0,0]}\t{BOX[0,1]}\n{BOX[1,0]}\t{BOX[1,1]}\n{BOX[2,0]}\t{BOX[2,1]}\n')
-        wfile.write('SHAPE:\n' + '\t'.join(str(s) for s in SHAPE) + '\n')
-        wfile.write('COORDINATES:\tX\tY\tZ\n')
-        wfile.write('\n'.join('\t'.join(
-            str(COORDS[i][j]) for i in range(3)
-        ) for j in range(len(COORDS[0]))))
-    
-    # generate dump files with random data
-    for i, step in enumerate(TIMESTEPS):
-        with open(
-                os.path.join(FIELDDIR, f'field_{step}.txt'),
-                "w",
-                encoding="utf-8"
-        ) as wfile:
-            wfile.write(f'TIMESTEP:\n{step}\nDATA:\t'+'\t'.join(FIELDS)+'\n')
-            wfile.write('\n'.join('\t'.join(
-                str(RNG.random()) for _ in FIELDS
-            ) for _ in range(len(COORDS[0]))))
-
-def create_lammps_data():
-    
-    # generate dump files with random data
-    for i, step in enumerate(TIMESTEPS):
-        with open(
-                os.path.join(LAMMPSDIR, f'dump{step}.txt'),
-                "w",
-                encoding="utf-8"
-        ) as wfile:
-            wfile.write(f'ITEM: TIMESTEP\n{step}\n')
-            wfile.write(f'ITEM: NUMBER OF ATOMS\n{NATOMS}\n')
-            wfile.write('ITEM: BOX BOUNDS pp pp pp\n')
-            wfile.write(f'{BOX[0,0]} {BOX[0,1]}\n{BOX[1,0]} {BOX[1,1]}\n{BOX[2,0]} {BOX[2,1]}\n')
-            wfile.write('ITEM: ATOMS '+' '.join(DATA)+'\n')
-            wfile.write('\n'.join('\t'.join(
-                str(RNG.random()) for _ in DATA
-            ) for _ in range(NATOMS)))
 
 # =============================================================================
 # PARTICLETRAJECTORY TESTS
@@ -116,7 +60,7 @@ def create_lammps_data():
 class TestParticleTrajectory(unittest.TestCase):
     
     @classmethod
-    def setUpClass(self):
+    def setUpClass(cls):
         """
         Basic setup. Generating test data.
 
@@ -126,22 +70,26 @@ def setUpClass(self):
 
         """
         # generate some test LAMMPS data
-        create_lammps_data()
         
         # load the data
-        self.traj = amep.load.traj(
+        cls.traj = amep.load.traj(
             LAMMPSDIR,
             mode='lammps',
             reload=True
         )
         
         # add some particle information
-        self.traj.add_particle_info(1, 'name', 'atom A')
-        self.traj.add_particle_info(1, 'type', 'carbon')
-    
-    def test_nojump(self):
-        pass
-    
+        cls.traj.add_particle_info(1, 'name', 'atom A')
+        cls.traj.add_particle_info(1, 'type', 'carbon')
+
+    def test_key_access(self):
+        self.assertIn(self.traj[0].coords(ptype=self.traj[0].ptypes[0]),
+                      self.traj[0].coords(ptype=self.traj[0].ptypes),
+                      "Something went wrong acessing data."
+                      )
+        self.assertWarns(UserWarning, self.traj[0].coords, ptype=2)
+        self.assertWarns(UserWarning, self.traj[0].coords, ptype="a")
+
     def test_add_particle_info(self):
         # add particle info
         self.traj.add_particle_info(1, 'test', 10)
@@ -167,7 +115,17 @@ def test_get_particle_info(self):
             f'''Invalid type. Got {type(self.traj.get_particle_info(1))}
             instead of dict.'''
         )
-    
+
+    def test_ptype_type(self):
+        self.traj.add_particle_info(1, "key", "value")
+        self.assertRaises(TypeError, self.traj.add_particle_info, "asdf", "key2", "value2")
+        self.assertTrue(
+            isinstance(self.traj.get_particle_info(1), dict),
+            f'''Invalid type. Got {type(self.traj.get_particle_info(1))}
+            instead of dict.'''
+        )
+        self.assertRaises(TypeError, self.traj.get_particle_info, "asdf")
+
     def test_delete_particle_info(self):
         # delete particle info
         self.traj.delete_particle_info(None)
@@ -188,7 +146,7 @@ def test_delete_particle_info(self):
 class TestFieldTrajectory(unittest.TestCase):
     
     @classmethod
-    def setUpClass(self):
+    def setUpClass(cls):
         """
         Basic setup. Generating test data.
 
@@ -197,22 +155,19 @@ def setUpClass(self):
         None.
 
         """
-        # generate some test field data
-        create_field_data()
-        
         # load the data
-        self.traj = amep.load.traj(
+        cls.traj = amep.load.traj(
             FIELDDIR,
             mode='field',
-            delimiter='\t',
+            delimiter=' ',
             dumps='field_*.txt',
             reload=True
         )
-        
+
         # add some field information
-        self.traj.add_field_info('c', 'name', 'chemicals')
-        self.traj.add_field_info('rho', 'name', 'bacterial density')
-    
+        cls.traj.add_field_info('c', 'name', 'chemicals')
+        cls.traj.add_field_info('rho', 'name', 'bacterial density')
+
     def test_add_field_info(self):
         # add field info
         self.traj.add_field_info('rho', 'test', 10)
@@ -252,6 +207,3 @@ def test_delete_field_info(self):
         # add deleted field information again
         self.traj.add_field_info('c', 'name', 'chemicals')
         self.traj.add_field_info('rho', 'name', 'bacterial density')
-
-if __name__ == '__main__':
-    unittest.main()