diff --git a/.gitignore b/.gitignore
index 7fc851d..4c915d1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -104,6 +104,3 @@ ENV/
 # IDE settings
 .vscode/
 .idea/
-
-# OS stuff
-.DS_Store
diff --git a/docs/conf.py b/docs/conf.py
index 63b769a..935162e 100755
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -38,6 +38,7 @@
     'sphinx.ext.napoleon',
     'sphinx.ext.autosummary',
     'sphinx.ext.mathjax',
+    'sphinx.ext.intersphinx',
     'autodocsumm']
 
 # show tocs for classes and functions of modules using the autodocsumm
@@ -99,6 +100,14 @@
 # default language for highlighting in source code
 highlight_language = "python3"
 
+# intersphinx mapping
+intersphinx_mapping = {
+    'numpy': ('https://numpy.org/doc/stable/', None),
+    'scipy': ('https://docs.scipy.org/doc/scipy/', None),
+    'matplotlib': ('https://matplotlib.org/stable/', None),
+    'pyfar': ('https://pyfar.readthedocs.io/en/stable/', None)
+    }
+
 # -- Options for HTML output -------------------------------------------
 
 # The theme to use for HTML and HTML Help pages.  See the documentation for
diff --git a/docs/modules.rst b/docs/modules.rst
index 52a15ec..657f0bd 100644
--- a/docs/modules.rst
+++ b/docs/modules.rst
@@ -7,4 +7,5 @@ according to their modules.
 .. toctree::
    :maxdepth: 1
 
-
+   modules/imkar.scattering
+   modules/imkar.utils
diff --git a/docs/modules/imkar.scattering.rst b/docs/modules/imkar.scattering.rst
new file mode 100644
index 0000000..69bb608
--- /dev/null
+++ b/docs/modules/imkar.scattering.rst
@@ -0,0 +1,7 @@
+imkar.scattering
+================
+
+.. automodule:: imkar.scattering
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/docs/modules/imkar.utils.rst b/docs/modules/imkar.utils.rst
new file mode 100644
index 0000000..89d6335
--- /dev/null
+++ b/docs/modules/imkar.utils.rst
@@ -0,0 +1,7 @@
+imkar.utils
+===========
+
+.. automodule:: imkar.utils
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/imkar/__init__.py b/imkar/__init__.py
index 2d1e55d..d221da6 100644
--- a/imkar/__init__.py
+++ b/imkar/__init__.py
@@ -3,3 +3,7 @@
 __author__ = """The pyfar developers"""
 __email__ = 'info@pyfar.org'
 __version__ = '0.1.0'
