diff --git a/pyLiSW/__pycache__/Atoms.cpython-312.pyc b/pyLiSW/__pycache__/Atoms.cpython-312.pyc
new file mode 100644
index 0000000..59b60e2
Binary files /dev/null and b/pyLiSW/__pycache__/Atoms.cpython-312.pyc differ
diff --git a/pyLiSW/__pycache__/Bonds.cpython-312.pyc b/pyLiSW/__pycache__/Bonds.cpython-312.pyc
new file mode 100644
index 0000000..99e839a
Binary files /dev/null and b/pyLiSW/__pycache__/Bonds.cpython-312.pyc differ
diff --git a/pyLiSW/__pycache__/LSWT.cpython-312.pyc b/pyLiSW/__pycache__/LSWT.cpython-312.pyc
new file mode 100644
index 0000000..5320cea
Binary files /dev/null and b/pyLiSW/__pycache__/LSWT.cpython-312.pyc differ
diff --git a/pyLiSW/__pycache__/QEspace.cpython-312.pyc b/pyLiSW/__pycache__/QEspace.cpython-312.pyc
new file mode 100644
index 0000000..d6c738d
Binary files /dev/null and b/pyLiSW/__pycache__/QEspace.cpython-312.pyc differ
diff --git a/pyLiSW/__pycache__/Sample.cpython-312.pyc b/pyLiSW/__pycache__/Sample.cpython-312.pyc
new file mode 100644
index 0000000..6823344
Binary files /dev/null and b/pyLiSW/__pycache__/Sample.cpython-312.pyc differ
diff --git a/pyLiSW/__pycache__/utils.cpython-312.pyc b/pyLiSW/__pycache__/utils.cpython-312.pyc
new file mode 100644
index 0000000..ee31715
Binary files /dev/null and b/pyLiSW/__pycache__/utils.cpython-312.pyc differ
diff --git a/requirements.txt b/requirements.txt
index 6ac4317..e0f82a8 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,25 +1,25 @@
-asteval==0.9.29
-black==22.3.0
-click==8.1.2
-cycler==0.11.0
-fonttools==4.33.3
-future==0.18.3
-kiwisolver==1.4.2
-lmfit==1.2.0
-matplotlib==3.5.1
-mpmath==1.2.1
-mypy-extensions==0.4.3
-numpy==1.22.3
-packaging==21.3
-pathspec==0.9.0
-periodictable==1.6.1
-Pillow==9.1.0
-platformdirs==2.5.2
-pyparsing==3.0.8
-python-dateutil==2.8.2
-scipy==1.8.0
-six==1.16.0
-sympy==1.10.1
-tabulate==0.9.0
-tomli==2.0.1
-uncertainties==3.1.7
+asteval>=0.9.29
+black>=22.3.0
+click>=8.1.2
+cycler>=0.11.0
+fonttools>=4.33.3
+future>=0.18.3
+kiwisolver>=1.4.2
+lmfit>=1.2.0
+matplotlib>=3.5.1
+mpmath>=1.2.1
+mypy-extensions>=0.4.3
+numpy>=1.22.3
+packaging>=21.3
+pathspec>=0.9.0
+periodictable>=1.6.1
+Pillow>=9.1.0
+platformdirs>=2.5.2
+pyparsing>=3.0.8
+python-dateutil>=2.8.2
+scipy>=1.8.0
+six>=1.16.0
+sympy>=1.10.1
+tabulate>=0.9.0
+tomli>=2.0.1
+uncertainties>=3.1.7
diff --git a/tutorials/17_YCrO3.py b/tutorials/17_YCrO3.py
new file mode 100644
index 0000000..a6f3dd2
--- /dev/null
+++ b/tutorials/17_YCrO3.py
@@ -0,0 +1,175 @@
+# -*- coding: utf-8 -*-
+import sys
+from pathlib import Path
+
+sys.path.append(str(Path(__file__).parent.parent) + "/pyLiSW")
+from Atoms import Atoms
+from Bonds import Bonds
+from Sample import Sample
+from LSWT import LSWT
+from utils import gamma_fnc
+from matplotlib import pyplot as plt
+import numpy as np
+
+# Information about YCrO3 can be found in the following reference
+# https://iopscience.iop.org/article/10.1088/1361-648X/abd781
+
+
+if __name__ == "__main__":
+    # lattice parameters in Angstrom
+    # a = 3
+    # # determin the effective lattice parameters
+    lat_params = [5.5225, 7.5474, 5.2521]
+    # propagation vector
+    tau = (0, 0, 0)
+    # vector perpendicular to the plane of rotation
+    n = (1, 0, 0)
+    # temperature
+    te = 0
+    ycro3 = Sample(lat_params, tau, n, te, gamma_fnc=gamma_fnc)
+
+    # -------------------------------------------------------------
+    # Add atoms
+    # -------------------------------------------------------------
+    s1 = 3 / 2
+    aniso = [[0, 0, 0], [0, 0, 0], [0, 0, -0.01]]  # force moments along c
+    # atom postions with effective lattice parameters
+    atoms = [
+        Atoms(
+            t=(0.5, 0.5, 0.0),
+            ion="Cr3",
+            spin=s1,
+            aniso=aniso,
+            theta=0,
+            n_Rp=(1, 0, 0),
+        ),  # moments along positive c-axis
+        Atoms(
+            t=(0.0, 0.5, 0.5),
+            ion="Cr3",
+            spin=s1,
+            aniso=aniso,
+            theta=np.pi,
+            n_Rp=(1, 0, 0),
+        ),  # moments along negative c-axis, rotation needed
+        Atoms(
+            t=(0.0, 0.0, 0.5),
+            ion="Cr3",
+            spin=s1,
+            aniso=aniso,
+            theta=0,
+            n_Rp=(1, 0, 0),
+        ),  # moments along positive c-axis
+        Atoms(
+            t=(0.5, 0.0, 0.0),
+            ion="Cr3",
+            spin=s1,
+            aniso=aniso,
+            theta=np.pi,
+            n_Rp=(1, 0, 0),
+        ),  # moments along negative c-axis, rotation needed
+    ]
+
+    ycro3.add_atoms(atoms)
+
+    # -------------------------------------------------------------
+    # Add bonds
+    # -------------------------------------------------------------
+    j1 = 2.725  # meV  # AFM is positive, FM is negative
+    # my_bonds = np.array([[1, 2], [0, 3], [1, 0], [2, 3]])
+    j2 = 0.1  # meV
+    # my_bonds = np.array([[0, 0], [1, 1], [2, 2], [3, 3], [3, 1], [0, 2]])
+    bonds = [
+        # j1_1 bonds
+        Bonds(ycro3, idx0=1, idx1=2, r0=(0, 0, 0), r1=(0, 0, 0), j=j1),
+        Bonds(ycro3, idx0=1, idx1=2, r0=(0, 0, 0), r1=(0, 1, 0), j=j1),
+        Bonds(ycro3, idx0=0, idx1=3, r0=(0, 0, 0), r1=(0, 0, 0), j=j1),
+        Bonds(ycro3, idx0=0, idx1=3, r0=(0, 0, 0), r1=(0, 1, 0), j=j1),
+        # j1_2 bonds
+        Bonds(ycro3, idx0=1, idx1=0, r0=(0, 0, 0), r1=(0, 0, 0), j=j1),
+        Bonds(ycro3, idx0=1, idx1=0, r0=(0, 0, 0), r1=(0, 0, 1), j=j1),
+        Bonds(ycro3, idx0=0, idx1=1, r0=(0, 0, 0), r1=(1, 0, 0), j=j1),
+        Bonds(ycro3, idx0=0, idx1=1, r0=(0, 0, 1), r1=(1, 0, 0), j=j1),
+        Bonds(ycro3, idx0=2, idx1=3, r0=(0, 0, 0), r1=(0, 0, 0), j=j1),
+        Bonds(ycro3, idx0=2, idx1=3, r0=(0, 0, 0), r1=(0, 0, 1), j=j1),
+        Bonds(ycro3, idx0=2, idx1=3, r0=(1, 0, 0), r1=(0, 0, 0), j=j1),
+        Bonds(ycro3, idx0=2, idx1=3, r0=(1, 0, 0), r1=(0, 0, 1), j=j1),
+        # j2_1 bonds
+        Bonds(ycro3, idx0=0, idx1=0, r0=(0, 0, 0), r1=(0, 0, 1), j=j2),
+        Bonds(ycro3, idx0=0, idx1=0, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
+        Bonds(ycro3, idx0=1, idx1=1, r0=(0, 0, 0), r1=(0, 0, 1), j=j2),
+        Bonds(ycro3, idx0=1, idx1=1, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
+        Bonds(ycro3, idx0=2, idx1=2, r0=(0, 0, 0), r1=(0, 0, 1), j=j2),
+        Bonds(ycro3, idx0=2, idx1=2, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
+        Bonds(ycro3, idx0=3, idx1=3, r0=(0, 0, 0), r1=(0, 0, 1), j=j2),
+        Bonds(ycro3, idx0=3, idx1=3, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
+        # j2_2 bonds
+        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 0, 0), r1=(0, 0, 0), j=j2),
+        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 1, 0), r1=(0, 0, 0), j=j2),
+        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 0, 1), r1=(0, 0, 0), j=j2),
+        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 1, 1), r1=(0, 0, 0), j=j2),
+        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
+        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 1, 0), r1=(1, 0, 0), j=j2),
+        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 0, 1), r1=(1, 0, 0), j=j2),
+        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 1, 1), r1=(1, 0, 0), j=j2),
+        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 0), r1=(0, 0, 0), j=j2),
+        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 1), r1=(0, 0, 0), j=j2),
+        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
+        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 1), r1=(1, 0, 0), j=j2),
+        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 0), r1=(0, 1, 0), j=j2),
+        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 1), r1=(0, 1, 0), j=j2),
+        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 0), r1=(1, 1, 0), j=j2),
+        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 1), r1=(1, 1, 0), j=j2),
+    ]
+    ycro3.add_bonds(bonds)
+
+    # -------------------------------------------------------------
+    # Simulate dispersion along (H00)
+    # -------------------------------------------------------------
+    proj = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
+    axes = ["(H,0,0)", "(0,K,0)", "(0,0,L)"]
+    qe_range = (
+        [0.00, 2.01, 0.01],
+        [0.00, 0.01, 0.01],
+        [0.00, 0.01, 0.01],
+        [-30, 30, 0.01],
+    )
+    sim_qespace = LSWT(qe_range, ycro3, proj_axes=proj, axes=axes)
+    sim_qespace.dispersion_calc()
+    sim_qespace.plot_disp("x")
+
+    slice_range = (
+        [0, 2.01, 0.01],
+        [0.00, 0.01],
+        [0.00, 0.01],
+        [-30, 30.01, 0.01],
+    )
+    sim_qespace.inten_calc()
+    sim_qespace.slice(slice_range, plot_axes=(0, 3), vmin=0, vmax=5)
+
+    # -------------------------------------------------------------
+    # Simulate dispersion along (HH0)
+    # -------------------------------------------------------------
+    proj = [[1, 1, 0], [-1, 1, 0], [0, 0, 1]]
+    axes = ["(H,H,0)", "(-K,K,0)", "(0,0,L)"]
+    qe_range_2 = (
+        [0.00, 2.01, 0.01],
+        [0.00, 0.01, 0.01],
+        [0.00, 0.01, 0.01],
+        [-30, 30, 0.01],
+    )
+    # Qlab  = {'[001]'  '\Gamma'   '[010]' '\Gamma' '[100]' };
+
+    sim_qespace_2 = LSWT(qe_range_2, ycro3, proj_axes=proj, axes=axes)
+    sim_qespace_2.dispersion_calc()
+    sim_qespace_2.plot_disp("x")
+
+    slice_range_2 = (
+        [0, 2.01, 0.01],
+        [0.00, 0.01],
+        [0.00, 0.01],
+        [-30, 30.01, 0.01],
+    )
+    sim_qespace_2.inten_calc()
+    sim_qespace_2.slice(slice_range_2, plot_axes=(0, 3), vmin=0, vmax=5)
+
+    plt.show()