diff --git a/poetry.lock b/poetry.lock
index daa2600..c5772bf 100644
--- a/poetry.lock
+++ b/poetry.lock
@@ -1261,6 +1261,33 @@ typing-extensions = {version = ">=4.0.0,<5.0", markers = "python_version < \"3.1
 [package.extras]
 jupyter = ["ipywidgets (>=7.5.1,<9)"]
 
+[[package]]
+name = "ruff"
+version = "0.7.0"
+description = "An extremely fast Python linter and code formatter, written in Rust."
+optional = false
+python-versions = ">=3.7"
+files = [
+    {file = "ruff-0.7.0-py3-none-linux_armv6l.whl", hash = "sha256:0cdf20c2b6ff98e37df47b2b0bd3a34aaa155f59a11182c1303cce79be715628"},
+    {file = "ruff-0.7.0-py3-none-macosx_10_12_x86_64.whl", hash = "sha256:496494d350c7fdeb36ca4ef1c9f21d80d182423718782222c29b3e72b3512737"},
+    {file = "ruff-0.7.0-py3-none-macosx_11_0_arm64.whl", hash = "sha256:214b88498684e20b6b2b8852c01d50f0651f3cc6118dfa113b4def9f14faaf06"},
+    {file = "ruff-0.7.0-py3-none-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:630fce3fefe9844e91ea5bbf7ceadab4f9981f42b704fae011bb8efcaf5d84be"},
+    {file = "ruff-0.7.0-py3-none-manylinux_2_17_armv7l.manylinux2014_armv7l.whl", hash = "sha256:211d877674e9373d4bb0f1c80f97a0201c61bcd1e9d045b6e9726adc42c156aa"},
+    {file = "ruff-0.7.0-py3-none-manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:194d6c46c98c73949a106425ed40a576f52291c12bc21399eb8f13a0f7073495"},
+    {file = "ruff-0.7.0-py3-none-manylinux_2_17_ppc64.manylinux2014_ppc64.whl", hash = "sha256:82c2579b82b9973a110fab281860403b397c08c403de92de19568f32f7178598"},
+    {file = "ruff-0.7.0-py3-none-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:9af971fe85dcd5eaed8f585ddbc6bdbe8c217fb8fcf510ea6bca5bdfff56040e"},
+    {file = "ruff-0.7.0-py3-none-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:b641c7f16939b7d24b7bfc0be4102c56562a18281f84f635604e8a6989948914"},
+    {file = "ruff-0.7.0-py3-none-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:d71672336e46b34e0c90a790afeac8a31954fd42872c1f6adaea1dff76fd44f9"},
+    {file = "ruff-0.7.0-py3-none-musllinux_1_2_aarch64.whl", hash = "sha256:ab7d98c7eed355166f367597e513a6c82408df4181a937628dbec79abb2a1fe4"},
+    {file = "ruff-0.7.0-py3-none-musllinux_1_2_armv7l.whl", hash = "sha256:1eb54986f770f49edb14f71d33312d79e00e629a57387382200b1ef12d6a4ef9"},
+    {file = "ruff-0.7.0-py3-none-musllinux_1_2_i686.whl", hash = "sha256:dc452ba6f2bb9cf8726a84aa877061a2462afe9ae0ea1d411c53d226661c601d"},
+    {file = "ruff-0.7.0-py3-none-musllinux_1_2_x86_64.whl", hash = "sha256:4b406c2dce5be9bad59f2de26139a86017a517e6bcd2688da515481c05a2cb11"},
+    {file = "ruff-0.7.0-py3-none-win32.whl", hash = "sha256:f6c968509f767776f524a8430426539587d5ec5c662f6addb6aa25bc2e8195ec"},
+    {file = "ruff-0.7.0-py3-none-win_amd64.whl", hash = "sha256:ff4aabfbaaba880e85d394603b9e75d32b0693152e16fa659a3064a85df7fce2"},
+    {file = "ruff-0.7.0-py3-none-win_arm64.whl", hash = "sha256:10842f69c245e78d6adec7e1db0a7d9ddc2fff0621d730e61657b64fa36f207e"},
+    {file = "ruff-0.7.0.tar.gz", hash = "sha256:47a86360cf62d9cd53ebfb0b5eb0e882193fc191c6d717e8bef4462bc3b9ea2b"},
+]
+
 [[package]]
 name = "six"
 version = "1.16.0"
@@ -1272,6 +1299,17 @@ files = [
     {file = "six-1.16.0.tar.gz", hash = "sha256:1e61c37477a1626458e36f7b1d82aa5c9b094fa4802892072e49de9c60c4c926"},
 ]
 
+[[package]]
+name = "style"
+version = "1.1.0"
+description = "🌈 Terminal string styling"
+optional = false
+python-versions = "*"
+files = [
+    {file = "style-1.1.0-py2.py3-none-any.whl", hash = "sha256:6485a4bcb84629341a5fd1587fe3ac4887daa4741f0c8a1d01b9c3c8a263afe7"},
+    {file = "style-1.1.0.tar.gz", hash = "sha256:8eb365fc15039b19b728bd4e6e85fb7daf24e7aeeec6a15a666f97484c564005"},
+]
+
 [[package]]
 name = "tenacity"
 version = "9.0.0"
@@ -1320,6 +1358,19 @@ files = [
     {file = "tzdata-2024.2.tar.gz", hash = "sha256:7d85cc416e9382e69095b7bdf4afd9e3880418a2413feec7069d533d6b4e31cc"},
 ]
 
+[[package]]
+name = "update"
+version = "0.0.1"
+description = ""
+optional = false
+python-versions = "*"
+files = [
+    {file = "update-0.0.1-py2.py3-none-any.whl", hash = "sha256:a25522b4bf60e3e3c1a3ff3ca3a4f5a328ac4b8ff400fdc9614483147e313323"},
+]
+
+[package.dependencies]
+style = "1.1.0"
+
 [[package]]
 name = "wcwidth"
 version = "0.2.13"
@@ -1356,4 +1407,4 @@ plotly = ["plotly"]
 [metadata]
 lock-version = "2.0"
 python-versions = ">=3.9,<4.0"
-content-hash = "f59d1f558858104329eb04fe0e9392843c2995c689ef2228a3c169d481f9d3d6"
+content-hash = "0f5462b3336d27bed4f4b51164d8f09a721ffb5c64f660103dd73d6e7bf0038c"
diff --git a/pyproject.toml b/pyproject.toml
index e74cdef..e60e049 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -49,6 +49,7 @@ exclude = ["test/**/*.py", "example/**", "ts-diagrams/**"]
         matplotlib = "^3.4.2"
         prettytable = ">=2.1,<4.0"
         plotly = { version = "^5.13.1", optional = true }
+        update = "^0.0.1"
 
         [tool.poetry.group.dev.dependencies]
         pytest = ">=7.2,<9.0"
@@ -58,6 +59,7 @@ exclude = ["test/**/*.py", "example/**", "ts-diagrams/**"]
         black = ">=22.12,<25.0"
         isort = "^5.11.4"
         pytest-console-scripts = "^1.3.1"
+        ruff = "^0.7.0"
 
 [build-system]
 requires = ["poetry-core>=1.0.0"]
@@ -68,3 +70,17 @@ tanabesugano = "tanabesugano.cmd:cmd_line"
 
 [tool.poetry.extras]
 plotly = ["plotly"]
+
+[tool.ruff]
+lint.select = ["ALL"]
+lint.ignore = ["N806"]
+target-version = "py38"
+src = ["tanabesugano"]
+
+[tool.ruff.lint.isort]
+known-first-party = ["umf"]
+force-single-line = true
+lines-between-types = 1
+lines-after-imports = 2
+known-third-party = ["poetry.core"]
+required-imports = ["from __future__ import annotations"]
diff --git a/tanabesugano/__init__.py b/tanabesugano/__init__.py
index a0e60d3..7a69a1d 100644
--- a/tanabesugano/__init__.py
+++ b/tanabesugano/__init__.py
@@ -1,2 +1,5 @@
 """tanabesugano: A Python package for Tanabe-Sugano diagrams."""
-__version__ = "1.4.3"
+from __future__ import annotations
+
+
+__version__ = "1.4.4"
diff --git a/tanabesugano/batch.py b/tanabesugano/batch.py
index a023b2b..3b3fcb1 100644
--- a/tanabesugano/batch.py
+++ b/tanabesugano/batch.py
@@ -1,8 +1,11 @@
+from __future__ import annotations
+
 from typing import List
 
 import numpy as np
 
-from tanabesugano import matrices, tools
+from tanabesugano import matrices
+from tanabesugano import tools
 
 
 class Batch:
@@ -22,26 +25,23 @@ def __init__(
             C = [3600.0, 4000, 10]
         if len(Dq) != 3:
             raise KeyError(
-                "The range of `Dq` is based on the three values: start, stop, steps!"
+                "The range of `Dq` is based on the three values: start, stop, steps!",
             )
-        else:
-            self.Dq = np.linspace(Dq[0], Dq[1], int(Dq[2]))  # Oh-crystalfield-splitting
+        self.Dq = np.linspace(Dq[0], Dq[1], int(Dq[2]))  # Oh-crystalfield-splitting
         if len(B) != 3:
             raise KeyError(
-                "The range of `B` is based on the three values: start, stop, steps!"
+                "The range of `B` is based on the three values: start, stop, steps!",
             )
-        else:
-            self.B = np.linspace(
-                B[0], B[1], int(B[2])
-            )  # Racah-Parameter B in wavenumbers
+        self.B = np.linspace(
+            B[0], B[1], int(B[2]),
+        )  # Racah-Parameter B in wavenumbers
         if len(C) != 3:
             raise KeyError(
-                "The range of `C` is based on the three values: start, stop, steps!"
+                "The range of `C` is based on the three values: start, stop, steps!",
             )
-        else:
-            self.C = np.linspace(
-                C[0], C[1], int(C[2])
-            )  # Racah-Parameter C in wavenumbers
+        self.C = np.linspace(
+            C[0], C[1], int(C[2]),
+        )  # Racah-Parameter C in wavenumbers
 
         if slater:
             # Transformin Racah to Slater-Condon
@@ -58,14 +58,11 @@ def __init__(
             self._size = 10
         self.result: List[dict] = []
 
-    def calculation(self):
-        """
-        Is filling the self.result with the iTS states of over-iterated energy range
-        """
+    def calculation(self) -> None:
+        """Is filling the self.result with the iTS states of over-iterated energy range."""
         for _Dq in self.Dq:
             for _B in self.B:
                 for _C in self.C:
-
                     if self.d_count == 2:  # d3
                         states = matrices.d2(Dq=_Dq, B=_B, C=_C).solver()
                         self.result.append(
@@ -75,7 +72,7 @@ def calculation(self):
                                 "B": _B,
                                 "C": _C,
                                 "states": states,
-                            }
+                            },
                         )
 
                     elif self.d_count == 3:  # d3
@@ -87,7 +84,7 @@ def calculation(self):
                                 "B": _B,
                                 "C": _C,
                                 "states": states,
-                            }
+                            },
                         )
                     elif self.d_count == 4:  # d4
                         states = matrices.d4(Dq=_Dq, B=_B, C=_C).solver()