+
+
+from . import scattering  # noqa
+from . import utils  # noqa
diff --git a/imkar/scattering/__init__.py b/imkar/scattering/__init__.py
new file mode 100644
index 0000000..b59d748
--- /dev/null
+++ b/imkar/scattering/__init__.py
@@ -0,0 +1,14 @@
+"""
+This module collects the functionality around sound scattering, such as,
+(random incident) scattering coefficients, and analytical solutions.
+"""
+
+from .scattering import (
+    freefield,
+    random,
+)
+
+__all__ = [
+    'freefield',
+    'random',
+    ]
diff --git a/imkar/scattering/scattering.py b/imkar/scattering/scattering.py
new file mode 100644
index 0000000..31d7867
--- /dev/null
+++ b/imkar/scattering/scattering.py
@@ -0,0 +1,155 @@
+import numpy as np
+import pyfar as pf
+from imkar import utils
+
+
+def freefield(sample_pressure, reference_pressure, microphone_weights):
+    r"""
+    Calculate the direction dependent free-field scattering coefficient.
+
+    Uses the Mommertz correlation method [1]_ to calculate the scattering
+    coefficient of the input data:
+
+    .. math::
+        s = 1 -
+            \frac{|\sum_w \underline{p}_{\text{sample}}(\vartheta,\varphi)
+            \cdot \underline{p}_{\text{reference}}^*(\vartheta,\varphi)
+            \cdot w(\vartheta,\varphi)|^2}
+            {\sum_w |\underline{p}_{\text{sample}}(\vartheta,\varphi)|^2
+            \cdot w(\vartheta,\varphi) \cdot \sum_w
+            |\underline{p}_{\text{reference}}(\vartheta,\varphi)|^2
+            \cdot w(\vartheta,\varphi) }
+
+    with the reflected sound pressure of the the sample under investigation
+    :math:`\underline{p}_{\text{sample}}`, the reflected sound pressure from
+    the reference sample (same dimension as the sample under investigation,
+    but with flat surface) :math:`\underline{p}_{\text{reference}}`, the
+    area weights of the sampling :math:`w`, and :math:`\vartheta` and
+    :math:`\varphi` are the ``colatitude``
+    angle and ``azimuth`` angles from the
+    :py:class:`~pyfar.classes.coordinates.Coordinates` object.
+    In other words, the test sample lies in the x-y-plane. See
+    :py:func:`random` to calculate the random incidence
+    scattering coefficient.
+
+    Parameters
+    ----------
+    sample_pressure : :py:class:`~pyfar.classes.audio.FrequencyData`
+        Reflected sound pressure or directivity of the test sample. Its cshape
+        needs to be (..., microphone_weights.size).
+    reference_pressure : :py:class:`~pyfar.classes.audio.FrequencyData`
+        Reflected sound pressure or directivity of the
+        reference sample. Needs to have the same cshape and frequencies as
+        `sample_pressure`.
+    microphone_weights : numpy.ndarray
+        Array containing the area weights for the microphone positions,
+        no normalization required.
+        Its shape needs to match the last dimension in the cshape of
+        `sample_pressure` and `reference_pressure`.
+
+    Returns
+    -------
+    scattering_coefficients : :py:class:`~pyfar.classes.audio.FrequencyData`
+        The scattering coefficient for each incident direction depending on
+        frequency.
+
+
+    References
+    ----------
+    .. [1]  E. Mommertz, „Determination of scattering coefficients from the
+            reflection directivity of architectural surfaces“, Applied
+            Acoustics, Bd. 60, Nr. 2, S. 201-203, June 2000,
+            doi: 10.1016/S0003-682X(99)00057-2.
+
+    """
+    # check inputs
+    if not isinstance(sample_pressure, pf.FrequencyData):
+        raise ValueError(
+            "sample_pressure has to be a pyfar.FrequencyData object")
+    if not isinstance(reference_pressure, pf.FrequencyData):
+        raise ValueError(
+            "reference_pressure has to be a pyfar.FrequencyData object")
+    microphone_weights = np.atleast_1d(
+        np.asarray(microphone_weights, dtype=float))
+    if sample_pressure.cshape != reference_pressure.cshape:
+        raise ValueError(
+            "sample_pressure and reference_pressure have to have the "
+            "same cshape.")
+    if microphone_weights.shape[0] != sample_pressure.cshape[-1]:
+        raise ValueError(
+            "the last dimension of sample_pressure needs be same as the "
+            "microphone_weights.shape.")
+    if not np.allclose(
+            sample_pressure.frequencies, reference_pressure.frequencies):
+        raise ValueError(
+            "sample_pressure and reference_pressure have to have the "
+            "same frequencies.")
+
+    # calculate according to mommertz correlation method Equation (5)
+    p_sample = np.moveaxis(sample_pressure.freq, -1, 0)
+    p_reference = np.moveaxis(reference_pressure.freq, -1, 0)
+    p_sample_sq = np.abs(p_sample)**2
+    p_reference_sq = np.abs(p_reference)**2
+    p_cross = p_sample * np.conj(p_reference)
+
+    p_sample_sum = np.sum(microphone_weights * p_sample_sq, axis=-1)
+    p_ref_sum = np.sum(microphone_weights * p_reference_sq, axis=-1)
+    p_cross_sum = np.sum(microphone_weights * p_cross, axis=-1)
+
+    data_scattering_coefficient \
+        = 1 - ((np.abs(p_cross_sum)**2)/(p_sample_sum*p_ref_sum))
+
+    scattering_coefficients = pf.FrequencyData(
+        np.moveaxis(data_scattering_coefficient, 0, -1),
+        sample_pressure.frequencies)
+
+    return scattering_coefficients
+
+
+def random(
+        scattering_coefficients, incident_directions):
+    r"""
+    Calculate the random-incidence scattering coefficient from free-field
+    data for several incident directions.
+
+    Uses the Paris formula [2]_.
+
+    .. math::
+        s_{rand} = \sum s(\vartheta,\varphi) \cdot cos(\vartheta) \cdot
+        w(\vartheta,\varphi)
+
+    with the scattering coefficients :math:`s(\vartheta,\varphi)`, the area
+    weights ``w`` taken from the `incident_directions.weights`,
+    and :math:`\vartheta` and :math:`\varphi` are the ``colatitude``
+    angle and ``azimuth`` angles from the
+    :py:class:`~pyfar.classes.coordinates.Coordinates` object.
+    Note that the incident directions should be
+    equally distributed to get a valid result. See
+    :py:func:`freefield` to calculate the free-field scattering coefficient.
+
+    Parameters
+    ----------
+    scattering_coefficients : :py:class:`~pyfar.classes.audio.FrequencyData`
+        Scattering coefficients for different incident directions. Its cshape
+        needs to be (..., incident_directions.csize)
+    incident_directions : :py:class:`~pyfar.classes.coordinates.Coordinates`
+        Defines the incidence directions of each `scattering_coefficients`
+        in a :py:class:`~pyfar.classes.coordinates.Coordinates` object.
+        Its cshape needs to match
+        the last dimension of `scattering_coefficients`.
+        Points contained in `incident_directions` must have the same radii.
+        The weights need to reflect the area `incident_directions.weights`.
+
+    Returns
+    -------
+    random_scattering : :py:class:`~pyfar.classes.audio.FrequencyData`
+        The random-incidence scattering coefficient depending on frequency.
+
+    References
+    ----------
+    .. [2]  H. Kuttruff, Room acoustics, Sixth edition. Boca Raton:
+            CRC Press/Taylor & Francis Group, 2017.
+    """
+    random_scattering = utils.paris_formula(
+        scattering_coefficients, incident_directions)
+    return random_scattering
diff --git a/imkar/utils/__init__.py b/imkar/utils/__init__.py
new file mode 100644
index 0000000..eaf3a04
--- /dev/null
+++ b/imkar/utils/__init__.py
@@ -0,0 +1,7 @@
+from .utils import (
+    paris_formula,
+)
+
+__all__ = [
+    'paris_formula',
+    ]
diff --git a/imkar/utils/utils.py b/imkar/utils/utils.py
new file mode 100644
index 0000000..ac4e05a
--- /dev/null
+++ b/imkar/utils/utils.py
@@ -0,0 +1,65 @@
+import numpy as np
+import pyfar as pf
+
+
+def paris_formula(coefficients, incident_directions):
+    r"""
+    Calculate the random-incidence coefficient from free-field
+    data for several incident directions.
+
+    Uses the Paris formula [2]_.
+
+    .. math::
+        c_{rand} = \sum c(\vartheta,\varphi) \cdot cos(\vartheta) \cdot
+        w(\vartheta,\varphi)
+
+    with the coefficients :math:`c(\vartheta,\varphi)`, the area
+    weights ``w`` taken from the `incident_directions.weights`,
+    and :math:`\vartheta` and :math:`\varphi` are the ``colatitude``
+    angle and ``azimuth`` angles from the
+    :py:class:`~pyfar.classes.coordinates.Coordinates` object.
+    Note that the incident directions should be
+    equally distributed to get a valid result.
+
+    Parameters
+    ----------
+    coefficients : :py:class:`~pyfar.classes.audio.FrequencyData`
+        Scattering coefficients for different incident directions. Its cshape
+        needs to be (..., `incident_directions.csize`)
+    incident_directions : :py:class:`~pyfar.classes.coordinates.Coordinates`
+        Defines the incidence directions of each `scattering_coefficients`
+        in a :py:class:`~pyfar.classes.coordinates.Coordinates` object.
+        Its cshape needs to match
+        the last dimension of `scattering_coefficients`.
+        Points contained in `incident_directions` must have the same radii.
+        The weights need to reflect the area ``incident_directions.weights``.
+
+    Returns
+    -------
+    random_coefficient : :py:class:`~pyfar.classes.audio.FrequencyData`
+        The random-incidence coefficient depending on frequency.
+
+    References
+    ----------
+    .. [2]  H. Kuttruff, Room acoustics, Sixth edition. Boca Raton:
+            CRC Press/Taylor & Francis Group, 2017.
+    """
+    if not isinstance(coefficients, pf.FrequencyData):
+        raise ValueError("coefficients has to be FrequencyData")
+    if not isinstance(incident_directions, pf.Coordinates):
+        raise ValueError("incident_directions have to be None or Coordinates")
+    if incident_directions.cshape[0] != coefficients.cshape[-1]:
+        raise ValueError(
+            "the last dimension of coefficients needs be same as "
+            "the incident_directions.cshape.")
+
+    theta = incident_directions.get_sph().T[1]
+    weight = np.cos(theta) * incident_directions.weights
+    norm = np.sum(weight)
+    coefficients_freq = np.swapaxes(coefficients.freq, -1, -2)
+    random_coefficient = pf.FrequencyData(
+        np.sum(coefficients_freq*weight/norm, axis=-1),
+        coefficients.frequencies,
+        comment='random-incidence coefficient'
+    )
+    return random_coefficient
diff --git a/setup.py b/setup.py
index 87153a2..705229b 100644
--- a/setup.py
+++ b/setup.py
@@ -10,7 +10,10 @@
 with open('HISTORY.rst') as history_file:
     history = history_file.read()
 
