Repository: amepproject/amep (language: python)

Issue:
Refactor tests
I changed the data requested by the tests so we only use our own example data.
Also I increased the test coverage by adding examples to our tests.
DOC: Trajectories not referencing to BaseFrame/BaseField
### Problem description:

In the documentation of the `ParticleTrajectory` and `FieldTrajectory` there does not exist a reference to `BaseFrame` or `BaseField` respectively. This makes it hard to understand how to get information from the trajectory objects.

### Suggested content improvement:

- [ ] Add clickable link to the individual frame types (`BaseFrame` or `BaseField`) to the docstring.
- [ ] Also, add a `__get_item__` documentation.
DOC: Documentation for fit functions should include mathematical expression
### Problem description:

The fit functions do fit mathematical functions to data. But there are not given any mathematic expression since we seem to expect that people know what a Maxwell-Boltzmann distribution is.

### Suggested content improvement:

For all predefined fit functions the mathematical expression of that function should be included in the documentation.
merge release 1.0.x into main
merge version 1.0.2 into the main branch


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/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..a4388fc 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
 
@@ -1840,16 +1841,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 
@@ -3191,15 +3198,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
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..2ca04a8 100644
--- a/amep/load.py
+++ b/amep/load.py
@@ -33,7 +33,7 @@
 import os
 import h5py
 
-from .reader import LammpsReader, H5amepReader, ContinuumReader
+from .reader import LammpsReader, H5amepReader, ContinuumReader, GSDReader
 from .trajectory import ParticleTrajectory,FieldTrajectory
 from .base import TRAJFILENAME, BaseEvalData, BaseDatabase, LOADMODES
 from .base import check_path, get_module_logger
@@ -221,6 +221,16 @@ def traj(
             **kwargs
         )
         return FieldTrajectory(reader)