@@ -98,7 +95,7 @@ def calculation(self):
                                 "B": _B,
                                 "C": _C,
                                 "states": states,
-                            }
+                            },
                         )
                     elif self.d_count == 5:  # d5
                         states = matrices.d5(Dq=_Dq, B=_B, C=_C).solver()
@@ -109,7 +106,7 @@ def calculation(self):
                                 "B": _B,
                                 "C": _C,
                                 "states": states,
-                            }
+                            },
                         )
                     elif self.d_count == 6:  # d6
                         states = matrices.d6(Dq=_Dq, B=_B, C=_C).solver()
@@ -120,7 +117,7 @@ def calculation(self):
                                 "B": _B,
                                 "C": _C,
                                 "states": states,
-                            }
+                            },
                         )
                     elif self.d_count == 7:  # d7
                         states = matrices.d7(Dq=_Dq, B=_B, C=_C).solver()
@@ -131,7 +128,7 @@ def calculation(self):
                                 "B": _B,
                                 "C": _C,
                                 "states": states,
-                            }
+                            },
                         )
                     elif self.d_count == 8:  # d8
                         states = matrices.d8(Dq=_Dq, B=_B, C=_C).solver()
@@ -142,7 +139,7 @@ def calculation(self):
                                 "B": _B,
                                 "C": _C,
                                 "states": states,
-                            }
+                            },
                         )
                     else:
                         raise ValueError("not a correct value!")
diff --git a/tanabesugano/cmd.py b/tanabesugano/cmd.py
index d5207ea..3f1af2f 100644
--- a/tanabesugano/cmd.py
+++ b/tanabesugano/cmd.py
@@ -1,21 +1,29 @@
 #!/usr/bin/env python
+
+from __future__ import annotations
+
 import argparse
 
+from pathlib import Path
+
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
+
 from prettytable import PrettyTable
-from pathlib import Path
+
 
 try:
     import plotly.express as px
 except ImportError:  # pragma: no cover
     px = None
 
-from tanabesugano import __version__, matrices, tools
+from tanabesugano import __version__
+from tanabesugano import matrices
+from tanabesugano import tools
 
 
-class CMDmain(object):
+class CMDmain:
     def __init__(
         self,
         Dq: float = 4000.0,
@@ -41,13 +49,13 @@ def __init__(
             Electron count, by default 5
         slater : bool, optional
              Transforming from Racah to Slater-Condon, by default False
+
         """
         self.Dq = Dq
         self.B = B
         self.C = C
 
         if slater:
-
             self.B, self.C = tools.racah(B, C)
         self.nroot = nroots
         energy = np.linspace(0.0, self.Dq, nroots)
@@ -62,7 +70,7 @@ def __init__(
         self.result = np.zeros((self._size + 1, nroots))
 
         self.df = pd.DataFrame(
-            {"Energy": energy, "delta_B": energy / self.B, "10Dq": energy * 10.0}
+            {"Energy": energy, "delta_B": energy / self.B, "10Dq": energy * 10.0},
         )
         self.title_TS = (
             f"TS-diagram_d{self.d_count}_10Dq_{int(self.Dq * 10.0)}_"
@@ -74,7 +82,6 @@ def __init__(
         )
 
     def plot(self) -> None:
-
         # Figure one for classical Tanabe-Sugano-Diagram with B-dependency
         plt.figure(1)
 
@@ -86,7 +93,7 @@ def plot(self) -> None:
             self.df.drop(["Energy", "delta_B", "10Dq"], axis=1).to_numpy() / self.B,
             ls="--",
         )
-        self.label_plot("Tanabe-Sugano-Diagram", "$E/B$", "$\Delta/B$")
+        self.label_plot("Tanabe-Sugano-Diagram", "$E/B$", r"$\Delta/B$")
         # Figure one for Energy-Correlation-Diagram Dq-Energy versus State-Energy
         plt.figure(2)
 
@@ -96,7 +103,7 @@ def plot(self) -> None:
             ls="--",
         )
         self.label_plot(
-            "DD excitations -Diagram", "$dd-state-energy\,(1/cm)$", "$10Dq\,(1/cm)$"
+            "DD excitations -Diagram", r"$dd-state-energy\,(1/cm)$", r"$10Dq\,(1/cm)$",
         )
         plt.show()
 
@@ -107,7 +114,6 @@ def label_plot(self, arg0: str, arg1: str, arg2: str) -> None:
         plt.xlabel(arg2)
 
     def savetxt(self) -> None:
-
         pd.concat(
             [
                 self.df["delta_B"],
@@ -122,55 +128,55 @@ def savetxt(self) -> None:
         ).to_csv(Path(f"{self.title_DD}.csv"), index=False)
 
     def calculation(self) -> None:
-        """
-        Is filling the self.result with the iTS states of over-iterated energy range
+        """Is filling the self.result with the iTS states of over-iterated energy range
         """
         result = []
         for dq in self.df["Energy"]:
             if self.d_count == 2:  # d2
                 result.append(
                     self.subsplit_states(
-                        matrices.d2(Dq=dq, B=self.B, C=self.C).solver()
-                    )
+                        matrices.d2(Dq=dq, B=self.B, C=self.C).solver(),
+                    ),
                 )
             elif self.d_count == 3:  # d3
                 result.append(
                     self.subsplit_states(
-                        matrices.d3(Dq=dq, B=self.B, C=self.C).solver()
-                    )
+                        matrices.d3(Dq=dq, B=self.B, C=self.C).solver(),
+                    ),
                 )
             elif self.d_count == 4:  # d4
                 result.append(
                     self.subsplit_states(
-                        matrices.d4(Dq=dq, B=self.B, C=self.C).solver()
-                    )
+                        matrices.d4(Dq=dq, B=self.B, C=self.C).solver(),
+                    ),
                 )
             elif self.d_count == 5:  # d5
                 result.append(
                     self.subsplit_states(
-                        matrices.d5(Dq=dq, B=self.B, C=self.C).solver()
-                    )
+                        matrices.d5(Dq=dq, B=self.B, C=self.C).solver(),
+                    ),
                 )
             elif self.d_count == 6:  # d6
                 result.append(
                     self.subsplit_states(
-                        matrices.d6(Dq=dq, B=self.B, C=self.C).solver()
-                    )
+                        matrices.d6(Dq=dq, B=self.B, C=self.C).solver(),
+                    ),
                 )
             elif self.d_count == 7:  # d7
                 result.append(
                     self.subsplit_states(
-                        matrices.d7(Dq=dq, B=self.B, C=self.C).solver()
-                    )
+                        matrices.d7(Dq=dq, B=self.B, C=self.C).solver(),
+                    ),
                 )
             elif self.d_count == 8:  # d8
                 result.append(
                     self.subsplit_states(
-                        matrices.d8(Dq=dq, B=self.B, C=self.C).solver()
-                    )
+                        matrices.d8(Dq=dq, B=self.B, C=self.C).solver(),
+                    ),
                 )
             else:
