diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 2736a60..ea5f814 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,35 +1,17 @@
-cmake_minimum_required(VERSION 3.18)
 
-set(CMAKE_CXX_STANDARD 17)
+add_subdirectory(times)
+add_subdirectory(motions)
+add_subdirectory(utils)
+
+add_library(simulation STATIC sims.cpp sims.h)
+target_link_libraries(simulation PRIVATE utils)
 
 
-add_executable(rwsim main.cpp
-        utils/functions.h
-        utils/functions.cpp
-        utils/io.cpp
-        utils/io.h
-        motions/base.cpp
-        motions/base.h
-        motions/random.cpp
-        motions/random.h
-        times/base.cpp
-        times/base.h
-        times/delta.cpp
-        times/delta.h
-        simulation/sims.cpp
-        simulation/sims.h
-        utils/ranges.cpp
-        utils/ranges.h
-        motions/tetrahedral.cpp
-        motions/tetrahedral.h
-        motions/isosmallangle.cpp
-        motions/isosmallangle.h
-        motions/coordinates.cpp
-        motions/coordinates.h
-        motions/bimodalangle.cpp
-        motions/bimodalangle.h
-        times/lognormal.cpp
-        times/lognormal.h
+add_executable(
+        rwsim
+        main.cpp
 )
 
+target_link_libraries(rwsim PUBLIC RWMotion RWTime utils simulation)
+
 target_compile_options(rwsim PUBLIC -Werror -Wall -Wextra -Wconversion -O2)
diff --git a/src/main.cpp b/src/main.cpp
index 5f782cb..ae699a0 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,10 +1,9 @@
 
+#include "sims.h"
 #include "utils/io.h"
-#include "simulation/sims.h"
 #include "motions/base.h"
 #include "times/base.h"
 
-
 #include <iostream>
 #include <unordered_map>
 #include <random>
@@ -37,8 +36,8 @@ int main (const int argc, char *argv[]) {
     std::random_device rd;
     std::mt19937_64 rng(rd());
 
-    Motion *motion = Motion::createFromInput(args.motion_type, rng);
-    Distribution *dist = Distribution::createFromInput(args.distribution_type, rng);
+    motions::BaseMotion *motion = motions::BaseMotion::createFromInput(args.motion_type, rng);
+    times::BaseDistribution *dist = times::BaseDistribution::createFromInput(args.distribution_type, rng);
     if (args.spectrum) {
        run_spectrum(parameter, args.optional, *motion, *dist);
     }
diff --git a/src/motions/CMakeLists.txt b/src/motions/CMakeLists.txt
new file mode 100644
index 0000000..76b5dad
--- /dev/null
+++ b/src/motions/CMakeLists.txt
@@ -0,0 +1,23 @@
+# Create a library target for motions
+add_library(
+        RWMotion STATIC
+        conewobble.cpp
+        conewobble.h
+        coordinates.cpp
+        coordinates.h
+        base.cpp
+        base.h
+        random.cpp
+        random.h
+        isosmallangle.cpp
+        isosmallangle.h
+        foursitejump.cpp
+        foursitejump.h
+        rjoac.cpp
+        rjoac.h
+        bimodalangle.cpp
+        bimodalangle.h
+        sixsitejump.cpp
+        sixsitejump.h
+)
+
diff --git a/src/motions/base.cpp b/src/motions/base.cpp
index 992d6f1..5a34531 100644
--- a/src/motions/base.cpp
+++ b/src/motions/base.cpp
@@ -4,68 +4,74 @@
 #include "bimodalangle.h"
 #include "isosmallangle.h"
 #include "random.h"
-#include "tetrahedral.h"
+#include "foursitejump.h"
 
 #include <stdexcept>
-#include <map>
 
-Motion::Motion(std::string name, const double delta, const double eta, std::mt19937_64& rng) : m_name(std::move(name)), m_delta(delta), m_eta(eta), m_rng(rng) {
-    m_uni_dist = std::uniform_real_distribution(0., 1.);
-}
+#include "sixsitejump.h"
 
-Motion::Motion(std::string name, std::mt19937_64& rng) : m_name(std::move(name)), m_rng(rng) {
-    m_uni_dist = std::uniform_real_distribution(0., 1.);
-}
+namespace motions {
+    BaseMotion::BaseMotion(std::string name, const double delta, const double eta, std::mt19937_64& rng) : m_name(std::move(name)), m_delta(delta), m_eta(eta), m_rng(rng) {
+        m_uni_dist = std::uniform_real_distribution(0., 1.);
+    }
 
-double Motion::omega_q(const double cos_theta, const double phi) const {
-    const double cos_theta_square = cos_theta * cos_theta;
-    const double sin_theta_square = 1. - cos_theta_square;
+    BaseMotion::BaseMotion(std::string name, std::mt19937_64& rng) : m_name(std::move(name)), m_rng(rng) {
+        m_uni_dist = std::uniform_real_distribution(0., 1.);
+    }
 
-    return M_PI * m_delta * (3. * cos_theta_square - 1. - m_eta * sin_theta_square * std::cos(2.*phi));
-}
+    double BaseMotion::omega_q(const double cos_theta, const double phi) const {
+        const double cos_theta_square = cos_theta * cos_theta;
+        const double sin_theta_square = 1. - cos_theta_square;
 
-double Motion::omega_q(const SphericalPos& pos) const {
-    auto [cos_theta, phi] = pos;
+        return M_PI * m_delta * (3. * cos_theta_square - 1. - m_eta * sin_theta_square * std::cos(2.*phi));
+    }
 
-    return omega_q(cos_theta, phi);
-}
+    double BaseMotion::omega_q(const coordinates::SphericalPos& pos) const {
+        auto [cos_theta, phi] = pos;
 
-SphericalPos Motion::draw_position() {
-    const double cos_theta = 1 - 2 * m_uni_dist(m_rng);
-    const double phi = 2.0 * M_PI * m_uni_dist(m_rng);
+        return omega_q(cos_theta, phi);
+    }
 
-    return {cos_theta, phi};
-}
+    coordinates::SphericalPos BaseMotion::draw_position() {
+        const double cos_theta = 1 - 2 * m_uni_dist(m_rng);
+        const double phi = 2.0 * M_PI * m_uni_dist(m_rng);
 
-Motion* Motion::createFromInput(const std::string& input, std::mt19937_64& rng) {
-    if (input == "TetrahedralJump")
-        return new TetrahedralJump(rng);
+        return {cos_theta, phi};
+    }
 
-    if (input == "IsotropicAngle")
-        return new SmallAngle(rng);
+    BaseMotion* BaseMotion::createFromInput(const std::string& input, std::mt19937_64& rng) {
+        if (input == "FourSiteTetrahedral")
+            return new FourSiteTetrahedron(rng);
 
-    if (input == "RandomJump")
-        return new RandomJump(rng);
+        if (input == "SixSiteOctahedral")
+            return new SixSiteOctahedron(rng);
 
-    if (input == "BimodalAngle")
-        return new BimodalAngle(rng);
+        if (input == "IsotropicAngle")
+            return new SmallAngle(rng);
 
-    throw std::invalid_argument("Invalid input " + input);
-}
+        if (input == "RandomJump")
+            return new RandomJump(rng);
 
-void Motion::setParameters(const std::unordered_map<std::string, double> &parameters) {
-    m_delta = parameters.at("delta");
-    m_eta = parameters.at("eta");
-}
+        if (input == "BimodalAngle")
+            return new BimodalAngle(rng);
 
-std::unordered_map<std::string, double> Motion::getParameters() const {
-    return std::unordered_map<std::string, double>{
-        {"delta", m_delta},
-        {"eta", m_eta}
-    };
-}
+        throw std::invalid_argument("Invalid input " + input);
+    }
 
-std::ostream& operator<<(std::ostream& os, const Motion& m) {
-    os << m.getName();
-    return os;
+    void BaseMotion::setParameters(const std::unordered_map<std::string, double> &parameters) {
+        m_delta = parameters.at("delta");
+        m_eta = parameters.at("eta");
+    }
+
+    std::unordered_map<std::string, double> BaseMotion::getParameters() const {
+        return std::unordered_map<std::string, double>{
+            {"delta", m_delta},
+            {"eta", m_eta}
+        };
+    }
+
+    std::ostream& operator<<(std::ostream& os, const BaseMotion& m) {
+        os << m.getName();
+        return os;
+    }
 }
\ No newline at end of file
diff --git a/src/motions/base.h b/src/motions/base.h
index fc5cc1a..8022b3d 100644
--- a/src/motions/base.h
+++ b/src/motions/base.h
@@ -6,42 +6,43 @@
 #include <random>
 #include <unordered_map>
 
-class Motion {
-public:
-    virtual ~Motion() = default;
+namespace motions {
+    class BaseMotion {
+    public:
+        virtual ~BaseMotion() = default;
 
-    Motion(std::string, double, double, std::mt19937_64&);
-    explicit Motion(std::string, std::mt19937_64&);
+        BaseMotion(std::string, double, double, std::mt19937_64&);
+        explicit BaseMotion(std::string, std::mt19937_64&);
 
-    SphericalPos draw_position();
-    [[nodiscard]] double omega_q(double, double) const;
-    [[nodiscard]] double omega_q(const SphericalPos&) const;
+        coordinates::SphericalPos draw_position();
+        [[nodiscard]] double omega_q(double, double) const;
+        [[nodiscard]] double omega_q(const coordinates::SphericalPos&) const;
 
-    virtual void initialize() = 0;
-    virtual double jump() = 0;
+        virtual void initialize() = 0;
+        virtual double jump() = 0;
 
-    virtual void setParameters(const std::unordered_map<std::string, double>&);
-    [[nodiscard]] virtual std::unordered_map<std::string, double> getParameters() const;
-    [[nodiscard]] double getDelta() const { return m_delta; }
-    void setDelta(const double delta) { m_delta = delta; }
-    [[nodiscard]] double getEta() const { return m_eta; }
-    void setEta(const double eta) { m_eta = eta; }
-    [[nodiscard]] std::string getName() const { return m_name; }
-    [[nodiscard]] double getInitOmega() const { return m_initial_omega; };
+        virtual void setParameters(const std::unordered_map<std::string, double>&);
+        [[nodiscard]] virtual std::unordered_map<std::string, double> getParameters() const;
+        [[nodiscard]] double getDelta() const { return m_delta; }
+        void setDelta(const double delta) { m_delta = delta; }
+        [[nodiscard]] double getEta() const { return m_eta; }
+        void setEta(const double eta) { m_eta = eta; }
+        [[nodiscard]] std::string getName() const { return m_name; }
+        [[nodiscard]] double getInitOmega() const { return m_initial_omega; };
 
-    [[nodiscard]] virtual std::string toString() const = 0;
+        [[nodiscard]] virtual std::string toString() const = 0;
 
-    static Motion* createFromInput(const std::string& input, std::mt19937_64& rng);
+        static BaseMotion* createFromInput(const std::string& input, std::mt19937_64& rng);
 
-protected:
-    std::string m_name{"BaseMotion"};
-    double m_delta{1.};
-    double m_eta{0.};
-    std::mt19937_64& m_rng;
-    std::uniform_real_distribution<> m_uni_dist;
-    double m_initial_omega{0.};
-};
-
-std::ostream& operator<<(std::ostream& os, const Motion& m);
+    protected:
+        std::string m_name{"BaseMotion"};
+        double m_delta{1.};
+        double m_eta{0.};
+        std::mt19937_64& m_rng;
+        std::uniform_real_distribution<> m_uni_dist;
+        double m_initial_omega{0.};
+    };
 
+    std::ostream& operator<<(std::ostream& os, const BaseMotion& m);
+}
 #endif //RWSIM_MOTIONBASE_H
diff --git a/src/motions/bimodalangle.cpp b/src/motions/bimodalangle.cpp
index f0f55b7..6712104 100644
--- a/src/motions/bimodalangle.cpp
+++ b/src/motions/bimodalangle.cpp
@@ -1,43 +1,46 @@
 
 #include "bimodalangle.h"
 #include "base.h"
+#include "coordinates.h"
 
-BimodalAngle::BimodalAngle(const double delta, const double eta, const double angle1, const double angle2, const double prob, std::mt19937_64 &rng) :
-    Motion(std::string("BimodalAngle"), delta, eta, rng),
-    m_angle1(angle1 * M_PI / 180.0),
-    m_angle2(angle2 * M_PI / 180.0),
-    m_prob(prob) {};
-BimodalAngle::BimodalAngle(std::mt19937_64 &rng) : Motion(std::string("BimodalAngle"), rng) {}
+namespace motions {
+    BimodalAngle::BimodalAngle(const double delta, const double eta, const double angle1, const double angle2, const double prob, std::mt19937_64 &rng) :
+        BaseMotion(std::string("BimodalAngle"), delta, eta, rng),
+        m_angle1(angle1 * M_PI / 180.0),
+        m_angle2(angle2 * M_PI / 180.0),
+        m_prob(prob) {};
+    BimodalAngle::BimodalAngle(std::mt19937_64 &rng) : BaseMotion(std::string("BimodalAngle"), rng) {}
 
-void BimodalAngle::initialize() {
-    m_prev_pos = draw_position();
-    m_initial_omega = omega_q(m_prev_pos);
-};
+    void BimodalAngle::initialize() {
+        m_prev_pos = draw_position();
+        m_initial_omega = omega_q(m_prev_pos);
+    };
 
-double BimodalAngle::jump() {
-    const double angle = m_uni_dist(m_rng) < m_prob ? m_angle1 : m_angle2;
-    const double gamma{2 * M_PI * m_uni_dist(m_rng)};
-    m_prev_pos = rotate(m_prev_pos, angle, gamma);
+    double BimodalAngle::jump() {
+        const double angle = m_uni_dist(m_rng) < m_prob ? m_angle1 : m_angle2;
+        const double gamma{2 * M_PI * m_uni_dist(m_rng)};
+        m_prev_pos = coordinates::rotate(m_prev_pos, angle, gamma);
 
-    return omega_q(m_prev_pos);
-}
-
-void BimodalAngle::setParameters(const std::unordered_map<std::string, double> &parameter) {
-    Motion::setParameters(parameter);
-
-    m_angle1 = parameter.at("angle1") * M_PI / 180.;
-    m_angle2 = parameter.at("angle2") * M_PI / 180.;
-    m_prob = parameter.at("probability1");
-}
-
-std::unordered_map<std::string, double> BimodalAngle::getParameters() const {
-    auto parameter = Motion::getParameters();
-    parameter["angle1"] = m_angle1 * 180 / M_PI;
-    parameter["angle2"] = m_angle2 * 180 / M_PI;
-    parameter["probality1"] = m_prob;
-    return parameter;
-}
-
-std::string BimodalAngle::toString() const {
-    return std::string{"BimodalAngle/angle1=" + std::to_string(m_angle1 * 180 / M_PI) + "/angle2=" + std::to_string(m_angle2 * 180 / M_PI) + "/probability1=" + std::to_string(m_prob)};
+        return omega_q(m_prev_pos);
+    }
+
+    void BimodalAngle::setParameters(const std::unordered_map<std::string, double> &parameter) {
+        BaseMotion::setParameters(parameter);
+
+        m_angle1 = parameter.at("angle1") * M_PI / 180.;
+        m_angle2 = parameter.at("angle2") * M_PI / 180.;
+        m_prob = parameter.at("probability1");
+    }
+
+    std::unordered_map<std::string, double> BimodalAngle::getParameters() const {
+        auto parameter = BaseMotion::getParameters();
+        parameter["angle1"] = m_angle1 * 180 / M_PI;
+        parameter["angle2"] = m_angle2 * 180 / M_PI;
+        parameter["probality1"] = m_prob;
+        return parameter;
+    }
+
+    std::string BimodalAngle::toString() const {
+        return std::string{"BimodalAngle/angle1=" + std::to_string(m_angle1 * 180 / M_PI) + "/angle2=" + std::to_string(m_angle2 * 180 / M_PI) + "/probability1=" + std::to_string(m_prob)};
+    }
 }
diff --git a/src/motions/bimodalangle.h b/src/motions/bimodalangle.h
index 8df624a..ac148ce 100644
--- a/src/motions/bimodalangle.h
+++ b/src/motions/bimodalangle.h
@@ -3,23 +3,25 @@
 #define BIMODALANGLE_H
 
 #include "base.h"
+#include "coordinates.h"
 
-class BimodalAngle : public Motion {
-public:
-    BimodalAngle(double, double, double, double, double, std::mt19937_64& );
-    explicit BimodalAngle(std::mt19937_64&);
+namespace motions {
+    class BimodalAngle : public BaseMotion {
+    public:
+        BimodalAngle(double, double, double, double, double, std::mt19937_64& );
+        explicit BimodalAngle(std::mt19937_64&);
 
-    void initialize() override;
-    double jump() override;
-    void setParameters(const std::unordered_map<std::string, double> &) override;
-    [[nodiscard]] std::unordered_map<std::string, double> getParameters() const override;
-    [[nodiscard]] std::string toString() const override;
-
-protected:
-    double m_angle1{0};
-    double m_angle2{0};
-    double m_prob{0};
-    SphericalPos m_prev_pos{0., 0.};
-};
+        void initialize() override;
+        double jump() override;
+        void setParameters(const std::unordered_map<std::string, double> &) override;
+        [[nodiscard]] std::unordered_map<std::string, double> getParameters() const override;
+        [[nodiscard]] std::string toString() const override;
 
+    protected:
+        double m_angle1{0};
+        double m_angle2{0};
+        double m_prob{0};
+        coordinates::SphericalPos m_prev_pos{0., 0.};
+    };
+}
 #endif //BIMODALANGLE_H
diff --git a/src/motions/coordinates.cpp b/src/motions/coordinates.cpp
index f7dbdf4..fe4b6c6 100644
--- a/src/motions/coordinates.cpp
+++ b/src/motions/coordinates.cpp
@@ -3,37 +3,39 @@
 #include <cmath>
 #include <iostream>
 
-SphericalPos rotate(const SphericalPos& old_pos, const double alpha, const double beta) {
-    const double sin_alpha{std::sin(alpha)};
-    const double cos_alpha{std::cos(alpha)};
+namespace coordinates {
+    SphericalPos rotate(const SphericalPos& old_pos, const double alpha, const double beta) {
+        const double sin_alpha{std::sin(alpha)};
+        const double cos_alpha{std::cos(alpha)};
 
-    const double sin_beta{std::sin(beta)};
-    const double cos_beta{std::cos(beta)};
+        const double sin_beta{std::sin(beta)};
+        const double cos_beta{std::cos(beta)};
 
-    const double cos_theta{old_pos.cos_theta};
+        const double cos_theta{old_pos.cos_theta};
 
-    if (std::abs(cos_theta) == 1) {
-        return xyz_to_spherical(CartesianPos{cos_alpha * cos_beta, cos_alpha * sin_beta, cos_alpha * cos_theta});
+        if (std::abs(cos_theta) == 1) {
+            return xyz_to_spherical(CartesianPos{cos_alpha * cos_beta, cos_alpha * sin_beta, cos_alpha * cos_theta});
+        }
+
+        const double norm{std::sqrt(1 - cos_theta * cos_theta) + 1e-15};
+
+        auto [x, y , z] = spherical_to_xyz(old_pos);
+
+        const auto new_pos = CartesianPos{
+            cos_alpha * x + sin_alpha * (-x * z * sin_beta - y * cos_beta) / norm,
+            cos_alpha * y + sin_alpha * (-y * z * sin_beta + x * cos_beta) / norm,
+            cos_alpha * z + sin_alpha * norm * sin_beta
+        };
+
+        return xyz_to_spherical(new_pos);
     }
 
-    const double norm{std::sqrt(1 - cos_theta * cos_theta) + 1e-15};
+    CartesianPos spherical_to_xyz(const SphericalPos& pos) {
+        const double sin_theta = std::sin(std::acos(pos.cos_theta));
+        return {sin_theta * std::cos(pos.phi), sin_theta * std::sin(pos.phi), pos.cos_theta};
+    }
 
-    auto [x, y , z] = spherical_to_xyz(old_pos);
-
-    const auto new_pos = CartesianPos{
-        cos_alpha * x + sin_alpha * (-x * z * sin_beta - y * cos_beta) / norm,
-        cos_alpha * y + sin_alpha * (-y * z * sin_beta + x * cos_beta) / norm,
-        cos_alpha * z + sin_alpha * norm * sin_beta
-    };
-
-    return xyz_to_spherical(new_pos);
-}
-
-CartesianPos spherical_to_xyz(const SphericalPos& pos) {
-    const double sin_theta = std::sin(std::acos(pos.cos_theta));
-    return {sin_theta * std::cos(pos.phi), sin_theta * std::sin(pos.phi), pos.cos_theta};
-}
-
-SphericalPos xyz_to_spherical(const CartesianPos& pos) {
-    return {pos.z, std::atan2(pos.y, pos.x)};
-}
+    SphericalPos xyz_to_spherical(const CartesianPos& pos) {
+        return {pos.z, std::atan2(pos.y, pos.x)};
+    }
+}
\ No newline at end of file
diff --git a/src/motions/coordinates.h b/src/motions/coordinates.h
index 56b64de..3e4fd2f 100644
--- a/src/motions/coordinates.h
+++ b/src/motions/coordinates.h
@@ -2,19 +2,21 @@
 #ifndef COORDINATES_H
 #define COORDINATES_H
 
-struct SphericalPos {
-    double cos_theta;
-    double phi;
-};
+namespace coordinates {
+    struct SphericalPos {
+        double cos_theta;
+        double phi;
+    };
 
-struct CartesianPos {
-    double x;
-    double y;
-    double z;
-};
+    struct CartesianPos {
+        double x;
+        double y;
+        double z;
+    };
 
-SphericalPos rotate(const SphericalPos&, double, double);
-CartesianPos spherical_to_xyz(const SphericalPos&);
-SphericalPos xyz_to_spherical(const CartesianPos&);
+    SphericalPos rotate(const SphericalPos&, double, double);
+    CartesianPos spherical_to_xyz(const SphericalPos&);
+    SphericalPos xyz_to_spherical(const CartesianPos&);
+}
 
 #endif //COORDINATES_H
diff --git a/src/motions/foursitejump.cpp b/src/motions/foursitejump.cpp
new file mode 100644
index 0000000..b35c374
--- /dev/null
+++ b/src/motions/foursitejump.cpp
@@ -0,0 +1,35 @@
+#include "foursitejump.h"
+#include "coordinates.h"
+
+#include <random>
+
+namespace motions {
+    FourSiteTetrahedron::FourSiteTetrahedron(const double delta, const double eta, std::mt19937_64& rng) :
+        BaseMotion(std::string{"FourSiteTetrahedral"}, delta, eta, rng) {}
+
+    FourSiteTetrahedron::FourSiteTetrahedron(std::mt19937_64& rng) : BaseMotion(std::string{"FourSiteTetrahedral"}, rng) {}
+
+    void FourSiteTetrahedron::initialize() {
+        const auto pos = draw_position();
+        m_corners[0] = omega_q(pos);
+
+        const double alpha = 2. * M_PI * m_uni_dist(m_rng);
+
+        for (int i = 1; i<4; i++) {
+            auto corner_pos = coordinates::rotate(pos, m_beta, alpha + (i-1) * 2*M_PI/3.);
+            m_corners[i] = omega_q(corner_pos);
+        }
+        m_initial_omega = FourSiteTetrahedron::jump();
+    }
+
+    double FourSiteTetrahedron::jump() {
+        m_corner_idx += m_chooser(m_rng);
+        m_corner_idx %= 4;
+
+        return m_corners[m_corner_idx];
+    }
+
+    std::string FourSiteTetrahedron::toString() const {
+        return {"FourSiteTetrahedral"};
+    }
+}
\ No newline at end of file
diff --git a/src/motions/foursitejump.h b/src/motions/foursitejump.h
new file mode 100644
index 0000000..7ae8538
--- /dev/null
+++ b/src/motions/foursitejump.h
@@ -0,0 +1,31 @@
+#ifndef RWSIM_MOTIONTETRAHEDRAL_H
+#define RWSIM_MOTIONTETRAHEDRAL_H
+
+#include "base.h"
+
+#include <random>
+#include <cmath>
+#include <array>
+
+namespace motions {
+    class FourSiteTetrahedron final : public BaseMotion {
+    public:
+        FourSiteTetrahedron(double, double, std::mt19937_64&);
+        explicit FourSiteTetrahedron(std::mt19937_64&);
+
+        void initialize() override;
+        double jump() override;
+
+        [[nodiscard]] std::string toString() const override;
+
+    private:
+        const double m_beta{std::acos(-1/3.)};
+
+        std::array<double, 4> m_corners{};
+        int m_corner_idx{0};
+
+        std::uniform_int_distribution<> m_chooser{1, 3};
+
+    };
+}
+#endif //RWSIM_MOTIONTETRAHEDRAL_H
diff --git a/src/motions/isosmallangle.cpp b/src/motions/isosmallangle.cpp
index 8fb815f..58f07b0 100644
--- a/src/motions/isosmallangle.cpp
+++ b/src/motions/isosmallangle.cpp
@@ -4,34 +4,35 @@
 #include <iostream>
 
 
+namespace motions {
+    SmallAngle::SmallAngle(const double delta, const double eta, const double chi, std::mt19937_64 &rng) :
+        BaseMotion(std::string("IsotropicAngle"), delta, eta, rng), m_chi(chi * M_PI / 180.0) {};
+    SmallAngle::SmallAngle(std::mt19937_64 &rng) : BaseMotion(std::string("IsotropicAngle"), rng) {}
 
-SmallAngle::SmallAngle(const double delta, const double eta, const double chi, std::mt19937_64 &rng) :
-    Motion(std::string("IsotropicAngle"), delta, eta, rng), m_chi(chi * M_PI / 180.0) {};
-SmallAngle::SmallAngle(std::mt19937_64 &rng) : Motion(std::string("IsotropicAngle"), rng) {}
+    void SmallAngle::initialize() {
+        m_prev_pos = draw_position();
+        m_initial_omega = omega_q(m_prev_pos);
+    };
 
-void SmallAngle::initialize() {
-    m_prev_pos = draw_position();
-    m_initial_omega = omega_q(m_prev_pos);
-};
+    double SmallAngle::jump() {
+        const double gamma{2 * M_PI * m_uni_dist(m_rng)};
+        m_prev_pos = coordinates::rotate(m_prev_pos, m_chi, gamma);
 
-double SmallAngle::jump() {
-    const double gamma{2 * M_PI * m_uni_dist(m_rng)};
-    m_prev_pos = rotate(m_prev_pos, m_chi, gamma);
+        return omega_q(m_prev_pos);
+    }
 
-    return omega_q(m_prev_pos);
-}
-
-void SmallAngle::setParameters(const std::unordered_map<std::string, double> &parameters) {
-    m_chi = parameters.at("angle") * M_PI / 180.0;
-    Motion::setParameters(parameters);
-}
-
-std::unordered_map<std::string, double> SmallAngle::getParameters() const {
-    auto parameter = Motion::getParameters();
-    parameter["angle"] = m_chi * 180 / M_PI;
-    return parameter;
-}
-
-std::string SmallAngle::toString() const {
-    return std::string{"IsotropicAngle/angle=" + std::to_string(m_chi * 180 / M_PI)};
+    void SmallAngle::setParameters(const std::unordered_map<std::string, double> &parameters) {
+        m_chi = parameters.at("angle") * M_PI / 180.0;
+        BaseMotion::setParameters(parameters);
+    }
+
+    std::unordered_map<std::string, double> SmallAngle::getParameters() const {
+        auto parameter = BaseMotion::getParameters();
+        parameter["angle"] = m_chi * 180 / M_PI;
+        return parameter;
+    }
+
+    std::string SmallAngle::toString() const {
+        return std::string{"IsotropicAngle/angle=" + std::to_string(m_chi * 180 / M_PI)};
+    }
 }
diff --git a/src/motions/isosmallangle.h b/src/motions/isosmallangle.h
index 66f805e..0e04f16 100644
--- a/src/motions/isosmallangle.h
+++ b/src/motions/isosmallangle.h
@@ -5,20 +5,21 @@
 #include "base.h"
 #include "coordinates.h"
 
-class SmallAngle final : public Motion {
-public:
-    SmallAngle(double, double, double, std::mt19937_64& );
-    explicit SmallAngle(std::mt19937_64&);
+namespace motions {
+    class SmallAngle final : public BaseMotion {
+    public:
+        SmallAngle(double, double, double, std::mt19937_64& );
+        explicit SmallAngle(std::mt19937_64&);
 
-    void initialize() override;
-    double jump() override;
-    void setParameters(const std::unordered_map<std::string, double> &) override;
-    [[nodiscard]] std::unordered_map<std::string, double> getParameters() const override;
-    [[nodiscard]] std::string toString() const override;
-
-private:
-    double m_chi{0};
-    SphericalPos m_prev_pos{0., 0.};
-};
+        void initialize() override;
+        double jump() override;
+        void setParameters(const std::unordered_map<std::string, double> &) override;
+        [[nodiscard]] std::unordered_map<std::string, double> getParameters() const override;
+        [[nodiscard]] std::string toString() const override;
 
+    private:
+        double m_chi{0};
+        coordinates::SphericalPos m_prev_pos{0., 0.};
+    };
+}
 #endif //RWSIM_MOTIONISOSMALLANGLE_H
diff --git a/src/motions/random.cpp b/src/motions/random.cpp
index 8a2f52f..eb34481 100644
--- a/src/motions/random.cpp
+++ b/src/motions/random.cpp
@@ -1,19 +1,20 @@
 
 #include "random.h"
 
+namespace motions {
+    RandomJump::RandomJump(const double delta, const double eta, std::mt19937_64 &rng) : BaseMotion(std::string("RandomJump"), delta, eta, rng) {}
 
-RandomJump::RandomJump(const double delta, const double eta, std::mt19937_64 &rng) : Motion(std::string("RandomJump"), delta, eta, rng) {}
+    RandomJump::RandomJump(std::mt19937_64 &rng) : BaseMotion(std::string("RandomJump"), rng) {}
 
-RandomJump::RandomJump(std::mt19937_64 &rng) : Motion(std::string("RandomJump"), rng) {}
+    std::string RandomJump::toString() const {
+        return {"RandomJump"};
+    }
 
-std::string RandomJump::toString() const {
-    return {"RandomJump"};
-}
-
-void RandomJump::initialize() {
-    m_initial_omega = RandomJump::jump();
-}
-
-double RandomJump::jump() {
-    return omega_q(draw_position());
+    void RandomJump::initialize() {
+        m_initial_omega = RandomJump::jump();
+    }
+
+    double RandomJump::jump() {
+        return omega_q(draw_position());
+    }
 }
diff --git a/src/motions/random.h b/src/motions/random.h
index 82e4c64..3eb465b 100644
--- a/src/motions/random.h
+++ b/src/motions/random.h
@@ -5,15 +5,17 @@
 #include "base.h"
 #include <random>
 
-class RandomJump final : public Motion {
-public:
-    RandomJump(double, double, std::mt19937_64&);
-    explicit RandomJump(std::mt19937_64&);
+namespace motions {
+    class RandomJump final : public BaseMotion {
+    public:
+        RandomJump(double, double, std::mt19937_64&);
+        explicit RandomJump(std::mt19937_64&);
 
-    [[nodiscard]] std::string toString() const override;
+        [[nodiscard]] std::string toString() const override;
 
-    void initialize() override;
-    double jump() override;
-};
+        void initialize() override;
+        double jump() override;
+    };
+}
 
 #endif //RWSIM_MOTIONRANDOMJUMP_H
diff --git a/src/motions/rjoac.cpp b/src/motions/rjoac.cpp
new file mode 100644
index 0000000..e64b3b3
--- /dev/null
+++ b/src/motions/rjoac.cpp
@@ -0,0 +1,32 @@
+
+#include "rjoac.h"
+
+namespace motions {
+    RandomJumpOnCone::RandomJumpOnCone(const double delta, const double eta, const double chi, std::mt19937_64 &rng) : BaseMotion("RJ on a Cone", delta, eta, rng), m_angle(chi) {}
+    RandomJumpOnCone::RandomJumpOnCone(std::mt19937_64 &rng) : BaseMotion("RJ on a Cone", rng) {}
+
+    void RandomJumpOnCone::initialize() {
+        m_axis = draw_position();
+    }
+
+    double RandomJumpOnCone::jump() {
+        const double phi = 2 * M_PI * m_uni_dist(m_rng);
+        return omega_q(coordinates::rotate(m_axis, m_angle, phi));
+    }
+
+    void RandomJumpOnCone::setParameters(const std::unordered_map<std::string, double> &parameters) {
+        BaseMotion::setParameters(parameters);
+        m_angle = parameters.at("angle");
+    }
+
+    std::unordered_map<std::string, double> RandomJumpOnCone::getParameters() const {
+        auto parameter =  BaseMotion::getParameters();
+        parameter["angle"] = m_angle;
+
+        return parameter;
+    }
+
+    std::string RandomJumpOnCone::toString() const {
+        return std::string("RandomJumpOnCone/angle=") + std::to_string(m_angle);
+    }
+}
\ No newline at end of file
diff --git a/src/motions/rjoac.h b/src/motions/rjoac.h
new file mode 100644
index 0000000..59eca9b
--- /dev/null
+++ b/src/motions/rjoac.h
@@ -0,0 +1,31 @@
+#ifndef RJOAC_H
+#define RJOAC_H
+
+#include "base.h"
+#include "coordinates.h"
+
+#include <random>
+
+namespace motions {
+    class RandomJumpOnCone: public BaseMotion {
+    public:
+        RandomJumpOnCone(double, double, double, std::mt19937_64&);
+        explicit RandomJumpOnCone(std::mt19937_64&);
+
+        void initialize() override;
+
+        double jump() override;
+
+        void setParameters(const std::unordered_map<std::string, double> &) override;
+
+        [[nodiscard]] std::unordered_map<std::string, double> getParameters() const override;
+
+        [[nodiscard]] std::string toString() const override;
+
+    private:
+        double m_angle{0};
+        coordinates::SphericalPos m_axis{1, 0};
+    };
+}
+
+#endif //RJOAC_H
diff --git a/src/motions/sixsitejump.cpp b/src/motions/sixsitejump.cpp
new file mode 100644
index 0000000..abf83c4
--- /dev/null
+++ b/src/motions/sixsitejump.cpp
@@ -0,0 +1,34 @@
+#include "sixsitejump.h"
+#include "coordinates.h"
+
+namespace motions {
+    SixSiteOctahedron::SixSiteOctahedron(const double delta, const double eta, std::mt19937_64& rng) :
+        BaseMotion(std::string{"SixSiteOctahedral"}, delta, eta, rng) {}
+
+    SixSiteOctahedron::SixSiteOctahedron(std::mt19937_64& rng) : BaseMotion(std::string{"SixSiteOctahedral"}, rng) {}
+
+    void SixSiteOctahedron::initialize() {
+        const auto pos = draw_position();
+        m_corners[0] = omega_q(pos);
+        m_corners[1] = omega_q(rotate(pos, M_PI, 0));
+
+        const double alpha = 2. * M_PI * m_uni_dist(m_rng);
+
+        for (int i = 2; i<6; i++) {
+            auto corner_pos = coordinates::rotate(pos, M_PI_2, alpha + (i-2) * M_PI_2);
+            m_corners[i] = omega_q(corner_pos);
+        }
+        m_initial_omega = SixSiteOctahedron::jump();
+    }
+
+    double SixSiteOctahedron::jump() {
+        m_corner_idx += m_chooser(m_rng);
+        m_corner_idx %= 6;
+
+        return m_corners[m_corner_idx];
+    }
+
+    std::string SixSiteOctahedron::toString() const {
+        return {"SixSiteOctahedral"};
+    }
+}
diff --git a/src/motions/sixsitejump.h b/src/motions/sixsitejump.h
new file mode 100644
index 0000000..32fe9d6
--- /dev/null
+++ b/src/motions/sixsitejump.h
@@ -0,0 +1,32 @@
+
+#ifndef SIXSITEJUMP_H
+#define SIXSITEJUMP_H
+
+#include "base.h"
+
+#include <random>
+#include <cmath>
+#include <array>
+
+namespace motions {
+    class SixSiteOctahedron final : public BaseMotion {
+    public:
+        SixSiteOctahedron(double, double, std::mt19937_64&);
+        explicit SixSiteOctahedron(std::mt19937_64&);
+
+        void initialize() override;
+        double jump() override;
+
+        [[nodiscard]] std::string toString() const override;
+
+    private:
+        const double m_beta{std::acos(-1/3.)};
+
+        std::array<double, 6> m_corners{};
+        int m_corner_idx{0};
+
+        std::uniform_int_distribution<> m_chooser{1, 5};
+
+    };
+}
+#endif //SIXSITEJUMP_H
diff --git a/src/motions/tetrahedral.cpp b/src/motions/tetrahedral.cpp
deleted file mode 100644
index 8dfea40..0000000
--- a/src/motions/tetrahedral.cpp
+++ /dev/null
@@ -1,34 +0,0 @@
-#include "tetrahedral.h"
-
-#include <random>
-
-#include "tetrahedral.h"
-
-TetrahedralJump::TetrahedralJump(const double delta, const double eta, std::mt19937_64& rng) :
-    Motion(std::string{"FourSiteTetrahedral"}, delta, eta, rng) {}
-
-TetrahedralJump::TetrahedralJump(std::mt19937_64& rng) : Motion(std::string{"FourSiteTetrahedral"}, rng) {}
-
-void TetrahedralJump::initialize() {
-    const auto pos = draw_position();
-    m_corners[0] = omega_q(pos);
-
-    const double alpha = 2. * M_PI * m_uni_dist(m_rng);
-
-    for (int i = 1; i<4; i++) {
-        auto corner_pos = rotate(pos, m_beta, alpha + (i-1) * 2*M_PI/3.);
-        m_corners[i] = omega_q(corner_pos);
-    }
-    m_initial_omega = TetrahedralJump::jump();
-}
-
-double TetrahedralJump::jump() {
-    m_corner_idx += m_chooser(m_rng);
-    m_corner_idx %= 4;
-
-    return m_corners[m_corner_idx];
-}
-
-std::string TetrahedralJump::toString() const {
-    return {"FourSiteTetrahedral"};
-}
diff --git a/src/motions/tetrahedral.h b/src/motions/tetrahedral.h
deleted file mode 100644
index 8c07c77..0000000
--- a/src/motions/tetrahedral.h
+++ /dev/null
@@ -1,29 +0,0 @@
-#ifndef RWSIM_MOTIONTETRAHEDRAL_H
-#define RWSIM_MOTIONTETRAHEDRAL_H
-
-#include "base.h"
-#include <random>
-#include <cmath>
-#include <array>
-
-class TetrahedralJump final : public Motion {
-public:
-    TetrahedralJump(double, double, std::mt19937_64&);
-    explicit TetrahedralJump(std::mt19937_64&);
-
-    void initialize() override;
-    double jump() override;
-
-    [[nodiscard]] std::string toString() const override;
-
-private:
-    const double m_beta{std::acos(-1/3.)};
-
-    std::array<double, 4> m_corners{};
-    int m_corner_idx{0};
-
-    std::uniform_int_distribution<> m_chooser{1, 3};
-
-};
-
-#endif //RWSIM_MOTIONTETRAHEDRAL_H
diff --git a/src/simulation/sims.cpp b/src/sims.cpp
similarity index 95%
rename from src/simulation/sims.cpp
rename to src/sims.cpp
index 77b3dc6..ec17575 100644
--- a/src/simulation/sims.cpp
+++ b/src/sims.cpp
@@ -1,9 +1,8 @@
-#include "sims.h"
-#include "../motions/base.h"
-#include "../times/base.h"
-#include "../utils/functions.h"
-#include "../utils/ranges.h"
-#include "../utils/io.h"
+#include "simulation/sims.h"
+#include "times/base.h"
+#include "utils/functions.h"
+#include "utils/ranges.h"
+#include "utils/io.h"
 
 #include <iostream>
 #include <algorithm>
@@ -15,11 +14,12 @@
 #include <chrono>
 
 
+
 void run_spectrum(
     std::unordered_map<std::string, double>& parameter,
     std::unordered_map<std::string, double> optional,
-    Motion& motion,
-    Distribution& dist
+    motions::BaseMotion& motion,
+    times::BaseDistribution& dist
     ) {
     const int num_walker = static_cast<int>(parameter["num_walker"]);
 
@@ -83,15 +83,15 @@ void run_spectrum(
 void run_ste(
     std::unordered_map<std::string, double>& parameter,
     std::unordered_map<std::string, double> optional,
-    Motion& motion,
-    Distribution& dist
+    motions::BaseMotion& motion,
+    times::BaseDistribution& dist
     ) {
     const int num_walker = static_cast<int>(parameter[std::string("num_walker")]);
 
     const int num_mix_times = static_cast<int>(parameter["tmix_steps"]);
     const std::vector<double> evolution_times = linspace(parameter["tevo_start"], parameter["tevo_stop"], static_cast<int>(parameter["tevo_steps"]));
     const std::vector<double> mixing_times = logspace(parameter["tmix_start"], parameter["tmix_stop"], num_mix_times);
-    const double tpulse4= parameter["tpulse4"];
+    const double tpulse4 = parameter["tpulse4"];
 
     // make ste decay vectors and set them to zero
     std::map<double, std::vector<double>> cc_dict;
@@ -175,8 +175,8 @@ void run_ste(
 
 
 void make_trajectory(
-    Motion& motion,
-    Distribution& dist,
+    motions::BaseMotion& motion,
+    times::BaseDistribution& dist,
     const double t_max,
     std::vector<double>& out_time,
     std::vector<double>& out_phase,
diff --git a/src/simulation/sims.h b/src/sims.h
similarity index 73%
rename from src/simulation/sims.h
rename to src/sims.h
index 08fbf5a..8414185 100644
--- a/src/simulation/sims.h
+++ b/src/sims.h
@@ -1,8 +1,8 @@
 #ifndef RWSIM_SIMS_H
 #define RWSIM_SIMS_H
 
-#include "../motions/base.h"
-#include "../times/base.h"
+#include "motions/base.h"
+#include "times/base.h"
 
 #include <unordered_map>
 #include <string>
@@ -16,7 +16,7 @@
  * @param motion Motion model
  * @param dist Distribution of correlation times
  */
-void run_spectrum(std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double> optional, Motion& motion, Distribution& dist);
+void run_spectrum(std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double> optional, motions::BaseMotion& motion, times::BaseDistribution& dist);
 
 /**
  * @brief Run simulation for stimulated echoes
@@ -26,7 +26,7 @@ void run_spectrum(std::unordered_map<std::string, double>& parameter, std::unord
  * @param motion Motion model
  * @param dist Distribution of correlation times
  */
-void run_ste(std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double> optional, Motion& motion, Distribution& dist);
+void run_ste(std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double> optional, motions::BaseMotion& motion, times::BaseDistribution& dist);
 
 /**
  * @brief Create trajectory of a single walker
@@ -38,7 +38,7 @@ void run_ste(std::unordered_map<std::string, double>& parameter, std::unordered_
  * @param out_phase Vector of phase between waiting times
  * @param out_omega Vector of omega at jump time
  */
-void make_trajectory(Motion& motion, Distribution& dist, double t_max, std::vector<double>& out_time, std::vector<double>& out_phase, std::vector<double>& out_omega);
+void make_trajectory(motions::BaseMotion& motion, times::BaseDistribution& dist, double t_max, std::vector<double>& out_time, std::vector<double>& out_phase, std::vector<double>& out_omega);
 
 std::chrono::system_clock::time_point printStart(std::unordered_map<std::string, double> &optional);
 void printEnd(std::chrono::system_clock::time_point start);
diff --git a/src/times/base.cpp b/src/times/base.cpp
index 1d36462..500e651 100644
--- a/src/times/base.cpp
+++ b/src/times/base.cpp
@@ -4,31 +4,33 @@
 
 #include <stdexcept>
 
-Distribution::Distribution(std::string name, const double tau, std::mt19937_64 &rng) : m_name(std::move(name)), m_tau(tau), m_tau_jump(tau), m_rng(rng) {}
+namespace times {
+    BaseDistribution::BaseDistribution(std::string name, const double tau, std::mt19937_64 &rng) : m_name(std::move(name)), m_tau(tau), m_tau_jump(tau), m_rng(rng) {}
 
-Distribution::Distribution(std::string name, std::mt19937_64 &rng) : m_name(std::move(name)), m_rng(rng) {}
+    BaseDistribution::BaseDistribution(std::string name, std::mt19937_64 &rng) : m_name(std::move(name)), m_rng(rng) {}
 
-double Distribution::tau_wait() const {
-    return std::exponential_distribution(1./m_tau_jump)(m_rng);
-}
+    double BaseDistribution::tau_wait() const {
+        return std::exponential_distribution(1./m_tau_jump)(m_rng);
+    }
 
-void Distribution::setParameters(const std::unordered_map<std::string, double> &parameters) {
-    m_tau = parameters.at("tau");
-}
+    void BaseDistribution::setParameters(const std::unordered_map<std::string, double> &parameters) {
+        m_tau = parameters.at("tau");
+    }
 
-std::unordered_map<std::string, double> Distribution::getParameters() const {
-    return std::unordered_map<std::string, double>{
-        {"tau", m_tau},
-    };
-}
+    std::unordered_map<std::string, double> BaseDistribution::getParameters() const {
+        return std::unordered_map<std::string, double>{
+            {"tau", m_tau},
+        };
+    }
 
 
-Distribution* Distribution::createFromInput(const std::string& input, std::mt19937_64& rng) {
-    if (input == "Delta")
-        return new DeltaDistribution(rng);
+    BaseDistribution* BaseDistribution::createFromInput(const std::string& input, std::mt19937_64& rng) {
+        if (input == "Delta")
+            return new DeltaDistribution(rng);
 
-    if (input == "LogNormal")
-        return new LogNormalDistribution(rng);
+        if (input == "LogNormal")
+            return new LogNormalDistribution(rng);
 
-    throw std::invalid_argument("Invalid input " + input);
+        throw std::invalid_argument("Invalid input " + input);
+    }
 }
\ No newline at end of file
diff --git a/src/times/base.h b/src/times/base.h
index b826798..afcf136 100644
--- a/src/times/base.h
+++ b/src/times/base.h
@@ -4,33 +4,35 @@
 #include <random>
 #include <unordered_map>
 
-class Distribution {
-public:
-    virtual ~Distribution() = default;
+namespace times {
+    class BaseDistribution {
+    public:
+        virtual ~BaseDistribution() = default;
 
-    Distribution(std::string, double, std::mt19937_64&);
-    explicit Distribution(std::string, std::mt19937_64&);
+        BaseDistribution(std::string, double, std::mt19937_64&);
+        explicit BaseDistribution(std::string, std::mt19937_64&);
 
-    [[nodiscard]] double getTau() const { return m_tau; }
-    void setTau(const double tau) { m_tau = tau; }
-    [[nodiscard]] std::string getName() const { return m_name; };
+        [[nodiscard]] double getTau() const { return m_tau; }
+        void setTau(const double tau) { m_tau = tau; }
+        [[nodiscard]] std::string getName() const { return m_name; };
 
-    virtual void setParameters(const std::unordered_map<std::string, double>&);
-    [[nodiscard]] virtual std::unordered_map<std::string, double> getParameters() const;
+        virtual void setParameters(const std::unordered_map<std::string, double>&);
+        [[nodiscard]] virtual std::unordered_map<std::string, double> getParameters() const;
 
-    virtual void initialize() = 0;
-    virtual void draw_tau() = 0;
-    [[nodiscard]] double tau_wait() const;
+        virtual void initialize() = 0;
+        virtual void draw_tau() = 0;
+        [[nodiscard]] double tau_wait() const;
 
-    [[nodiscard]] virtual std::string toString() const = 0;
+        [[nodiscard]] virtual std::string toString() const = 0;
 
-    static Distribution* createFromInput(const std::string& input, std::mt19937_64& rng);
+        static BaseDistribution* createFromInput(const std::string& input, std::mt19937_64& rng);
 
-protected:
-    std::string m_name{"BaseDistribution"};
-    double m_tau{1.};
-    double m_tau_jump{1.};
-    std::mt19937_64& m_rng;
-};
+    protected:
+        std::string m_name{"BaseDistribution"};
+        double m_tau{1.};
+        double m_tau_jump{1.};
+        std::mt19937_64& m_rng;
+    };
+}
 
 #endif //RWSIM_TIMESBASE_H
diff --git a/src/times/delta.cpp b/src/times/delta.cpp
index f8201e8..8f59508 100644
--- a/src/times/delta.cpp
+++ b/src/times/delta.cpp
@@ -1,15 +1,16 @@
 #include "delta.h"
 
-DeltaDistribution::DeltaDistribution(const double tau, std::mt19937_64& rng) : Distribution(std::string("Delta"), tau, rng) {}
-DeltaDistribution::DeltaDistribution(std::mt19937_64& rng) : Distribution(std::string("Delta"), rng) {}
+namespace times {
+    DeltaDistribution::DeltaDistribution(const double tau, std::mt19937_64& rng) : BaseDistribution(std::string("Delta"), tau, rng) {}
+    DeltaDistribution::DeltaDistribution(std::mt19937_64& rng) : BaseDistribution(std::string("Delta"), rng) {}
 
-void DeltaDistribution::initialize() {
-    m_tau_jump = m_tau;
+    void DeltaDistribution::initialize() {
+        m_tau_jump = m_tau;
+    }
+
+    void DeltaDistribution::draw_tau() {}
+
+    std::string DeltaDistribution::toString() const {
+        return {"Delta/tau=" + std::to_string(m_tau)};
+    }
 }
-
-void DeltaDistribution::draw_tau() {}
-
-std::string DeltaDistribution::toString() const {
-    return {"Delta/tau=" + std::to_string(m_tau)};
-}
-
diff --git a/src/times/delta.h b/src/times/delta.h
index 4b0b35a..33b527f 100644
--- a/src/times/delta.h
+++ b/src/times/delta.h
@@ -3,15 +3,16 @@
 
 #include "base.h"
 
-class DeltaDistribution final : public Distribution {
-public:
-    DeltaDistribution(double, std::mt19937_64&);
-    explicit DeltaDistribution(std::mt19937_64 &rng);
+namespace times {
+    class DeltaDistribution final : public BaseDistribution {
+    public:
+        DeltaDistribution(double, std::mt19937_64&);
+        explicit DeltaDistribution(std::mt19937_64 &rng);
 
-    void initialize() override;
-    void draw_tau() override;
-
-    [[nodiscard]] std::string toString() const override;
-};
+        void initialize() override;
+        void draw_tau() override;
 
+        [[nodiscard]] std::string toString() const override;
+    };
+}
 #endif //RWSIM_TIMESDELTA_H
diff --git a/src/times/lognormal.cpp b/src/times/lognormal.cpp
index 43fa0d2..023e1ea 100644
--- a/src/times/lognormal.cpp
+++ b/src/times/lognormal.cpp
@@ -1,29 +1,31 @@
 #include "lognormal.h"
 #include <cmath>
 
-LogNormalDistribution::LogNormalDistribution(const double tau, const double sigma, std::mt19937_64& rng) : Distribution(std::string("LogNormal"), tau, rng), m_sigma(sigma), m_distribution(std::log(tau), sigma) {}
-LogNormalDistribution::LogNormalDistribution(std::mt19937_64& rng) : Distribution(std::string("LogNormal"), rng) {}
+namespace times {
+    LogNormalDistribution::LogNormalDistribution(const double tau, const double sigma, std::mt19937_64& rng) : BaseDistribution(std::string("LogNormal"), tau, rng), m_sigma(sigma), m_distribution(std::log(tau), sigma) {}
+    LogNormalDistribution::LogNormalDistribution(std::mt19937_64& rng) : BaseDistribution(std::string("LogNormal"), rng) {}
 
-void LogNormalDistribution::setParameters(const std::unordered_map<std::string, double> &parameters) {
-    m_sigma = parameters.at("sigma");
-    Distribution::setParameters(parameters);
-}
+    void LogNormalDistribution::setParameters(const std::unordered_map<std::string, double> &parameters) {
+        m_sigma = parameters.at("sigma");
+        BaseDistribution::setParameters(parameters);
+    }
 
-std::unordered_map<std::string, double> LogNormalDistribution::getParameters() const {
-    auto parameter = Distribution::getParameters();
-    parameter["sigma"] = m_sigma;
-    return parameter;
-}
+    std::unordered_map<std::string, double> LogNormalDistribution::getParameters() const {
+        auto parameter = BaseDistribution::getParameters();
+        parameter["sigma"] = m_sigma;
+        return parameter;
+    }
 
-void LogNormalDistribution::initialize() {
-    m_distribution = std::lognormal_distribution(std::log(m_tau), m_sigma);
-    m_tau_jump = m_distribution(m_rng);
-}
+    void LogNormalDistribution::initialize() {
+        m_distribution = std::lognormal_distribution(std::log(m_tau), m_sigma);
+        m_tau_jump = m_distribution(m_rng);
+    }
 
-void LogNormalDistribution::draw_tau() {
-    m_tau_jump = m_distribution(m_rng);
-}
+    void LogNormalDistribution::draw_tau() {
+        m_tau_jump = m_distribution(m_rng);
+    }
 
-std::string LogNormalDistribution::toString() const {
-    return {"LogNormal/tau=" + std::to_string(m_tau) + "/sigma=" + std::to_string(m_sigma)};
+    std::string LogNormalDistribution::toString() const {
+        return {"LogNormal/tau=" + std::to_string(m_tau) + "/sigma=" + std::to_string(m_sigma)};
+    }
 }
diff --git a/src/times/lognormal.h b/src/times/lognormal.h
index e41e984..b2991ac 100644
--- a/src/times/lognormal.h
+++ b/src/times/lognormal.h
@@ -4,22 +4,24 @@
 #include "base.h"
 #include <random>
 
-class LogNormalDistribution final : public Distribution {
-public:
-    LogNormalDistribution(double, double, std::mt19937_64&);
-    explicit LogNormalDistribution(std::mt19937_64 &rng);
+namespace times {
+    class LogNormalDistribution final : public BaseDistribution {
+    public:
+        LogNormalDistribution(double, double, std::mt19937_64&);
+        explicit LogNormalDistribution(std::mt19937_64 &rng);
 
-    void setParameters(const std::unordered_map<std::string, double> &) override;
-    [[nodiscard]] std::unordered_map<std::string, double> getParameters() const override;
+        void setParameters(const std::unordered_map<std::string, double> &) override;
+        [[nodiscard]] std::unordered_map<std::string, double> getParameters() const override;
 
-    [[nodiscard]] std::string toString() const override;
+        [[nodiscard]] std::string toString() const override;
 
-    void initialize() override;
-    void draw_tau() override;
+        void initialize() override;
+        void draw_tau() override;
 
-private:
-    double m_sigma{1};
-    std::lognormal_distribution<> m_distribution;
-};
+    private:
+        double m_sigma{1};
+        std::lognormal_distribution<> m_distribution;
+    };
+}
 
 #endif //LOGNORMAL_H
diff --git a/src/utils/functions.cpp b/src/utils/functions.cpp
index 0399869..beed0f8 100644
--- a/src/utils/functions.cpp
+++ b/src/utils/functions.cpp
@@ -1,9 +1,6 @@
 #include <vector>
-#include <array>
 #include <chrono>
 #include <iostream>
-#include <cmath>
-#include <utility>
 
 
 int nearest_index(const std::vector<double> &x_ref, const double x, int start=0) {
diff --git a/src/utils/io.cpp b/src/utils/io.cpp
index 0dd5740..3664e57 100644
--- a/src/utils/io.cpp
+++ b/src/utils/io.cpp
@@ -129,8 +129,8 @@ std::unordered_map<std::string, double> read_parameter(const std::filesystem::pa
 
 
 std::string make_directory(
-    const Motion& motion,
-    const Distribution& distribution
+    const motions::BaseMotion& motion,
+    const times::BaseDistribution& distribution
     ) {
     std::ostringstream path_name;
     path_name << motion.toString() << "/" << distribution.toString();
diff --git a/src/utils/io.h b/src/utils/io.h
index 772a52e..e2493b1 100644
--- a/src/utils/io.h
+++ b/src/utils/io.h
@@ -23,7 +23,7 @@ Arguments parse_args(int argc, char* argv[]);
 std::pair<std::string, double> get_optional_parameter(std::vector<std::string>::const_iterator &it);
 std::unordered_map<std::string, double> read_parameter(const std::filesystem::path&);
 
-std::string make_directory(const Motion&, const Distribution&);
+std::string make_directory(const motions::BaseMotion&, const times::BaseDistribution&);
 
 void save_parameter_to_file(const std::string&, const std::string&, const std::unordered_map<std::string, double>&, const std::unordered_map<std::string, double>&);
 void save_data_to_file(const std::string&, const std::string&, const std::vector<double>&, const std::map<double, std::vector<double>>&, const std::unordered_map<std::string, double>&);