diff --git a/.github/workflows/workflow.yml b/.github/workflows/workflow.yml
index e2d0a90..6dd3a39 100644
--- a/.github/workflows/workflow.yml
+++ b/.github/workflows/workflow.yml
@@ -24,9 +24,6 @@ jobs:
     - name: Check formatting with black
       run: |
         poetry run black --check .
-    - name: Docformatter
-      run: |
-        poetry run docformatter --check -r blobmodel/
 
   pytest:
 
diff --git a/blobmodel/blobs.py b/blobmodel/blobs.py
index 5152ef2..0f1177b 100644
--- a/blobmodel/blobs.py
+++ b/blobmodel/blobs.py
@@ -1,3 +1,5 @@
+"""This module defines a Blob class and related functions for discretizing and manipulating blobs."""
+
 import warnings
 from typing import Tuple, Union
 from nptyping import NDArray
@@ -7,7 +9,7 @@
 
 
 class Blob:
-    """A single blob."""
+    """Define a single blob."""
 
     def __init__(
         self,
@@ -16,8 +18,8 @@ def __init__(
         amplitude: float,
         width_prop: float,
         width_perp: float,
-        velocity_x: float,
-        velocity_y: float,
+        v_x: float,
+        v_y: float,
         pos_x: float,
         pos_y: float,
         t_init: float,
@@ -26,14 +28,51 @@ def __init__(
         perp_shape_parameters: dict = None,
         blob_alignment: bool = True,
     ) -> None:
+
+        """
+        Initialize a single blob.
+
+        Parameters
+        ----------
+        blob_id : int
+            Identifier for the blob.
+        blob_shape : AbstractBlobShape
+            Shape of the blob.
+        amplitude : float
+            Amplitude of the blob.
+        width_prop : float
+            Width of the blob in the propagation direction.
+        width_perp : float
+            Width of the blob in the perpendicular direction.
+        v_x : float
+            Velocity of the blob in the x-direction.
+        v_y : float
+            Velocity of the blob in the y-direction.
+        pos_x : float
+            Initial position of the blob in the x-direction.
+        pos_y : float
+            Initial position of the blob in the y-direction.
+        t_init : float
+            Initial time of the blob.
+        t_drain : Union[float, NDArray]
+            Time scale for the blob to drain.
+        prop_shape_parameters : dict
+            Additional shape parameters for the propagation direction.
+        perp_shape_parameters : dict
+            Additional shape parameters for the perpendicular direction.
+        blob_alignment : bool, optional
+            If blob_aligment == True, the blob shapes are rotated in the propagation direction of the blob
+            If blob_aligment == False, the blob shapes are independent of the propagation direction
+
+        """
         self.int = int
         self.blob_id = blob_id
         self.blob_shape = blob_shape
         self.amplitude = amplitude
         self.width_prop = width_prop
         self.width_perp = width_perp
-        self.v_x = velocity_x
-        self.v_y = velocity_y
+        self.v_x = v_x
+        self.v_y = v_y
         self.pos_x = pos_x
         self.pos_y = pos_y
         self.t_init = t_init
@@ -45,7 +84,7 @@ def __init__(
             {} if perp_shape_parameters is None else perp_shape_parameters
         )
         self.blob_alignment = blob_alignment
-        self.theta = cmath.phase(self.v_x + self.v_y * 1j) if blob_alignment else 0.0
+        self._theta = cmath.phase(self.v_x + self.v_y * 1j) if blob_alignment else 0.0
 
     def discretize_blob(
         self,
@@ -58,13 +97,27 @@ def discretize_blob(
     ) -> NDArray:
         """
         Discretize blob on grid. If one_dimensional the perpendicular pulse shape is ignored.
-        The following blob shapes are implemented:
-                gauss: 2D gaussian function
-                exp: one sided exponential in x and gaussian in y
 
-                Returns
-                -------
-                discretized blob on 3d array with dimensions x,y and t : np.array
+        Parameters
+        ----------
+        x : NDArray
+            Grid coordinates in the x-direction.
+        y : NDArray
+            Grid coordinates in the y-direction.
+        t : NDArray
+            Time coordinates.
+        Ly : float
+            Length of domain in the y-direction.
+        periodic_y : bool, optional
+            Flag indicating periodicity in the y-direction (default: False).
+        one_dimensional : bool, optional
+            Flag indicating a one-dimensional blob (default: False).
+
+        Returns
+        -------
+        discretized_blob : NDArray
+            Discretized blob on a 3D array with dimensions (x, y, t).
+
         """
         # If one_dimensional, then Ly should be 0.
         assert (one_dimensional and Ly == 0) or not one_dimensional
@@ -73,17 +126,17 @@ def discretize_blob(
             warnings.warn("blob width big compared to Ly")
 
         x_perp, y_perp = self._rotate(
-            origin=(self.pos_x, self.pos_y), x=x, y=y, angle=-self.theta
+            origin=(self.pos_x, self.pos_y), x=x, y=y, angle=-self._theta
         )
         if not periodic_y or one_dimensional:
             return self._single_blob(
                 x_perp, y_perp, t, Ly, periodic_y, one_dimensional=one_dimensional
             )
-        if self.theta == 0:
+        if self._theta == 0:
             number_of_y_propagations = 0
         else:
-            x_border = (Ly - self.pos_y) / np.sin(self.theta)
-            adjusted_Ly = Ly / np.sin(self.theta)
+            x_border = (Ly - self.pos_y) / np.sin(self._theta)
+            adjusted_Ly = Ly / np.sin(self._theta)
             prop_dir = (
                 self._prop_dir_blob_position(t)
                 if type(t) in [int, float]  # t has dimensionality = 0, used for testing
@@ -103,8 +156,8 @@ def discretize_blob(
                 Ly,
                 periodic_y,
                 number_of_y_propagations,
-                x_offset=Ly * np.sin(self.theta),
-                y_offset=Ly * np.cos(self.theta),
+                x_offset=Ly * np.sin(self._theta),
+                y_offset=Ly * np.cos(self._theta),
                 one_dimensional=one_dimensional,
             )
             + self._single_blob(
@@ -114,15 +167,15 @@ def discretize_blob(
                 Ly,
                 periodic_y,
                 number_of_y_propagations,
-                x_offset=-Ly * np.sin(self.theta),
-                y_offset=-Ly * np.cos(self.theta),
+                x_offset=-Ly * np.sin(self._theta),
+                y_offset=-Ly * np.cos(self._theta),
                 one_dimensional=one_dimensional,
             )
         )
 
     def _single_blob(
         self,
-        x_perp: NDArray,
+        x_prop: NDArray,
         y_perp: NDArray,
         t: NDArray,
         Ly: float,
@@ -132,11 +185,41 @@ def _single_blob(
         y_offset: NDArray = 0,
         one_dimensional: bool = False,
     ) -> NDArray:
+        """
+        Calculate the discretized blob for a single blob instance.
+
+        Parameters
+        ----------
+        x_prop : NDArray
+            Propagation direction coordinates.
+        y_perp : NDArray
+            Perpendicular coordinates.
+        t : NDArray
+            Time coordinates.
+        Ly : float
+            Length of domain in the y-direction.
+        periodic_y : bool
+            Flag indicating periodicity in the y-direction.
+        number_of_y_propagations : NDArray, optional
+            Number of times the blob propagates through the domain in y-direction (default: 0).
+        x_offset : NDArray, optional
+            Offset in the x-direction (default: 0).
+        y_offset : NDArray, optional
+            Offset in the y-direction (default: 0).
+        one_dimensional : bool, optional
+            Flag indicating a one-dimensional blob (default: False).
+
+        Returns
+        -------
+        blob : NDArray
+            Discretized blob.
+
+        """
         return (
             self.amplitude
             * self._drain(t)
             * self._propagation_direction_shape(
-                x_perp + x_offset,
+                x_prop + x_offset,
                 t,
                 Ly,
                 periodic_y,
@@ -156,6 +239,20 @@ def _single_blob(
         )
 
     def _drain(self, t: NDArray) -> NDArray:
+        """
+        Calculate the drain factor for the blob.
+
+        Parameters
+        ----------
+        t : NDArray
+            Time coordinates.
+
+        Returns
+        -------
+        drain_factor : NDArray
+            Drain factor.
+
+        """
         if isinstance(self.t_drain, (int, float)):
             return np.exp(-(t - self.t_init) / float(self.t_drain))
         return np.exp(-(t - self.t_init) / self.t_drain[np.newaxis, :, np.newaxis])
@@ -168,11 +265,33 @@ def _propagation_direction_shape(
         periodic_y: bool,
         number_of_y_propagations: NDArray,
     ) -> NDArray:
+        """
+        Calculate the shape in the propagation direction.
+
+        Parameters
+        ----------
+        x : NDArray
+            Coordinates in the x-direction.
+        t : NDArray
+            Time coordinates.
+        Ly : float
+            Length of domain in the y-direction.
+        periodic_y : bool
+            Flag indicating periodicity in the y-direction.
+        number_of_y_propagations : NDArray
+            Number of times the blob propagates through the domain in y-direction.
+
+        Returns
+        -------
+        shape : NDArray
+            Shape in the propagation direction.
+
+        """
         if periodic_y:
             x_diffs = (
                 x
                 - self._prop_dir_blob_position(t)
-                + number_of_y_propagations * Ly * np.sin(self.theta)
+                + number_of_y_propagations * Ly * np.sin(self._theta)
             )
         else:
             x_diffs = x - self._prop_dir_blob_position(t)
@@ -189,11 +308,31 @@ def _perpendicular_direction_shape(
         periodic_y: bool,
         number_of_y_propagations: NDArray,
     ) -> NDArray:
+        """
+        Calculate the shape in the perpendicular direction.
+
+        Parameters
+        ----------
+        y : NDArray
+            Coordinates in the y-direction.
+        Ly : float
+            Length of domain in the y-direction.
+        periodic_y : bool
+            Flag indicating periodicity in the y-direction.
+        number_of_y_propagations : NDArray
+            Number of times the blob propagates through the domain in y-direction.
+
+        Returns
+        -------
+        shape : NDArray
+            Blob shape in the perpendicular direction.
+
+        """
         if periodic_y:
             y_diffs = (
                 y
                 - self._perp_dir_blob_position(t)
-                + number_of_y_propagations * Ly * np.cos(self.theta)
+                + number_of_y_propagations * Ly * np.cos(self._theta)
             )
         else:
             y_diffs = y - self._perp_dir_blob_position(t)
@@ -203,6 +342,20 @@ def _perpendicular_direction_shape(
         )
 
     def _prop_dir_blob_position(self, t: NDArray) -> NDArray:
+        """
+        Calculate the position of the blob in the propagation direction.
+
+        Parameters
+        ----------
+        t : NDArray
+            Time coordinates.
+
+        Returns
+        -------
+        position : NDArray
+            Position of the blob in the propagation direction.
+
+        """
         return (
             self.pos_x + (self.v_x**2 + self.v_y**2) ** 0.5 * (t - self.t_init)
             if self.blob_alignment
@@ -210,6 +363,20 @@ def _prop_dir_blob_position(self, t: NDArray) -> NDArray:
         )
 
     def _perp_dir_blob_position(self, t: NDArray) -> float:
+        """
+        Return the position of the blob in the perpendicular direction.
+
+        Parameters
+        ----------
+        t : NDArray
+            Time coordinates.
+
+        Returns
+        -------
+        position : float
+            Position of the blob in the perpendicular direction.
+
+        """
         return (
             self.pos_y
             if self.blob_alignment
@@ -219,6 +386,26 @@ def _perp_dir_blob_position(self, t: NDArray) -> float:
     def _rotate(
         self, origin: Tuple[float, float], x: NDArray, y: NDArray, angle: float
     ) -> Tuple[NDArray, NDArray]:
+        """
+        Rotate the coordinates around a given origin point.
+
+        Parameters
+        ----------
+        origin : Tuple[float, float]
+            Origin point of rotation.
+        x : NDArray
+            Coordinates in the x-direction.
+        y : NDArray
+            Coordinates in the y-direction.
+        angle : float
+            Rotation angle.
+
+        Returns
+        -------
+        rotated_coordinates : Tuple[NDArray, NDArray]
+            Rotated coordinates.
+
+        """
         ox, oy = origin
         px, py = x, y
 
diff --git a/blobmodel/geometry.py b/blobmodel/geometry.py
index 73c2be5..eb91f0e 100644
--- a/blobmodel/geometry.py
+++ b/blobmodel/geometry.py
@@ -1,3 +1,5 @@
+"""This module defines the Geometry class for creating a grid for the Model."""
+
 from typing import Any
 from nptyping import NDArray, Float
 import numpy as np
@@ -17,16 +19,24 @@ def __init__(
         periodic_y: bool,
     ) -> None:
         """
-        Attributes
+        Initialize a Geometry object.
+
+        Parameters
         ----------
-        Nx: int, grid points in x
-        Ny: int, grid points in y
-        Lx: float, length of grid in x
-        Ly: float, length of grid in y
-        dt: float, time step
-        T: float, time length
-        periodic_y: bool, optional
-            allow periodicity in y-direction
+        Nx : int
+            Number of grid points in the x-direction.
+        Ny : int
+            Number of grid points in the y-direction.
+        Lx : float
+            Length of domain in the x-direction.
+        Ly : float
+            Length of domain in the y-direction.
+        dt : float
+            Time step.
+        T : float
+            Time length.
+        periodic_y : bool
+            Flag indicating whether periodicity is allowed in the y-direction.
         """
         self.Nx = Nx
         self.Ny = Ny
@@ -48,7 +58,14 @@ def __init__(
         )
 
     def __str__(self) -> str:
-        """string representation of Geometry."""
+        """
+        Return a string representation of the Geometry object.
+
+        Returns
+        -------
+        str
+            String representation of the Geometry object.
+        """
         return (
             f"Geometry parameters:  Nx:{self.Nx},  Ny:{self.Ny}, Lx:{self.Lx}, Ly:{self.Ly}, "
             + f"dt:{self.dt}, T:{self.T}, y-periodicity:{self.periodic_y}"
diff --git a/blobmodel/model.py b/blobmodel/model.py
index 8e9fd03..2e3154c 100644
--- a/blobmodel/model.py
+++ b/blobmodel/model.py
@@ -1,3 +1,5 @@
+"""This module defines a 2D model of propagating blobs."""
+
 import numpy as np
 import xarray as xr
 from tqdm import tqdm
@@ -31,37 +33,60 @@ def __init__(
         one_dimensional: bool = False,
     ) -> None:
         """
-        Attributes
-        ----------
-        Nx: int, grid points in x
-        Ny: int, grid points in y
-        Lx: float, length of grid in x
-        Ly: float, length of grid in y
-        dt: float, time step
-        T: float, time length
-        periodic_y: bool, optional
-            allow periodicity in y-direction
+        Initialize the 2D Model of propagating blobs.
 
+        Parameters
+        ----------
+        Nx : int, optional
+            Number of grid points in x.
+        Ny : int, optional
+            Number of grid points in y.
+        Lx : float, optional
+            Length of the domain in x.
+        Ly : float, optional
+            Length of the domain in y.
+        dt : float, optional
+            Time step.
+        T : float, optional
+            Time length.
+        periodic_y : bool, optional
+            Allow periodicity in the y-direction.
             Important: only good approximation for Ly >> blob width
-        num_blobs:
-            number of blobs
-        blob_shape: AbstractBlobShape or str, optional
-            see AbstractBlobShape and BlobShapeImpl dataclass for available shapes
-        t_drain: float or array of length Nx, optional
-            drain time for blobs
-        blob_factory: BlobFactory, optional
-            sets distributions of blob parameters
-        labels: str, optional
+        blob_shape : AbstractBlobShape or str, optional
+            Shape of the blobs. Can be an instance of AbstractBlobShape or a string
+            specifying the shape.
+        num_blobs : int, optional
+            Number of blobs.
+        t_drain : float or array-like, optional
+            Drain time for the blobs. Can be a single float value or an array-like
+            of length Nx.
+        blob_factory : BlobFactory, optional
+            BlobFactory instance for setting blob parameter distributions.
+        labels : str, optional
+            Blob label setting. Possible values: "off", "same", "individual".
             "off": no blob labels returned
             "same": regions where blobs are present are set to label 1
             "individual": different blobs return individual labels
-            used for creating training data for supervised machine learning algorithms
-        label_border: float, optional
-            defines region of blob as region where density >= label_border * amplitude of Blob
-            only used if labels = "same" or "individual"
-        one_dimensional: bool, optional
-            If True, the perpendicular shape of the blobs will be discarded. Parameters
-            for the y-component (Ny and Ly) will be overwritten to Ny=1, Ly=0.
+            Used for creating training data for supervised machine learning algorithms
+        label_border : float, optional
+            Defines region of blob as region where density >= label_border * amplitude of Blob
+            Only used if labels = "same" or "individual"
+        one_dimensional : bool, optional
+            If True, the perpendicular shape of the blobs will be discarded.
+            Parameters for the y-component (Ny and Ly) will be overwritten to Ny=1, Ly=0.
+
+        Raises
+        ------
+        AssertionError
+            If t_drain is not a single value or an array-like of length Nx.
+
+        Warns
+        -----
+        UserWarning
+            If the model is one-dimensional and Ny and Ly are not appropriate.
+
+        UserWarning
+            If the model is one-dimensional and the blob factory is not one-dimensional.
         """
         self._one_dimensional = one_dimensional
         if self._one_dimensional and (Ny != 1 or Ly != 0):
@@ -99,17 +124,32 @@ def __init__(
         self._reset_fields()
 
     def __str__(self) -> str:
-        """string representation of Model."""
+        """
+        Return a string representation of the Model.
+
+        Returns
+        -------
+        str
+            String representation of the Model.
+        """
         return (
             f"2d Blob Model with"
             + f" num_blobs:{self.num_blobs} and t_drain:{self.t_drain}"
         )
 
     def get_blobs(self) -> List[Blob]:
-        """Returns blobs list.
+        """
+        Return the list of blobs.
+
+        Returns
+        -------
+        List[Blob]
+            List of Blob objects.
+
+        Notes
+        -----
+        - Note that if Model.sample_blobs has not been called, the list will be empty
 
-        Note that if Model.sample_blobs has not been called, the list
-        will be empty
         """
         return self._blobs
 
@@ -119,24 +159,27 @@ def make_realization(
         speed_up: bool = False,
         error: float = 1e-10,
     ) -> xr.Dataset:
-        """Integrate Model over time and write out data as xarray dataset.
+        """
+        Integrate the Model over time and write out data as an xarray dataset.
 
         Parameters
         ----------
-        file_name: str, optional
-            file name for .nc file containing data as xarray dataset
-        speed_up: bool, optional
-            speeding up code by discretizing each single blob at smaller time window
+        file_name : str, optional
+            File name for the .nc file containing data as an xarray dataset.
+        speed_up : bool, optional
+            Flag for speeding up the code by discretizing each single blob at smaller time window
             when blob values fall under given error value the blob gets discarded
-            !!!  this is only a good approximation for blob_shape='exp' !!!
-
-        error: float, optional
-            numerical error at x = Lx when blob gets truncated
-            only used if speed_up = True
+        error : float, optional
+            Numerical error at x = Lx when the blob gets truncated.
 
         Returns
-        ----------
-            xarray dataset with result data
+        -------
+        xr.Dataset
+            xarray dataset with the resulting data.
+
+        Notes
+        -----
+        - speed_up is only a good approximation for blob_shape="exp"
         """
 
         # Reset density field
@@ -161,6 +204,14 @@ def make_realization(
         return dataset
 
     def _create_xr_dataset(self) -> xr.Dataset:
+        """
+        Create an xarray dataset from the density field.
+
+        Returns
+        -------
+        xr.Dataset
+            xarray dataset with the density field data.
+        """
         if self._geometry.Ly == 0:
             dataset = xr.Dataset(
                 data_vars=dict(
@@ -195,6 +246,18 @@ def _sum_up_blobs(
         speed_up: bool,
         error: float,
     ):
+        """
+        Sum up the contribution of a single blob to the density field.
+
+        Parameters
+        ----------
+        blob : Blob
+            Blob object.
+        speed_up : bool
+            Flag for speeding up the code by discretizing each single blob at a smaller time window.
+        error : float
+            Numerical error when the blob gets truncated.
+        """
         _start, _stop = self._compute_start_stop(blob, speed_up, error)
         _single_blob = blob.discretize_blob(
             x=self._geometry.x_matrix[:, :, _start:_stop],
@@ -221,6 +284,23 @@ def _sum_up_blobs(
             ] = (blob.blob_id + 1)
 
     def _compute_start_stop(self, blob: Blob, speed_up: bool, error: float):
+        """
+        Compute the start and stop indices for summing up the contribution of a single blob.
+
+        Parameters
+        ----------
+        blob : Blob
+            Blob object.
+        speed_up : bool
+            Flag for speeding up the code by discretizing each single blob at a smaller time window.
+        error : float
+            Numerical error when the blob gets truncated.
+
+        Returns
+        -------
+        Tuple[int, int]
+            Start and stop indices.
+        """
         if not speed_up or blob.v_x == 0:
             return 0, self._geometry.t.size
         start = np.maximum(
@@ -250,6 +330,7 @@ def _compute_start_stop(self, blob: Blob, speed_up: bool, error: float):
         return start, stop
 
     def _reset_fields(self):
+        """Reset the density and labels fields."""
         self._density = np.zeros(
             shape=(self._geometry.Ny, self._geometry.Nx, self._geometry.t.size)
         )
diff --git a/blobmodel/plotting.py b/blobmodel/plotting.py
index 58060d0..f26a7e8 100644
--- a/blobmodel/plotting.py
+++ b/blobmodel/plotting.py
@@ -1,3 +1,5 @@
+"""This module provides functions to create and display animations of model output."""
+
 import matplotlib.pyplot as plt
 from mpl_toolkits.axes_grid1 import make_axes_locatable
 import numpy as np
@@ -13,22 +15,32 @@ def show_model(
     gif_name: str = "blobs.gif",
     fps: int = 10,
 ) -> None:
-    """Show animation of Model output.
+    """
+    Creates an animation that shows the evolution of a specific variable over time.
 
     Parameters
     ----------
-    dataset: xarray Dataset,
-        Model data
-    variable: str, optional
-        variable to be animated
-    interval: int, optional
-        time interval between frames in ms
-    save: bool, optional
-        if True save animation as gif
-    gif_name: str, optional
-        set name for gif
-    fps: int, optional
-        set fps for gif
+    dataset : xr.Dataset
+        Model data.
+    variable : str, optional
+        Variable to be animated (default: "n").
+    interval : int, optional
+        Time interval between frames in milliseconds (default: 100).
+    save : bool, optional
+        If True, save the animation as a GIF (default: False).
+    gif_name : str, optional
+        Set the name for the saved GIF (default: "blobs.gif").
+    fps : int, optional
+        Set the frames per second for the saved GIF (default: 10).
+
+    Returns
+    -------
+    None
+
+    Notes
+    -----
+    - This function chooses between a 1D and 2D visualizations based on the dimensionality of the dataset.
+
     """
     fig = plt.figure()
 
@@ -41,12 +53,38 @@ def show_model(
         frames.append(frame)
 
     def animate_1d(i: int) -> None:
+        """
+        Create the 1D plot for each frame of the animation.
+
+        Parameters
+        ----------
+        i : int
+            Frame index.
+
+        Returns
+        -------
+        None
+
+        """
         x = dataset.x
         y = frames[i]
         line.set_data(x, y)
         plt.title(f"t = {i*dt:.2f}")
 
     def animate_2d(i: int) -> None:
+        """
+        Create the 2D plot for each frame of the animation.
+
+        Parameters
+        ----------
+        i : int
+            Frame index.
+
+        Returns
+        -------
+        None
+
+        """
         arr = frames[i]
         vmax = np.max(arr)
         vmin = np.min(arr)
@@ -71,6 +109,24 @@ def animate_2d(i: int) -> None:
 
 
 def _setup_1d_plot(dataset, variable):
+    """
+    Set up a 1D plot for the animation.
+
+    Parameters
+    ----------
+    dataset : xr.Dataset
+        Model data.
+    variable : str
+        Variable to be animated.
+
+    Returns
+    -------
+    line : matplotlib.lines.Line2D
+        Line object representing the plot.
+    tx : matplotlib.text.Text
+        Text object for the plot title.
+
+    """
     ax = plt.axes(xlim=(0, dataset.x[-1]), ylim=(0, dataset[variable].max()))
     tx = ax.set_title(r"$t = 0$")
     (line,) = ax.plot([], [], lw=2)
@@ -80,6 +136,24 @@ def _setup_1d_plot(dataset, variable):
 
 
 def _setup_2d_plot(fig, cv0):
+    """
+    Set up a 2D plot for the animation.
+
+    Parameters
+    ----------
+    fig : matplotlib.figure.Figure
+        Figure object for the plot.
+    cv0 : numpy.ndarray
+        Initial 2D array for the plot.
+
+    Returns
+    -------
+    im : matplotlib.image.AxesImage
+        Image object representing the plot.
+    tx : matplotlib.text.Text
+        Text object for the plot title.
+
+    """
     ax = fig.add_subplot(111)
     tx = ax.set_title("t = 0")
     div = make_axes_locatable(ax)
diff --git a/blobmodel/pulse_shape.py b/blobmodel/pulse_shape.py
index f825cc4..ca05c64 100644
--- a/blobmodel/pulse_shape.py
+++ b/blobmodel/pulse_shape.py
@@ -1,3 +1,5 @@
+"""This module defines classes for blob pulse shapes used in two-dimensional simulations."""
+
 from abc import ABC, abstractmethod
 import numpy as np
 
@@ -22,6 +24,15 @@ class BlobShapeImpl(AbstractBlobShape):
     """Implementation of the AbstractPulseShape class."""
 
     def __init__(self, pulse_shape_prop="gauss", pulse_shape_perp="gauss"):
+        """Initialize the BlobShapeImpl object.
+
+        Attributes
+        ----------
+        pulse_shape_prop : str, optional
+            Type of pulse shape in the propagation direction, by default "gauss".
+        pulse_shape_perp : str, optional
+            Type of pulse shape perpendicular to the propagation direction, by default "gauss".
+        """
         if (
             pulse_shape_prop not in BlobShapeImpl.__GENERATORS.keys()
             or pulse_shape_perp not in BlobShapeImpl.__GENERATORS.keys()
@@ -33,21 +44,93 @@ def __init__(self, pulse_shape_prop="gauss", pulse_shape_perp="gauss"):
         self.get_pulse_shape_perp = BlobShapeImpl.__GENERATORS.get(pulse_shape_perp)
 
     def get_pulse_shape_prop(self, theta: np.ndarray, **kwargs) -> np.ndarray:
+        """Compute the pulse shape in the propagation direction.
+
+        Parameters
+        ----------
+        theta : np.ndarray
+            Array of theta values.
+        **kwargs
+            Additional keyword arguments.
+
+        Returns
+        -------
+        np.ndarray
+            Array representing the pulse shape in the propagation direction.
+        """
         raise NotImplementedError
 
     def get_pulse_shape_perp(self, theta: np.ndarray, **kwargs) -> np.ndarray:
+        """Compute the pulse shape perpendicular to the propagation direction.
+
+        Parameters
+        ----------
+        theta : np.ndarray
+            Array of theta values.
+        **kwargs
+            Additional keyword arguments.
+
+        Returns
+        -------
+        np.ndarray
+            Array representing the pulse shape perpendicular to the propagation direction.
+        """
         raise NotImplementedError
 
     @staticmethod
     def _get_exponential_shape(theta: np.ndarray, **kwargs) -> np.ndarray:
+        """Compute the exponential pulse shape.
+
+        Parameters
+        ----------
+        theta : np.ndarray
+            Array of theta values.
+        **kwargs
+            Additional keyword arguments.
+
+        Returns
+        -------
+        np.ndarray
+            Array representing the exponential pulse shape.
+        """
         return np.exp(theta) * np.heaviside(-1.0 * theta, 1)
 
     @staticmethod
     def _get_lorentz_shape(theta: np.ndarray, **kwargs) -> np.ndarray:
+        """Compute the Lorentzian pulse shape.
+
+        Parameters
+        ----------
+        theta : np.ndarray
+            Array of theta values.
+        **kwargs
+            Additional keyword arguments.
+
+        Returns
+        -------
+        np.ndarray
+            Array representing the Lorentzian pulse shape.
+        """
         return 1 / (np.pi * (1 + theta**2))
 
     @staticmethod
     def _get_double_exponential_shape(theta: np.ndarray, **kwargs) -> np.ndarray:
+        """Compute the double-exponential pulse shape.
+
+        Parameters
+        ----------
+        theta : np.ndarray
+            Array of theta values.
+        **kwargs
+            Additional keyword arguments.
+            lam : float
+                Asymmetry parameter controlling the shape.
+
+        Returns
+        -------
+        np.ndarray
+            Array representing the double-exponential pulse shape.
+        """
         lam = kwargs["lam"]
         assert (lam > 0.0) & (lam < 1.0)
         kern = np.zeros(shape=np.shape(theta))
@@ -57,10 +140,38 @@ def _get_double_exponential_shape(theta: np.ndarray, **kwargs) -> np.ndarray:
 
     @staticmethod
     def _get_gaussian_shape(theta: np.ndarray, **kwargs) -> np.ndarray:
+        """Compute the Gaussian pulse shape.
+
+        Parameters
+        ----------
+        theta : np.ndarray
+            Array of theta values.
+        **kwargs
+            Additional keyword arguments.
+
+        Returns
+        -------
+        np.ndarray
+            Array representing the Gaussian pulse shape.
+        """
         return 1 / np.sqrt(np.pi) * np.exp(-(theta**2))
 
     @staticmethod
     def _get_secant_shape(theta: np.ndarray, **kwargs) -> np.ndarray:
+        """Compute the secant pulse shape.
+
+        Parameters
+        ----------
+        theta : np.ndarray
+            Array of theta values.
+        **kwargs
+            Additional keyword arguments.
+
+        Returns
+        -------
+        np.ndarray
+            Array representing the secant pulse shape.
+        """
         return 2 / np.pi / (np.exp(theta) + np.exp(-theta))
 
     __GENERATORS = {
diff --git a/blobmodel/stochasticality.py b/blobmodel/stochasticality.py
index ff37cfa..fb23fab 100644
--- a/blobmodel/stochasticality.py
+++ b/blobmodel/stochasticality.py
@@ -1,3 +1,5 @@
+"""This module defines a class for generating blob parameters."""
+
 from abc import ABC, abstractmethod
 import numpy as np
 from nptyping import NDArray, Float
@@ -54,19 +56,58 @@ def __init__(
         shape_param_y_parameter: float = 0.5,
         blob_alignment: bool = True,
     ) -> None:
-        """The following distributions are implemented:
-
-        exp: exponential distribution with mean 1
-        gamma: gamma distribution with `free_parameter` as shape parameter and mean 1
-        normal: normal distribution with zero mean and `free_parameter` as scale parameter
-        uniform: uniorm distribution with mean 1 and `free_parameter` as width
-        ray: rayleight distribution with mean 1
-        deg: degenerate distribution at `free_parameter`
-        zeros: array of zeros
-
-        blob_alignment: bool = True, optional, If True, the blobs are aligned with their velocity.
-            If blob_aligment == True, the blob shapes are rotated in the propagation direction of the blob.
-            If blob_aligment == False, the blob shapes are independent of the propagation direction.
+
+        """
+        Default implementation of BlobFactory.
+
+        Generates blob parameters for different possible random distributions.
+        All random variables are independent from each other.
+
+        Parameters
+        ----------
+        A_dist : str, optional
+            Distribution type for amplitude, by default "exp"
+        wx_dist : str, optional
+            Distribution type for width in the x-direction, by default "deg"
+        wy_dist : str, optional
+            Distribution type for width in the y-direction, by default "deg"
+        vx_dist : str, optional
+            Distribution type for velocity in the x-direction, by default "deg"
+        vy_dist : str, optional
+            Distribution type for velocity in the y-direction, by default "normal"
+        spx_dist : str, optional
+            Distribution type for shape parameter in the x-direction, by default "deg"
+        spy_dist : str, optional
+            Distribution type for shape parameter in the y-direction, by default "deg"
+        A_parameter : float, optional
+            Free parameter for the amplitude distribution, by default 1.0
+        wx_parameter : float, optional
+            Free parameter for the width distribution in the x-direction, by default 1.0
+        wy_parameter : float, optional
+            Free parameter for the width distribution in the y-direction, by default 1.0
+        vx_parameter : float, optional
+            Free parameter for the velocity distribution in the x-direction, by default 1.0
+        vy_parameter : float, optional
+            Free parameter for the velocity distribution in the y-direction, by default 1.0
+        shape_param_x_parameter : float, optional
+            Free parameter for the shape parameter distribution in the x-direction, by default 0.5
+        shape_param_y_parameter : float, optional
+            Free parameter for the shape parameter distribution in the y-direction, by default 0.5
+        blob_alignment : bool, optional
+            If blob_aligment == True, the blob shapes are rotated in the propagation direction of the blob
+            If blob_aligment == False, the blob shapes are independent of the propagation direction
+
+        Notes
+        -----
+        - The following distributions are implemented:
+            - exp: exponential distribution with mean 1
+            - gamma: gamma distribution with `free_parameter` as shape parameter and mean 1
+            - normal: normal distribution with zero mean and `free_parameter` as scale parameter
+            - uniform: uniorm distribution with mean 1 and `free_parameter` as width
+            - ray: rayleight distribution with mean 1
+            - deg: degenerate distribution at `free_parameter`
+            - zeros: array of zeros
+
         """
         self.amplitude_dist = A_dist
         self.width_x_dist = wx_dist
@@ -90,7 +131,23 @@ def _draw_random_variables(
         free_parameter: float,
         num_blobs: int,
     ) -> NDArray[Any, Float[64]]:
+        """
+        Draws random variables from a specified distribution.
 
+        Parameters
+        ----------
+        dist_type : str
+            Type of distribution.
+        free_parameter : float
+            Free parameter for the distribution.
+        num_blobs : int
+            Number of random variables to draw.
+
+        Returns
+        -------
+        NDArray[Any, Float[64]]
+            Array of random variables drawn from the specified distribution.
+        """
         if dist_type == "exp":
             return np.random.exponential(scale=1, size=num_blobs)
         elif dist_type == "gamma":
@@ -122,6 +179,27 @@ def sample_blobs(
         blob_shape: AbstractBlobShape,
         t_drain: float,
     ) -> List[Blob]:
+        """
+        Creates a list of Blobs used in the Model.
+
+        Parameters
+        ----------
+        Ly : float
+            Size of the domain in the y-direction.
+        T : float
+            Total time duration.
+        num_blobs : int
+            Number of blobs to generate.
+        blob_shape : AbstractBlobShape
+            Object representing the shape of the blobs.
+        t_drain : float
+            Time at which the blobs start draining.
+
+        Returns
+        -------
+        List[Blob]
+            List of Blob objects generated for the Model.
+        """
         amps = self._draw_random_variables(
             dist_type=self.amplitude_dist,
             free_parameter=self.amplitude_parameter,
@@ -159,8 +237,8 @@ def sample_blobs(
                 amplitude=amps[i],
                 width_prop=wxs[i],
                 width_perp=wys[i],
-                velocity_x=vxs[i],
-                velocity_y=vys[i],
+                v_x=vxs[i],
+                v_y=vys[i],
                 pos_x=posxs[i],
                 pos_y=posys[i],
                 t_init=t_inits[i],
@@ -176,5 +254,18 @@ def sample_blobs(
         return np.array(blobs)[np.argsort(amps)]
 
     def is_one_dimensional(self) -> bool:
-        # Perpendicular width parameters are irrelevant since perp shape should be ignored by the Bolb class.
+        """
+        Returns True if the BlobFactory is compatible with a one-dimensional model.
+
+        Returns
+        -------
+        bool
+            True if the BlobFactory is compatible with a one-dimensional model,
+            False otherwise.
+
+        Notes
+        -----
+        - Perpendicular width parameters are irrelevant since perp shape should be ignored by the Bolb class.
+
+        """
         return self.velocity_y_dist == "zeros"
diff --git a/examples/custom_blobfactory.py b/examples/custom_blobfactory.py
index de457b7..d75edfa 100644
--- a/examples/custom_blobfactory.py
+++ b/examples/custom_blobfactory.py
@@ -42,8 +42,8 @@ def sample_blobs(
                 amplitude=amp[i],
                 width_prop=width[i],
                 width_perp=width[i],
-                velocity_x=vx[i],
-                velocity_y=vy[i],
+                v_x=vx[i],
+                v_y=vy[i],
                 pos_x=posx[i],
                 pos_y=posy[i],
                 t_init=t_init[i],
diff --git a/examples/single_blob.py b/examples/single_blob.py
index e6a5e63..421393b 100644
--- a/examples/single_blob.py
+++ b/examples/single_blob.py
@@ -32,8 +32,8 @@ def sample_blobs(
                 amplitude=amp[i],
                 width_prop=width[i],
                 width_perp=width[i],
-                velocity_x=vx[i],
-                velocity_y=vy[i],
+                v_x=vx[i],
+                v_y=vy[i],
                 pos_x=posx[i],
                 pos_y=posy[i],
                 t_init=t_init[i],
diff --git a/tests/test_blob.py b/tests/test_blob.py
index 594658f..8741290 100644
--- a/tests/test_blob.py
+++ b/tests/test_blob.py
@@ -8,8 +8,8 @@
     amplitude=1,
     width_prop=1,
     width_perp=1,
-    velocity_x=1,
-    velocity_y=5,
+    v_x=1,
+    v_y=5,
     pos_x=0,
     pos_y=0,
     t_init=0,
@@ -38,8 +38,8 @@ def test_blob_non_alignment():
         amplitude=1,
         width_prop=1,
         width_perp=1,
-        velocity_x=1,
-        velocity_y=1,
+        v_x=1,
+        v_y=1,
         pos_x=0,
         pos_y=0,
         t_init=0,
@@ -81,8 +81,8 @@ def test_single_point():
         amplitude=1,
         width_prop=1,
         width_perp=1,
-        velocity_x=1,
-        velocity_y=1,
+        v_x=1,
+        v_y=1,
         pos_x=0,
         pos_y=6,
         t_init=0,
@@ -111,8 +111,8 @@ def test_negative_radial_velocity():
         amplitude=1,
         width_prop=1,
         width_perp=1,
-        velocity_x=vx,
-        velocity_y=1,
+        v_x=vx,
+        v_y=1,
         pos_x=0,
         pos_y=6,
         t_init=0,
@@ -144,8 +144,8 @@ def test_theta_0():
         amplitude=1,
         width_prop=1,
         width_perp=1,
-        velocity_x=1,
-        velocity_y=0,
+        v_x=1,
+        v_y=0,
         pos_x=0,
         pos_y=6,
         t_init=0,
@@ -182,8 +182,8 @@ def test_kwargs():
         amplitude=1,
         width_prop=1,
         width_perp=1,
-        velocity_x=1,
-        velocity_y=0,
+        v_x=1,
+        v_y=0,
         pos_x=0,
         pos_y=0,
         t_init=0,
diff --git a/tests/test_drain.py b/tests/test_drain.py
index bb5323c..9d25e1c 100644
--- a/tests/test_drain.py
+++ b/tests/test_drain.py
@@ -9,8 +9,8 @@ def test_high_t_drain():
         amplitude=1,
         width_prop=1,
         width_perp=1,
-        velocity_x=1,
-        velocity_y=1,
+        v_x=1,
+        v_y=1,
         pos_x=0,
         pos_y=6,
         t_init=0,
diff --git a/tests/test_labels.py b/tests/test_labels.py
index 90444bc..3e340c5 100644
--- a/tests/test_labels.py
+++ b/tests/test_labels.py
@@ -34,8 +34,8 @@ def sample_blobs(
                 amplitude=_amp[i],
                 width_prop=_width[i],
                 width_perp=_width[i],
-                velocity_x=_vx[i],
-                velocity_y=_vy[i],
+                v_x=_vx[i],
+                v_y=_vy[i],
                 pos_x=_posx[i],
                 pos_y=_posy[i],
                 t_init=_t_init[i],
diff --git a/tests/test_vx_or_vy=0.py b/tests/test_vx_or_vy=0.py
index 4a189f7..12c9d6c 100644
--- a/tests/test_vx_or_vy=0.py
+++ b/tests/test_vx_or_vy=0.py
@@ -33,8 +33,8 @@ def sample_blobs(
                 amplitude=_amp[i],
                 width_prop=_width[i],
                 width_perp=_width[i],
-                velocity_x=_vx[i],
-                velocity_y=_vy[i],
+                v_x=_vx[i],
+                v_y=_vy[i],
                 pos_x=_posx[i],
                 pos_y=_posy[i],
                 t_init=_t_init[i],
@@ -77,8 +77,8 @@ def sample_blobs(
                 amplitude=_amp[i],
                 width_prop=_width[i],
                 width_perp=_width[i],
-                velocity_x=_vx[i],
-                velocity_y=_vy[i],
+                v_x=_vx[i],
+                v_y=_vy[i],
                 pos_x=_posx[i],
                 pos_y=_posy[i],
                 t_init=_t_init[i],