-                raise ValueError("`d_count` must be in {2,3,4,5,6,7,8}")
+                msg = "The number of unpaired electrons should be between 2 and 8."
+                raise ValueError(msg)
 
         # Transform list of dictionaries to dictionary of arrays
         result = {
@@ -191,41 +197,33 @@ def subsplit_states(states: dict) -> dict:
         return rearranged_states
 
     def ci_cut(self, dq_ci: float = None) -> None:
-        """
-        Extracting the atomic-termsymbols for a specific dq depending on the oxidation state
+        """Extracting the atomic-termsymbols for a specific dq depending on the oxidation state
         """
         if self.d_count == 2:  # d2
-
             states = matrices.d2(Dq=dq_ci / 10.0, B=self.B, C=self.C).solver()
             self.ts_print(states, dq_ci=dq_ci)
 
         elif self.d_count == 3:  # d3
-
             states = matrices.d3(Dq=dq_ci / 10.0, B=self.B, C=self.C).solver()
             self.ts_print(states, dq_ci=dq_ci)
 
         elif self.d_count == 4:  # d4
-
             states = matrices.d4(Dq=dq_ci / 10.0, B=self.B, C=self.C).solver()
             self.ts_print(states, dq_ci=dq_ci)
 
         elif self.d_count == 5:  # d5
-
             states = matrices.d5(Dq=dq_ci / 10.0, B=self.B, C=self.C).solver()
             self.ts_print(states, dq_ci=dq_ci)
 
         elif self.d_count == 6:  # d6
-
             states = matrices.d6(Dq=dq_ci / 10.0, B=self.B, C=self.C).solver()
             self.ts_print(states, dq_ci=dq_ci)
 
         elif self.d_count == 7:  # d7
-
             states = matrices.d7(Dq=dq_ci / 10.0, B=self.B, C=self.C).solver()
             self.ts_print(states, dq_ci=dq_ci)
 
         elif self.d_count == 8:  # d8
-
             states = matrices.d8(Dq=dq_ci / 10.0, B=self.B, C=self.C).solver()
             self.ts_print(states, dq_ci=dq_ci)
 
@@ -241,11 +239,12 @@ def ts_print(self, states: dict, dq_ci: float = None) -> None:
             List of atomic-termsymbols for a specific oxidation state
         dq_ci : float, optional
             Specific crystalfield-splitting in Dq, by default None
+
         """
         count = 0
         cut = np.zeros(
             self._size + 1,
-            dtype=[("state", np.unicode_, 7), ("cm", int), ("eV", float)],
+            dtype=[("state", np.str_, 7), ("cm", int), ("eV", float)],
         )
         for irreducible in states:
             for energy in states[irreducible]:
@@ -288,7 +287,7 @@ def interactive_plot(self) -> None:
         if px is None:
             raise ImportError(
                 "Plotly is not installed. Please install plotly "
-                "with 'pip install tanabesugano[plotly]'!"
+                "with 'pip install tanabesugano[plotly]'!",
             )
 
         _col = self.df.drop(["Energy", "delta_B", "10Dq"], axis=1).columns
@@ -367,7 +366,7 @@ def cmd_line() -> None:
 
     parser = argparse.ArgumentParser(description=description)
     parser.add_argument(
-        "-d", type=int, default=6, help="Number of unpaired electrons (default d5)"
+        "-d", type=int, default=6, help="Number of unpaired electrons (default d5)",
     )
     parser.add_argument(
         "-Dq",
@@ -399,7 +398,7 @@ def cmd_line() -> None:
         "1.)",
     )
     parser.add_argument(
-        "-n", type=int, default=500, help="Number of roots (default nroots = 500)"
+        "-n", type=int, default=500, help="Number of roots (default nroots = 500)",
     )
     parser.add_argument(
         "-ndisp",
diff --git a/tanabesugano/matrices.py b/tanabesugano/matrices.py
index a14ae18..406801c 100644
--- a/tanabesugano/matrices.py
+++ b/tanabesugano/matrices.py
@@ -1,8 +1,14 @@
-from __future__ import print_function
+from __future__ import annotations
+
+from typing import Dict
+# from typing import TypeAlias
 
 import numpy as np
+
+from numpy._typing._array_like import NDArray
 from numpy.linalg import eigh
 
+
 _sqrt2 = np.sqrt(2.0)
 _sqrt3 = np.sqrt(3.0)
 _sqrt6 = np.sqrt(6.0)
@@ -13,105 +19,131 @@
 _3sqrt3 = _sqrt3 * 3.0
 _3sqrt6 = _sqrt6 * 3.0
 
+# Float64Array: TypeAlias = NDArray[np.float64]
+Float64Array = NDArray[np.float64]
+
+
+class LigandFieldTheory:
+    """Parent class for ligand field theory configurations."""
+
+    def __init__(self, Dq: float, B: float, C: float) -> None:
+        """Initializes the configuration with given parameters.
+
+        Args:
+            dq (float): Crystal field splitting in wavenumbers (cm-1).
+            b (float): Racah parameter B in wavenumbers (cm-1).
+            c (float): Racah parameter C in wavenumbers (cm-1).
 
-class d2(object):
-    def __init__(self, Dq: float = 0.0, B: float = 860.0, C: float = 3801.0):
-        """
-        parameter
-        ---------
-        All parameters in wavenumbers (cm-)
-        Dq: float
-                Crystalfield-Splitting
-        B: float
-                Racah-Parameter
-        C: float
-                Racah-Parameter
-        returns
-        -------
-        dictionary with elements of:
-                * Atomic-Termsymbols: str
-                * Eigen-Energies: float numpy-array
-                        Eigen-Energies of the atomic states depending on the crystalfield
         """
         self.Dq = np.float64(Dq)
         self.B = np.float64(B)
         self.C = np.float64(C)
 
-    def A_1_1_states(self):
-        # -  diagonal elements
+    def eigensolver(self, matrix: Float64Array) -> Float64Array:
+        """Solve for the eigenvalues of the given matrix.
+
+        Args:
+            matrix (Float64Array): 2-dimensional square array representing the TS matrix of the ligand field Hamiltonian.
+
+        Returns:
+            Float64Array: 1-dimensional array of eigenvalues of the diagonalized ligand field Hamiltonian.
+
+        """
+        return eigh(matrix)[0]
+
+    def solver(self) -> Dict[str, Float64Array]:
+        """Solve for all states and return a dictionary of results.
+
+        Returns:
+            Dict[str, Float64Array]: Dictionary with atomic term symbols as keys and eigenvalues as values.
+
+        """
+        msg = "Subclasses should implement this method."
+        raise NotImplementedError(msg)
+
+
+class d2(LigandFieldTheory):
+    """Class representing the d2 configuration in ligand field theory.
+    """
+
+    def __init__(self, Dq: float = 0.0, B: float = 860.0, C: float = 3801.0) -> None:
+        """Initializes the d2 configuration with given parameters.
+
+        Args:
+            Dq (float): Crystal field splitting in wavenumbers (cm-1).
+            B (float): Racah parameter B in wavenumbers (cm-1).
+            C (float): Racah parameter C in wavenumbers (cm-1).
+
+        """
+        super().__init__(Dq, B, C)
 
+    def A_1_1_states(self) -> Float64Array:
+        """Calculate the A_1_1 states."""
+        # Diagonal elements
         AA = -8 * self.Dq + 10 * self.B + 5 * self.C
         BB = +12 * self.Dq + 8 * self.B + 4 * self.C
 
-        # non diagonal elements
-
+        # Non-diagonal elements
         AB = BA = _sqrt6 * (2 * self.B + self.C)
 
         states = np.array([[AA, AB], [BA, BB]])
 
         return self.eigensolver(states)
 
-    def E_1_states(self):
-        # -  diagonal elements
-
+    def E_1_states(self) -> Float64Array:
+        """Calculate the E_1 states."""
+        # Diagonal elements
         AA = -8 * self.Dq + self.B + 2 * self.C
         BB = +12 * self.Dq + 2 * self.C
 
-        # non diagonal elements
-
+        # Non-diagonal elements
         AB = BA = -_2sqrt3 * self.B
 
         states = np.array([[AA, AB], [BA, BB]])
 
         return self.eigensolver(states)
 
-    def T_1_2_states(self):
-        # -  diagonal elements
-
+    def T_1_2_states(self) -> Float64Array:
+        """Calculate the T_1_2 states."""
+        # Diagonal elements
         AA = -8 * self.Dq + self.B + 2 * self.C
         BB = +2 * self.Dq + 2 * self.C
 
-        # non diagonal elements
-
+        # Non-diagonal elements
         AB = BA = +_2sqrt3 * self.B
 
         states = np.array([[AA, AB], [BA, BB]])
 
         return self.eigensolver(states)
 
-    def T_3_1_states(self):
-        # -  diagonal elements
-
+    def T_3_1_states(self) -> Float64Array:
+        """Calculate the T_3_1 states."""
+        # Diagonal elements
         AA = -8 * self.Dq - 5 * self.B
         BB = +2 * self.Dq + 4 * self.B
 
-        # non diagonal elements
-
+        # Non-diagonal elements
         AB = BA = 6 * self.B
 
         states = np.array([[AA, AB], [BA, BB]])
 
         return self.eigensolver(states)
 
-    def eigensolver(self, M):
-        """
-        :param M: 2 dimensional square array == TS matrics of Ligand field Hamiltonian
-        :return: 1 dimensiona                                 l array == eigenvalues of the diagonalized Ligand field Hamiltonian
-        """
+    def solver(self) -> Dict[str, Float64Array]:
+        """Solve for all states and return a dictionary of results.
 
-        return eigh(M)[0]
+        Returns:
+            Dict[str, Float64Array]: Dictionary with atomic term symbols as keys and eigenvalues as values.
 
-    def solver(self):
+        """
         # Ligand field independent states
-
-        # Ligendfield single depentent states
-
         GS = self.T_3_1_states()[0]
 
         T_1_1 = np.array([+2 * self.Dq + 4 * self.B + 2 * self.C]) - GS
         T_3_2 = np.array([+2 * self.Dq - 8 * self.B]) - GS
         A_3_2 = np.array([12 * self.Dq - 8 * self.B]) - GS
-        # Ligandfield dependent
+
+        # Ligand field dependent states
         A_1_1 = self.A_1_1_states() - GS
         E_1 = self.E_1_states() - GS
         T_1_2 = self.T_1_2_states() - GS
@@ -128,41 +160,30 @@ def solver(self):
         }
 
 
-class d3(object):
-    def __init__(self, Dq: float = 0.0, B: float = 918.0, C: float = 4133.0):
-        """
-        parameter
-        ---------
-        All parameters in wavenumbers (cm-)
-        Dq: float
-                Crystalfield-Splitting
-        B: float
-                Racah-Parameter
-        C: float
-                Racah-Parameter
-
-        returns
-        -------
-        dictionary with elements of:
-                * Atomic-Termsymbols: str
-                * Eigen-Energies: float numpy-array
-                        Eigen-Energies of the atomic states depending on the crystalfield
-        """
-        self.Dq = np.float64(Dq)
-        self.B = np.float64(B)
-        self.C = np.float64(C)
+class d3(LigandFieldTheory):
+    """Class representing the d3 configuration in ligand field theory."""
 
-    def T_2_2_states(self):
-        # -  diagonal elements
+    def __init__(self, Dq: float = 0.0, B: float = 918.0, C: float = 4133.0) -> None:
+        """Initializes the d3 configuration with given parameters.
+
+        Args:
+            Dq (float): Crystal field splitting in wavenumbers (cm-1).
+            B (float): Racah parameter B in wavenumbers (cm-1).
+            C (float): Racah parameter C in wavenumbers (cm-1).
 
+        """
+        super().__init__(Dq, B, C)
+
+    def T_2_2_states(self) -> Float64Array:
+        """Calculate the T_2_2 states."""
+        # Diagonal elements
         AA = -12 * self.Dq + 5 * self.C
         BB = -2 * self.Dq - 6 * self.B + 3 * self.C
         CC = -2 * self.Dq + 4 * self.B + 3 * self.C
         DD = +8 * self.Dq + 6 * self.B + 5 * self.C
         EE = +8 * self.Dq - 2 * self.B + 3 * self.C
 
-        # non diagonal elements
-
+        # Non-diagonal elements
         AB = BA = -_3sqrt3 * self.B
         AC = CA = -5 * _sqrt3 * self.B
         AD = DA = 4 * self.B + 2 * self.C
@@ -184,22 +205,21 @@ def T_2_2_states(self):
                 [CA, CB, CC, CD, CE],
                 [DA, DB, DC, DD, DE],
                 [EA, EB, EC, ED, EE],
-            ]
+            ],
         )
 
         return self.eigensolver(states)
 
-    def T_2_1_states(self):
-        # -  diagonal elements
-
+    def T_2_1_states(self) -> Float64Array:
+        """Calculate the T_2_1 states."""
+        # Diagonal elements
         AA = -12 * self.Dq - 6 * self.B + 3 * self.C
         BB = -2 * self.Dq + 3 * self.C
         CC = -2 * self.Dq - 6 * self.B + 3 * self.C
         DD = +8 * self.Dq - 6 * self.B + 3 * self.C
         EE = +8 * self.Dq - 2 * self.B + 3 * self.C
 
-        # non diagonal elements
-
+        # Non-diagonal elements
         AB = BA = -3 * self.B
         AC = CA = +3 * self.B
         AD = DA = 0.0
@@ -221,21 +241,20 @@ def T_2_1_states(self):
                 [CA, CB, CC, CD, CE],
                 [DA, DB, DC, DD, DE],
                 [EA, EB, EC, ED, EE],
-            ]
+            ],
         )
 
         return self.eigensolver(states)
 
-    def E_2_states(self):
-        # -  diagonal elements
-
+    def E_2_states(self) -> Float64Array:
+        """Calculate the E_2 states."""
+        # Diagonal elements
         AA = -12 * self.Dq - 6 * self.B + 3 * self.C
         BB = -2 * self.Dq + 8 * self.B + 6 * self.C
         CC = -2 * self.Dq - 1 * self.B + 3 * self.C
         DD = +18 * self.Dq - 8 * self.B + 4 * self.C
 
-        # non diagonal elements
-
+        # Non-diagonal elements
         AB = BA = -6 * _sqrt2 * self.B
         AC = CA = -_3sqrt2 * self.B
         AD = DA = 0.0
@@ -246,38 +265,32 @@ def E_2_states(self):
         CD = DC = _2sqrt3 * self.B
 
         states = np.array(
-            [[AA, AB, AC, AD], [BA, BB, BC, BD], [CA, CB, CC, CD], [DA, DB, DC, DD]]
+            [[AA, AB, AC, AD], [BA, BB, BC, BD], [CA, CB, CC, CD], [DA, DB, DC, DD]],
         )
 
         return self.eigensolver(states)
 
-    def T_4_1_states(self):
-        # -  diagonal elements
-
+    def T_4_1_states(self) -> Float64Array:
+        """Calculate the T_4_1 states."""
+        # Diagonal elements
         AA = -2 * self.Dq - 3 * self.B
         BB = +8 * self.Dq - 12 * self.B
 
-        # non diagonal elements
-
+        # Non-diagonal elements
         AB = BA = 6 * self.B
 
         states = np.array([[AA, AB], [BA, BB]])
 
         return self.eigensolver(states)
 
-    def eigensolver(self, M):
-        """
-        :param M: 2 dimensional square array == TS matrics of Ligand field Hamiltonian
-        :return: 1 dimensiona                                 l array == eigenvalues of the diagonalized Ligand field Hamiltonian
-        """
+    def solver(self) -> Dict[str, Float64Array]:
+        """Solve for all states and return a dictionary of results.
 
-        return eigh(M)[0]
+        Returns:
+            Dict[str, Float64Array]: Dictionary with atomic term symbols as keys and eigenvalues as values.
 
-    def solver(self):
+        """
         # Ligand field independent states
-
-        # Ligendfield single depentent states
-
         GS = np.array([-12 * self.Dq - 15 * self.B])
 
         A_4_2 = np.array([0], dtype=np.float64)
@@ -286,7 +299,7 @@ def solver(self):
         A_2_1 = np.array([-2 * self.Dq - 11 * self.B + 3 * self.C]) - GS
         A_2_2 = np.array([-2 * self.Dq + 9 * self.B + 3 * self.C]) - GS
 
-        # Ligandfield dependent
+        # Ligand field dependent states
         T_2_2 = self.T_2_2_states() - GS
         T_2_1 = self.T_2_1_states() - GS
         E_2 = self.E_2_states() - GS
@@ -304,32 +317,23 @@ def solver(self):
         }
 
 
-class d4(object):
-    def __init__(self, Dq: float = 0.0, B: float = 965.0, C: float = 4449.0):
-        """
-        parameter
-        ---------
-        All parameters in wavenumbers (cm-)
-        Dq: float
-                Crystalfield-Splitting
-        B: float
-                Racah-Parameter
-        C: float
-                Racah-Parameter
-        returns
-        -------
-        dictionary with elements of:
-                * Atomic-Termsymbols: str
-                * Eigen-Energies: float numpy-array
-                * Eigen-Energies of the atomic states depending on the crystalfield
-        """
-        self.Dq = np.float64(Dq)
-        self.B = np.float64(B)
-        self.C = np.float64(C)
+class d4(LigandFieldTheory):
+    """Class representing the d4 configuration in ligand field theory."""
 
-    def T_3_1_states(self):
-        # -  diagonal elements
+    def __init__(self, Dq: float = 0.0, B: float = 965.0, C: float = 4449.0) -> None:
+        """Initializes the d4 configuration with given parameters.
 
+        Args:
+            Dq (float): Crystal field splitting in wavenumbers (cm-1).
+            B (float): Racah parameter B in wavenumbers (cm-1).
+            C (float): Racah parameter C in wavenumbers (cm-1).
+
+        """
+        super().__init__(Dq, B, C)
+
+    def T_3_1_states(self) -> Float64Array:
+        """Calculate the T_3_1 states."""
+        # Diagonal elements
         AA = -16 * self.Dq - 15 * self.B + 5 * self.C
         BB = -6 * self.Dq - 11 * self.B + 4 * self.C
         CC = -6 * self.Dq - 3 * self.B + 6 * self.C
@@ -338,8 +342,7 @@ def T_3_1_states(self):
         FF = +4 * self.Dq - 11 * self.B + 4 * self.C
         GG = +14 * self.Dq - 16 * self.B + 5 * self.C
 
-        # non diagonal elements
-
+        # Non-diagonal elements
         AB = BA = -_sqrt6 * self.B
         AC = CA = -_3sqrt2 * self.B
         AD = DA = _sqrt2 * (2 * self.B + self.C)
@@ -376,14 +379,14 @@ def T_3_1_states(self):
                 [EA, EB, EC, ED, EE, EF, EG],
                 [FA, FB, FC, FD, FE, FF, FG],
                 [GA, GB, GC, GD, GE, GF, GG],
-            ]
+            ],
         )
 
         return self.eigensolver(states)
 
-    def T_1_2_states(self):
-        # diagonal elements
-
+    def T_1_2_states(self) -> Float64Array:
+        """Calculate the T_1_2 states."""
+        # Diagonal elements
         AA = -16 * self.Dq - 9 * self.B + 7 * self.C
         BB = -6 * self.Dq - 9 * self.B + 6 * self.C
         CC = -6 * self.Dq + 3 * self.B + 8 * self.C
@@ -392,8 +395,7 @@ def T_1_2_states(self):
         FF = +4 * self.Dq + 5 * self.B + 8 * self.C
         GG = +14 * self.Dq + 7 * self.C
 
-        # non diagonal elements
-
+        # Non-diagonal elements
         AB = BA = _3sqrt2 * self.B
         AC = CA = -5 * _sqrt6 * self.B
         AD = DA = 0.0
@@ -430,22 +432,21 @@ def T_1_2_states(self):
                 [EA, EB, EC, ED, EE, EF, EG],
                 [FA, FB, FC, FD, FE, FF, FG],
                 [GA, GB, GC, GD, GE, GF, GG],
-            ]
+            ],
         )
 
         return self.eigensolver(states)
 
-    def A_1_1_states(self):
-        # diagonal elements
-
+    def A_1_1_states(self) -> Float64Array:
+        """Calculate the A_1_1 states."""
+        # Diagonal elements
         AA = -16 * self.Dq + 10 * self.C
         BB = -6 * self.Dq + 6 * self.C
         CC = +4 * self.Dq + 14 * self.B + 11 * self.C
         DD = +4 * self.Dq - 3 * self.B + 6 * self.C
         EE = +24 * self.Dq - 16 * self.B + 8 * self.C
 
-        # non diagonal elements
-
+        # Non-diagonal elements
         AB = BA = -12 * _sqrt2 * self.B
         AC = CA = _sqrt2 * (4 * self.B + 2 * self.C)
         AD = DA = _2sqrt2 * self.B
@@ -467,22 +468,21 @@ def A_1_1_states(self):
                 [CA, CB, CC, CD, CE],
                 [DA, DB, DC, DD, DE],
                 [EA, EB, EC, ED, EE],
-            ]
+            ],
         )
 
         return self.eigensolver(states)
 
-    def E_1_1_states(self):
-        # diagonal elements
-
+    def E_1_1_states(self) -> Float64Array:
+        """Calculate the E_1_1 states."""
+        # Diagonal elements
         AA = -16 * self.Dq - 9 * self.B + 7 * self.C
         BB = -6 * self.Dq - 6 * self.B + 6 * self.C
         CC = +4 * self.Dq + 5 * self.B + 8 * self.C
         DD = +4 * self.Dq + 6 * self.B + 9 * self.C
         EE = +4 * self.Dq - 3 * self.B + 6 * self.C
 
-        # non diagonal elements
-
+        # Non-diagonal elements
         AB = BA = 6 * self.B
         AC = CA = _sqrt2 * (2 * self.B + self.C)
         AD = DA = -2 * self.B
@@ -504,22 +504,21 @@ def E_1_1_states(self):
                 [CA, CB, CC, CD, CE],
                 [DA, DB, DC, DD, DE],
                 [EA, EB, EC, ED, EE],
-            ]
+            ],
         )
 
         return self.eigensolver(states)
 
-    def T_3_2_states(self):
-        # diagonal elements
-
+    def T_3_2_states(self) -> Float64Array:
+        """Calculate the T_3_2 states."""
+        # Diagonal elements
         AA = -6 * self.Dq - 9 * self.B + 4 * self.C
         BB = -6 * self.Dq - 5 * self.B + 6 * self.C
         CC = +4 * self.Dq - 13 * self.B + 4 * self.C
         DD = +4 * self.Dq - 9 * self.B + 4 * self.C
         EE = +14 * self.Dq - 8 * self.B + 5 * self.C
 
-        # non diagonal elements
-
+        # Non-diagonal elements
         AB = BA = -5 * _sqrt3 * self.B
         AC = CA = _sqrt6 * self.B
         AD = DA = _sqrt3 * self.B
@@ -541,21 +540,20 @@ def T_3_2_states(self):
                 [CA, CB, CC, CD, CE],
                 [DA, DB, DC, DD, DE],
                 [EA, EB, EC, ED, EE],
-            ]
+            ],
         )
 
         return self.eigensolver(states)
 
-    def T_1_1_states(self):
-        # diagonal elements
-
+    def T_1_1_states(self) -> Float64Array:
+        """Calculate the T_1_1 states."""
+        # Diagonal elements
         AA = -6 * self.Dq - 3 * self.B + 6 * self.C
         BB = -6 * self.Dq - 3 * self.B + 8 * self.C
         CC = +4 * self.Dq - 3 * self.B + 6 * self.C
         DD = +14 * self.Dq - 16 * self.B + 7 * self.C
 
-        # non diagonal elements
-
+        # Non-diagonal elements
         AB = BA = 5 * _sqrt3 * self.B
         AC = CA = 3 * self.B
         AD = DA = _sqrt6 * self.B
@@ -566,20 +564,19 @@ def T_1_1_states(self):
         CD = DC = -_sqrt6 * self.B
 
         states = np.array(
-            [[AA, AB, AC, AD], [BA, BB, BC, BD], [CA, CB, CC, CD], [DA, DB, DC, DD]]
+            [[AA, AB, AC, AD], [BA, BB, BC, BD], [CA, CB, CC, CD], [DA, DB, DC, DD]],
         )
 
         return self.eigensolver(states)
 
-    def E_3_1_states(self):
-        # diagonal elements
-
+    def E_3_1_states(self) -> Float64Array:
+        """Calculate the E_3_1 states."""
+        # Diagonal elements
         AA = -6 * self.Dq - 13 * self.B + 4 * self.C
         BB = -6 * self.Dq - 10 * self.B + 4 * self.C
         CC = +4 * self.Dq - 11 * self.B + 4 * self.C
 
-        # non diagonal elements
-
+        # Non-diagonal elements
         AB = BA = -4 * self.B
         AC = CA = 0.0
 
@@ -589,52 +586,40 @@ def E_3_1_states(self):
 
         return self.eigensolver(states)
 
-    def A_3_2_states(self):
-        # diagonal elements
-
+    def A_3_2_states(self) -> Float64Array:
+        """Calculate the A_3_2 states."""
+        # Diagonal elements
         AA = -6 * self.Dq - 8 * self.B + 4 * self.C
         BB = +4 * self.Dq - 2 * self.B + 7 * self.C
 
-        # non diagonal elements
-
+        # Non-diagonal elements
         AB = BA = -12 * self.B
 
         states = np.array([[AA, AB], [BA, BB]])
 
         return self.eigensolver(states)
 
-    def A_1_2_states(self):
-        # diagonal elements
-
+    def A_1_2_states(self) -> Float64Array:
+        """Calculate the A_1_2 states."""
+        # Diagonal elements
         AA = -6 * self.Dq - 12 * self.B + 6 * self.C
         BB = +4 * self.Dq - 3 * self.B + 6 * self.C
 
-        # non diagonal elements
-
+        # Non-diagonal elements
         AB = BA = 6 * self.B
 
         states = np.array([[AA, AB], [BA, BB]])
 
         return self.eigensolver(states)
 
-    def eigensolver(self, M):
-        """
-        :param M: 2 dimensional square array == TS matrics of Ligand field Hamiltonian
-        :return: 1 dimensional array == eigenvalues of the diagonalized Ligand field Hamiltonian
-        """
+    def solver(self) -> Dict[str, Float64Array]:
+        """Solve for all states and return a dictionary of results.
 
-        return eigh(M)[0]
+        Returns:
+            Dict[str, Float64Array]: Dictionary with atomic term symbols as keys and eigenvalues as values.
 
-    def solver(self):
+        """
         # Ligand field independent states
-
-        # A_6_1 = np.array( [ 0 ] )  # Starting value is -35. * B, but has to set to zero per definition
-        # E_4 = self.E_4_states( ) + 35 * self.B
-        # A_4_1 = np.array( [ -25 * self.B + 5 * self.C ] ) + 35 * self.B
-        # A_4_2 = np.array( [ -13 * self.B + 7 * self.C ] ) + 35 * self.B
-
-        # Ligendfield single depentent states
-
         GS = np.array([-6 * self.Dq - 21 * self.B])
 
         E_5_1 = np.array([0], dtype=np.float64)
@@ -642,7 +627,7 @@ def solver(self):
 
         A_3_1 = np.array([-6 * self.Dq - 12 * self.B + 4 * self.C]) - GS
 
-        # Ligandfield dependent
+        # Ligand field dependent states
         T_1_2 = self.T_1_2_states() - GS
         T_3_1 = self.T_3_1_states() - GS
         A_1_1 = self.A_1_1_states() - GS
@@ -666,6 +651,7 @@ def solver(self):
             T_5_2 -= T_3_1[0]
             A_3_1 -= T_3_1[0]
             T_3_1 -= T_3_1[0]
+
         return {
             "3_T_1": T_3_1,
             "1_T_2": T_1_2,
@@ -682,31 +668,22 @@ def solver(self):
         }
 
 
-class d5(object):
-    def __init__(self, Dq: float = 0.0, B: float = 860.0, C: float = 3850.0):
-        """
-        parameter
-        ---------
-        All parameters in wavenumbers (cm-)
-        Dq: float
-                Crystalfield-Splitting
-        B: float
-                Racah-Parameter
-        C: float
-                Racah-Parameter
-
-        returns
-        -------
-        dictionary with elements of:
-                * Atomic-Termsymbols: str
-                * Eigen-Energies: float numpy-array
-                * Eigen-Energies of the atomic states depending on the crystalfield
+class d5(LigandFieldTheory):
+    """Class representing the d5 configuration in ligand field theory."""
+
+    def __init__(self, Dq: float = 0.0, B: float = 860.0, C: float = 3850.0) -> None:
+        """Initializes the d5 configuration with given parameters.
+
+        Args:
+            Dq (float): Crystal field splitting in wavenumbers (cm-1).
+            B (float): Racah parameter B in wavenumbers (cm-1).
+            C (float): Racah parameter C in wavenumbers (cm-1).
+
         """
-        self.Dq = np.float64(Dq)
-        self.B = np.float64(B)
-        self.C = np.float64(C)
+        super().__init__(Dq, B, C)
 
-    def T_2_2_states(self):
+    def T_2_2_states(self) -> Float64Array:
+        """Calculate the T_2_2 states."""
         # diagonal elements
 
         AA = -20 * self.Dq - 20 * self.B + 10 * self.C
@@ -788,12 +765,13 @@ def T_2_2_states(self):
                 [HA, HB, HC, HD, HE, HF, HG, HH, HI, HJ],
                 [IA, IB, IC, ID, IE, IF, IG, IH, II, IJ],
                 [JA, JB, JC, JD, JE, JF, JG, JH, JI, JJ],
-            ]
+            ],
         )
 
         return self.eigensolver(states)
 