-requirements = [ ]
+requirements = [
+    'numpy>=1.23.0',
+    'pyfar',
+]
 
 test_requirements = ['pytest>=3', ]
 
diff --git a/tests/conftest.py b/tests/conftest.py
new file mode 100644
index 0000000..a98287c
--- /dev/null
+++ b/tests/conftest.py
@@ -0,0 +1,79 @@
+import pytest
+import pyfar as pf
+import numpy as np
+
+
+@pytest.fixture
+def half_sphere():
+    """return 42th order gaussian sampling for the half sphere and radius 1.
+
+    Returns
+    -------
+    pf.Coordinates
+        half sphere sampling grid
+    """
+    mics = pf.samplings.sph_gaussian(42)
+    # delete lower part of sphere
+    return mics[mics.get_sph().T[1] <= np.pi/2]
+
+
+@pytest.fixture
+def quarter_half_sphere():
+    """return 10th order gaussian sampling for the quarter half sphere
+    and radius 1.
+
+    Returns
+    -------
+    pf.Coordinates
+        quarter half sphere sampling grid
+    """
+    incident_directions = pf.samplings.sph_gaussian(10)
+    incident_directions = incident_directions[
+        incident_directions.get_sph().T[1] <= np.pi/2]
+    return incident_directions[
+        incident_directions.get_sph().T[0] <= np.pi/2]
+
+
+@pytest.fixture
+def pressure_data_mics(half_sphere):
+    """returns a sound pressure data example, with sound pressure 0 and
+    two frequency bins
+
+    Parameters
+    ----------
+    half_sphere : pf.Coordinates
+        half sphere sampling grid for mics
+
+    Returns
+    -------
+    pyfar.FrequencyData
+        output sound pressure data
+    """
+    frequencies = [200, 300]
+    shape_new = np.append(half_sphere.cshape, len(frequencies))
+    return pf.FrequencyData(np.zeros(shape_new), frequencies)
+
+
+@pytest.fixture
+def pressure_data_mics_incident_directions(
+        half_sphere, quarter_half_sphere):
+    """returns a sound pressure data example, with sound pressure 0 and
+    two frequency bins
+
+    Parameters
+    ----------
+    half_sphere : pf.Coordinates
+        half sphere sampling grid for mics
+    quarter_half_sphere : pf.Coordinates
+        quarter half sphere sampling grid for incident directions
+
+    Returns
+    -------
+    pyfar.FrequencyData
+        output sound pressure data
+    """
+    frequencies = [200, 300]
+    shape_new = np.append(
+        quarter_half_sphere.cshape, half_sphere.cshape)
+    shape_new = np.append(shape_new, len(frequencies))
+    return pf.FrequencyData(np.zeros(shape_new), frequencies)
diff --git a/tests/test_imkar.py b/tests/test_imkar.py
index e6d07f7..f801097 100644
--- a/tests/test_imkar.py
+++ b/tests/test_imkar.py
@@ -1,5 +1,24 @@
-def test_import_imkar():
-    try:
-        import imkar           # noqa
-    except ImportError:
-        assert False
+#!/usr/bin/env python
+
+"""Tests for `imkar` package."""
+
+import pytest
+
+
+from imkar import imkar
+
+
+@pytest.fixture
+def response():
+    """Sample pytest fixture.
+
+    See more at: http://doc.pytest.org/en/latest/fixture.html
+    """
+    # import requests
+    # return requests.get('https://github.com/mberz/cookiecutter-pypackage')
+
+
+def test_content(response):
+    """Sample pytest test function with the pytest fixture as an argument."""
+    # from bs4 import BeautifulSoup
+    # assert 'GitHub' in BeautifulSoup(response.content).title.string
diff --git a/tests/test_scattering.py b/tests/test_scattering.py
new file mode 100644
index 0000000..3120d22
--- /dev/null
+++ b/tests/test_scattering.py
@@ -0,0 +1,118 @@
+import pytest
+import numpy as np
+
+from imkar import scattering
+
+
+def test_freefield_1(
+        half_sphere, pressure_data_mics):
+    mics = half_sphere
+    p_sample = pressure_data_mics.copy()
+    p_sample.freq.fill(1)
+    p_reference = pressure_data_mics.copy()
+    p_sample.freq[5, :] = 0
+    p_reference.freq[5, :] = np.sum(p_sample.freq.flatten())/2
+    s = scattering.freefield(p_sample, p_reference, mics.weights)
+    np.testing.assert_allclose(s.freq, 1)
+
+
+def test_freefield_wrong_input(
+        half_sphere, pressure_data_mics):
+    mics = half_sphere
+    p_sample = pressure_data_mics.copy()
+    p_reference = pressure_data_mics.copy()
+
+    with pytest.raises(ValueError, match='sample_pressure'):
+        scattering.freefield(1, p_reference, mics.weights)
+    with pytest.raises(ValueError, match='reference_pressure'):
+        scattering.freefield(p_sample, 1, mics.weights)
+    with pytest.raises(ValueError, match='microphone_weights'):
+        scattering.freefield(p_sample, p_reference, 1)
+    with pytest.raises(ValueError, match='cshape'):
+        scattering.freefield(
+            p_sample[:-2, ...], p_reference, mics.weights)
+    with pytest.raises(ValueError, match='microphone_weights'):
+        scattering.freefield(
+            p_sample, p_reference, mics.weights[:10])
+    with pytest.raises(ValueError, match='same frequencies'):
+        p_sample.frequencies[0] = 1
+        scattering.freefield(p_sample, p_reference, mics.weights)
+
+
+def test_freefield_05(
+        half_sphere, pressure_data_mics):
+    mics = half_sphere
+    p_sample = pressure_data_mics.copy()
+    p_reference = pressure_data_mics.copy()
+    p_sample.freq[7, :] = 1
+    p_sample.freq[28, :] = 1
+    p_reference.freq[7, :] = 1
+    s = scattering.freefield(p_sample, p_reference, mics.weights)
+    np.testing.assert_allclose(s.freq, 0.5)
+
+
+def test_freefield_0(
+        half_sphere, pressure_data_mics):
+    mics = half_sphere
+    p_sample = pressure_data_mics.copy()
+    p_reference = pressure_data_mics.copy()
+    p_reference.freq[5, :] = 1
+    p_sample.freq[5, :] = 1
+    s = scattering.freefield(p_sample, p_reference, mics.weights)
+    np.testing.assert_allclose(s.freq, 0)
+    assert s.freq.shape[-1] == p_sample.n_bins
+
+
+def test_freefield_0_with_incident(
+        half_sphere, quarter_half_sphere,
+        pressure_data_mics_incident_directions):
+    mics = half_sphere
+    incident_directions = quarter_half_sphere
+    p_sample = pressure_data_mics_incident_directions.copy()
+    p_reference = pressure_data_mics_incident_directions.copy()
+    p_reference.freq[:, 2, :] = 1
+    p_sample.freq[:, 2, :] = 1
+    s = scattering.freefield(p_sample, p_reference, mics.weights)
+    s_rand = scattering.random(s, incident_directions)
+    np.testing.assert_allclose(s.freq, 0)
+    np.testing.assert_allclose(s_rand.freq, 0)
+    assert s.freq.shape[-1] == p_sample.n_bins
+    assert s.cshape == incident_directions.cshape
+    assert s_rand.freq.shape[-1] == p_sample.n_bins
+
+
+def test_freefield_1_with_incidence(
+        half_sphere, quarter_half_sphere,
+        pressure_data_mics_incident_directions):
+    mics = half_sphere
+    incident_directions = quarter_half_sphere
+    p_sample = pressure_data_mics_incident_directions.copy()
+    p_reference = pressure_data_mics_incident_directions.copy()
+    p_reference.freq[:, 2, :] = 1
+    p_sample.freq[:, 3, :] = 1
+    s = scattering.freefield(p_sample, p_reference, mics.weights)
+    s_rand = scattering.random(s, incident_directions)
+    np.testing.assert_allclose(s.freq, 1)
+    np.testing.assert_allclose(s_rand.freq, 1)
+    assert s.freq.shape[-1] == p_sample.n_bins
+    assert s.cshape == incident_directions.cshape
+    assert s_rand.freq.shape[-1] == p_sample.n_bins
+
+
+def test_freefield_05_with_inci(
+        half_sphere, quarter_half_sphere,
+        pressure_data_mics_incident_directions):
+    mics = half_sphere
+    incident_directions = quarter_half_sphere
+    p_sample = pressure_data_mics_incident_directions.copy()
+    p_reference = pressure_data_mics_incident_directions.copy()
+    p_sample.freq[:, 7, :] = 1
+    p_sample.freq[:, 28, :] = 1
+    p_reference.freq[:, 7, :] = 1
+    s = scattering.freefield(p_sample, p_reference, mics.weights)
+    s_rand = scattering.random(s, incident_directions)
+    np.testing.assert_allclose(s.freq, 0.5)
+    np.testing.assert_allclose(s_rand.freq, 0.5)
+    assert s.freq.shape[-1] == p_sample.n_bins
+    assert s.cshape == incident_directions.cshape
+    assert s_rand.freq.shape[-1] == p_sample.n_bins
diff --git a/tests/test_utils.py b/tests/test_utils.py
new file mode 100644
index 0000000..1440c0b
--- /dev/null
+++ b/tests/test_utils.py
@@ -0,0 +1,34 @@
+import pytest
+import numpy as np
+import pyfar as pf
+
+from imkar import utils
+
+
+@pytest.mark.parametrize("c_value",  [
+    (0), (0.2), (0.5), (0.8), (1)])
+@pytest.mark.parametrize("frequencies",  [
+    ([100, 200]), ([100])])
+def test_random_constant_coefficient(
+        c_value, frequencies, half_sphere):
+    incident_directions = half_sphere
+    shape = np.append(half_sphere.cshape, len(frequencies))
+    coefficient = pf.FrequencyData(np.zeros(shape)+c_value, frequencies)
+    c_rand = utils.paris_formula(coefficient, incident_directions)
+    np.testing.assert_allclose(c_rand.freq, c_value)
+    assert c_rand.comment == 'random-incidence coefficient'
+
+
+def test_random_non_constant_coefficient():
+    data = pf.samplings.sph_gaussian(10)
+    incident_directions = data[data.get_sph().T[1] <= np.pi/2]
+    incident_cshape = incident_directions.cshape
+    s_value = np.arange(
+        incident_cshape[0]).reshape(incident_cshape) / incident_cshape[0]
+    theta = incident_directions.get_sph().T[1]
+    actual_weight = np.cos(theta) * incident_directions.weights
+    actual_weight /= np.sum(actual_weight)
+    coefficient = pf.FrequencyData(s_value.reshape((50, 1)), [100])
+    c_rand = utils.paris_formula(coefficient, incident_directions)
+    desired = np.sum(s_value*actual_weight)
+    np.testing.assert_allclose(c_rand.freq, desired)