diff --git a/.coveragerc b/.coveragerc
new file mode 100644
index 0000000..98111fa
--- /dev/null
+++ b/.coveragerc
@@ -0,0 +1,16 @@
+[run]
+branch = True
+source = wiimatch
+
+[report]
+exclude_lines =
+    if self.debug:
+    pragma: no cover
+    raise NotImplementedError
+    if __name__ == .__main__.:
+ignore_errors = True
+omit =
+    conftest.py
+    setup.py
+    tests/*
+    wiimatch/tests/*
diff --git a/wiimatch/tests/__init__.py b/wiimatch/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/wiimatch/tests/test_match.py b/wiimatch/tests/test_match.py
new file mode 100644
index 0000000..28319fc
--- /dev/null
+++ b/wiimatch/tests/test_match.py
@@ -0,0 +1,78 @@
+"""
+A module containing unit tests for the `wcsutil` module.
+
+Licensed under a 3-clause BSD style license - see LICENSE.rst
+
+"""
+import pytest
+import numpy as np
+
+from wiimatch import match
+
+
+@pytest.mark.parametrize('solver', ['RLU', 'PINV'])
+def test_match_lsq_solver(solver):
+    # simulate background image data:
+    c = [1.32, 0.15, 0.62, 0, 0.74, 0, 0, 0]
+    im1 = np.zeros((5, 5, 4), dtype=np.float)
+    cbg = c[0] * np.ones_like(im1)  # constand background level image
+
+    # add slope:
+    ind = np.indices(im1.shape, dtype=np.float)
+    im3 = cbg + c[1] * ind[0] + c[2] * ind[1] + c[4] * ind[2]
+
+    mask = np.ones_like(im1, dtype=np.int8)
+    sigma = np.ones_like(im1, dtype=np.float)
+
+    p = match.match_lsq(
+        [im1, im3], [mask, mask], [sigma, sigma],
+        degree=(1, 1, 1), center=(0, 0, 0), solver=solver
+    )
+
+    assert np.allclose(-p[0], p[1], rtol=1.e-8, atol=1.e-12)
+    assert np.allclose(c, 2 * np.abs(p[0]), rtol=1.e-8, atol=1.e-12)
+
+
+def test_match_lsq_extended_return():
+    # simulate background image data:
+    c = [1.32, 0.15, 0.62, 0, 0.74, 0, 0, 0]
+    im1 = np.zeros((5, 5, 4), dtype=np.float)
+    cbg = c[0] * np.ones_like(im1)  # constand background level image
+
+    # add slope:
+    ind = np.indices(im1.shape, dtype=np.float)
+    im3 = cbg + c[1] * ind[0] + c[2] * ind[1] + c[4] * ind[2]
+
+    mask = np.ones_like(im1, dtype=np.int8)
+    sigma = np.ones_like(im1, dtype=np.float)
+
+    p, a, b, coord_arrays, eff_center, coord_system = match.match_lsq(
+        [im1, im3], [mask, mask], [sigma, sigma],
+        degree=1, center=(0, 0, 0), ext_return=True
+    )
+
+    assert np.allclose(-p[0], p[1], rtol=1.e-8, atol=1.e-12)
+    assert np.allclose(c, 2 * np.abs(p[0]), rtol=1.e-8, atol=1.e-12)
+
+
+@pytest.mark.parametrize('degree', [1, (1, 1, 1)])
+def test_match_lsq_num_degree(degree):
+    # simulate background image data:
+    c = [1.32, 0.15, 0.62, 0, 0.74, 0, 0, 0]
+    im1 = np.zeros((5, 5, 4), dtype=np.float)
+    cbg = c[0] * np.ones_like(im1)  # constand background level image
+
+    # add slope:
+    ind = np.indices(im1.shape, dtype=np.float)
+    im3 = cbg + c[1] * ind[0] + c[2] * ind[1] + c[4] * ind[2]
+
+    mask = np.ones_like(im1, dtype=np.int8)
+    sigma = np.ones_like(im1, dtype=np.float)
+
+    p = match.match_lsq(
+        [im1, im3], [mask, mask], [sigma, sigma],
+        degree=degree, center=(0, 0, 0), solver='RLU', ext_return=False
+    )
+
+    assert np.allclose(-p[0], p[1], rtol=1.e-8, atol=1.e-12)
+    assert np.allclose(c, 2 * np.abs(p[0]), rtol=1.e-8, atol=1.e-12)
diff --git a/wiimatch/tests/test_utils.py b/wiimatch/tests/test_utils.py
new file mode 100644
index 0000000..9c7d08e
--- /dev/null
+++ b/wiimatch/tests/test_utils.py
@@ -0,0 +1,32 @@
+"""
+A module containing unit tests for the `wcsutil` module.
+
+Licensed under a 3-clause BSD style license - see LICENSE.rst
+
+"""
+import pytest
+import numpy as np
+
+from wiimatch import utils
+
+
+def test_utils_coordinates():
+    image_shape = (3, 5, 4)
+    center = (0, 0, 0)
+    c = utils.create_coordinate_arrays(image_shape, center=center)
+    ind = np.indices(image_shape, dtype=np.float)[::-1]
+
+    assert np.allclose(c[0], ind, rtol=1.e-8, atol=1.e-12)
+    assert np.allclose(c[1], center, rtol=1.e-8, atol=1.e-12)
+
+
+def test_utils_coordinates_no_center():
+    image_shape = (3, 5, 4)
+    c = utils.create_coordinate_arrays(image_shape, center=None)
+    ind = np.indices(image_shape, dtype=np.float)[::-1]
+
+    center = tuple(i // 2 for i in image_shape)
+    for orig, cc, i in zip(center, c[0], ind):
+        assert(i - orig, cc)
+
+    assert np.allclose(c[1], center, rtol=1.e-8, atol=1.e-12)