-    def T_2_1_states(self):
+    def T_2_1_states(self) -> Float64Array:
+        """Calculate the T_2_1 states."""
         # diagonal elements
 
         AA = -10 * self.Dq - 22 * self.B + 9 * self.C
@@ -852,12 +830,13 @@ def T_2_1_states(self):
                 [FA, FB, FC, FD, FE, FF, FG, FH],
                 [GA, GB, GC, GD, GE, GF, GG, GH],
                 [HA, HB, HC, HD, HE, HF, HG, HH],
-            ]
+            ],
         )
 
         return self.eigensolver(states)
 
-    def E_2_states(self):
+    def E_2_states(self) -> Float64Array:
+        """Calculate the E_2 states."""
         # diagonal elements
 
         AA = -10 * self.Dq - 4 * self.B + 12 * self.C
@@ -906,12 +885,13 @@ def E_2_states(self):
                 [EA, EB, EC, ED, EE, EF, EG],
                 [FA, FB, FC, FD, FE, FF, FG],
                 [GA, GB, GC, GD, GE, GF, GG],
-            ]
+            ],
         )
 
         return self.eigensolver(states)
 
-    def A_2_1_states(self):
+    def A_2_1_states(self) -> Float64Array:
+        """Calculate the A_2_1 states."""
         # diagonal elements
 
         AA = -10 * self.Dq - 3 * self.B + 9 * self.C