+    elif mode == 'gsd':
+        reader = GSDReader(
+            directory,
+            savedir,
+            trajfile = trajfile,
+            deleteold = deleteold,
+            verbose = verbose,
+            **kwargs
+        )
+        return ParticleTrajectory(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/reader.py b/amep/reader.py
index 363115a..5d09e5e 100644
--- a/amep/reader.py
+++ b/amep/reader.py
@@ -37,7 +37,7 @@
 import warnings
 import h5py
 import string
-
+from gsd import hoomd
 import numpy as np
 
 from tqdm.autonotebook import tqdm
@@ -125,6 +125,7 @@ def version(self):
         '''
         return self.__version
 
+
 class LammpsReader(BaseReader):
     '''Reads LAMMPS simulation data and writes it to an hdf5 file.
     '''
@@ -562,6 +563,7 @@ def __init__(
         finally:
             pass
 
+
     def __sorter(self, item):
         r'''
         Returns the time step of a dump file that is given
@@ -580,7 +582,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.
@@ -1084,3 +1085,470 @@ def __detect_delimiter(
         if delimiter == []:
             return None
         return delimiter
+
+
+class GSDReader(BaseReader):
+    '''Reads GSD files and writes the containing data to an hdf5 file.
+    '''
+
+    def __init__(
+            self, directory: str, savedir: str, start: float = 0.0,
+            stop: float = 1.0, nth: int = 1, filename: str = 'trajectory.gsd',
+            trajfile: str = TRAJFILENAME, deleteold: bool = False,
+            verbose: bool = False) -> None:
+        r'''
+        Reader for simulation data in the gsd file format (hoomd-blue).
+
+        TODO in progress
+        - real-time not in gsd file
+
+        Parameters
+        ----------
+        directory : str
+            Simulation directory.
+        savedir : str
+            Directory in which the trajectory file will be stored.
+        start : float, optional
+            Start reading data from this fraction of the whole trajectory.
+            The default is None.
+        stop : float, optional
+            Stop reading data from this fraction of the whole trajectory.
+            The default is None.
+        nth : int, optional
+            Read each nth dump file. The default is None.
+        filename : str, optional
+            File name of the gsd trajectory file. The default is 'trajectory.gsd'.
+        trajfile : str, optional
+            Name of the hdf5 trajectory file that is created when an object of
+            this class is initialized. The default is TRAJFILENAME.
+        deleteold : bool, optional
+            If True, an existing old h5amep trajectory file #<trajfile>
+            will be removed. The default is False.
+        verbose : bool, optional
+            If True, runtime information is printed. The default is False.
+
+        Returns
+        -------
+        None.
+
+        '''
+        # init class logger
+        self.__log = get_class_logger(__name__, self.__class__.__name__)
+        
+        # back up trajectory file
+        if os.path.exists(os.path.join(savedir, trajfile)):
+            # rename trajectory to backup filename "#+<trajfile>"
+            if os.path.exists(os.path.join(savedir, "#"+trajfile)):
+                if verbose:
+                    self.__log.info(
+                        f"Existing old file #{trajfile} will be replaced."
+                    )
+                os.remove(os.path.join(savedir, "#"+trajfile))
+            os.rename(
+                os.path.join(savedir, trajfile),
+                os.path.join(savedir, "#"+trajfile)
+            )
+            if verbose:
+                self.__log.info(
+                    f"Renamed {trajfile} to #{trajfile}."
+                )
+        # super(LammpsReader, self).__init__(directory, start, stop, nth, "#temp#"+trajfile)
+        super().__init__(os.path.abspath(savedir), start, stop, nth, "#temp#"+trajfile)
+        self.directory = directory
+
+        try:
+            # self.__filename = filename # important: self.filename ≠ self.__filename
+
+            # load gsd file
+            gsd_file = hoomd.open(os.path.join(directory,filename))
+
+            first = int(self.start*len(gsd_file))  # first frame index
+            last  = int(self.stop*len(gsd_file))   # last frame index
+            
+            # select part of gsd.hoomd.HOOMDTrajectory ( = gsd.hoomd._HOOMDTrajectoryView )
+            gsd_file = gsd_file[first:last:self.nth]
+
+            # total number of files
+            self.__nframes = len(gsd_file)
+
+            # array of number of time steps
+            steps = np.zeros(self.__nframes, dtype=int)
+
+            # load data from dump files
+            with h5py.File(os.path.join(self.savedir, "#temp#"+trajfile), 'a') as root:
+                # Set type of h5amep file to particle
+                root.attrs["type"] = "particle"
+
+                # If no ID is specified only one warning is issued to the user.
+                # Trigger idwarning is set to false:
+                idwarning=False
+
+                # loop through all dump files
+                for n, gsd_frame in enumerate(tqdm(gsd_file)):
+
+                    # get number of time steps
+                    step = gsd_frame.configuration.step
+
+                    # get total number of atoms
+                    N = gsd_frame.particles.N
+
+                    # create new group for the given step in the 'frames' group
+                    if str(step) not in root['frames'].keys():
+                        frame = root['frames'].create_group(str(step))
+                    else:
+                        frame = root['frames'][str(step)]
+
+                    # get boundaries of the simulation box
+                    gsd_box = gsd_frame.configuration.box
+                    if np.any(gsd_box[3:] != 0):
+                        raise Exception(f"GSDReader: box tilt factors {gsd_box[3:]} must be 0.")
+                    box=np.zeros((3,2))
+                    box[:,0]=-gsd_box[:3]/2
+                    box[:,1]=gsd_box[:3]/2 # hoomd-gsd box always centered around 0
+
+                    # add box to hdf5 file
+                    if 'box' not in frame.keys():
+                        frame.create_dataset('box',
+                                             (3, 2),
+                                             data=box,
+                                             dtype=DTYPE,
+                                             compression=COMPRESSION,
+                                             shuffle=SHUFFLE,
+                                             fletcher32=FLETCHER)
+                    else:
+                        frame['box'][:] = box
+
+                    # get center of the simulation box
+                    center = (box[:, 1]+box[:, 0])/2
+
+                    # add center to hdf5 file
+                    if 'center' not in frame.keys():
+                        frame.create_dataset('center',
+                                             (3,),
+                                             data=center,
+                                             dtype=DTYPE,
+                                             compression=COMPRESSION,
+                                             shuffle=SHUFFLE,
+                                             fletcher32=FLETCHER)
+                    else:
+                        frame['center'][:] = center
+
+                    print(center)
+                    continue
+
+                    # extract the data
+                    data = np.zeros((N, nparams), dtype=DTYPE)
+                    # ignore other lines (N+9 and larger)
+                    # (LAMMPS sometimes adds an empty line at the end
+                    #  which would cause problems)
+                    for l, line in enumerate(lines[9:N+9]):
+                        data[l] = np.fromstring(line, sep=' ')
+
+                    # add data to hdf5 file
+                    coords       = np.zeros((N, 3), dtype=DTYPE)
+                    uwcoords     = np.zeros((N, 3), dtype=DTYPE)  # unwrapped
+                    velocities   = np.zeros((N, 3), dtype=DTYPE)
+                    forces       = np.zeros((N, 3), dtype=DTYPE)
+                    orientations = np.zeros((N, 3), dtype=DTYPE)
+                    omegas       = np.zeros((N, 3), dtype=DTYPE)
+                    torque       = np.zeros((N, 3), dtype=DTYPE)
+                    angmom       = np.zeros((N, 3), dtype=DTYPE)
+
+                    # sort data by id if available
+                    if "id" in keys:
+                        sorting = np.argsort(data[:, keys.index("id")])
+                        data = data[sorting, :]
+                    else:
+                        # Warn user if data could not be sorted.
+                        idwarning = True
+                        # side-note from the documentation:
+                        # """Repetitions of a particular warning for the same
+                        # source location are typically suppressed."""
+
+                    unwrapped_available = False
+                    d = 0
+                    for i, key in enumerate(keys):
+                        if key == 'x':
+                            coords[:, 0] = data[:, i]
+                            if not np.all((data[:, i] == 0.0)):
+                                d += 1
+                        elif key == 'y':
+                            coords[:, 1] = data[:, i]
+                            if not np.all((data[:, i] == 0.0)):
+                                d += 1
+                        elif key == 'z':
+                            coords[:, 2] = data[:, i]
+                            if not np.all((data[:, i] == 0.0)):
+                                d += 1
+                        elif key == 'xu':
+                            uwcoords[:, 0] = data[:, i]
+                            unwrapped_available = True
+                        elif key == 'yu':
+                            uwcoords[:, 1] = data[:, i]
+                            unwrapped_available = True
+                        elif key == 'zu':
+                            uwcoords[:, 2] = data[:, i]
+                            unwrapped_available = True
+                        elif key == 'vx':
+                            velocities[:, 0] = data[:, i]
+                        elif key == 'vy':
+                            velocities[:, 1] = data[:, i]
+                        elif key == 'vz':
+                            velocities[:, 2] = data[:, i]
+                        elif key == 'fx':
+                            forces[:, 0] = data[:, i]
+                        elif key == 'fy':
+                            forces[:, 1] = data[:, i]
+                        elif key == 'fz':
+                            forces[:, 2] = data[:, i]
+                        elif key == 'mux':
+                            orientations[:, 0] = data[:, i]
+                        elif key == 'muy':
+                            orientations[:, 1] = data[:, i]
+                        elif key == 'muz':
+                            orientations[:, 2] = data[:, i]
+                        elif key == 'omegax':
+                            omegas[:, 0] = data[:, i]
+                        elif key == 'omegay':
+                            omegas[:, 1] = data[:, i]
+                        elif key == 'omegaz':
+                            omegas[:, 2] = data[:, i]
+                        elif key == 'tqx':
+                            torque[:, 0] = data[:, i]
+                        elif key == 'tqy':
+                            torque[:, 1] = data[:, i]
+                        elif key == 'tqz':
+                            torque[:, 2] = data[:, i]
+                        elif key == 'angmomx':
+                            angmom[:, 0] = data[:, i]
+                        elif key == 'angmomy':
+                            angmom[:, 1] = data[:, i]
+                        elif key == 'angmomz':
+                            angmom[:, 2] = data[:, i]
+                        elif key == 'type':
+                            if key not in frame.keys():
+                                frame.create_dataset(key,
+                                                     (N,),
+                                                     data=data[:, i],
+                                                     dtype=int,
+                                                     compression=COMPRESSION,
+                                                     shuffle=SHUFFLE,
+                                                     fletcher32=FLETCHER)
+                            else:
+                                frame[key][:] = data[:, i]
+                        elif key == 'id':
+                            if key not in frame.keys():
+                                frame.create_dataset(key,
+                                                     (N,),
+                                                     data=data[:, i],
+                                                     dtype=int,
+                                                     compression=COMPRESSION,
+                                                     shuffle=SHUFFLE,
+                                                     fletcher32=FLETCHER)
+                            else:
+                                frame[key][:] = data[:, i]
+                        else:
+                            if key not in frame.keys():
+                                frame.create_dataset(key,
+                                                     (N,),
+                                                     data=data[:, i],
+                                                     dtype=DTYPE,
+                                                     compression=COMPRESSION,
+                                                     shuffle=SHUFFLE,
+                                                     fletcher32=FLETCHER)
+                            else:
+                                frame[key][:] = data[:, i]
+
+                    if 'coords' not in frame.keys():
+                        frame.create_dataset('coords',
+                                             (N, 3),
+                                             data=coords,
+                                             dtype=DTYPE,
+                                             compression=COMPRESSION,
+                                             shuffle=SHUFFLE,
+                                             fletcher32=FLETCHER)
+                    else:
+                        frame['coords'][:] = coords
+
+                    if 'uwcoords' not in frame.keys() and unwrapped_available:
+                        frame.create_dataset('uwcoords',
+                                             (N, 3),
+                                             data=uwcoords,
+                                             dtype=DTYPE,
+                                             compression=COMPRESSION,
+                                             shuffle=SHUFFLE,
+                                             fletcher32=FLETCHER)
+                    elif 'uwcoords' in frame.keys() and unwrapped_available:
+                        frame['uwcoords'][:] = uwcoords
+
+                    if 'velocities' not in frame.keys():
+                        frame.create_dataset('velocities',
+                                             (N, 3),
+                                             data=velocities,
+                                             dtype=DTYPE,
+                                             compression=COMPRESSION,
+                                             shuffle=SHUFFLE,
+                                             fletcher32=FLETCHER)
+                    else:
+                        frame['velocities'][:] = velocities
+
+                    if 'forces' not in frame.keys():
+                        frame.create_dataset('forces',
+                                             (N, 3),
+                                             data=forces,
+                                             dtype=DTYPE,
+                                             compression=COMPRESSION,
+                                             shuffle=SHUFFLE,
+                                             fletcher32=FLETCHER)
+                    else:
+                        frame['forces'][:] = forces
+
+                    if 'orientations' not in frame.keys():
+                        frame.create_dataset('orientations',
+                                             (N, 3),
+                                             data=orientations,
+                                             dtype=DTYPE,
+                                             compression=COMPRESSION,
+                                             shuffle=SHUFFLE,
+                                             fletcher32=FLETCHER)
+                    else:
+                        frame['orientations'][:] = orientations
+
+                    if 'omegas' not in frame.keys():
+                        frame.create_dataset('omegas',
+                                             (N, 3),
+                                             data=omegas,
+                                             dtype=DTYPE,
+                                             compression=COMPRESSION,
+                                             shuffle=SHUFFLE,
+                                             fletcher32=FLETCHER)
+                    else:
+                        frame['omegas'][:] = omegas
+
+                    if 'torque' not in frame.keys():
+                        frame.create_dataset('torque',
+                                             (N, 3),
+                                             data=torque,
+                                             dtype=DTYPE,
+                                             compression=COMPRESSION,
+                                             shuffle=SHUFFLE,
+                                             fletcher32=FLETCHER)
+                    else:
+                        frame['torque'][:] = torque
+
+                    # spatial dimension
+                    frame.attrs['d'] = d
+
+                    # add type if not specified
+                    if 'type' not in frame.keys():
+                        frame.create_dataset('type',
+                                             (N,),
+                                             data=np.ones(N),
+                                             dtype=int,
+                                             compression=COMPRESSION,
+                                             shuffle=SHUFFLE,
+                                             fletcher32=FLETCHER)
+
+                    steps[n] = step
+
+                    del lines
+
+                if idwarning:
+                    # Warn user if data could not be sorted.
+                    self.__log.warning(
+                        "No IDs specified. Data may be unsorted between "\
+                        "individual frames."
+                    )
+
+            self.steps = steps
+            # get time step from log file (this also sets self.times)
+            self.dt = self.__get_timestep_from_logfile()
+            self.d = d
+        except Exception as e:
+            os.remove(os.path.join(savedir, "#temp#"+trajfile))
+            # print("LammpsReader: loading dump files failed.")
+            if "LammpsReader: no dump files in this directory." in e.args[0]:
+                # print(e)
+                # we should not raise an error in this case.
+                pass
+            # print(f"LammpsReader: error while loading dump files. {e}")
+            raise
+        else:
+            # if no exception. rename #temp#<trajfile> and delete #<trajfile>.
+            if os.path.exists(os.path.join(savedir, trajfile)) and verbose:
+                # This case should not occur if AMEP is used properly.
+                # Possible scenario:
+                # two instances of AMEP reading and writing in the same directory.
+                self.__log.info(
+                    f"File {trajfile} already exists. "\
+                    "Overwriting existing file."
+                )
+            os.rename(os.path.join(savedir, "#temp#"+trajfile), os.path.join(savedir, trajfile))
+            self.filename = trajfile
+
+            # delete old trajectory file #<trajfile> if specified by user.
+            if deleteold and os.path.exists(os.path.join(savedir, "#"+trajfile)):
+                os.remove(os.path.join(savedir, "#"+trajfile))
+                if verbose:
+                    self.__log.info(
+                        f"Deleted old trajectory file #{trajfile}"
+                    )
+        finally:
+            pass
+
+
+    def __sorter(self, item):
+        r'''
+        Returns the time step of a dump file that is given
+        in the filename of the dump file.
+
+        INPUT:
+            item: dump-file name (str)
+
+        OUTPUT:
+            time step (float)
+        '''
+        key = self.__dumps.split('*')[0]
+        basedir, basename = os.path.split(item)
+        if key=='':
+            return float(basename.split('.')[0])
+        return float(basename.split(key)[1].split('.')[0])
+
+
+    def __get_timestep_from_logfile(self):
+        r'''
+        Reads the time step from the log.lammps file.
+
+        Returns
+        -------
+        dt : float
+            Time step.
+        '''
+        if os.path.exists(os.path.join(self.directory, 'log.lammps')):
+            with open(os.path.join(
+                    self.directory, 'log.lammps'
+                ), encoding="utf-8") as f:
+                lines = f.readlines()
+                relevant = [line for line in lines if line.startswith('timestep')]
+                if len(relevant)>=1:
+                    dt = float(relevant[-1].split('timestep')[-1])
+                else:
+                    dt = 1
+                    self.__log.warning(
+                        "No timestep mentioned in logfile. Using dt=1. "\
+                        "The timestep can be set manually by traj.dt=<dt>."
+                    )
+                if len(relevant)>1:
+                    self.__log.warning(
+                        "More than one timestep found. Using last mention, "\
+                        f"dt={dt}. "\
+                        "The timestep can be set manually by traj.dt=<dt>."
+                    )
+        else:
+            self.__log.warning(
+                "No log file found. Using dt=1. "\
+                "The timestep can be set manually by traj.dt=<dt>."
+            )
+            dt = 1
+
+        return dt
diff --git a/amep/trajectory.py b/amep/trajectory.py
index 3b8379d..84625b2 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,7 +192,7 @@ def add_particle_info(
         >>> print(traj.get_particle_info())
         {1: {'name': 'active'}, 2: {'name': 'passive'}}
         >>> traj.delete_particle_info(None)
-        >>> 
+        >>>
 
         '''
         with h5py.File(os.path.join(self.reader.savedir, self.reader.filename), 'a') as root:
@@ -213,7 +215,7 @@ def get_particle_info(self, ptype: int | None = None) -> dict:
         -------
         p : dict
             Parameters.
-        
+
         Examples
         --------
         >>> import amep
@@ -223,14 +225,14 @@ 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)
-        >>> 
+        >>>
 
         '''
         with h5py.File(os.path.join(self.reader.savedir, self.reader.filename), 'r') as root:
             if ptype is None:
                 p = {}
                 for t in list(root['particles'].keys()):
-                    p[int(t)] = dict(a for a in root['particles'][t].attrs.items())
+                    p[t] = dict(a for a in root['particles'][t].attrs.items())
             elif str(ptype) in list(root['particles'].keys()):
                 p = dict(a for a in root['particles'][str(ptype)].attrs.items())
             else:
@@ -320,6 +322,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..8d8fb31 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
 
@@ -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/datastructures/index.rst b/doc/source/datastructures/index.rst
index 0a3fb5c..1af028d 100644
--- a/doc/source/datastructures/index.rst
+++ b/doc/source/datastructures/index.rst
@@ -11,3 +11,7 @@ the sections below.
 .. toctree::
    
    field_data
+
+.. toctree::
+
+    particle_data_reader
diff --git a/doc/source/datastructures/particle_data_reader.md b/doc/source/datastructures/particle_data_reader.md
new file mode 100644
index 0000000..e77ffe4
--- /dev/null
+++ b/doc/source/datastructures/particle_data_reader.md
@@ -0,0 +1,20 @@
+# Creating a new data reader
+
+If you want to implement a reader for a particle based simulation data format.
+You need to read out the following data:
+
+
+For each timestep:
+1. Timestep
+2. Simulation time
+3. Simulation Box
+4. Center of the Box
+5. particle coordinates
+6. Unwrapped particle coordinates
+7. particle velocities
+8. forces
+9. orientations
+10. rotation speeds
+11. torques
+12. Angular momenta
+
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..0996a6d 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -33,10 +33,11 @@ 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"
+    "tqdm >= 4.65.0",
+    "gsd >= 3.2.0"
 ]
 
 [project.urls]


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..5ba16b6 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=2000, nybins=2000,
+                      njobs=4, skip=0.9
+                      )
+        pcf2d.save(RESULT_DIR/"pcf2d.h5")
+        pcfangle = PCFangle(
+            traj, nav=2, ndbins=1000, nabins=1000,
+            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=5, 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..7a3a05b 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")
+        plot.animate_trajectory(p_trajectory, out_particles)
 
     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..f582cb4 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,21 @@ 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.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.assertTrue(
+            isinstance(self.traj.get_particle_info("asdf"), dict),
+            f'''Invalid type. Got {type(self.traj.get_particle_info(1))}
+            instead of dict.'''
+        )
+
     def test_delete_particle_info(self):
         # delete particle info
         self.traj.delete_particle_info(None)
@@ -188,7 +150,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 +159,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 +211,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()