@@ -931,12 +911,13 @@ def A_2_1_states(self):
         CD = DC = 0.0
 
         states = np.array(
-            [[AA, AB, AC, AD], [BA, BB, BC, BD], [CA, CB, CC, CD], [DA, DB, DC, DD]]
+            [[AA, AB, AC, AD], [BA, BB, BC, BD], [CA, CB, CC, CD], [DA, DB, DC, DD]],
         )
 
         return self.eigensolver(states)
 
-    def A_2_2_states(self):
+    def A_2_2_states(self) -> Float64Array:
+        """Calculate the A_2_2 states."""
         # diagonal elements
 
         AA = -10 * self.Dq - 23 * self.B + 9 * self.C
@@ -954,7 +935,8 @@ def A_2_2_states(self):
 
         return self.eigensolver(states)
 
-    def T_4_1_states(self):
+    def T_4_1_states(self) -> Float64Array:
+        """Calculate the T_4_1 states."""
         # diagonal elements
 
         AA = -10 * self.Dq - 25 * self.B + 6 * self.C
@@ -972,7 +954,8 @@ def T_4_1_states(self):
 
         return self.eigensolver(states)
 
-    def T_4_2_states(self):
+    def T_4_2_states(self) -> Float64Array:
+        """Calculate the T_4_2 states."""
         # diagonal elements
 
         AA = -10 * self.Dq - 17 * self.B + 6 * self.C
@@ -991,7 +974,8 @@ def T_4_2_states(self):
 
         return self.eigensolver(states)
 
-    def E_4_states(self):
+    def E_4_states(self) -> Float64Array:
+        """Calculate the E_4 states."""
         # diagonal elements
 
         AA = -22 * self.B + 5 * self.C
@@ -1005,20 +989,18 @@ def E_4_states(self):
 
         return self.eigensolver(states)
 
-    def eigensolver(self, M):
-        """
-        :param M: 2 dimensional square array == TS matrics of Ligand field Hamiltonian
-        :return: 1 dimensional array == eigenvalues of the diagonalized Ligand field Hamiltonian
-        """
+    def solver(self) -> Dict[str, Float64Array]:
+        """Solve for all states and return a dictionary of results.
 
-        return eigh(M)[0]
+        Returns:
+            Dict[str, Float64Array]: Dictionary with atomic term symbols as keys and eigenvalues as values.
 
-    def solver(self):
+        """
         # Ligand field independent states
         GS = -35 * self.B
 
         A_6_1 = np.array(
-            [0.0], dtype=float
+            [0.0], dtype=float,
         )  # Starting value is -35. * B, but has to set to zero per definition
         E_4 = self.E_4_states() - GS
         A_4_1 = np.array([-25 * self.B + 5 * self.C]) - GS
@@ -1062,31 +1044,22 @@ def solver(self):
         }
 
 
-class d6(object):
+class d6(LigandFieldTheory):
+    """Class representing the d6 configuration in ligand field theory."""
+
     def __init__(self, Dq: float = 0.0, B: float = 1065.0, C: float = 5120.0):
+        """Initializes the d6 configuration with given parameters.
+
+        Args:
+            Dq (float): Crystal field splitting in wavenumbers (cm-1).
+            B (float): Racah parameter B in wavenumbers (cm-1).
+            C (float): Racah parameter C in wavenumbers (cm-1).
+
         """
-        parameter
-        ---------
-        All parameters in wavenumbers (cm-)
-        Dq: float
-                Crystalfield-Splitting
-        B: float
-                Racah-Parameter
-        C: float
-                Racah-Parameter
-
-        returns
-        -------
-        dictionary with elements of:
-                * Atomic-Termsymbols: str
-                * Eigen-Energies: float numpy-array
-                        Eigen-Energies of the atomic states depending on the crystalfield
-        """
-        self.Dq = np.float64(Dq)
-        self.B = np.float64(B)
-        self.C = np.float64(C)
+        super().__init__(Dq, B, C)
 
-    def T_3_1_states(self):
+    def T_3_1_states(self) -> Float64Array:
+        """Calculate the T_3_1 states."""
         # -  diagonal elements
 
         AA = +16 * self.Dq - 15 * self.B + 5 * self.C
@@ -1135,12 +1108,13 @@ def T_3_1_states(self):
                 [EA, EB, EC, ED, EE, EF, EG],
                 [FA, FB, FC, FD, FE, FF, FG],
                 [GA, GB, GC, GD, GE, GF, GG],
-            ]
+            ],
         )
 
         return self.eigensolver(states)
 
-    def T_1_2_states(self):
+    def T_1_2_states(self) -> Float64Array:
+        """Calculate the T_1_2 states."""
         # diagonal elements
 
         AA = +16 * self.Dq - 9 * self.B + 7 * self.C
@@ -1189,12 +1163,13 @@ def T_1_2_states(self):
                 [EA, EB, EC, ED, EE, EF, EG],
                 [FA, FB, FC, FD, FE, FF, FG],
                 [GA, GB, GC, GD, GE, GF, GG],
-            ]
+            ],
         )
 
         return self.eigensolver(states)
 
-    def A_1_1_states(self):
+    def A_1_1_states(self) -> Float64Array:
+        """Calculate the A_1_1 states."""
         # diagonal elements
 
         AA = +16 * self.Dq + 10 * self.C
@@ -1226,12 +1201,13 @@ def A_1_1_states(self):
                 [CA, CB, CC, CD, CE],
                 [DA, DB, DC, DD, DE],
                 [EA, EB, EC, ED, EE],
-            ]
+            ],
         )
 
         return self.eigensolver(states)
 
-    def E_1_1_states(self):
+    def E_1_1_states(self) -> Float64Array:
+        """Calculate the E_1_1 states."""
         # diagonal elements
 
         AA = +16 * self.Dq - 9 * self.B + 7 * self.C
@@ -1263,12 +1239,13 @@ def E_1_1_states(self):
                 [CA, CB, CC, CD, CE],
                 [DA, DB, DC, DD, DE],
                 [EA, EB, EC, ED, EE],
-            ]
+            ],
         )
 
         return self.eigensolver(states)
 
-    def T_3_2_states(self):
+    def T_3_2_states(self) -> Float64Array:
+        """Calculate the T_3_2 states."""
         # diagonal elements
 
         AA = +6 * self.Dq - 9 * self.B + 4 * self.C
@@ -1300,12 +1277,13 @@ def T_3_2_states(self):
                 [CA, CB, CC, CD, CE],
                 [DA, DB, DC, DD, DE],
                 [EA, EB, EC, ED, EE],
-            ]
+            ],
         )
 
         return self.eigensolver(states)
 
-    def T_1_1_states(self):
+    def T_1_1_states(self) -> Float64Array:
+        """Calculate the T_1_1 states."""
         # diagonal elements
 
         AA = +6 * self.Dq - 3 * self.B + 6 * self.C
@@ -1325,12 +1303,13 @@ def T_1_1_states(self):
         CD = DC = -_sqrt6 * self.B
 
         states = np.array(
-            [[AA, AB, AC, AD], [BA, BB, BC, BD], [CA, CB, CC, CD], [DA, DB, DC, DD]]
+            [[AA, AB, AC, AD], [BA, BB, BC, BD], [CA, CB, CC, CD], [DA, DB, DC, DD]],
         )
 
         return self.eigensolver(states)
 
-    def E_3_1_states(self):
+    def E_3_1_states(self) -> Float64Array:
+        """Calculate the E_3_1 states."""
         # diagonal elements
 
         AA = +6 * self.Dq - 13 * self.B + 4 * self.C
@@ -1348,7 +1327,8 @@ def E_3_1_states(self):
 
         return self.eigensolver(states)
 
-    def A_3_2_states(self):
+    def A_3_2_states(self) -> Float64Array:
+        """Calculate the A_3_2 states."""
         # diagonal elements
 
         AA = +6 * self.Dq - 8 * self.B + 4 * self.C
@@ -1362,7 +1342,8 @@ def A_3_2_states(self):
 
         return self.eigensolver(states)
 
-    def A_1_2_states(self):
+    def A_1_2_states(self) -> Float64Array:
+        """Calculate the A_1_2 states."""
         # diagonal elements
 
         AA = +6 * self.Dq - 12 * self.B + 6 * self.C
@@ -1376,24 +1357,15 @@ def A_1_2_states(self):
 
         return self.eigensolver(states)
 
-    def eigensolver(self, M):
-        """
-        :param M: 2 dimensional square array == TS matrics of Ligand field Hamiltonian
-        :return: 1 dimensional array == eigenvalues of the diagonalized Ligand field Hamiltonian
-        """
-
-        return eigh(M)[0]
 
-    def solver(self):
-        # Ligand field independent states
 
-        # A_6_1 = np.array( [ 0 ] )  # Starting value is -35. * B, but has to set to zero per definition
-        # E_4 = self.E_4_states( ) + 35 * self.B
-        # A_4_1 = np.array( [ -25 * self.B + 5 * self.C ] ) + 35 * self.B
-        # A_4_2 = np.array( [ -13 * self.B + 7 * self.C ] ) + 35 * self.B
+    def solver(self) -> Dict[str, Float64Array]:
+        """Solve for all states and return a dictionary of results.
 
-        # Ligendfield single depentent states
+        Returns:
+            Dict[str, Float64Array]: Dictionary with atomic term symbols as keys and eigenvalues as values.
 
+        """
         GS = np.array([-4 * self.Dq - 21 * self.B])
 
         T_5_2 = np.array([0], dtype=np.float64)
@@ -1445,31 +1417,22 @@ def solver(self):
         }
 
 
-class d7(object):
+class d7(LigandFieldTheory):
+    """Class for d7 configuration."""
+
     def __init__(self, Dq: float = 0.0, B: float = 971.0, C: float = 4499.0):
+        """Initializes the d7 configuration with given parameters.
+
+        Args:
+            Dq (float): Crystal field splitting in wavenumbers (cm-1).
+            B (float): Racah parameter B in wavenumbers (cm-1).
+            C (float): Racah parameter C in wavenumbers (cm-1).
+
         """
-        parameter
-        ---------
-        All parameters in wavenumbers (cm-)
-        Dq: float
-                Crystalfield-Splitting
-        B: float
-                Racah-Parameter
-        C: float
-                Racah-Parameter
-
-        returns
-        -------
-        dictionary with elements of:
-                * Atomic-Termsymbols: str
-                * Eigen-Energies: float numpy-array
-                * Eigen-Energies of the atomic states depending on the crystalfield
-        """
-        self.Dq = np.float64(Dq)
-        self.B = np.float64(B)
-        self.C = np.float64(C)
+        super().__init__(Dq, B, C)
 
     def T_2_2_states(self):
+        """Calculate the T_2_2 states."""
         # -  diagonal elements
 
         AA = +12 * self.Dq + 5 * self.C
@@ -1501,12 +1464,13 @@ def T_2_2_states(self):
                 [CA, CB, CC, CD, CE],
                 [DA, DB, DC, DD, DE],
                 [EA, EB, EC, ED, EE],
-            ]
+            ],
         )
 
         return self.eigensolver(states)
 
     def T_2_1_states(self):
+        """Calculate the T_2_1 states."""
         # -  diagonal elements
 
         AA = +12 * self.Dq - 6 * self.B + 3 * self.C
@@ -1538,12 +1502,13 @@ def T_2_1_states(self):
                 [CA, CB, CC, CD, CE],
                 [DA, DB, DC, DD, DE],
                 [EA, EB, EC, ED, EE],
-            ]
+            ],
         )
 
         return self.eigensolver(states)
 
     def E_2_states(self):
+        """Calculate the E_2 states."""
         # -  diagonal elements
 
         AA = +12 * self.Dq - 6 * self.B + 3 * self.C
@@ -1563,12 +1528,13 @@ def E_2_states(self):
         CD = DC = _2sqrt3 * self.B
 
         states = np.array(
-            [[AA, AB, AC, AD], [BA, BB, BC, BD], [CA, CB, CC, CD], [DA, DB, DC, DD]]
+            [[AA, AB, AC, AD], [BA, BB, BC, BD], [CA, CB, CC, CD], [DA, DB, DC, DD]],
         )
 
         return self.eigensolver(states)
 
     def T_4_1_states(self):
+        """Calculate the T_4_1 states."""
         # -  diagonal elements
 
         AA = +2 * self.Dq - 3 * self.B
@@ -1582,15 +1548,15 @@ def T_4_1_states(self):
 
         return self.eigensolver(states)
 
-    def eigensolver(self, M: np.ndarray):
-        """
-        :param M: 2 dimensional square array == TS matrics of Ligand field Hamiltonian
-        :return: 1 dimensiona                                 l array == eigenvalues of the diagonalized Ligand field Hamiltonian
-        """
 
-        return eigh(M)[0]
 
-    def solver(self):
+    def solver(self) -> Dict[str, np.ndarray]:
+        """Solve for all states and return a dictionary of results.
+
+        Returns:
+            Dict[str, np.ndarray]: Dictionary with atomic term symbols as keys and eigenvalues as values.
+
+        """
         # Ligand field independent states
 
         # Ligandfield multi depnedent state become GS
@@ -1636,31 +1602,20 @@ def solver(self):
         }
 
 
-class d8(object):
+class d8(LigandFieldTheory):
     def __init__(self, Dq: float = 0.0, B: float = 1030.0, C: float = 4850.0):
+        """Initializes the d8 configuration with given parameters.
+
+        Args:
+            Dq (float): Crystal field splitting in wavenumbers (cm-1).
+            B (float): Racah parameter B in wavenumbers (cm-1).
+            C (float): Racah parameter C in wavenumbers (cm-1).
+
         """
-        parameter
-        ---------
-        All parameters in wavenumbers (cm-)
-        Dq: float
-                Crystalfield-Splitting
-        B: float
-                Racah-Parameter
-        C: float
-                Racah-Parameter
-
-        returns
-        -------
-        dictionary with elements of:
-                * Atomic-Termsymbols: str
-                * Eigen-Energies: float numpy-array
-                * Eigen-Energies of the atomic states depending on the crystalfield
-        """
-        self.Dq = np.float64(Dq)
-        self.B = np.float64(B)
-        self.C = np.float64(C)
+        super().__init__(Dq, B, C)
 
     def A_1_1_states(self):
+        """Calculate the A_1_1 states."""
         # -  diagonal elements
 
         AA = +8 * self.Dq + 10 * self.B + 5 * self.C
@@ -1675,6 +1630,7 @@ def A_1_1_states(self):
         return self.eigensolver(states)
 
     def E_1_states(self):
+        """Calculate the E_1 states."""
         # -  diagonal elements
 
         AA = +8 * self.Dq + self.B + 2 * self.C
@@ -1689,6 +1645,7 @@ def E_1_states(self):
         return self.eigensolver(states)
 
     def T_1_2_states(self):
+        """Calculate the T_1_2 states."""
         # -  diagonal elements
 
         AA = +8 * self.Dq + self.B + 2 * self.C
@@ -1703,6 +1660,7 @@ def T_1_2_states(self):
         return self.eigensolver(states)
 
     def T_3_1_states(self):
+        """Calculate the T_3_1 states."""
         # -  diagonal elements
 
         AA = +8 * self.Dq - 5 * self.B
@@ -1716,15 +1674,15 @@ def T_3_1_states(self):
 
         return self.eigensolver(states)
 
-    def eigensolver(self, M):
-        """
-        :param M: 2 dimensional square array == TS matrics of Ligand field Hamiltonian
-        :return: 1 dimensiona                                 l array == eigenvalues of the diagonalized Ligand field Hamiltonian
-        """
 
-        return eigh(M)[0]
 
-    def solver(self):
+    def solver(self) -> Dict[str, Float64Array]:
+        """Solve for all states and return a dictionary of results.
+
+        Returns:
+            Dict[str, Float64Array]: Dictionary with atomic term symbols as keys and eigenvalues as values.
+
+        """
         # Ligand field independent states
 
         # Ligendfield single depentent states
diff --git a/tanabesugano/test/test_batch.py b/tanabesugano/test/test_batch.py
index 51375ff..9463848 100644
--- a/tanabesugano/test/test_batch.py
+++ b/tanabesugano/test/test_batch.py
@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 from tanabesugano.batch import Batch
 
 
diff --git a/tanabesugano/test/test_front.py b/tanabesugano/test/test_front.py
index ddf0935..3698143 100644
--- a/tanabesugano/test/test_front.py
+++ b/tanabesugano/test/test_front.py
@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 from typing import Any
 
 from tanabesugano import cmd as frontapp
@@ -5,7 +7,7 @@
 
 def test_frontapp():
     return frontapp.CMDmain(
-        Dq=4000.0, B=400.0, C=3600.0, nroots=100, d_count=5
+        Dq=4000.0, B=400.0, C=3600.0, nroots=100, d_count=5,
     ).calculation()
 
 
diff --git a/tanabesugano/test/test_num.py b/tanabesugano/test/test_num.py
index 3efbd04..994aebf 100644
--- a/tanabesugano/test/test_num.py
+++ b/tanabesugano/test/test_num.py
@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import numpy as np
 
 from tanabesugano import matrices
@@ -8,16 +10,16 @@ def state_check(x):
         if x == 3:
             states = matrices.d3(Dq=i).solver()
             return len(states)
-        elif x == 4:
+        if x == 4:
             states = matrices.d4(Dq=i).solver()
             return len(states)
-        elif x == 5:
+        if x == 5:
             states = matrices.d5(Dq=i).solver()
             return len(states)
-        elif x == 6:
+        if x == 6:
             states = matrices.d6(Dq=i).solver()
             return len(states)
-        elif x == 7:
+        if x == 7:
             states = matrices.d7(Dq=i).solver()
             return len(states)
 
diff --git a/tanabesugano/tools.py b/tanabesugano/tools.py
index b885e29..f62c1a9 100644
--- a/tanabesugano/tools.py
+++ b/tanabesugano/tools.py
@@ -1,10 +1,13 @@
-from typing import Tuple, Union
+from __future__ import annotations
+
+from typing import Tuple
+from typing import Union
 
 import numpy as np
 
 
 def racah(
-    F2: Union[float, np.array], F4: Union[float, np.array]
+    F2: Union[float, np.array], F4: Union[float, np.array],
 ) -> Union[Tuple[float, float], Tuple[np.array, np.array]]:
     """Transform the Slater-Condon-Parameter to Racah-Parameter.
 
@@ -17,6 +20,7 @@ def racah(
 
     Returns:
         Union[Tuple[float, float], Tuple[np.array, np.array]]: Return the racah parameters.
+
     """
     eVcm = 8065.54
     B = eVcm * (F2 / 49.0 - 5 / 441.0 * F4)