Comments

Dynamic scheduling of threads
Reuse trajectories for multiple experiments
2026-03-08 15:33:56 +01:00 · 2026-03-08 15:02:47 +01:00 · 2026-03-08 14:28:33 +01:00 · 2026-03-08 14:01:37 +01:00 · 2026-03-08 13:43:51 +01:00 · 2026-03-08 13:26:05 +01:00
68 changed files with 2197 additions and 1051 deletions
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -3,6 +3,9 @@ project(RWSim VERSION 1.0)

 set(CMAKE_CXX_STANDARD 17)

+find_package(OpenMP REQUIRED)
+
 add_subdirectory(src)

 target_compile_options(rwsim PUBLIC -Werror -Wall -Wextra -Wconversion -O2)
+target_link_libraries(rwsim PUBLIC OpenMP::OpenMP_CXX)
--- a/README.md
+++ b/README.md
@@ -1,32 +1,37 @@
-# Download
+# NMR Random Walk Simulation

-1. Clone repository with `git clone https://gitea.pkm.physik.tu-darmstadt.de/NMRRWSims/cpp.git` 
-2. After download, change permissions of build.sh in terminal with `chmod a+x build.sh`
+Simulate solid-state NMR spectra and stimulated echo (STE) decays using random walk models.

 # Build

-# Running
-
-## Command line
-
-To run a random walk simulation Execute `rwsim` in the command line with 
+Requires CMake 3.18+, a C++17 compiler, and OpenMP.

 ```bash
-  ./rwsim /path/to/config.txt MotionModel DistributionModel --ste --spectrum -ARG1 1 -ARG2 2
+cmake -B build
+cmake --build build
 ```

-It needs three positional arguments: the path to your basic configuration file (see below) and the names of the motional model and of the distribution of correlation times.
-Set the optional arguments `--ste` and/or `--spectrum` to create Stimulated Echos or normal echo spectra, respectively.
-You can overwrite any parameter given in the configuration file by adding it as optional argument with a numerical value, e.g. `-TAU 1e-3` for a correlation time of 1 ms.
+The executable is `build/src/rwsim`.

+# Running

-## Configuration file `config.txt`
+```bash
+./build/src/rwsim /path/to/config.txt MotionModel DistributionModel --ste --spectrum -ARG1 1 -ARG2 2
+```

-Change the values of delta, eta, mixing times, echo times,... in this file. If a paramter is defined multiple times, only the last one is used.  
+Three positional arguments are required: the path to a configuration file, the motion model name, and the distribution model name.
+Set `--ste` and/or `--spectrum` to run stimulated echo or spectrum simulations.
+Any parameter from the configuration file can be overridden on the command line, e.g. `-tau 1e-3`.

-### General parameter 
+Use `-seed 42` for reproducible results. Without a seed, a random one is used.

-This list of general parameter are necessary for all simulations:
+**Output files are overwritten if a simulation with the same parameters is run again.**
+
+## Configuration file
+
+Change delta, eta, mixing times, echo times, etc. in this file. If a parameter appears multiple times, only the last value is used.
+
+### General parameters

 | Parameter  | Description                |
 |------------|----------------------------|
@@ -36,48 +41,88 @@ This list of general parameter are necessary for all simulations:

 ### Distribution of correlation times

-Two distributions are available: A delta distribution `Delta`, i.e. the same correlation time for every walker, and a log-normal distribution `LogNormal`.
-
-+ Parameters for delta distribution `Delta`
-
-  | Parameter | Description           |
-  |-----------|-----------------------|
-  | tau       | Jump time in s        |
-
-+ Parameters for log-normal distribution `LogNormal`
-
-  | Parameter | Description                                |
-  |-----------|--------------------------------------------|
-  | tau       | Maximum jump time of the distribution in s |
-  | sigma     | Standard deviation of the distribution     |
+| Model | CLI name | Parameters |
+|-------|----------|------------|
+| Delta distribution (same tau for all walkers) | `Delta` | `tau` (jump time in s) |
+| Log-normal distribution | `LogNormal` | `tau` (peak time in s), `sigma` (width) |

 ### Motion models

-Four different jump models are implemented: Isotropic random jump `RandomJump`, isotropic jump with a certain angle, isotropic jump with a bimodal distribution of angles, and a tetrahedral jump `TetrahedralJump`.
+| Model | CLI name | Parameters |
+|-------|----------|------------|
+| Isotropic random jump | `RandomJump` | — |
+| Isotropic jump by fixed angle | `IsotropicAngle` | `angle` (degrees) |
+| Bimodal angle distribution | `BimodalAngle` | `angle1`, `angle2` (degrees), `probability1` (0–1) |
+| Tetrahedral 4-site jump | `FourSiteTetrahedral` | — |
+| Octahedral 6-site jump (C3) | `SixSiteOctahedralC3` | — |
+| N-site jump on cone | `NSiteConeJump` | `num_sites`, `chi` (cone angle, degrees) |
+| Random jump on cone | `RandomJumpOnCone` | `angle` (cone angle) |
+| Wobbling in cone | `ConeWobble` | `angle` (cone angle) |

-+ Isotropic random jump `RandomJump` does not have additional parameters. 
-+ Tetrahedral jump `TetrahedralJump` does not have additional parameters.
-+ Parameters for isotropic jump by an angle `IsotropicAngle`
+### Spectrum parameters

-  | Parameter | Description          |
-  |-----------|----------------------|
-  | angle     | Jump angle in degree |
+| Parameter    | Description                     |
+|--------------|---------------------------------|
+| dwell_time   | Acquisition dwell time in s     |
+| num_acq      | Number of acquisition points    |
+| techo_start  | First echo time in s            |
+| techo_stop   | Last echo time in s             |
+| techo_steps  | Number of echo times            |

-+ Parameters for isotropic jump with bimodal angle distribution `BimodalAngle`
+### STE parameters

-  | Parameter    | Description                          |
-  |--------------|--------------------------------------|
-  | angle1       | First jump angle in degree           |
-  | angle2       | Second jump angle in degree          |
-  | probability1 | Probality that jump has angle1 (0-1) |  
+| Parameter   | Description                     |
+|-------------|---------------------------------|
+| tevo_start  | First evolution time in s       |
+| tevo_stop   | Last evolution time in s        |
+| tevo_steps  | Number of evolution times       |
+| tmix_start  | First mixing time in s          |
+| tmix_stop   | Last mixing time in s           |
+| tmix_steps  | Number of mixing times (log-spaced) |
+| tpulse4     | Delay of 4th pulse in s         |

+# Architecture

+The simulation is built around three abstractions:

+## Dynamics

+`Dynamics` (`src/dynamics/base.h`) produces trajectory steps — pairs of `{dt, omega}` representing a waiting time and the NMR frequency after a jump. This is the core interface for defining physical models.

+**`DecoupledDynamics`** wraps a separate motion model and time distribution, calling them independently. This covers all existing models.

-# Running
+To implement a **coupled model** where motion and time are interdependent (e.g. jump angle affects waiting time, or position-dependent dynamics), implement the `Dynamics` interface directly:

-**Because filenames are always the same, previous simulations and results are overwritten!**
+```cpp
+class MyModel : public Dynamics {
+    void initialize(std::mt19937_64& rng) override { /* set initial state */ }
+    Step next(std::mt19937_64& rng) override {
+        // draw angle, compute rate from angle, draw wait time from rate
+        return {dt, omega};
+    }
+    // ... clone(), setParameters(), etc.
+};
+```

-1. Execute `build.sh` (in terminal with `./build.sh`). It compiles the source code and executes test.py
+## Experiment
+
+`Experiment` (`src/experiments/base.h`) defines how trajectory data is turned into observables:
+
+- `setup(params)` — configure time axes, allocate accumulators
+- `accumulate(trajectory)` — process one walker's trajectory
+- `save(directory)` — write results to files
+- `clone()` / `merge()` — for thread-parallel accumulation
+
+Existing experiments: `SpectrumExperiment`, `STEExperiment`.
+
+## Simulation runner
+
+`run_simulation()` (`src/sims.h`) connects a `Dynamics` and an `Experiment`:
+
+1. Sets parameters on both
+2. Generates trajectories in parallel (OpenMP, one clone per thread)
+3. Each thread accumulates into its own experiment clone
+4. Merges thread-local results
+5. Saves output
+
+The walker loop is parallelized with OpenMP using static scheduling for deterministic reproducibility with a given seed and thread count.
--- a/angle_dependence.py
+++ b/angle_dependence.py
@@ -0,0 +1,61 @@
+import pathlib
+import re
+
+import numpy as np
+from matplotlib import pyplot as plt
+
+from python.helpers import read_parameter_file
+
+angles = []
+tau_cc = []
+tau_ss = []
+tau2 = []
+
+tevo = [5e-6, 10e-6]
+
+
+for fit_files in pathlib.Path('.').glob(f'IsotropicAngle/angle=*/Delta/tau=*/ste_fit_*.dat'):
+    folder = fit_files.parent
+
+    file_base = fit_files.stem.replace('ste_fit_', '')
+    parameter = read_parameter_file(folder / ('ste_' + file_base + '_parameter.txt'))
+    angles.append(parameter['angle'])
+
+    with fit_files.open('r') as f:
+        # tau of F2 is hidden in the second header line
+        for _ in range(2):
+            line = f.readline()
+        tau2.append(float(re.search('tau=(.+?)\s', line).group(1)))
+
+    fit_values = np.loadtxt(fit_files)
+    x = fit_values[:, 0]
+
+    # get indexes for given evolution times
+    nearest_idx = [np.searchsorted(x, tt) for tt in tevo]
+    tau_cc.append(fit_values[nearest_idx, 1])
+    tau_ss.append(fit_values[nearest_idx, 4])
+
+angles = np.array(angles)
+tau_cc = np.array(tau_cc)
+tau_ss = np.array(tau_ss)
+tau2 = np.array(tau2)
+
+
+sortidx = np.argsort(angles)
+angles = angles[sortidx]
+tau_cc = tau_cc[sortidx]
+tau_ss = tau_ss[sortidx]
+tau2 = tau2[sortidx]
+
+
+plt.semilogy(angles, tau_cc, '-.')
+plt.semilogy(angles, tau_ss, '--')
+plt.semilogy(angles, tau2)
+
+np.savetxt(
+    'tau_angles.dat',
+    np.c_[angles, tau_cc, tau_ss, tau2],
+    header=f"Angle dependence of correlation times\nfor evolution times {tevo}\nangle->tau_cc->tauss->tau2"
+)
+
+plt.show()
--- a/config.txt
+++ b/config.txt
@@ -5,13 +5,10 @@ delta=126e3
 eta=0.0
 # Distribution part
 tau=1e-3
-angle1=2
-angle2=30
-probability1=0
 # Spectrum part
 dwell_time=1e-8
 num_acq=4096
-techo_start=0e-6
+techo_start=1e-6
 techo_stop=40e-6
 techo_steps=5
 # STE part
@@ -20,4 +17,5 @@ tevo_stop=120e-6
 tevo_steps=121
 tmix_start=1e-5
 tmix_stop=1e1
-tmix_steps=61
+tmix_steps=31
+tpulse4=10e-6
--- a/process_spectrum.py
+++ b/process_spectrum.py
@@ -0,0 +1,41 @@
+import pathlib
+
+import matplotlib.pyplot as plt
+
+from python.spectrum import *
+
+tpulse = 3e-6
+gb = 4e3
+
+fig_spec, ax_spec = plt.subplots()
+
+# for conf_file in pathlib.Path('.').glob(f'SixSiteOctahedral/angle=59/Delta/tau=*/timesignal_*_parameter.txt'):
+for conf_file in pathlib.Path('.').glob(f'RandomJumpOnCone/angle=128.*/Delta/tau=*/timesignal_*_parameter.txt'):
+    vary_string = post_process_spectrum(conf_file, tpulse=tpulse, apod=gb)
+    print(conf_file)
+
+#     ax_spec
+#
+#     ax_tau_cc.semilogy(tau_cc[:, 0], tau_cc[:, 1], label=vary_string)
+#     ax_tau_cc.axhline(tau_2[:, 1], color='k', linestyle='--')
+#     ax_beta_cc.plot(*beta_cc.T, label=vary_string)
+#     ax_finfty_cc.plot(*finfty_cc.T, label=vary_string)
+#     ax_tau_ss.semilogy(tau_ss[:, 0], tau_ss[:, 1], label=vary_string)
+#     ax_tau_ss.axhline(tau_2[:, 1], color='k', linestyle='--')
+#     ax_beta_ss.plot(*beta_ss.T, label=vary_string)
+#     ax_finfty_ss.plot(*finfty_ss.T, label=vary_string)
+#
+#     np.savetxt(
+#         conf_file.with_name(f'ste_fit_{vary_string}.dat'),
+#         np.c_[
+#             tau_cc, beta_cc[:, 1], finfty_cc[:, 1],
+#             tau_ss[:, 1], beta_ss[:, 1], finfty_ss[:, 1],
+#         ],
+#         header=f'Fit STE {vary_string}\n'
+#                f'F2: tau={tau_2[0, 1]} beta={beta_2[0, 1]} finfty={finfty_2[0, 1]}\n'
+#                f'tevo\ttaucc\tbetacc\tfinftycc\ttauss\tbetass\tfinftyss',
+#     )
+#
+# for ax in [ax_tau_cc, ax_beta_cc, ax_finfty_cc, ax_tau_ss, ax_beta_ss, ax_finfty_ss]:
+#     ax.legend()
+# plt.show()
--- a/process_ste.py
+++ b/process_ste.py
@@ -1,17 +1,4 @@
-import matplotlib.pyplot as plt
-
 from python.ste import *
-from python.helpers import *
-
-# Simulation parameter
-motion = 'IsotropicAngle'
-distribution = 'Delta'
-# parameter = {}
-parameter = {
-    "angle": [10, 109.47],
-}
-
-parameter = prepare_rw_parameter(parameter)

 fig_tau_cc, ax_tau_cc = plt.subplots()
 ax_tau_cc.set_title('tau_cc')
@@ -31,15 +18,11 @@ ax_beta_ss.set_title('beta_ss')
 fig_finfty_ss, ax_finfty_ss = plt.subplots()
 ax_finfty_ss.set_title('f_infty_ss')

-for variation in parameter:
-    print(f"\nRun RW for {motion}/{distribution} with arguments {variation}\n")
-    run_sims(motion, distribution, ste=True, spectrum=False, **variation)
-
-    conf_file = find_config_file(motion, distribution, variation)
-
-    vary_string, tau_cc, beta_cc, finfty_cc = fit_and_save_ste(conf_file, 'coscos', plot_decays=False, verbose=False)
-    _, tau_ss, beta_ss, finfty_ss = fit_and_save_ste(conf_file, 'sinsin', plot_decays=False, verbose=False)
-    _, tau_2, beta_2, finfty_2 = fit_and_save_ste(conf_file, 'f2', plot_decays=True, verbose=True)
+for conf_file in pathlib.Path('.').glob(f'IsotropicAngle/angle=*/Delta/tau=*/*_parameter.txt'):
+    print(conf_file)
+    vary_string, tau_cc, beta_cc, finfty_cc = fit_ste(conf_file, f'coscos', plot_decays=False, verbose=False)
+    _, tau_ss, beta_ss, finfty_ss = fit_ste(conf_file, f'sinsin', plot_decays=False, verbose=False)
+    _, tau_2, beta_2, finfty_2 = fit_ste(conf_file, f'f2', plot_decays=True, verbose=True)

    ax_tau_cc.semilogy(tau_cc[:, 0], tau_cc[:, 1], label=vary_string)
    ax_tau_cc.axhline(tau_2[:, 1], color='k', linestyle='--')
@@ -50,6 +33,17 @@ for variation in parameter:
    ax_beta_ss.plot(*beta_ss.T, label=vary_string)
    ax_finfty_ss.plot(*finfty_ss.T, label=vary_string)

+    np.savetxt(
+        conf_file.with_name(f'ste_fit_{vary_string}.dat'),
+        np.c_[
+            tau_cc, beta_cc[:, 1], finfty_cc[:, 1],
+            tau_ss[:, 1], beta_ss[:, 1], finfty_ss[:, 1],
+        ],
+        header=f'Fit STE {vary_string}\n'
+               f'F2: tau={tau_2[0, 1]} beta={beta_2[0, 1]} finfty={finfty_2[0, 1]}\n'
+               f'tevo\ttaucc\tbetacc\tfinftycc\ttauss\tbetass\tfinftyss',
+    )
+
 for ax in [ax_tau_cc, ax_beta_cc, ax_finfty_cc, ax_tau_ss, ax_beta_ss, ax_finfty_ss]:
    ax.legend()
-plt.show()
+plt.show()
--- a/python/helpers.py
+++ b/python/helpers.py
@@ -1,10 +1,10 @@
 from __future__ import annotations

 import pathlib
-import re
 import subprocess
 from itertools import product

+
 def prepare_rw_parameter(parameter: dict) -> list:
    """
    Converts a dictionary of iterables to list of dictionaries
@@ -53,23 +53,11 @@ def run_sims(
    subprocess.run(arguments)


-def find_config_file(motion: str, distribution: str, var_params: dict) -> pathlib.Path:
-    # TODO handle situation if multiple files fit
-    p_file = None
-    if var_params:
-        var_string = '|'.join(([f'{k}={v:1.6e}' for (k, v) in var_params.items()])).replace('.', '\.').replace('+', '\+')
-        pattern = re.compile(var_string)
-        for p_file in pathlib.Path('.').glob('*_parameter.txt'):
-            if len(re.findall(pattern, str(p_file))) == len(var_params) and re.search(f'{motion}_{distribution}', str(p_file)):
-                return p_file
-        raise ValueError(f'No parameter file found for {motion}, {distribution}, {var_params}')
-    else:
-        for p_file in pathlib.Path('.').glob('*_parameter.txt'):
-            if re.search(f'{motion}_{distribution}', str(p_file)):
-                return p_file
-        raise ValueError(f'No parameter file found for {motion}, {distribution}, {var_params}')
-
+def find_config_file(config_path: str | pathlib.Path, varied_params: dict[str, float]) -> dict[str, float]:
+    parameter = read_parameter_file(config_path)
+    parameter.update(varied_params)

+    return parameter


 def read_parameter_file(path: str | pathlib.Path) -> dict[str, float]:
@@ -80,8 +68,9 @@ def read_parameter_file(path: str | pathlib.Path) -> dict[str, float]:
    parameter_dict = {}
    with path.open('r') as f:
        for line in f.readlines():
+            if line.startswith('#'):
+                continue
            k, v = line.split('=')
            parameter_dict[k] = float(v)

-            k, v = line.split('=')
    return parameter_dict
--- a/python/spectrum.py
+++ b/python/spectrum.py
@@ -1,8 +1,7 @@
-import numpy
 import numpy as np
-from matplotlib import pyplot
-
+import matplotlib.pyplot as plt

+from python.helpers import read_parameter_file

 # parameter for spectrum simulations
 lb = 2e3
@@ -41,30 +40,46 @@ def pulse_attn(freq: np.ndarray, t_pulse: float) -> np.ndarray:
    return np.pi * numerator / denominator / 2


-def post_process_spectrum(taus, apod, tpulse):
-    reduction_factor = np.zeros((taus.size, 5)) # hard-coded t_echo :(
+def post_process_spectrum(parameter_file, apod, tpulse):
+    parameter = read_parameter_file(parameter_file)

-    for i, tau in enumerate(taus):
-        try:
-            raw_data = np.loadtxt(f'fid_tau={tau:.6e}.dat')
-        except OSError:
-            continue
+    # files have form ste_arg=0.000000e+01_parameter.txt, first remove ste part then parameter.txt to get variables
+    varied_string = str(parameter_file).split('/')[-1].partition('_')[-1].rpartition('_')[0]
+    print(varied_string)

-        t = raw_data[:, 0]
-        timesignal = raw_data[:, 1:]
+    # make evolution times
+    tevo = np.linspace(parameter['techo_start'], parameter['techo_stop'], num=int(parameter['techo_steps']))

-        timesignal *= dampening(t, apod)[:, None]
-        timesignal[0, :] /= 2
+    if varied_string:
+        raw_data = np.loadtxt(parameter_file.with_name(f'timesignal_{varied_string}.dat'))
+    else:
+        raw_data = np.loadtxt(parameter_file.with_name(f'timesignal.dat'))

-        # FT to spectrum
-        freq = np.fft.fftshift(np.fft.fftfreq(t.size, d=1e-6))
-        spec = np.fft.fftshift(np.fft.fft(timesignal, axis=0), axes=0).real
-        spec *= pulse_attn(freq, t_pulse=tpulse)[:, None]
+    t = raw_data[:, 0]
+    timesignal = raw_data[:, 1:]

-        reduction_factor[i, :] = 2*timesignal[0, :]
+    timesignal *= dampening(t, apod)[:, None]
+    timesignal[0, :] /= 2

-        plt.plot(freq, spec)
-        plt.show()
+    # FT to spectrum
+    freq = np.fft.fftshift(np.fft.fftfreq(t.size, d=parameter['dwell_time']))
+    spec = np.fft.fftshift(np.fft.fft(timesignal, axis=0), axes=0).real
+    spec *= pulse_attn(freq, t_pulse=tpulse)[:, None]

-    plt.semilogx(taus, reduction_factor, '.')
+    if varied_string:
+        np.savetxt(parameter_file.with_name(f"spectrum_{varied_string}.dat"), np.c_[freq, spec])
+    else:
+        np.savetxt(parameter_file.with_name(f"spectrum.dat"), np.c_[freq, spec])
+
+
+    #
+    #
+    # reduction_factor[i, :] = 2*timesignal[0, :]
+
+
+    plt.plot(freq, spec)
+    plt.gca().set_title(varied_string)
    plt.show()
+    #
+    # plt.semilogx(taus, reduction_factor, '.')
+    # plt.show()
--- a/python/ste.py
+++ b/python/ste.py
@@ -1,6 +1,5 @@
 import pathlib

-import numpy
 import numpy as np
 from matplotlib import pyplot as plt
 from scipy.optimize import curve_fit
@@ -59,12 +58,12 @@ def fit_decay(x: np.ndarray, y: np.ndarray, tevo: np.ndarray, verbose: bool = Tr
    return tau_fit, beta_fit, finfty_fit


-def fit_and_save_ste(
+def fit_ste(
        parameter_file: pathlib.Path,
        prefix: str,
        plot_decays: bool = True,
        verbose: bool = True,
-) -> tuple[str, np.ndarray, np.ndarray, np.ndarray]:
+    ) -> tuple[str, np.ndarray, np.ndarray, np.ndarray]:
    # read simulation parameters
    parameter = read_parameter_file(parameter_file)

@@ -74,7 +73,10 @@ def fit_and_save_ste(
    # make evolution times
    tevo = np.linspace(parameter['tevo_start'], parameter['tevo_stop'], num=int(parameter['tevo_steps']))

-    raw_data = np.loadtxt(f'{prefix}_{varied_string}.dat')
+    if varied_string:
+        raw_data = np.loadtxt(parameter_file.with_name(f'{prefix}_{varied_string}.dat'))
+    else:
+        raw_data = np.loadtxt(parameter_file.with_name(f'{prefix}.dat'))

    t_mix = raw_data[:, 0]
    decay = raw_data[:, 1:]
@@ -89,9 +91,4 @@ def fit_and_save_ste(
    print(f'Fit {prefix}')
    tau, beta, finfty = fit_decay(t_mix, decay, tevo, verbose=verbose)

-
-    np.savetxt(f'tau_{prefix}_{varied_string}.dat', tau)
-    np.savetxt(f'beta_{prefix}_{varied_string}.dat', beta)
-    np.savetxt(f'finfty_{prefix}_{varied_string}.dat', finfty)
-
    return varied_string, tau, beta, finfty
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,35 +1,28 @@
-cmake_minimum_required(VERSION 3.18)

-set(CMAKE_CXX_STANDARD 17)
+add_subdirectory(times)
+add_subdirectory(motions)
+add_subdirectory(utils)
+add_subdirectory(experiments)
+add_subdirectory(dynamics)

-
-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_library(simulation STATIC 
+  sims.cpp sims.h
+)
+target_link_libraries(simulation PRIVATE 
+  utils 
+  experiments 
+  dynamics 
+  OpenMP::OpenMP_CXX
+)
+
+add_executable(rwsim main.cpp)
+
+target_link_libraries(rwsim PUBLIC 
+  RWMotion 
+  RWTime 
+  utils 
+  experiments 
+  dynamics 
+  simulation
 )

-target_compile_options(rwsim PUBLIC -Werror -Wall -Wextra -Wconversion -O2)
--- a/src/dynamics/CMakeLists.txt
+++ b/src/dynamics/CMakeLists.txt
@@ -0,0 +1,4 @@
+add_library(dynamics OBJECT
+  base.cpp base.h
+  decoupled.cpp decoupled.h
+)
--- a/src/dynamics/base.cpp
+++ b/src/dynamics/base.cpp
@@ -0,0 +1,6 @@
+#include "base.h"
+#include "../utils/registry.h"
+
+std::unique_ptr<Dynamics> Dynamics::createFromInput(const std::string &input) {
+  return Registry<Dynamics>::create(input);
+}
--- a/src/dynamics/base.h
+++ b/src/dynamics/base.h
@@ -0,0 +1,49 @@
+#ifndef RWSIM_DYNAMICSBASE_H
+#define RWSIM_DYNAMICSBASE_H
+
+#include <memory>
+#include <random>
+#include <string>
+#include <unordered_map>
+
+/// A single step in the dynamics: waiting time dt and new frequency omega.
+struct Step {
+  double dt;
+  double omega;
+};
+
+/// Interface for trajectory dynamics (motion + timing).
+///
+/// Produces a sequence of (dt, omega) steps that define a walker's trajectory.
+/// DecoupledDynamics wraps independent motion and time models.
+/// Coupled models (e.g., position-dependent jump rates) implement this directly.
+///
+/// Use 2 positional CLI args for coupled: ./rwsim config.txt MyDynamics
+/// Use 3 positional CLI args for decoupled: ./rwsim config.txt Motion Dist
+class Dynamics {
+public:
+  virtual ~Dynamics() = default;
+
+  virtual void
+  setParameters(const std::unordered_map<std::string, double> &parameters) = 0;
+
+  /// Set up initial state for a new walker (e.g., draw random orientation).
+  virtual void initialize(std::mt19937_64 &rng) = 0;
+
+  /// Generate the next step: waiting time and new NMR frequency.
+  virtual Step next(std::mt19937_64 &rng) = 0;
+
+  /// Return the NMR frequency at the walker's initial position.
+  [[nodiscard]] virtual double getInitOmega() const = 0;
+
+  /// Deep copy for per-thread cloning in parallel simulation.
+  [[nodiscard]] virtual std::unique_ptr<Dynamics> clone() const = 0;
+
+  /// String identifier used for output directory naming.
+  [[nodiscard]] virtual std::string toString() const = 0;
+
+  /// Factory: create a Dynamics instance by registered name.
+  static std::unique_ptr<Dynamics> createFromInput(const std::string &input);
+};
+
+#endif // RWSIM_DYNAMICSBASE_H
--- a/src/dynamics/decoupled.cpp
+++ b/src/dynamics/decoupled.cpp
@@ -0,0 +1,36 @@
+#include "decoupled.h"
+
+DecoupledDynamics::DecoupledDynamics(
+    std::unique_ptr<motions::BaseMotion> motion,
+    std::unique_ptr<times::BaseDistribution> dist)
+    : m_motion(std::move(motion)), m_dist(std::move(dist)) {}
+
+void DecoupledDynamics::setParameters(
+    const std::unordered_map<std::string, double> &parameters) {
+  m_motion->setParameters(parameters);
+  m_dist->setParameters(parameters);
+}
+
+void DecoupledDynamics::initialize(std::mt19937_64 &rng) {
+  m_motion->initialize(rng);
+  m_dist->initialize(rng);
+}
+
+Step DecoupledDynamics::next(std::mt19937_64 &rng) {
+  double dt = m_dist->tau_wait(rng);
+  double omega = m_motion->jump(rng);
+  return {dt, omega};
+}
+
+double DecoupledDynamics::getInitOmega() const {
+  return m_motion->getInitOmega();
+}
+
+std::unique_ptr<Dynamics> DecoupledDynamics::clone() const {
+  return std::make_unique<DecoupledDynamics>(m_motion->clone(),
+                                             m_dist->clone());
+}
+
+std::string DecoupledDynamics::toString() const {
+  return m_motion->toString() + "/" + m_dist->toString();
+}
--- a/src/dynamics/decoupled.h
+++ b/src/dynamics/decoupled.h
@@ -0,0 +1,32 @@
+#ifndef RWSIM_DECOUPLED_H
+#define RWSIM_DECOUPLED_H
+
+#include "../motions/base.h"
+#include "../times/base.h"
+#include "base.h"
+
+#include <memory>
+
+/// Dynamics where motion and waiting times are independent.
+/// Wraps a BaseMotion (geometry) and BaseDistribution (timing) and delegates
+/// to each independently. This is the default mode when two positional CLI
+/// args are given (motion + distribution).
+class DecoupledDynamics final : public Dynamics {
+public:
+  DecoupledDynamics(std::unique_ptr<motions::BaseMotion> motion,
+                    std::unique_ptr<times::BaseDistribution> dist);
+
+  void setParameters(
+      const std::unordered_map<std::string, double> &parameters) override;
+  void initialize(std::mt19937_64 &rng) override;
+  Step next(std::mt19937_64 &rng) override;
+  [[nodiscard]] double getInitOmega() const override;
+  [[nodiscard]] std::unique_ptr<Dynamics> clone() const override;
+  [[nodiscard]] std::string toString() const override;
+
+private:
+  std::unique_ptr<motions::BaseMotion> m_motion;
+  std::unique_ptr<times::BaseDistribution> m_dist;
+};
+
+#endif // RWSIM_DECOUPLED_H
--- a/src/experiments/CMakeLists.txt
+++ b/src/experiments/CMakeLists.txt
@@ -0,0 +1,7 @@
+add_library(experiments STATIC
+  base.h 
+  spectrum.cpp spectrum.h
+  ste.cpp ste.h
+)
+
+target_link_libraries(experiments PRIVATE utils)
--- a/src/experiments/base.h
+++ b/src/experiments/base.h
@@ -0,0 +1,47 @@
+#ifndef RWSIM_EXPERIMENTBASE_H
+#define RWSIM_EXPERIMENTBASE_H
+
+#include <memory>
+#include <string>
+#include <unordered_map>
+#include <vector>
+
+/// Time series of a single walker: NMR frequency, accumulated phase, and time.
+struct Trajectory {
+  std::vector<double> time;   ///< Cumulative time at each step
+  std::vector<double> phase;  ///< Accumulated NMR phase (integral of omega*dt)
+  std::vector<double> omega;  ///< Instantaneous NMR frequency at each step
+};
+
+/// Interface for NMR experiments that process walker trajectories.
+///
+/// Lifecycle: setup() -> [accumulate() per walker] -> merge() -> save()
+/// Each thread gets a clone(); results are merged after the parallel loop.
+class Experiment {
+public:
+  virtual ~Experiment() = default;
+
+  /// Configure time axes and allocate output buffers from simulation parameters.
+  virtual void
+  setup(const std::unordered_map<std::string, double> &parameter,
+        const std::unordered_map<std::string, double> &optional) = 0;
+
+  /// Maximum trajectory time needed by this experiment.
+  [[nodiscard]] virtual double tmax() const = 0;
+
+  /// Process one walker's trajectory, adding its contribution to the result.
+  /// Results are normalized by num_walkers during accumulation.
+  virtual void accumulate(const Trajectory &traj, double init_omega,
+                          int num_walkers) = 0;
+
+  /// Write results to files in the given directory.
+  virtual void save(const std::string &directory) = 0;
+
+  /// Deep copy for per-thread cloning in parallel simulation.
+  [[nodiscard]] virtual std::unique_ptr<Experiment> clone() const = 0;
+
+  /// Combine results from another (per-thread) experiment instance.
+  virtual void merge(const Experiment &other) = 0;
+};
+
+#endif // RWSIM_EXPERIMENTBASE_H
--- a/src/experiments/spectrum.cpp
+++ b/src/experiments/spectrum.cpp
@@ -0,0 +1,69 @@
+#include "spectrum.h"
+#include "../utils/functions.h"
+#include "../utils/io.h"
+#include "../utils/ranges.h"
+
+#include <algorithm>
+#include <cmath>
+
+void SpectrumExperiment::setup(
+    const std::unordered_map<std::string, double> &parameter,
+    const std::unordered_map<std::string, double> &optional) {
+  m_parameter = parameter;
+  m_optional = optional;
+  m_num_acq = static_cast<int>(parameter.at("num_acq"));
+  m_t_fid = arange(m_num_acq, parameter.at("dwell_time"));
+  m_echo_times =
+      linspace(parameter.at("techo_start"), parameter.at("techo_stop"),
+               static_cast<int>(parameter.at("techo_steps")));
+
+  m_fid_dict.clear();
+  for (auto t_echo_i : m_echo_times) {
+    m_fid_dict[t_echo_i] = std::vector<double>(m_num_acq, 0.);
+  }
+
+  m_tmax = *std::max_element(m_echo_times.begin(), m_echo_times.end()) * 2 +
+           m_t_fid.back();
+}
+
+double SpectrumExperiment::tmax() const { return m_tmax; }
+
+void SpectrumExperiment::accumulate(const Trajectory &traj, double,
+                                    int num_walkers) {
+  for (auto &[t_echo_j, fid_j] : m_fid_dict) {
+    int current_pos = nearest_index(traj.time, t_echo_j, 0);
+    const double phase_techo =
+        lerp(traj.time, traj.phase, t_echo_j, current_pos);
+
+    for (int acq_idx = 0; acq_idx < m_num_acq; acq_idx++) {
+      const double real_time = m_t_fid[acq_idx] + 2 * t_echo_j;
+
+      current_pos = nearest_index(traj.time, real_time, current_pos);
+      const double phase_acq =
+          lerp(traj.time, traj.phase, real_time, current_pos);
+
+      fid_j[acq_idx] += std::cos(phase_acq - 2 * phase_techo) / num_walkers;
+    }
+  }
+}
+
+void SpectrumExperiment::save(const std::string &directory) {
+  save_parameter_to_file(std::string("timesignal"), directory, m_parameter,
+                         m_optional);
+  save_data_to_file(std::string("timesignal"), directory, m_t_fid, m_fid_dict,
+                    m_optional);
+}
+
+std::unique_ptr<Experiment> SpectrumExperiment::clone() const {
+  return std::make_unique<SpectrumExperiment>(*this);
+}
+
+void SpectrumExperiment::merge(const Experiment &other) {
+  const auto &o = dynamic_cast<const SpectrumExperiment &>(other);
+  for (auto &[t_echo, fid] : m_fid_dict) {
+    const auto &other_fid = o.m_fid_dict.at(t_echo);
+    for (size_t i = 0; i < fid.size(); i++) {
+      fid[i] += other_fid[i];
+    }
+  }
+}
--- a/src/experiments/spectrum.h
+++ b/src/experiments/spectrum.h
@@ -0,0 +1,39 @@
+#ifndef RWSIM_SPECTRUM_H
+#define RWSIM_SPECTRUM_H
+
+#include "base.h"
+
+#include <map>
+#include <unordered_map>
+#include <vector>
+
+/// Solid-echo NMR spectrum experiment.
+///
+/// Computes the free induction decay (FID) after a solid echo at various echo
+/// times. For each echo time t_echo, the FID is sampled at t_fid points after
+/// refocusing at 2*t_echo. The signal is cos(phase_acq - 2*phase_echo),
+/// averaged over all walkers.
+///
+/// Output: timesignal files with FID data for each echo time.
+class SpectrumExperiment final : public Experiment {
+public:
+  void setup(const std::unordered_map<std::string, double> &parameter,
+             const std::unordered_map<std::string, double> &optional) override;
+  [[nodiscard]] double tmax() const override;
+  void accumulate(const Trajectory &traj, double init_omega,
+                  int num_walkers) override;
+  void save(const std::string &directory) override;
+  [[nodiscard]] std::unique_ptr<Experiment> clone() const override;
+  void merge(const Experiment &other) override;
+
+private:
+  int m_num_acq{0};                              ///< Number of acquisition points per FID
+  std::vector<double> m_t_fid;                   ///< FID time axis (dwell_time spacing)
+  std::vector<double> m_echo_times;              ///< Echo delay times
+  std::map<double, std::vector<double>> m_fid_dict; ///< FID data per echo time
+  std::unordered_map<std::string, double> m_parameter;
+  std::unordered_map<std::string, double> m_optional;
+  double m_tmax{0};
+};
+
+#endif // RWSIM_SPECTRUM_H
--- a/src/experiments/ste.cpp
+++ b/src/experiments/ste.cpp
@@ -0,0 +1,108 @@
+#include "ste.h"
+#include "../utils/functions.h"
+#include "../utils/io.h"
+#include "../utils/ranges.h"
+
+#include <algorithm>
+#include <cmath>
+
+void STEExperiment::setup(
+    const std::unordered_map<std::string, double> &parameter,
+    const std::unordered_map<std::string, double> &optional) {
+  m_parameter = parameter;
+  m_optional = optional;
+  m_num_mix_times = static_cast<int>(parameter.at("tmix_steps"));
+  m_evolution_times =
+      linspace(parameter.at("tevo_start"), parameter.at("tevo_stop"),
+               static_cast<int>(parameter.at("tevo_steps")));
+  m_mixing_times = logspace(parameter.at("tmix_start"),
+                            parameter.at("tmix_stop"), m_num_mix_times);
+  m_tpulse4 = parameter.at("tpulse4");
+
+  m_cc_dict.clear();
+  m_ss_dict.clear();
+  for (auto t_evo_i : m_evolution_times) {
+    m_cc_dict[t_evo_i] = std::vector<double>(m_num_mix_times, 0.);
+    m_ss_dict[t_evo_i] = std::vector<double>(m_num_mix_times, 0.);
+  }
+  m_f2.assign(m_num_mix_times, 0.);
+
+  m_tmax =
+      *std::max_element(m_evolution_times.begin(), m_evolution_times.end()) *
+          2 +
+      *std::max_element(m_mixing_times.begin(), m_mixing_times.end()) +
+      2 * m_tpulse4;
+}
+
+double STEExperiment::tmax() const { return m_tmax; }
+
+void STEExperiment::accumulate(const Trajectory &traj, double init_omega,
+                               int num_walkers) {
+  int f2_pos = 0;
+  for (int f2_idx = 0; f2_idx < m_num_mix_times; f2_idx++) {
+    const double t_mix_f2 = m_mixing_times[f2_idx];
+    f2_pos = nearest_index(traj.time, t_mix_f2, f2_pos);
+    m_f2[f2_idx] += traj.omega[f2_pos] * init_omega / num_walkers;
+  }
+
+  for (auto &[t_evo_j, _] : m_cc_dict) {
+    auto &cc_j = m_cc_dict[t_evo_j];
+    auto &ss_j = m_ss_dict[t_evo_j];
+
+    int current_pos = nearest_index(traj.time, t_evo_j, 0);
+    const double dephased = lerp(traj.time, traj.phase, t_evo_j, current_pos);
+    const double cc_tevo = std::cos(dephased);
+    const double ss_tevo = std::sin(dephased);
+
+    for (int mix_idx = 0; mix_idx < m_num_mix_times; mix_idx++) {
+      const double time_end_mix = m_mixing_times[mix_idx] + t_evo_j;
+      current_pos = nearest_index(traj.time, time_end_mix, current_pos);
+      const double phase_mix_end =
+          lerp(traj.time, traj.phase, time_end_mix, current_pos);
+
+      const double time_pulse4 = time_end_mix + m_tpulse4;
+      current_pos = nearest_index(traj.time, time_pulse4, current_pos);
+      const double phase_4pulse =
+          lerp(traj.time, traj.phase, time_pulse4, current_pos);
+
+      const double time_echo = time_pulse4 + m_tpulse4 + t_evo_j;
+      current_pos = nearest_index(traj.time, time_echo, current_pos);
+      double rephased = lerp(traj.time, traj.phase, time_echo, current_pos) +
+                        phase_mix_end - 2 * phase_4pulse;
+
+      cc_j[mix_idx] += cc_tevo * std::cos(rephased) / num_walkers;
+      ss_j[mix_idx] += ss_tevo * std::sin(rephased) / num_walkers;
+    }
+  }
+}
+
+void STEExperiment::save(const std::string &directory) {
+  save_parameter_to_file(std::string("ste"), directory, m_parameter,
+                         m_optional);
+  save_data_to_file(std::string("coscos"), directory, m_mixing_times, m_cc_dict,
+                    m_optional);
+  save_data_to_file(std::string("sinsin"), directory, m_mixing_times, m_ss_dict,
+                    m_optional);
+  save_data_to_file(std::string("f2"), directory, m_mixing_times, m_f2,
+                    m_optional);
+}
+
+std::unique_ptr<Experiment> STEExperiment::clone() const {
+  return std::make_unique<STEExperiment>(*this);
+}
+
+void STEExperiment::merge(const Experiment &other) {
+  const auto &o = dynamic_cast<const STEExperiment &>(other);
+  for (auto &[t_evo, cc] : m_cc_dict) {
+    const auto &other_cc = o.m_cc_dict.at(t_evo);
+    const auto &other_ss = o.m_ss_dict.at(t_evo);
+    auto &ss = m_ss_dict[t_evo];
+    for (size_t i = 0; i < cc.size(); i++) {
+      cc[i] += other_cc[i];
+      ss[i] += other_ss[i];
+    }
+  }
+  for (size_t i = 0; i < m_f2.size(); i++) {
+    m_f2[i] += o.m_f2[i];
+  }
+}
--- a/src/experiments/ste.h
+++ b/src/experiments/ste.h
@@ -0,0 +1,44 @@
+#ifndef RWSIM_STE_H
+#define RWSIM_STE_H
+
+#include "base.h"
+
+#include <map>
+#include <unordered_map>
+#include <vector>
+
+/// Stimulated echo (STE) experiment.
+///
+/// Computes the stimulated echo decay as a function of mixing time for various
+/// evolution times. The STE probes slow molecular reorientation by comparing
+/// frequencies before and after a mixing period.
+///
+/// Output files:
+///   coscos — cos-cos correlation function F_cc(t_evo, t_mix)
+///   sinsin — sin-sin correlation function F_ss(t_evo, t_mix)
+///   f2     — two-time correlation function <omega(0)*omega(t_mix)>
+class STEExperiment final : public Experiment {
+public:
+  void setup(const std::unordered_map<std::string, double> &parameter,
+             const std::unordered_map<std::string, double> &optional) override;
+  [[nodiscard]] double tmax() const override;
+  void accumulate(const Trajectory &traj, double init_omega,
+                  int num_walkers) override;
+  void save(const std::string &directory) override;
+  [[nodiscard]] std::unique_ptr<Experiment> clone() const override;
+  void merge(const Experiment &other) override;
+
+private:
+  int m_num_mix_times{0};
+  double m_tpulse4{0};                            ///< Duration of the 4th pulse
+  std::vector<double> m_evolution_times;           ///< t_evo values (linear)
+  std::vector<double> m_mixing_times;              ///< t_mix values (logarithmic)
+  std::map<double, std::vector<double>> m_cc_dict; ///< cos-cos per t_evo
+  std::map<double, std::vector<double>> m_ss_dict; ///< sin-sin per t_evo
+  std::vector<double> m_f2;                        ///< Two-time correlator
+  std::unordered_map<std::string, double> m_parameter;
+  std::unordered_map<std::string, double> m_optional;
+  double m_tmax{0};
+};
+
+#endif // RWSIM_STE_H
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,50 +1,70 @@

-#include "utils/io.h"
-#include "simulation/sims.h"
+#include "dynamics/base.h"
+#include "dynamics/decoupled.h"
+#include "experiments/spectrum.h"
+#include "experiments/ste.h"
 #include "motions/base.h"
+#include "sims.h"
 #include "times/base.h"
-
+#include "utils/io.h"

 #include <iostream>
-#include <unordered_map>
 #include <random>
+#include <unordered_map>

+int main(const int argc, char *argv[]) {
+  Arguments args;
+  try {
+    args = parse_args(argc, argv);
+  } catch (std::exception &error) {
+    std::cerr << error.what() << std::endl;
+    return 1;
+  }

+  std::unordered_map parameter{read_parameter(args.parameter_file)};

+  for (const auto &[key, value] : args.optional) {
+    parameter[key] = value;
+  }

-int main (const int argc, char *argv[]) {
-    Arguments args;
-    try {
-        args = parse_args(argc, argv);
-    } catch (std::invalid_argument& error) {
-        std::cerr << error.what() << std::endl;
-        return 1;
-    }
-
-    std::unordered_map parameter { read_parameter(args.parameter_file) };
-
-    for (const auto& [key, value]: args.optional) {
-        parameter[key] = value;
-    }
-
-    // print parameter of simulation to inform user
-    std::cout << "Found parameter\n";
-    for (const auto& [key, value]: parameter) {
-        std::cout << key << ": " << std::to_string(value) << "\n";
-    }
-    std::cout << std::endl;
+  // print parameter of simulation to inform user
+  std::cout << "Found parameter\n";
+  for (const auto &[key, value] : parameter) {
+    std::cout << key << ": " << std::to_string(value) << "\n";
+  }
+  std::cout << std::endl;

+  std::mt19937_64 rng;
+  if (parameter.count("seed")) {
+    rng.seed(static_cast<uint64_t>(parameter.at("seed")));
+  } else {
    std::random_device rd;
-    std::mt19937_64 rng(rd());
+    rng.seed(rd());
+  }

-    Motion *motion = Motion::createFromInput(args.motion_type, rng);
-    Distribution *dist = Distribution::createFromInput(args.distribution_type, rng);
-    if (args.spectrum) {
-       run_spectrum(parameter, args.optional, *motion, *dist);
-    }
-    if (args.ste) {
-        run_ste(parameter, args.optional, *motion, *dist);
-    }
+  std::unique_ptr<Dynamics> dynamics;
+  if (!args.dynamics_type.empty()) {
+    dynamics = Dynamics::createFromInput(args.dynamics_type);
+  } else {
+    dynamics = std::make_unique<DecoupledDynamics>(
+        motions::BaseMotion::createFromInput(args.motion_type),
+        times::BaseDistribution::createFromInput(args.distribution_type));
+  }

-    return 0;
+  SpectrumExperiment spectrum_exp;
+  STEExperiment ste_exp;
+
+  std::vector<Experiment *> experiments;
+  if (args.spectrum) {
+    experiments.push_back(&spectrum_exp);
+  }
+  if (args.ste) {
+    experiments.push_back(&ste_exp);
+  }
+
+  if (!experiments.empty()) {
+    run_simulation(experiments, parameter, args.optional, *dynamics, rng);
+  }
+
+  return 0;
 }
--- a/src/motions/CMakeLists.txt
+++ b/src/motions/CMakeLists.txt
@@ -0,0 +1,15 @@
+add_library(RWMotion OBJECT
+  base.cpp base.h 
+  bimodalangle.cpp bimodalangle.h 
+  conewobble.cpp conewobble.h 
+  conemotion.cpp conemotion.h 
+  coordinates.cpp coordinates.h 
+  diffusivemotion.cpp diffusivemotion.h 
+  foursitejump.cpp foursitejump.h 
+  isosmallangle.cpp isosmallangle.h 
+  nsiteconejump.cpp nsiteconejump.h
+  random.cpp random.h 
+  rjoac.cpp rjoac.h 
+  sixsitejump.cpp sixsitejump.h 
+)
+
--- a/src/motions/base.cpp
+++ b/src/motions/base.cpp
@@ -1,64 +1,57 @@

 #include "base.h"
-#include "coordinates.h"
-#include "bimodalangle.h"
-#include "isosmallangle.h"
-#include "random.h"
-#include "tetrahedral.h"
+#include "../utils/registry.h"

-#include <stdexcept>
+#include <iostream>

-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.);
+namespace motions {
+BaseMotion::BaseMotion(std::string name, const double delta, const double eta)
+    : m_name(std::move(name)), m_delta(delta), m_eta(eta) {}
+
+BaseMotion::BaseMotion(std::string name) : m_name(std::move(name)) {}
+
+// Quadrupolar NMR frequency: omega_q(theta, phi) =
+//   pi * delta * (3*cos^2(theta) - 1 - eta*sin^2(theta)*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;
+
+  return M_PI * m_delta *
+         (3. * cos_theta_square - 1. -
+          m_eta * sin_theta_square * std::cos(2. * phi));
 }

-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.);
+double BaseMotion::omega_q(const coordinates::SphericalPos &pos) const {
+  auto [cos_theta, phi] = pos;
+
+  return omega_q(cos_theta, phi);
 }

-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;
+coordinates::SphericalPos BaseMotion::draw_position(std::mt19937_64 &rng) {
+  const double cos_theta = 1 - 2 * m_uni_dist(rng);
+  const double phi = 2.0 * M_PI * m_uni_dist(rng);

-    return M_PI * m_delta * (3. * cos_theta_square - 1. - m_eta * sin_theta_square * std::cos(2.*phi));
+  return {cos_theta, phi};
 }

-double Motion::omega_q(const SphericalPos& pos) const {
-    auto [cos_theta, phi] = pos;
-
-    return omega_q(cos_theta, phi);
+std::unique_ptr<BaseMotion>
+BaseMotion::createFromInput(const std::string &input) {
+  return Registry<BaseMotion>::create(input);
 }

-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 {cos_theta, phi};
+void BaseMotion::setParameters(
+    const std::unordered_map<std::string, double> &parameters) {
+  m_delta = parameters.at("delta");
+  m_eta = parameters.at("eta");
 }

-Motion* Motion::createFromInput(const std::string& input, std::mt19937_64& rng) {
-    if (input == "TetrahedralJump")
-        return new TetrahedralJump(rng);
-
-    if (input == "IsotropicAngle")
-        return new SmallAngle(rng);
-
-    if (input == "RandomJump")
-        return new RandomJump(rng);
-
-    if (input == "BimodalAngle")
-        return new BimodalAngle(rng);
-
-    throw std::invalid_argument("Invalid input " + input);
+std::unordered_map<std::string, double> BaseMotion::getParameters() const {
+  return std::unordered_map<std::string, double>{{"delta", m_delta},
+                                                 {"eta", m_eta}};
 }

-void Motion::setParameters(const std::unordered_map<std::string, double> &parameters) {
-    m_delta = parameters.at("delta");
-    m_eta = parameters.at("eta");
+std::ostream &operator<<(std::ostream &os, const BaseMotion &m) {
+  os << m.getName();
+  return os;
 }
-
-
-std::ostream& operator<<(std::ostream& os, const Motion& m) {
-    os << m.getName();
-    return os;
-}
+} // namespace motions
--- a/src/motions/base.h
+++ b/src/motions/base.h
@@ -3,43 +3,68 @@

 #include "coordinates.h"

+#include <memory>
 #include <random>
 #include <unordered_map>

-class Motion {
+namespace motions {
+
+/// Base class for NMR motion models.
+///
+/// Each model defines how a molecular orientation evolves over discrete jumps.
+/// The quadrupolar frequency omega_q depends on orientation relative to the
+/// static magnetic field via the second-order Legendre polynomial:
+///   omega_q = pi * delta * (3*cos^2(theta) - 1 - eta*sin^2(theta)*cos(2*phi))
+///
+/// Intermediate base classes:
+///   ConeMotion    — models with a fixed symmetry axis (wobble, random jump on cone)
+///   DiffusiveMotion — models with continuous rotational diffusion (small angle, bimodal)
+class BaseMotion {
 public:
-    virtual ~Motion() = default;
+  virtual ~BaseMotion() = default;

-    Motion(std::string, double, double, std::mt19937_64&);
-    explicit Motion(std::string, std::mt19937_64&);
+  BaseMotion(std::string, double, double);
+  explicit BaseMotion(std::string);

-    SphericalPos draw_position();
-    [[nodiscard]] double omega_q(double, double) const;
-    [[nodiscard]] double omega_q(const SphericalPos&) const;
+  /// Draw a uniform random orientation on the unit sphere.
+  coordinates::SphericalPos draw_position(std::mt19937_64 &rng);

-    virtual void initialize() = 0;
-    virtual double jump() = 0;
+  /// Compute quadrupolar frequency from orientation angles.
+  [[nodiscard]] double omega_q(double, double) const;
+  [[nodiscard]] double omega_q(const coordinates::SphericalPos &) const;

-    virtual void setParameters(const std::unordered_map<std::string, double>&);
+  /// Set up initial state for a new walker (called once per walker).
+  virtual void initialize(std::mt19937_64 &rng) = 0;

-    [[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; };
+  /// Perform one jump, returning the new omega_q.
+  virtual double jump(std::mt19937_64 &rng) = 0;

-    static Motion* createFromInput(const std::string& input, std::mt19937_64& rng);
+  /// Deep copy for per-thread cloning.
+  [[nodiscard]] virtual std::unique_ptr<BaseMotion> clone() const = 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; }
+
+  [[nodiscard]] virtual std::string toString() const = 0;
+
+  /// Factory: create a motion model by registered name.
+  static std::unique_ptr<BaseMotion> createFromInput(const std::string &input);

 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::string m_name{"BaseMotion"};
+  double m_delta{1.};  ///< Quadrupolar coupling constant (Hz)
+  double m_eta{0.};    ///< Asymmetry parameter [0, 1]
+  std::uniform_real_distribution<> m_uni_dist{0., 1.};
+  double m_initial_omega{0.}; ///< Frequency at initial position
 };

-std::ostream& operator<<(std::ostream& os, const Motion& m);
-
-#endif //RWSIM_MOTIONBASE_H
+std::ostream &operator<<(std::ostream &os, const BaseMotion &m);
+} // namespace motions
+#endif // RWSIM_MOTIONBASE_H
--- a/src/motions/bimodalangle.cpp
+++ b/src/motions/bimodalangle.cpp
@@ -1,31 +1,53 @@

 #include "bimodalangle.h"
-#include "base.h"
+#include "../utils/registry.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) {}
+static AutoRegister<motions::BaseMotion, motions::BimodalAngle>
+    reg("BimodalAngle");

-void BimodalAngle::initialize() {
-    m_prev_pos = draw_position();
-    m_initial_omega = omega_q(m_prev_pos);
-};
+namespace motions {
+BimodalAngle::BimodalAngle(const double delta, const double eta,
+                           const double angle1, const double angle2,
+                           const double prob)
+    : DiffusiveMotion(std::string("BimodalAngle"), delta, eta),
+      m_angle1(angle1 * M_PI / 180.0), m_angle2(angle2 * M_PI / 180.0),
+      m_prob(prob) {}
+BimodalAngle::BimodalAngle() : DiffusiveMotion(std::string("BimodalAngle")) {}

-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(std::mt19937_64 &rng) {
+  const double angle = m_uni_dist(rng) < m_prob ? m_angle1 : m_angle2;
+  const double gamma{2 * M_PI * m_uni_dist(rng)};
+  m_prev_pos = rotate(m_prev_pos, angle, gamma);

-    return omega_q(m_prev_pos);
+  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::unique_ptr<BaseMotion> BimodalAngle::clone() const {
+  return std::make_unique<BimodalAngle>(*this);
 }
+
+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["probability1"] = 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)};
+}
+} // namespace motions
--- a/src/motions/bimodalangle.h
+++ b/src/motions/bimodalangle.h
@@ -2,22 +2,25 @@
 #ifndef BIMODALANGLE_H
 #define BIMODALANGLE_H

-#include "base.h"
+#include "diffusivemotion.h"

-class BimodalAngle : public Motion {
+namespace motions {
+class BimodalAngle final : public DiffusiveMotion {
 public:
-    BimodalAngle(double, double, double, double, double, std::mt19937_64& );
-    explicit BimodalAngle(std::mt19937_64&);
+  BimodalAngle(double, double, double, double, double);
+  BimodalAngle();

-    void initialize() override;
-    double jump() override;
-    void setParameters(const std::unordered_map<std::string, double> &) override;
+  double jump(std::mt19937_64 &rng) override;
+  [[nodiscard]] std::unique_ptr<BaseMotion> clone() 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;

 protected:
-    double m_angle1{0};
-    double m_angle2{0};
-    double m_prob{0};
-    SphericalPos m_prev_pos{0., 0.};
+  double m_angle1{0};
+  double m_angle2{0};
+  double m_prob{0};
 };
-
-#endif //BIMODALANGLE_H
+} // namespace motions
+#endif // BIMODALANGLE_H
--- a/src/motions/conemotion.cpp
+++ b/src/motions/conemotion.cpp
@@ -0,0 +1,19 @@
+#include "conemotion.h"
+
+namespace motions {
+void ConeMotion::initialize(std::mt19937_64 &rng) {
+  m_axis = draw_position(rng);
+}
+
+void ConeMotion::setParameters(
+    const std::unordered_map<std::string, double> &parameters) {
+  BaseMotion::setParameters(parameters);
+  m_angle = parameters.at("angle");
+}
+
+std::unordered_map<std::string, double> ConeMotion::getParameters() const {
+  auto parameter = BaseMotion::getParameters();
+  parameter["angle"] = m_angle;
+  return parameter;
+}
+} // namespace motions
--- a/src/motions/conemotion.h
+++ b/src/motions/conemotion.h
@@ -0,0 +1,29 @@
+#ifndef CONEMOTION_H
+#define CONEMOTION_H
+
+#include "base.h"
+#include "coordinates.h"
+
+namespace motions {
+
+/// Intermediate base for motions with a fixed symmetry axis and cone angle.
+/// Provides shared axis initialization and angle parameter handling.
+/// Used by WobbleCone and RandomJumpOnCone.
+class ConeMotion : public BaseMotion {
+public:
+  using BaseMotion::BaseMotion;
+
+  /// Draws a random symmetry axis on the unit sphere.
+  void initialize(std::mt19937_64 &rng) override;
+
+  void setParameters(const std::unordered_map<std::string, double> &) override;
+  [[nodiscard]] std::unordered_map<std::string, double>
+  getParameters() const override;
+
+protected:
+  double m_angle{0};                       ///< Half-opening angle of the cone (radians)
+  coordinates::SphericalPos m_axis{1, 0};  ///< Symmetry axis orientation
+};
+} // namespace motions
+
+#endif // CONEMOTION_H
--- a/src/motions/conewobble.cpp
+++ b/src/motions/conewobble.cpp
@@ -0,0 +1,28 @@
+
+#include "conewobble.h"
+#include "../utils/registry.h"
+#include "coordinates.h"
+
+static AutoRegister<motions::BaseMotion, motions::WobbleCone> reg("ConeWobble");
+
+namespace motions {
+WobbleCone::WobbleCone(const double delta, const double eta, const double chi)
+    : ConeMotion("Wobble in Cone", delta, eta) {
+  m_angle = chi;
+}
+WobbleCone::WobbleCone() : ConeMotion("Wobble in Cone") {}
+
+double WobbleCone::jump(std::mt19937_64 &rng) {
+  const double real_angle = m_uni_dist(rng) * m_angle;
+  const double phi = 2 * M_PI * m_uni_dist(rng);
+  return omega_q(rotate(m_axis, real_angle, phi));
+}
+
+std::unique_ptr<BaseMotion> WobbleCone::clone() const {
+  return std::make_unique<WobbleCone>(*this);
+}
+
+std::string WobbleCone::toString() const {
+  return std::string("ConeWobble/angle=") + std::to_string(m_angle);
+}
+} // namespace motions
--- a/src/motions/conewobble.h
+++ b/src/motions/conewobble.h
@@ -0,0 +1,17 @@
+#ifndef CONEWOBBLE_H
+#define CONEWOBBLE_H
+
+#include "conemotion.h"
+
+namespace motions {
+class WobbleCone final : public ConeMotion {
+public:
+  WobbleCone(double, double, double);
+  WobbleCone();
+
+  double jump(std::mt19937_64 &rng) override;
+  [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
+  [[nodiscard]] std::string toString() const override;
+};
+} // namespace motions
+#endif // CONEWOBBLE_H
--- a/src/motions/coordinates.cpp
+++ b/src/motions/coordinates.cpp
@@ -1,39 +1,43 @@
 #include "coordinates.h"

 #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};
+  const double norm{std::sqrt(1 - cos_theta * cos_theta) + 1e-15};

-    auto [x, y , z] = spherical_to_xyz(old_pos);
+  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
-    };
+  const auto new_pos = CartesianPos{
+      cos_alpha * x + sin_alpha * (-cos_beta * y - sin_beta * x * z) / norm,
+      cos_alpha * y + sin_alpha * (cos_beta * x - sin_beta * y * z) / norm,
+      cos_alpha * z + sin_alpha * sin_beta * norm,
+  };

-    return xyz_to_spherical(new_pos);
+  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};
+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)};
 }
+} // namespace coordinates
--- a/src/motions/coordinates.h
+++ b/src/motions/coordinates.h
@@ -2,19 +2,27 @@
 #ifndef COORDINATES_H
 #define COORDINATES_H

+/// Coordinate types and rotations for molecular orientations on the unit sphere.
+namespace coordinates {
+
+/// Spherical position stored as (cos_theta, phi) to avoid repeated cos/acos.
 struct SphericalPos {
-    double cos_theta;
-    double phi;
+  double cos_theta;
+  double phi;
 };

 struct CartesianPos {
-    double x;
-    double y;
-    double z;
+  double x;
+  double y;
+  double z;
 };

-SphericalPos rotate(const SphericalPos&, double, double);
-CartesianPos spherical_to_xyz(const SphericalPos&);
-SphericalPos xyz_to_spherical(const CartesianPos&);
+/// Rotate a spherical position by polar angle alpha around an axis defined
+/// by azimuthal angle beta (rotation about the original position's z-axis).
+SphericalPos rotate(const SphericalPos &, double alpha, double beta);

-#endif //COORDINATES_H
+CartesianPos spherical_to_xyz(const SphericalPos &);
+SphericalPos xyz_to_spherical(const CartesianPos &);
+} // namespace coordinates
+
+#endif // COORDINATES_H
--- a/src/motions/diffusivemotion.cpp
+++ b/src/motions/diffusivemotion.cpp
@@ -0,0 +1,8 @@
+#include "diffusivemotion.h"
+
+namespace motions {
+void DiffusiveMotion::initialize(std::mt19937_64 &rng) {
+  m_prev_pos = draw_position(rng);
+  m_initial_omega = omega_q(m_prev_pos);
+}
+} // namespace motions
--- a/src/motions/diffusivemotion.h
+++ b/src/motions/diffusivemotion.h
@@ -0,0 +1,24 @@
+#ifndef DIFFUSIVEMOTION_H
+#define DIFFUSIVEMOTION_H
+
+#include "base.h"
+#include "coordinates.h"
+
+namespace motions {
+
+/// Intermediate base for motions with continuous rotational diffusion.
+/// Tracks the previous position and provides shared initialization.
+/// Used by SmallAngle and BimodalAngle.
+class DiffusiveMotion : public BaseMotion {
+public:
+  using BaseMotion::BaseMotion;
+
+  /// Draws a random initial position and sets the initial frequency.
+  void initialize(std::mt19937_64 &rng) override;
+
+protected:
+  coordinates::SphericalPos m_prev_pos{0., 0.}; ///< Current orientation
+};
+} // namespace motions
+
+#endif // DIFFUSIVEMOTION_H
--- a/src/motions/foursitejump.cpp
+++ b/src/motions/foursitejump.cpp
@@ -0,0 +1,42 @@
+#include "foursitejump.h"
+#include "../utils/registry.h"
+#include "coordinates.h"
+
+static AutoRegister<motions::BaseMotion, motions::FourSiteTetrahedron>
+    reg("FourSiteTetrahedral");
+
+namespace motions {
+FourSiteTetrahedron::FourSiteTetrahedron(const double delta, const double eta)
+    : BaseMotion(std::string{"FourSiteTetrahedral"}, delta, eta) {}
+
+FourSiteTetrahedron::FourSiteTetrahedron()
+    : BaseMotion(std::string{"FourSiteTetrahedral"}) {}
+
+void FourSiteTetrahedron::initialize(std::mt19937_64 &rng) {
+  const auto pos = draw_position(rng);
+  m_corners[0] = omega_q(pos);
+
+  const double alpha = 2. * M_PI * m_uni_dist(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 = FourSiteTetrahedron::jump(rng);
+}
+
+double FourSiteTetrahedron::jump(std::mt19937_64 &rng) {
+  m_corner_idx += m_chooser(rng);
+  m_corner_idx %= 4;
+
+  return m_corners[m_corner_idx];
+}
+
+std::unique_ptr<BaseMotion> FourSiteTetrahedron::clone() const {
+  return std::make_unique<FourSiteTetrahedron>(*this);
+}
+
+std::string FourSiteTetrahedron::toString() const {
+  return {"FourSiteTetrahedral"};
+}
+} // namespace motions
--- a/src/motions/foursitejump.h
+++ b/src/motions/foursitejump.h
@@ -0,0 +1,30 @@
+#ifndef RWSIM_MOTIONTETRAHEDRAL_H
+#define RWSIM_MOTIONTETRAHEDRAL_H
+
+#include "base.h"
+
+#include <array>
+#include <cmath>
+
+namespace motions {
+class FourSiteTetrahedron final : public BaseMotion {
+public:
+  FourSiteTetrahedron(double, double);
+  FourSiteTetrahedron();
+
+  void initialize(std::mt19937_64 &rng) override;
+  double jump(std::mt19937_64 &rng) override;
+  [[nodiscard]] std::unique_ptr<BaseMotion> clone() const 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};
+};
+} // namespace motions
+#endif // RWSIM_MOTIONTETRAHEDRAL_H
--- a/src/motions/isosmallangle.cpp
+++ b/src/motions/isosmallangle.cpp
@@ -1,27 +1,41 @@
 #include "isosmallangle.h"
+#include "../utils/registry.h"
 #include "coordinates.h"

-#include <iostream>
+static AutoRegister<motions::BaseMotion, motions::SmallAngle>
+    reg("IsotropicAngle");

+namespace motions {
+SmallAngle::SmallAngle(const double delta, const double eta, const double chi)
+    : DiffusiveMotion(std::string("IsotropicAngle"), delta, eta),
+      m_chi(chi * M_PI / 180.0) {}
+SmallAngle::SmallAngle() : DiffusiveMotion(std::string("IsotropicAngle")) {}

+double SmallAngle::jump(std::mt19937_64 &rng) {
+  const double gamma{2 * M_PI * m_uni_dist(rng)};
+  m_prev_pos = rotate(m_prev_pos, m_chi, gamma);

-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);
-};
-
-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::unique_ptr<BaseMotion> SmallAngle::clone() const {
+  return std::make_unique<SmallAngle>(*this);
 }
+
+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)};
+}
+} // namespace motions
--- a/src/motions/isosmallangle.h
+++ b/src/motions/isosmallangle.h
@@ -2,21 +2,23 @@
 #ifndef RWSIM_MOTIONISOSMALLANGLE_H
 #define RWSIM_MOTIONISOSMALLANGLE_H

-#include "base.h"
-#include "coordinates.h"
+#include "diffusivemotion.h"

-class SmallAngle final : public Motion {
+namespace motions {
+class SmallAngle final : public DiffusiveMotion {
 public:
-    SmallAngle(double, double, double, std::mt19937_64& );
-    explicit SmallAngle(std::mt19937_64&);
+  SmallAngle(double, double, double);
+  SmallAngle();

-    void initialize() override;
-    double jump() override;
-    void setParameters(const std::unordered_map<std::string, double> &) override;
+  double jump(std::mt19937_64 &rng) override;
+  [[nodiscard]] std::unique_ptr<BaseMotion> clone() 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;

 private:
-    double m_chi{0};
-    SphericalPos m_prev_pos{0., 0.};
+  double m_chi{0};
 };
-
-#endif //RWSIM_MOTIONISOSMALLANGLE_H
+} // namespace motions
+#endif // RWSIM_MOTIONISOSMALLANGLE_H
--- a/src/motions/nsiteconejump.cpp
+++ b/src/motions/nsiteconejump.cpp
@@ -0,0 +1,63 @@
+#include "nsiteconejump.h"
+#include "../utils/registry.h"
+#include "coordinates.h"
+
+#include <cmath>
+#include <vector>
+
+static AutoRegister<motions::BaseMotion, motions::NSiteJumpOnCone>
+    reg("NSiteConeJump");
+
+namespace motions {
+NSiteJumpOnCone::NSiteJumpOnCone(const double delta, const double eta,
+                                 const int num_sites, const double chi)
+    : BaseMotion("NSiteJumpOnCone", delta, eta), m_num_sites(num_sites),
+      m_chi(chi * M_PI / 180.) {}
+
+NSiteJumpOnCone::NSiteJumpOnCone() : BaseMotion("NSiteJumpOnCone") {}
+
+void NSiteJumpOnCone::initialize(std::mt19937_64 &rng) {
+  m_sites = std::vector<double>(m_num_sites);
+  m_chooser = std::uniform_int_distribution<>{1, m_num_sites - 1};
+
+  m_axis = draw_position(rng);
+
+  const double alpha = m_uni_dist(rng) * 2 * M_PI;
+  const double steps = 2 * M_PI / m_num_sites;
+  for (int i = 0; i < m_num_sites; i++) {
+    m_sites[i] = omega_q(rotate(m_axis, m_chi, i * steps + alpha));
+  }
+}
+
+double NSiteJumpOnCone::jump(std::mt19937_64 &rng) {
+  m_cone_idx += m_chooser(rng);
+  m_cone_idx %= m_num_sites;
+
+  return m_sites[m_cone_idx];
+}
+
+std::unique_ptr<BaseMotion> NSiteJumpOnCone::clone() const {
+  return std::make_unique<NSiteJumpOnCone>(*this);
+}
+
+void NSiteJumpOnCone::setParameters(
+    const std::unordered_map<std::string, double> &parameters) {
+  BaseMotion::setParameters(parameters);
+  m_num_sites = static_cast<int>(parameters.at("num_sites"));
+  m_chi = parameters.at("chi") * M_PI / 180.;
+}
+
+std::unordered_map<std::string, double> NSiteJumpOnCone::getParameters() const {
+  auto parameter = BaseMotion::getParameters();
+  parameter["num_sites"] = m_num_sites;
+  parameter["chi"] = m_chi * 180. / M_PI;
+
+  return parameter;
+}
+
+std::string NSiteJumpOnCone::toString() const {
+  return std::to_string(m_num_sites) +
+         "SiteJumpOnCone/angle=" + std::to_string(m_chi * 180 / M_PI);
+}
+
+} // namespace motions
--- a/src/motions/nsiteconejump.h
+++ b/src/motions/nsiteconejump.h
@@ -0,0 +1,34 @@
+#ifndef NSITEJUMPONCONE_H
+#define NSITEJUMPONCONE_H
+
+#include "base.h"
+
+#include <vector>
+
+namespace motions {
+class NSiteJumpOnCone final : public BaseMotion {
+public:
+  NSiteJumpOnCone(double, double, int, double);
+  NSiteJumpOnCone();
+
+  void initialize(std::mt19937_64 &rng) override;
+  double jump(std::mt19937_64 &rng) override;
+  [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
+
+  [[nodiscard]] std::string toString() const override;
+
+  void setParameters(const std::unordered_map<std::string, double> &) override;
+  [[nodiscard]] std::unordered_map<std::string, double>
+  getParameters() const override;
+
+private:
+  int m_num_sites{2};
+  std::vector<double> m_sites{};
+  double m_chi{1. / 2.};
+  int m_cone_idx = 0;
+  coordinates::SphericalPos m_axis{1, 0};
+  std::uniform_int_distribution<> m_chooser;
+};
+} // namespace motions
+
+#endif // NSITEJUMPONCONE_H
--- a/src/motions/random.cpp
+++ b/src/motions/random.cpp
@@ -1,15 +1,25 @@
-
 #include "random.h"
+#include "../utils/registry.h"

+static AutoRegister<motions::BaseMotion, motions::RandomJump> reg("RandomJump");

-RandomJump::RandomJump(const double delta, const double eta, std::mt19937_64 &rng) : Motion(std::string("RandomJump"), delta, eta, rng) {}
+namespace motions {
+RandomJump::RandomJump(const double delta, const double eta)
+    : BaseMotion(std::string("RandomJump"), delta, eta) {}

-RandomJump::RandomJump(std::mt19937_64 &rng) : Motion(std::string("RandomJump"), rng) {}
+RandomJump::RandomJump() : BaseMotion(std::string("RandomJump")) {}

-void RandomJump::initialize() {
-    m_initial_omega = RandomJump::jump();
+std::string RandomJump::toString() const { return {"RandomJump"}; }
+
+std::unique_ptr<BaseMotion> RandomJump::clone() const {
+  return std::make_unique<RandomJump>(*this);
 }

-double RandomJump::jump() {
-    return omega_q(draw_position());
+void RandomJump::initialize(std::mt19937_64 &rng) {
+  m_initial_omega = RandomJump::jump(rng);
 }
+
+double RandomJump::jump(std::mt19937_64 &rng) {
+  return omega_q(draw_position(rng));
+}
+} // namespace motions
--- a/src/motions/random.h
+++ b/src/motions/random.h
@@ -1,17 +1,20 @@
 #ifndef RWSIM_MOTIONRANDOMJUMP_H
 #define RWSIM_MOTIONRANDOMJUMP_H

-
 #include "base.h"
-#include <random>

-class RandomJump final : public Motion {
+namespace motions {
+class RandomJump final : public BaseMotion {
 public:
-    RandomJump(double, double, std::mt19937_64&);
-    explicit RandomJump(std::mt19937_64&);
+  RandomJump(double, double);
+  RandomJump();

-    void initialize() override;
-    double jump() override;
+  [[nodiscard]] std::string toString() const override;
+  [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
+
+  void initialize(std::mt19937_64 &rng) override;
+  double jump(std::mt19937_64 &rng) override;
 };
+} // namespace motions

-#endif //RWSIM_MOTIONRANDOMJUMP_H
+#endif // RWSIM_MOTIONRANDOMJUMP_H
--- a/src/motions/rjoac.cpp
+++ b/src/motions/rjoac.cpp
@@ -0,0 +1,28 @@
+#include "rjoac.h"
+#include "../utils/registry.h"
+#include "coordinates.h"
+
+static AutoRegister<motions::BaseMotion, motions::RandomJumpOnCone>
+    reg("RandomJumpOnCone");
+
+namespace motions {
+RandomJumpOnCone::RandomJumpOnCone(const double delta, const double eta,
+                                   const double chi)
+    : ConeMotion("RJ on a Cone", delta, eta) {
+  m_angle = chi;
+}
+RandomJumpOnCone::RandomJumpOnCone() : ConeMotion("RJ on a Cone") {}
+
+double RandomJumpOnCone::jump(std::mt19937_64 &rng) {
+  const double phi = 2 * M_PI * m_uni_dist(rng);
+  return omega_q(rotate(m_axis, m_angle, phi));
+}
+
+std::unique_ptr<BaseMotion> RandomJumpOnCone::clone() const {
+  return std::make_unique<RandomJumpOnCone>(*this);
+}
+
+std::string RandomJumpOnCone::toString() const {
+  return std::string("RandomJumpOnCone/angle=") + std::to_string(m_angle);
+}
+} // namespace motions
--- a/src/motions/rjoac.h
+++ b/src/motions/rjoac.h
@@ -0,0 +1,18 @@
+#ifndef RJOAC_H
+#define RJOAC_H
+
+#include "conemotion.h"
+
+namespace motions {
+class RandomJumpOnCone final : public ConeMotion {
+public:
+  RandomJumpOnCone(double, double, double);
+  RandomJumpOnCone();
+
+  double jump(std::mt19937_64 &rng) override;
+  [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
+  [[nodiscard]] std::string toString() const override;
+};
+} // namespace motions
+
+#endif // RJOAC_H
--- a/src/motions/sixsitejump.cpp
+++ b/src/motions/sixsitejump.cpp
@@ -0,0 +1,47 @@
+#include "sixsitejump.h"
+#include "../utils/registry.h"
+#include "coordinates.h"
+
+static AutoRegister<motions::BaseMotion, motions::SixSiteOctahedronC3>
+    reg("SixSiteOctahedralC3");
+
+namespace motions {
+SixSiteOctahedronC3::SixSiteOctahedronC3(const double delta, const double eta,
+                                         const double chi)
+    : BaseMotion(std::string{"SixSiteOctahedral"}, delta, eta),
+      m_chi{chi * M_PI / 180.} {}
+
+SixSiteOctahedronC3::SixSiteOctahedronC3()
+    : BaseMotion(std::string{"SixSiteOctahedralC3"}) {}
+
+void SixSiteOctahedronC3::initialize(std::mt19937_64 &rng) {
+  const coordinates::SphericalPos c3_axis = draw_position(rng);
+  const double alpha = 2. * M_PI * m_uni_dist(rng);
+
+  const double m_chi_opposite = M_PI - m_chi;
+
+  for (int i = 0; i < 3; i++) {
+    m_corners[2 * i] =
+        omega_q(rotate(c3_axis, m_chi, alpha + i * 2. / 3. * M_PI));
+    m_corners[2 * i + 1] = omega_q(rotate(
+        c3_axis, m_chi_opposite, alpha + i * 2. / 3. * M_PI + M_PI / 3.));
+  }
+
+  m_initial_omega = SixSiteOctahedronC3::jump(rng);
+}
+
+double SixSiteOctahedronC3::jump(std::mt19937_64 &rng) {
+  m_corner_idx += m_chooser(rng);
+  m_corner_idx %= 6;
+
+  return m_corners[m_corner_idx];
+}
+
+std::unique_ptr<BaseMotion> SixSiteOctahedronC3::clone() const {
+  return std::make_unique<SixSiteOctahedronC3>(*this);
+}
+
+std::string SixSiteOctahedronC3::toString() const {
+  return {"SixSiteOctahedral/angle=" + std::to_string(m_chi / M_PI * 180.)};
+}
+} // namespace motions
--- a/src/motions/sixsitejump.h
+++ b/src/motions/sixsitejump.h
@@ -0,0 +1,30 @@
+#ifndef SIXSITEJUMP_H
+#define SIXSITEJUMP_H
+
+#include "base.h"
+
+#include <array>
+#include <cmath>
+
+namespace motions {
+class SixSiteOctahedronC3 final : public BaseMotion {
+public:
+  SixSiteOctahedronC3(double, double, double);
+  SixSiteOctahedronC3();
+
+  void initialize(std::mt19937_64 &rng) override;
+  double jump(std::mt19937_64 &rng) override;
+  [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
+
+  [[nodiscard]] std::string toString() const override;
+
+private:
+  double m_chi{std::acos(-1.0 / 3.0)}; // 54.74 deg
+
+  std::array<double, 6> m_corners{};
+  int m_corner_idx{0};
+
+  std::uniform_int_distribution<> m_chooser{1, 5};
+};
+} // namespace motions
+#endif // SIXSITEJUMP_H
--- a/src/motions/tetrahedral.cpp
+++ b/src/motions/tetrahedral.cpp
@@ -1,30 +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];
-}
--- a/src/motions/tetrahedral.h
+++ b/src/motions/tetrahedral.h
@@ -1,28 +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;
-
-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
--- a/src/sims.cpp
+++ b/src/sims.cpp
@@ -0,0 +1,151 @@
+#include "sims.h"
+#include "utils/functions.h"
+#include "utils/io.h"
+
+#include <algorithm>
+#include <chrono>
+#include <iostream>
+#include <omp.h>
+
+void run_simulation(std::vector<Experiment *> experiments,
+                    const std::unordered_map<std::string, double> &parameter,
+                    const std::unordered_map<std::string, double> &optional,
+                    Dynamics &dynamics, std::mt19937_64 &rng) {
+  const int num_walker = static_cast<int>(parameter.at("num_walker"));
+
+  dynamics.setParameters(parameter);
+
+  double tmax = 0;
+  for (auto *exp : experiments) {
+    exp->setup(parameter, optional);
+    tmax = std::max(tmax, exp->tmax());
+  }
+
+  const auto start = printStart(optional);
+
+  const int num_threads = omp_get_max_threads();
+
+  // Create per-thread RNGs seeded deterministically from the main RNG
+  std::vector<std::mt19937_64> thread_rngs;
+  thread_rngs.reserve(num_threads);
+  for (int i = 0; i < num_threads; i++) {
+    thread_rngs.emplace_back(rng());
+  }
+
+  // Create per-thread clones of dynamics and experiments
+  std::vector<std::unique_ptr<Dynamics>> thread_dynamics;
+  std::vector<std::vector<std::unique_ptr<Experiment>>> thread_experiments(
+      num_threads);
+  for (int i = 0; i < num_threads; i++) {
+    thread_dynamics.push_back(dynamics.clone());
+    for (auto *exp : experiments) {
+      thread_experiments[i].push_back(exp->clone());
+    }
+  }
+
+  int steps_done = 0;
+  auto last_print_out = std::chrono::system_clock::now();
+
+#pragma omp parallel
+  {
+    const int tid = omp_get_thread_num();
+    auto &local_rng = thread_rngs[tid];
+    auto &local_dynamics = *thread_dynamics[tid];
+    auto &local_experiments = thread_experiments[tid];
+
+    Trajectory traj;
+
+#pragma omp for schedule(dynamic, 16)
+    for (int mol_i = 0; mol_i < num_walker; mol_i++) {
+      make_trajectory(traj, local_dynamics, tmax, local_rng);
+      const double init_omega = local_dynamics.getInitOmega();
+
+      for (auto &exp : local_experiments) {
+        exp->accumulate(traj, init_omega, num_walker);
+      }
+
+      if (tid == 0) {
+#pragma omp atomic
+        steps_done++;
+        last_print_out =
+            printSteps(last_print_out, start, num_walker, steps_done);
+      } else {
+#pragma omp atomic
+        steps_done++;
+      }
+    }
+  }
+
+  // Merge per-thread results and save
+  const auto directory = make_directory(dynamics.toString());
+  for (size_t e = 0; e < experiments.size(); e++) {
+    for (int i = 0; i < num_threads; i++) {
+      experiments[e]->merge(*thread_experiments[i][e]);
+    }
+    experiments[e]->save(directory);
+  }
+
+  printEnd(start);
+}
+
+void make_trajectory(Trajectory &traj, Dynamics &dynamics, const double t_max,
+                     std::mt19937_64 &rng) {
+  traj.time.clear();
+  traj.phase.clear();
+  traj.omega.clear();
+
+  // Reserve is only useful on first call per thread; after that
+  // clear() preserves capacity from the previous walker
+  if (traj.time.capacity() == 0) {
+    const size_t estimate = 1024;
+    traj.time.reserve(estimate);
+    traj.phase.reserve(estimate);
+    traj.omega.reserve(estimate);
+  }
+
+  double t_passed = 0;
+  double phase = 0;
+
+  dynamics.initialize(rng);
+
+  double omega = dynamics.getInitOmega();
+
+  traj.time.emplace_back(t_passed);
+  traj.phase.emplace_back(phase);
+  traj.omega.emplace_back(omega);
+
+  while (t_passed < t_max) {
+    auto [dt, new_omega] = dynamics.next(rng);
+    phase += omega * dt;
+    t_passed += dt;
+    omega = new_omega;
+
+    traj.time.emplace_back(t_passed);
+    traj.phase.emplace_back(phase);
+    traj.omega.emplace_back(omega);
+  }
+}
+
+std::chrono::system_clock::time_point
+printStart(const std::unordered_map<std::string, double> &optional) {
+  const auto start = std::chrono::system_clock::now();
+  const time_t start_time = std::chrono::system_clock::to_time_t(start);
+
+  std::cout << "Random walk for ";
+  for (const auto &[key, value] : optional) {
+    std::cout << key << " = " << value << "; ";
+  }
+  std::cout << std::endl;
+  std::cout << "Start: " << ctime(&start_time);
+
+  return start;
+}
+
+void printEnd(const std::chrono::system_clock::time_point start) {
+  const auto end = std::chrono::system_clock::now();
+
+  const std::chrono::duration<float> duration = end - start;
+  const time_t end_time = std::chrono::system_clock::to_time_t(end);
+  std::cout << "End:   " << ctime(&end_time);
+  std::cout << "Duration: " << duration.count() << "s\n" << std::endl;
+}
--- a/src/sims.h
+++ b/src/sims.h
@@ -0,0 +1,32 @@
+#ifndef RWSIM_SIMS_H
+#define RWSIM_SIMS_H
+
+#include "dynamics/base.h"
+#include "experiments/base.h"
+
+#include <chrono>
+#include <random>
+#include <string>
+#include <unordered_map>
+#include <vector>
+
+/// Run a parallel random walk simulation over multiple experiments.
+///
+/// Generates one trajectory per walker (using the longest tmax across all
+/// experiments) and accumulates results into each experiment. Parallelized
+/// with OpenMP using per-thread clones of dynamics and experiments.
+void run_simulation(std::vector<Experiment *> experiments,
+                    const std::unordered_map<std::string, double> &parameter,
+                    const std::unordered_map<std::string, double> &optional,
+                    Dynamics &dynamics, std::mt19937_64 &rng);
+
+/// Generate a trajectory by repeatedly calling dynamics.next() until t_max.
+/// Reuses the trajectory's allocated memory across calls (clear preserves capacity).
+void make_trajectory(Trajectory &traj, Dynamics &dynamics, double t_max,
+                     std::mt19937_64 &rng);
+
+std::chrono::system_clock::time_point
+printStart(const std::unordered_map<std::string, double> &optional);
+void printEnd(std::chrono::system_clock::time_point start);
+
+#endif // RWSIM_SIMS_H
--- a/src/simulation/sims.cpp
+++ b/src/simulation/sims.cpp
@@ -1,225 +0,0 @@
-#include "sims.h"
-#include "../motions/base.h"
-#include "../times/base.h"
-#include "../utils/functions.h"
-#include "../utils/ranges.h"
-#include "../utils/io.h"
-
-#include <iostream>
-#include <algorithm>
-#include <unordered_map>
-#include <map>
-#include <string>
-#include <vector>
-#include <cmath>
-#include <chrono>
-
-
-void run_spectrum(
-    std::unordered_map<std::string, double>& parameter,
-    std::unordered_map<std::string, double> optional,
-    Motion& motion,
-    Distribution& dist
-    ) {
-    const int num_walker = static_cast<int>(parameter["num_walker"]);
-
-    // time axis for all time signals
-    const int num_acq = static_cast<int>(parameter["num_acq"]);
-    const std::vector<double> t_fid = arange(num_acq, parameter["dwell_time"]);
-    const std::vector<double> echo_times = linspace(parameter["techo_start"], parameter["techo_stop"], static_cast<int>(parameter["techo_steps"]));
-
-    // make timesignal vectors and set them to zero
-    std::map<double, std::vector<double>> fid_dict;
-    for (auto t_echo_i: echo_times) {
-        fid_dict[t_echo_i] = std::vector<double>(num_acq);
-        std::fill(fid_dict[t_echo_i].begin(), fid_dict[t_echo_i].end(), 0.);
-    }
-
-    // calculate min length of a trajectory
-    const double tmax = *std::max_element(echo_times.begin(), echo_times.end()) * 2 + t_fid.back();
-
-    // set parameter in distribution and motion model
-    dist.setParameters(parameter);
-    motion.setParameters(parameter);
-
-    const auto start = printStart(optional);
-    auto last_print_out = std::chrono::system_clock::now();
-
-    // let the walker walk
-    for (int mol_i = 0; mol_i < num_walker; mol_i++){
-        std::vector<double> traj_time{};
-        std::vector<double> traj_phase{};
-        std::vector<double> traj_omega{};
-
-        make_trajectory(motion, dist, tmax, traj_time, traj_phase, traj_omega);
-
-        for (auto& [t_echo_j, fid_j] : fid_dict) {
-            // get phase at echo pulse
-            int current_pos = nearest_index(traj_time, t_echo_j, 0);
-            const double phase_techo = lerp(traj_time, traj_phase, t_echo_j, current_pos);
-
-
-            for (int acq_idx = 0; acq_idx < num_acq; acq_idx++) {
-                const double real_time = t_fid[acq_idx] + 2 * t_echo_j;
-
-                current_pos = nearest_index(traj_time, real_time, current_pos);
-                const double phase_acq = lerp(traj_time, traj_phase, real_time, current_pos);
-
-                fid_j[acq_idx] += std::cos(phase_acq - 2 * phase_techo) / num_walker;
-            }
-            last_print_out = printSteps(last_print_out, start, num_walker, mol_i);
-        }
-    }
-
-    // write fid to files
-    save_parameter_to_file("timesignal", motion.getName(), dist.getName(), parameter, optional);
-    save_data_to_file("timesignal", motion.getName(), dist.getName(), t_fid, fid_dict, optional);
-
-    printEnd(start);
-}
-
-
-void run_ste(
-    std::unordered_map<std::string, double>& parameter,
-    std::unordered_map<std::string, double> optional,
-    Motion& motion,
-    Distribution& dist
-    ) {
-    const int num_walker = static_cast<int>(parameter[std::string("num_walker")]);
-
-    const int num_mix_times = static_cast<int>(parameter[std::string("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);
-
-    // make ste decay vectors and set them to zero
-    std::map<double, std::vector<double>> cc_dict;
-    std::map<double, std::vector<double>> ss_dict;
-    for (auto t_evo_i: evolution_times) {
-        cc_dict[t_evo_i] = std::vector<double>(num_mix_times);
-        ss_dict[t_evo_i] = std::vector<double>(num_mix_times);
-        std::fill(cc_dict[t_evo_i].begin(), cc_dict[t_evo_i].end(), 0.);
-        std::fill(ss_dict[t_evo_i].begin(), ss_dict[t_evo_i].end(), 0.);
-    }
-    std::vector<double> f2(num_mix_times);
-
-    // each trajectory must have a duration of at least tmax
-    const double tmax = *std::max_element(evolution_times.begin(), evolution_times.end()) * 2 + *std::max_element(mixing_times.begin(), mixing_times.end());
-
-    // set parameter in distribution and motion model
-    dist.setParameters(parameter);
-    motion.setParameters(parameter);
-
-    const auto start = printStart(optional);
-    auto last_print_out = std::chrono::system_clock::now();
-
-    // let the walker walk
-    for (int mol_i = 0; mol_i < num_walker; mol_i++){
-        std::vector<double> traj_time{};
-        std::vector<double> traj_phase{};
-        std::vector<double> traj_omega{};
-
-        make_trajectory(motion, dist, tmax, traj_time, traj_phase, traj_omega);
-
-        int f2_pos = 0;
-        for (int f2_idx=0; f2_idx < num_mix_times; f2_idx++) {
-            const double t_mix_f2 = mixing_times[f2_idx];
-            f2_pos = nearest_index(traj_time, t_mix_f2, f2_pos);
-            f2[f2_idx] += traj_omega[f2_pos] * motion.getInitOmega() / num_walker;
-        }
-
-        for (auto& [t_evo_j, _] : cc_dict) {
-            auto& cc_j = cc_dict[t_evo_j];
-            auto& ss_j = ss_dict[t_evo_j];
-
-            // get phase at beginning of mixing time
-            int current_pos = nearest_index(traj_time, t_evo_j, 0);
-            const double dephased = lerp(traj_time, traj_phase, t_evo_j, current_pos);
-            const double cc_tevo = std::cos(dephased);
-            const double ss_tevo = std::sin(dephased);
-
-            for (int mix_idx = 0; mix_idx < num_mix_times; mix_idx++) {
-
-                // get phase at end of mixing time
-                const double time_end_mix = mixing_times[mix_idx] + t_evo_j;
-                current_pos = nearest_index(traj_time, time_end_mix, current_pos);
-                const double phase_mix_end = lerp(traj_time, traj_phase, time_end_mix, current_pos);
-
-                // get phase at echo position
-                const double time_echo = mixing_times[mix_idx] + 2 * t_evo_j;
-                current_pos = nearest_index(traj_time, time_echo, current_pos);
-                const double rephased = lerp(traj_time, traj_phase, time_echo, current_pos) - phase_mix_end;
-
-                cc_j[mix_idx] += cc_tevo * std::cos(rephased) / num_walker;
-                ss_j[mix_idx] += ss_tevo * std::sin(rephased) / num_walker;
-            }
-        }
-        last_print_out = printSteps(last_print_out, start, num_walker, mol_i);
-    }
-
-    // write to files
-    save_parameter_to_file("ste", motion.getName(), dist.getName(), parameter, optional);
-    save_data_to_file("coscos", motion.getName(), dist.getName(), mixing_times, cc_dict, optional);
-    save_data_to_file("sinsin", motion.getName(), dist.getName(), mixing_times, ss_dict, optional);
-    save_data_to_file("f2", motion.getName(), dist.getName(), mixing_times, f2, optional);
-
-    printEnd(start);
-}
-
-
-void make_trajectory(
-    Motion& motion,
-    Distribution& dist,
-    const double t_max,
-    std::vector<double>& out_time,
-    std::vector<double>& out_phase,
-    std::vector<double>& out_omega
-    ) {
-    // Starting position
-    double t_passed = 0;
-    double phase = 0;
-
-    motion.initialize();
-    dist.initialize();
-
-    double omega = motion.getInitOmega();
-
-    out_time.emplace_back(t_passed);
-    out_phase.emplace_back(phase);
-    out_omega.emplace_back(omega);
-
-    while (t_passed < t_max) {
-        const double t = dist.tau_wait();
-        t_passed += t;
-        omega = motion.jump();
-        phase += omega * t;
-
-        out_time.emplace_back(t_passed);
-        out_phase.emplace_back(phase);
-        out_omega.emplace_back(omega);
-    }
-}
-
-
-std::chrono::system_clock::time_point printStart(std::unordered_map<std::string, double> &optional) {
-    const auto start = std::chrono::system_clock::now();
-    const time_t start_time = std::chrono::system_clock::to_time_t(start);
-
-    std::cout << "Random walk for ";
-    for (const auto& [key, value] : optional) {
-        std::cout << key << " = " << value << "; ";
-    }
-    std::cout << std::endl;
-    std::cout << "Start: " << ctime(&start_time);
-
-    return start;
-}
-
-
-void printEnd(const std::chrono::system_clock::time_point start) {
-    const auto end = std::chrono::system_clock::now();
-
-    const std::chrono::duration<float> duration = end - start;
-    const time_t end_time = std::chrono::system_clock::to_time_t(end);
-    std::cout << "End:   " << ctime(&end_time);
-    std::cout << "Duration: " << duration.count() << "s\n" << std::endl;
-}
--- a/src/simulation/sims.h
+++ b/src/simulation/sims.h
@@ -1,46 +0,0 @@
-#ifndef RWSIM_SIMS_H
-#define RWSIM_SIMS_H
-
-#include "../motions/base.h"
-#include "../times/base.h"
-
-#include <unordered_map>
-#include <string>
-#include <chrono>
-
-/**
- * @brief Run simulation for spectra
- *
- * @param parameter Dictionary of parameter for simulation
- * @param optional Dictionary of parameter set via command line
- * @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);
-
-/**
- * @brief Run simulation for stimulated echoes
- *
- * @param parameter Dictionary of parameter for simulation
- * @param optional Dictionary of parameter set via command line
- * @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);
-
-/**
- * @brief Create trajectory of a single walker
- *
- * @param motion Motion model
- * @param dist Distribution of correlation times
- * @param t_max Double that defines maximum time of trajectory
- * @param out_time Vector of waiting times
- * @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);
-
-std::chrono::system_clock::time_point printStart(std::unordered_map<std::string, double> &optional);
-void printEnd(std::chrono::system_clock::time_point start);
-
-#endif //RWSIM_SIMS_H
--- a/src/times/CMakeLists.txt
+++ b/src/times/CMakeLists.txt
@@ -0,0 +1,5 @@
+add_library(RWTime OBJECT
+  base.cpp base.h 
+  delta.cpp delta.h 
+  lognormal.cpp lognormal.h
+)
--- a/src/times/base.cpp
+++ b/src/times/base.cpp
@@ -1,27 +1,31 @@
 #include "base.h"
-#include "delta.h"
-#include "lognormal.h"
+#include "../utils/registry.h"

-#include <stdexcept>
+namespace times {
+BaseDistribution::BaseDistribution(std::string name, const double tau)
+    : m_name(std::move(name)), m_tau(tau), m_tau_jump(tau) {}

-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) {}
+BaseDistribution::BaseDistribution(std::string name)
+    : m_name(std::move(name)) {}

-Distribution::Distribution(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(std::mt19937_64 &rng) const {
+  return std::exponential_distribution(1. / m_tau_jump)(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");
 }

-Distribution* Distribution::createFromInput(const std::string& input, std::mt19937_64& rng) {
-    if (input == "Delta")
-        return new DeltaDistribution(rng);
+std::unordered_map<std::string, double>
+BaseDistribution::getParameters() const {
+  return std::unordered_map<std::string, double>{
+      {"tau", m_tau},
+  };
+}

-    if (input == "LogNormal")
-        return new LogNormalDistribution(rng);
-
-    throw std::invalid_argument("Invalid input " + input);
-}
+std::unique_ptr<BaseDistribution>
+BaseDistribution::createFromInput(const std::string &input) {
+  return Registry<BaseDistribution>::create(input);
+}
+} // namespace times
--- a/src/times/base.h
+++ b/src/times/base.h
@@ -1,33 +1,54 @@
 #ifndef RWSIM_TIMESBASE_H
 #define RWSIM_TIMESBASE_H

+#include <memory>
 #include <random>
 #include <unordered_map>

-class Distribution {
+namespace times {
+
+/// Base class for correlation time distributions.
+///
+/// Models the distribution of correlation times tau for individual walkers.
+/// Each walker draws its own tau from the distribution during initialize(),
+/// then waiting times between jumps are drawn from an exponential distribution
+/// with mean tau_jump (Poisson process).
+///
+/// Delta: all walkers have the same tau (no distribution).
+/// LogNormal: tau is drawn from a log-normal distribution with width sigma.
+class BaseDistribution {
 public:
-    virtual ~Distribution() = default;
+  virtual ~BaseDistribution() = default;

-    Distribution(std::string, double, std::mt19937_64&);
-    explicit Distribution(std::string, std::mt19937_64&);
+  BaseDistribution(std::string, double);
+  explicit BaseDistribution(std::string);

-    [[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>&);
+  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;
+  /// Draw a walker-specific correlation time from the distribution.
+  virtual void initialize(std::mt19937_64 &rng) = 0;

-    static Distribution* createFromInput(const std::string& input, std::mt19937_64& rng);
+  /// Draw a waiting time from the exponential distribution with mean m_tau_jump.
+  [[nodiscard]] double tau_wait(std::mt19937_64 &rng) const;
+
+  [[nodiscard]] virtual std::unique_ptr<BaseDistribution> clone() const = 0;
+  [[nodiscard]] virtual std::string toString() const = 0;
+
+  /// Factory: create a distribution by registered name.
+  static std::unique_ptr<BaseDistribution>
+  createFromInput(const std::string &input);

 protected:
-    std::string m_name{"BaseDistribution"};
-    double m_tau{1.};
-    double m_tau_jump{1.};
-    std::mt19937_64& m_rng;
+  std::string m_name{"BaseDistribution"};
+  double m_tau{1.};      ///< Mean correlation time (center of distribution)
+  double m_tau_jump{1.}; ///< Walker-specific correlation time (drawn in initialize)
 };
+} // namespace times

-#endif //RWSIM_TIMESBASE_H
+#endif // RWSIM_TIMESBASE_H
--- a/src/times/delta.cpp
+++ b/src/times/delta.cpp
@@ -1,11 +1,22 @@
 #include "delta.h"
+#include "../utils/registry.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) {}
+static AutoRegister<times::BaseDistribution, times::DeltaDistribution>
+    reg("Delta");

-void DeltaDistribution::initialize() {
-    m_tau_jump = m_tau;
+namespace times {
+DeltaDistribution::DeltaDistribution(const double tau)
+    : BaseDistribution(std::string("Delta"), tau) {}
+DeltaDistribution::DeltaDistribution()
+    : BaseDistribution(std::string("Delta")) {}
+
+void DeltaDistribution::initialize(std::mt19937_64 &) { m_tau_jump = m_tau; }
+
+std::unique_ptr<BaseDistribution> DeltaDistribution::clone() const {
+  return std::make_unique<DeltaDistribution>(*this);
 }

-void DeltaDistribution::draw_tau() {}
-
+std::string DeltaDistribution::toString() const {
+  return {"Delta/tau=" + std::to_string(m_tau)};
+}
+} // namespace times
--- a/src/times/delta.h
+++ b/src/times/delta.h
@@ -3,13 +3,16 @@

 #include "base.h"

-class DeltaDistribution final : public Distribution {
+namespace times {
+class DeltaDistribution final : public BaseDistribution {
 public:
-    DeltaDistribution(double, std::mt19937_64&);
-    explicit DeltaDistribution(std::mt19937_64 &rng);
+  explicit DeltaDistribution(double);
+  DeltaDistribution();

-    void initialize() override;
-    void draw_tau() override;
+  void initialize(std::mt19937_64 &rng) override;
+  [[nodiscard]] std::unique_ptr<BaseDistribution> clone() const override;
+
+  [[nodiscard]] std::string toString() const override;
 };
-
-#endif //RWSIM_TIMESDELTA_H
+} // namespace times
+#endif // RWSIM_TIMESDELTA_H
--- a/src/times/lognormal.cpp
+++ b/src/times/lognormal.cpp
@@ -1,19 +1,42 @@
 #include "lognormal.h"
+#include "../utils/registry.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) {}
+static AutoRegister<times::BaseDistribution, times::LogNormalDistribution>
+    reg("LogNormal");

-void LogNormalDistribution::setParameters(const std::unordered_map<std::string, double> &parameters) {
-    m_sigma = parameters.at("sigma");
-    Distribution::setParameters(parameters);
+namespace times {
+LogNormalDistribution::LogNormalDistribution(const double tau,
+                                             const double sigma)
+    : BaseDistribution(std::string("LogNormal"), tau), m_sigma(sigma),
+      m_distribution(std::log(tau), sigma) {}
+LogNormalDistribution::LogNormalDistribution()
+    : BaseDistribution(std::string("LogNormal")) {}
+
+void LogNormalDistribution::setParameters(
+    const std::unordered_map<std::string, double> &parameters) {
+  m_sigma = parameters.at("sigma");
+  BaseDistribution::setParameters(parameters);
 }

-void LogNormalDistribution::initialize() {
-    m_distribution = std::lognormal_distribution(std::log(m_tau), m_sigma);
-    m_tau_jump = m_distribution(m_rng);
+std::unordered_map<std::string, double>
+LogNormalDistribution::getParameters() const {
+  auto parameter = BaseDistribution::getParameters();
+  parameter["sigma"] = m_sigma;
+  return parameter;
 }

-void LogNormalDistribution::draw_tau() {
-    m_tau_jump = m_distribution(m_rng);
+void LogNormalDistribution::initialize(std::mt19937_64 &rng) {
+  m_distribution = std::lognormal_distribution(std::log(m_tau), m_sigma);
+  m_tau_jump = m_distribution(rng);
 }
+
+std::unique_ptr<BaseDistribution> LogNormalDistribution::clone() const {
+  return std::make_unique<LogNormalDistribution>(*this);
+}
+
+std::string LogNormalDistribution::toString() const {
+  return {"LogNormal/tau=" + std::to_string(m_tau) +
+          "/sigma=" + std::to_string(m_sigma)};
+}
+} // namespace times
--- a/src/times/lognormal.h
+++ b/src/times/lognormal.h
@@ -3,21 +3,26 @@

 #include "base.h"
 #include <random>
-#include <set>

-class LogNormalDistribution final : public Distribution {
+namespace times {
+class LogNormalDistribution final : public BaseDistribution {
 public:
-    LogNormalDistribution(double, double, std::mt19937_64&);
-    explicit LogNormalDistribution(std::mt19937_64 &rng);
+  LogNormalDistribution(double, double);
+  LogNormalDistribution();

-    void setParameters(const std::unordered_map<std::string, double> &) override;
+  void setParameters(const std::unordered_map<std::string, double> &) override;
+  [[nodiscard]] std::unordered_map<std::string, double>
+  getParameters() const override;

-    void initialize() override;
-    void draw_tau() override;
+  [[nodiscard]] std::string toString() const override;
+  [[nodiscard]] std::unique_ptr<BaseDistribution> clone() const override;
+
+  void initialize(std::mt19937_64 &rng) override;

 private:
-    double m_sigma{1};
-    std::lognormal_distribution<> m_distribution;
+  double m_sigma{1};
+  std::lognormal_distribution<> m_distribution;
 };
+} // namespace times

-#endif //LOGNORMAL_H
+#endif // LOGNORMAL_H
--- a/src/utils/CMakeLists.txt
+++ b/src/utils/CMakeLists.txt
@@ -0,0 +1,9 @@
+add_library(
+        utils STATIC
+        functions.cpp
+        functions.h
+        io.cpp
+        io.h
+        ranges.cpp
+        ranges.h
+)
--- a/src/utils/functions.cpp
+++ b/src/utils/functions.cpp
@@ -1,53 +1,62 @@
-#include <vector>
-#include <array>
+#include <algorithm>
 #include <chrono>
 #include <iostream>
-#include <cmath>
-#include <utility>
+#include <vector>

-
-int nearest_index(const std::vector<double> &x_ref, const double x, int start=0) {
-    while (x > x_ref[start+1]) {
-        start++;
-    }
-    return start;
+/// Find index i such that x_ref[i] <= x < x_ref[i+1].
+/// When start > 0, scans linearly from the hint (efficient for sequential lookups).
+/// When start == 0, uses binary search on large vectors.
+int nearest_index(const std::vector<double> &x_ref, const double x,
+                  int start = 0) {
+  const int last = static_cast<int>(x_ref.size()) - 2;
+  if (start == 0 && last > 32) {
+    // Binary search for cold starts on large trajectories
+    auto it = std::upper_bound(x_ref.begin(), x_ref.end(), x);
+    int idx = static_cast<int>(it - x_ref.begin()) - 1;
+    return std::clamp(idx, 0, last);
+  }
+  while (start < last && x > x_ref[start + 1]) {
+    start++;
+  }
+  return start;
 }

-double lerp(const std::vector<double>& x_ref, const std::vector<double>& y_ref, const double x, const int i) {
-    /*
-     * Linear interpolation between two
-     */
-    const double x_left = x_ref[i];
-    const double y_left = y_ref[i];
-    const double x_right = x_ref[i+1];
-    const double y_right = y_ref[i+1];
+double lerp(const std::vector<double> &x_ref, const std::vector<double> &y_ref,
+            const double x, const int i) {
+  // Linear interpolation between points (x_ref[i], y_ref[i]) and (x_ref[i+1], y_ref[i+1]).
+  const double x_left = x_ref[i];
+  const double y_left = y_ref[i];
+  const double x_right = x_ref[i + 1];
+  const double y_right = y_ref[i + 1];

-    const double dydx = (y_right - y_left) / ( x_right - x_left );
+  const double dydx = (y_right - y_left) / (x_right - x_left);

-    return y_left + dydx * (x - x_left);
+  return y_left + dydx * (x - x_left);
 }

-
 std::chrono::time_point<std::chrono::system_clock> printSteps(
    /*
-     * Prints roughly every 10 seconds how many runs were done and gives a time estimation
+     * Prints roughly every 10 seconds how many runs were done and gives a time
+     * estimation
     */
-    const std::chrono::time_point<std::chrono::system_clock> last_print_out,
+    const std::chrono::time_point<std::chrono::system_clock>
+        last_print_out,
    const std::chrono::time_point<std::chrono::system_clock> start,
-    const int total,
-    const int steps
-    ) {
-    const auto now = std::chrono::high_resolution_clock::now();
+    const int total, const int steps) {
+  const auto now = std::chrono::high_resolution_clock::now();

-    if (const std::chrono::duration<float> duration = now - last_print_out; duration.count() < 10.) {
-        return last_print_out;
-    }
+  if (const std::chrono::duration<float> duration = now - last_print_out;
+      duration.count() < 10.) {
+    return last_print_out;
+  }

-    const std::chrono::duration<float> duration = now - start;
-    const auto passed = duration.count();
+  const std::chrono::duration<float> duration = now - start;
+  const auto passed = duration.count();

+  std::cout << steps << " of " << total << " steps: " << passed << "s passed; ~"
+            << passed * static_cast<float>(total - steps) /
+                   static_cast<float>(steps)
+            << "s remaining\n";

-    std::cout << steps << " of " << total << " steps: " << passed << "s passed; ~" << passed * static_cast<float>(total-steps) / static_cast<float>(steps) << "s remaining\n";
-
-    return now;
-}
+  return now;
+}
--- a/src/utils/functions.h
+++ b/src/utils/functions.h
@@ -1,14 +1,16 @@
 #ifndef RWSIM_FUNCTIONS_H
 #define RWSIM_FUNCTIONS_H

-#include <vector>
 #include <chrono>
+#include <vector>

+int nearest_index(const std::vector<double> &, double, int);

-int nearest_index(const std::vector<double>&, double, int);
+double lerp(const std::vector<double> &, const std::vector<double> &, double,
+            int);

-double lerp(const std::vector<double>&, const std::vector<double>&, double, int);
-
-std::chrono::time_point<std::chrono::system_clock> printSteps(std::chrono::time_point<std::chrono::system_clock>, std::chrono::time_point<std::chrono::system_clock>, int, int);
+std::chrono::time_point<std::chrono::system_clock>
+printSteps(std::chrono::time_point<std::chrono::system_clock>,
+           std::chrono::time_point<std::chrono::system_clock>, int, int);

 #endif
--- a/src/utils/io.cpp
+++ b/src/utils/io.cpp
@@ -1,29 +1,28 @@
 #include "io.h"

-#include <sstream>
-#include <fstream>
-#include <iostream>
 #include <algorithm>
-#include <complex>
-#include <vector>
-#include <iomanip>
-#include <unordered_map>
-#include <map>
-#include <string>
 #include <filesystem>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <map>
+#include <sstream>
+#include <string>
+#include <vector>

+std::pair<std::string, double>
+get_optional_parameter(std::vector<std::string>::const_iterator &it) {
+  std::string stripped_arg;
+  if (it->size() > 2 && it->at(0) == '-' && it->at(1) == '-') {
+    stripped_arg = it->substr(2, it->size());
+  } else if (it->size() > 1 && it->at(0) == '-') {
+    stripped_arg = it->substr(1, it->size());
+  }
+  std::transform(stripped_arg.begin(), stripped_arg.end(), stripped_arg.begin(),
+                 [](unsigned char c) { return std::tolower(c); });
+  const auto stripped_value = std::stod(*++it, nullptr);

-std::pair<std::string, double> get_optional_parameter(std::vector<std::string>::const_iterator &it) {
-    std::string stripped_arg;
-    if (it->size() > 2 && it->at(0) == '-' && it->at(1) == '-') {
-        stripped_arg = it->substr(2, it->size());
-    } else if (it->size() > 1 && it->at(0) == '-') {
-        stripped_arg = it->substr(1, it->size());
-    }
-    std::transform(stripped_arg.begin(), stripped_arg.end(), stripped_arg.begin(), [](unsigned char c) { return std::tolower(c); });
-    const auto stripped_value = std::stod(*(++it), nullptr);
-
-    return std::make_pair(stripped_arg, stripped_value);
+  return std::make_pair(stripped_arg, stripped_value);
 }

 /**
@@ -32,194 +31,194 @@ std::pair<std::string, double> get_optional_parameter(std::vector<std::string>::
 * @param argv List of command-line arguments
 * @return Arguments
 */
-Arguments parse_args(const int argc, char* argv[]) {
-    if (argc < 3) {
-        throw std::runtime_error("Not enough arguments: missing parameter file");
+Arguments parse_args(const int argc, char *argv[]) {
+  if (argc < 3) {
+    throw std::runtime_error("Not enough arguments: missing parameter file");
+  }
+
+  Arguments args;
+
+  // convert to vector to use iterator for loop
+  const std::vector<std::string> input_args(argv + 1, argv + argc);
+
+  // collect positional arguments separately
+  std::vector<std::string> positional;
+
+  for (auto it = input_args.begin(); it != input_args.end(); ++it) {
+    // check for optional parameter
+    if (it->at(0) == '-') {
+
+      // only --spectrum and --ste are parameter that are predefined
+      if (*it == "--spectrum") {
+        args.spectrum = true;
+        continue;
+      }
+
+      if (*it == "--ste") {
+        args.ste = true;
+        continue;
+      }
+
+      // all other arguments are optional parameter
+      auto [option_name, option_value] = get_optional_parameter(it);
+      args.optional[option_name] = option_value;
+      continue;
    }

-    Arguments args;
+    positional.push_back(*it);
+  }

-    // convert to vector to use iterator for loop
-    const std::vector<std::string> input_args(argv + 1, argv + argc);
+  if (positional.size() == 2) {
+    // config_file DynamicsModel
+    args.parameter_file = positional[0];
+    args.dynamics_type = positional[1];
+  } else if (positional.size() == 3) {
+    // config_file MotionModel DistributionModel
+    args.parameter_file = positional[0];
+    args.motion_type = positional[1];
+    args.distribution_type = positional[2];
+  } else {
+    throw std::invalid_argument(
+        "Expected 2 positional args (config dynamics) or "
+        "3 positional args (config motion distribution), got " +
+        std::to_string(positional.size()));
+  }

-    for (auto it = input_args.begin();  it != input_args.end(); ++it) {
-        // check for optional parameter
-        if (it->at(0) == '-') {
-
-            // only --spectrum and --ste are parameter that are predefined
-            if (*it == "--spectrum") {
-                args.spectrum = true;
-                continue;
-            }
-
-            if (*it == "--ste") {
-                args.ste = true;
-                continue;
-            }
-
-            // all other arguments are optional parameter
-            auto [option_name, option_value] = get_optional_parameter(it);
-            args.optional[option_name] = option_value;
-            continue;
-        }
-
-        // Two positional parameters are defined: 1. Location of config file; 2. Name of motion model
-        if (args.parameter_file.empty()) {
-            args.parameter_file = *it;
-            continue;
-        }
-
-        if (args.motion_type.empty()) {
-            args.motion_type = *it;
-            continue;
-        }
-
-        if (args.distribution_type.empty()) {
-            args.distribution_type = *it;
-            continue;
-        }
-
-        throw std::invalid_argument("too many positional arguments");
-    }
-
-    if (args.parameter_file.empty() || args.motion_type.empty() || args.distribution_type.empty()) {
-        throw std::invalid_argument("Missing argument");
-    }
-
-    return args;
+  return args;
 }

+std::unordered_map<std::string, double>
+read_parameter(const std::filesystem::path &infile) {
+  if (!exists(infile)) {
+    std::cerr << "File " << infile << " does not exist" << std::endl;
+    exit(1);
+  }

-std::unordered_map<std::string, double> read_parameter(const std::filesystem::path& infile) {
-    if (!std::filesystem::exists(infile)) {
-        std::cerr << "File " << infile << " does not exist" << std::endl;
-        exit(1);
-    }
+  std::ifstream instream(infile);

-    std::ifstream instream(infile);
+  std::unordered_map<std::string, double> parameter;

-    std::unordered_map<std::string, double> parameter;
+  std::string line;
+  std::string delim = "=";
+  std::string key;
+  std::string value;
+  size_t delim_pos;

-    std::string line;
-    std::string delim = "=";
-    std::string key;
-    std::string value;
-    size_t delim_pos;
+  while (std::getline(instream, line)) {
+    // skip comment lines starting with #, and empty lines
+    if (line[0] == '#' || line.length() == 1)
+      continue;

-    while (std::getline(instream, line)) {
-        // skip comment lines starting with #, and empty lines
-        if (line[0] == '#' || line.length() == 1) continue;
+    // strip spaces from line to have always key=value
+    line.erase(std::remove(line.begin(), line.end(), ' '), line.end());

-        // strip spaces from line to have always key=value
-        line.erase(std::remove(line.begin(), line.end(), ' '), line.end());
+    // split at '=' character and add to map
+    delim_pos = line.find('=');
+    key = line.substr(0, delim_pos);
+    value = line.substr(delim_pos + 1);
+    parameter[key] = std::stod(value);
+  }

-        // split at '=' character and add to map
-        delim_pos = line.find('=');
-        key = line.substr(0, delim_pos);
-        value = line.substr(delim_pos+1);
-        parameter[key] = std::stod(value);
-    }
-
-    return parameter;
+  return parameter;
 }

+std::string make_directory(const std::string &path_name) {
+  if (!std::filesystem::create_directories(path_name)) {
+    std::cout << "Created directory " << path_name << std::endl;
+  }
+
+  return path_name;
+}

 void save_parameter_to_file(
-    const std::string& resulttype,
-    const std::string& motiontype,
-    const std::string& disttype,
-    std::unordered_map<std::string, double>& parameter,
-    std::unordered_map<std::string, double>& optional
-    ) {
+    const std::string &resulttype, const std::string &directory,
+    const std::unordered_map<std::string, double> &parameter,
+    const std::unordered_map<std::string, double> &optional) {
+  std::ostringstream parameter_filename;
+  parameter_filename << directory << "/" << resulttype;

-    std::ostringstream parameter_filename;
-    parameter_filename << resulttype << "_" << motiontype << "_" << disttype;
+  parameter_filename << std::setprecision(6) << std::scientific;
+  for (const auto &[key, value] : optional) {
+    parameter_filename << "_" << key << "=" << value;
+  }
+  parameter_filename << "_parameter.txt";

-    parameter_filename << std::setprecision(6) << std::scientific;
-    for (const auto& [key, value]: optional) {
-        parameter_filename << "_" << key << "=" << value;
-    }
-    parameter_filename << "_parameter.txt";
-
-    {
-        // write data to file, columns are secondary axis (echo delay, evolution times)
-        std::string datafile = parameter_filename.str();
-        std::ofstream filestream(datafile, std::ios::out);
-
-        for (const auto& [key, value]: parameter) {
-            filestream << key << "=" << value << "\n";
-        }
-    }
-}
-
-
-void save_data_to_file(
-    const std::string& resulttype,
-    const std::string& motiontype,
-    const std::string& disttype,
-    const std::vector<double>& x,
-    const std::map<double, std::vector<double>>& y,
-    std::unordered_map<std::string, double>& optional
-    ) {
-    // make file name
-    std::ostringstream datafile_name;
-    datafile_name << resulttype << "_" << motiontype << "_" << disttype;
-    datafile_name << std::setprecision(6) << std::scientific;
-    for (const auto& [key, value]: optional) {
-        datafile_name << "_" << key << "=" << value;
-    }
-    datafile_name << ".dat";
-
-    {
-        // write data to file, columns are secondary axis (echo delay, evolution times)
-        std::string datafile = datafile_name.str();
-        std::ofstream filestream(datafile, std::ios::out);
-
-        // first line are values of secondary axis
-        filestream << "#";
-        for (const auto& [t_echo_j, _] : y) {
-            filestream << t_echo_j << "\t";
-        }
-        filestream << std::endl;
-
-        // write values to file
-        auto size = x.size();
-        for (unsigned int i = 0; i < size; i++) {
-            filestream << x[i];
-            for (const auto& [_, fid_j] : y) {
-                filestream << "\t" << fid_j[i];
-            }
-            filestream << "\n";
-        }
+  {
+    // write data to file, columns are secondary axis (echo delay, evolution
+    // times)
+    std::string datafile = parameter_filename.str();
+    std::ofstream filestream(datafile, std::ios::out);
+
+    for (const auto &[key, value] : parameter) {
+      filestream << key << "=" << value << "\n";
    }
+  }
 }

 void save_data_to_file(
-    const std::string& resulttype,
-    const std::string& motiontype,
-    const std::string& disttype,
-    const std::vector<double>& x,
-    const std::vector<double>& y,
-    std::unordered_map<std::string, double>& optional
-    ) {
-    // make file name
-    std::ostringstream datafile_name;
-    datafile_name << resulttype << "_" << motiontype << "_" << disttype;
-    datafile_name << std::setprecision(6) << std::scientific;
-    for (const auto& [key, value]: optional) {
-        datafile_name << "_" << key << "=" << value;
-    }
-    datafile_name << ".dat";
+    const std::string &resulttype, const std::string &directory,
+    const std::vector<double> &x,
+    const std::map<double, std::vector<double>> &y,
+    const std::unordered_map<std::string, double> &optional) {
+  // make file name
+  std::ostringstream datafile_name;
+  datafile_name << directory << "/" << resulttype;

-    {
-        // write data to file, columns are secondary axis (echo delay, evolution times)
-        std::string datafile = datafile_name.str();
-        std::ofstream filestream(datafile, std::ios::out);
+  datafile_name << std::setprecision(6) << std::scientific;
+  for (const auto &[key, value] : optional) {
+    datafile_name << "_" << key << "=" << value;
+  }
+  datafile_name << ".dat";

-        // write values to file
-        auto size = x.size();
-        for (unsigned int i = 0; i < size; i++) {
-            filestream << x[i] << "\t" << y[i] << "\n";
-        }
+  {
+    // write data to file, columns are secondary axis (echo delay, evolution
+    // times)
+    std::string datafile = datafile_name.str();
+    std::ofstream filestream(datafile, std::ios::out);
+
+    // first line are values of secondary axis
+    filestream << "#";
+    for (const auto &[t_echo_j, _] : y) {
+      filestream << t_echo_j << "\t";
    }
+    filestream << std::endl;
+
+    // write values to file
+    auto size = x.size();
+    for (unsigned int i = 0; i < size; i++) {
+      filestream << x[i];
+      for (const auto &[_, fid_j] : y) {
+        filestream << "\t" << fid_j[i];
+      }
+      filestream << "\n";
+    }
+  }
+}
+
+void save_data_to_file(
+    const std::string &resulttype, const std::string &directory,
+    const std::vector<double> &x, const std::vector<double> &y,
+    const std::unordered_map<std::string, double> &optional) {
+  // make file name
+  std::ostringstream datafile_name;
+  datafile_name << directory << "/" << resulttype;
+
+  datafile_name << std::setprecision(6) << std::scientific;
+  for (const auto &[key, value] : optional) {
+    datafile_name << "_" << key << "=" << value;
+  }
+  datafile_name << ".dat";
+
+  {
+    // write data to file, columns are secondary axis (echo delay, evolution
+    // times)
+    std::string datafile = datafile_name.str();
+    std::ofstream filestream(datafile, std::ios::out);
+
+    // write values to file
+    auto size = x.size();
+    for (unsigned int i = 0; i < size; i++) {
+      filestream << x[i] << "\t" << y[i] << "\n";
+    }
+  }
 }
--- a/src/utils/io.h
+++ b/src/utils/io.h
@@ -1,29 +1,39 @@
 #ifndef RWSIM_IO_H
 #define RWSIM_IO_H

-#include <unordered_map>
+#include <filesystem>
 #include <map>
 #include <string>
-#include <filesystem>
+#include <unordered_map>
 #include <vector>

 struct Arguments {
-    std::string parameter_file{};
-    bool ste = false;
-    bool spectrum = false;
-    std::string motion_type{};
-    std::string distribution_type{};
-    std::unordered_map<std::string, double> optional;
+  std::string parameter_file{};
+  bool ste = false;
+  bool spectrum = false;
+  std::string motion_type{};
+  std::string distribution_type{};
+  std::string dynamics_type{};
+  std::unordered_map<std::string, double> optional;
 };

-Arguments parse_args(int argc, char* argv[]);
-std::pair<std::string, double> get_optional_parameter(std::vector<std::string>::const_iterator &it);
+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::unordered_map<std::string, double> read_parameter(const std::filesystem::path&);
+std::string make_directory(const std::string &path_name);

-
-void save_parameter_to_file(const std::string&, const std::string&, const std::string&, std::unordered_map<std::string, double>&, std::unordered_map<std::string, double>&);
-void save_data_to_file(const std::string&, const std::string&, const std::string&, const std::vector<double>&, const std::map<double, std::vector<double>>&, std::unordered_map<std::string, double>&);
-void save_data_to_file(const std::string&, const std::string&, const std::string&, const std::vector<double>&, const std::vector<double>&, std::unordered_map<std::string, double>&);
+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> &);
+void save_data_to_file(const std::string &, const std::string &,
+                       const std::vector<double> &, const std::vector<double> &,
+                       const std::unordered_map<std::string, double> &);

 #endif
--- a/src/utils/ranges.cpp
+++ b/src/utils/ranges.cpp
@@ -1,55 +1,58 @@
-#include <vector>
 #include <algorithm>
 #include <cmath>
+#include <vector>

 #include "ranges.h"

+std::vector<double> arange(const int size, const double spacing = 1.) {
+  std::vector<double> out(size);
+  std::generate(out.begin(), out.end(),
+                [n = 0, spacing]() mutable { return n++ * spacing; });

-std::vector<double> arange(const int size, const double spacing=1.) {
-    std::vector<double> out(size);
-    std::generate(out.begin(), out.end(), [n = 0, spacing]() mutable { return n++ * spacing; });
-
-    return out;
+  return out;
 }

-std::vector<double> linspace(const double start, const double stop, const int steps) {
-    std::vector<double> range;
+std::vector<double> linspace(const double start, const double stop,
+                             const int steps) {
+  std::vector<double> range;
+  range.reserve(steps);

-    if (steps == 0) {
-        return range;
-    }
-
-    if (steps == 1) {
-        range.push_back(start);
-        return range;
-    }
-
-    const double stepsize = (stop-start) / (steps-1);
-    for (int i=0; i<steps; i++) {
-        range.push_back(start + stepsize * i);
-    }
+  if (steps == 0) {
    return range;
+  }
+
+  if (steps == 1) {
+    range.push_back(start);
+    return range;
+  }
+
+  const double stepsize = (stop - start) / (steps - 1);
+  for (int i = 0; i < steps; i++) {
+    range.push_back(start + stepsize * i);
+  }
+  return range;
 }

+std::vector<double> logspace(const double start, const double stop,
+                             const int steps) {
+  std::vector<double> range;
+  range.reserve(steps);

-std::vector<double> logspace(const double start, const double stop, const int steps) {
-    std::vector<double> range;
-
-    if (steps == 0) {
-        return range;
-    }
-
-    if (steps == 1) {
-        range.push_back(start);
-        return range;
-    }
-
-    const double logstart = std::log10(start);
-    const double logstop = std::log10(stop);
-
-    const double stepsize = (logstop-logstart) / (steps-1);
-    for (int i=0; i<steps; i++) {
-        range.push_back(pow(10, logstart + stepsize * i));
-    }
+  if (steps == 0) {
    return range;
-}
+  }
+
+  if (steps == 1) {
+    range.push_back(start);
+    return range;
+  }
+
+  const double logstart = std::log10(start);
+  const double logstop = std::log10(stop);
+
+  const double stepsize = (logstop - logstart) / (steps - 1);
+  for (int i = 0; i < steps; i++) {
+    range.push_back(std::pow(10, logstart + stepsize * i));
+  }
+  return range;
+}
--- a/src/utils/registry.h
+++ b/src/utils/registry.h
@@ -0,0 +1,59 @@
+#ifndef RWSIM_REGISTRY_H
+#define RWSIM_REGISTRY_H
+
+#include <functional>
+#include <memory>
+#include <stdexcept>
+#include <string>
+#include <unordered_map>
+#include <vector>
+
+/// Self-registering factory for polymorphic types.
+///
+/// Models register themselves via static AutoRegister objects in their .cpp files:
+///   static AutoRegister<BaseMotion, RandomJump> reg("RandomJump");
+///
+/// Important: libraries containing AutoRegister statics must be OBJECT libraries
+/// (not STATIC) in CMake, otherwise the linker strips unreferenced translation
+/// units and the registrations never execute.
+template <typename Base> class Registry {
+public:
+  using Creator = std::function<std::unique_ptr<Base>()>;
+
+  static std::unordered_map<std::string, Creator> &entries() {
+    static std::unordered_map<std::string, Creator> map;
+    return map;
+  }
+
+  static void add(const std::string &name, Creator creator) {
+    entries()[name] = std::move(creator);
+  }
+
+  /// Create an instance by name. Throws invalid_argument with available names
+  /// if the name is not found.
+  static std::unique_ptr<Base> create(const std::string &name) {
+    auto &map = entries();
+    auto it = map.find(name);
+    if (it == map.end()) {
+      std::string msg = "Unknown model '" + name + "'. Available: ";
+      for (const auto &[key, _] : map) {
+        msg += key + ", ";
+      }
+      if (!map.empty()) {
+        msg.pop_back();
+        msg.pop_back();
+      }
+      throw std::invalid_argument(msg);
+    }
+    return it->second();
+  }
+};
+
+/// Place a static instance in each model's .cpp to register it with the factory.
+template <typename Base, typename Derived> struct AutoRegister {
+  explicit AutoRegister(const std::string &name) {
+    Registry<Base>::add(name, []() { return std::make_unique<Derived>(); });
+  }
+};
+
+#endif // RWSIM_REGISTRY_H
--- a/start_sims.py
+++ b/start_sims.py
@@ -0,0 +1,16 @@
+import numpy as np
+
+from python.helpers import *
+
+#Simulation parameter
+motion = 'IsotropicAngle'
+distribution = 'Delta'
+parameter = {
+    "angle": [10, 109.47],
+}
+
+parameter = prepare_rw_parameter(parameter)
+
+for variation in parameter:
+    print(f"\nRun RW for {motion}/{distribution} with arguments {variation}\n")
+    run_sims(motion, distribution, ste=True, spectrum=False, **variation)
Author	SHA1	Message	Date
Dominik Demuth	5acbaaa5f8	Comments	2026-03-08 15:33:56 +01:00
Dominik Demuth	072e211e90	Dynamic scheduling of threads	2026-03-08 15:02:47 +01:00
Dominik Demuth	e178e9dd21	Reuse trajectories for multiple experiments	2026-03-08 14:28:33 +01:00
Dominik Demuth	285c78bed5	Formatting	2026-03-08 14:01:37 +01:00
Dominik Demuth	6579bb028b	Smaller fixes	2026-03-08 13:43:51 +01:00
Dominik Demuth	86ec6b220a	Registry to make easier additions of new models	2026-03-08 13:26:05 +01:00
Dominik Demuth	c4485aac6f	Create dynamics interface to combine motion and dist	2026-03-08 13:04:40 +01:00
Dominik Demuth	6b42051e45	Parallelization with OpenMP, 2-3x faster than before	2026-03-08 12:32:30 +01:00
Dominik Demuth	e45ef8162d	Create Experiment class	2026-03-08 11:43:47 +01:00
Dominik Demuth	ef3cd31be6	Use unique_ptr instead of raw pointer	2026-03-08 11:28:02 +01:00
Dominik Demuth	d44e762bf7	add seed option	2026-03-08 10:47:39 +01:00
Dominik Demuth	e242b6b11f	Add missing CMakeLists.txt	2026-03-08 10:28:13 +01:00
Dominik Demuth	a488c47713	save spectrum in post_process_spectrum.py	2025-02-10 18:23:15 +01:00
Dominik Demuth	bfb1cb314c	add more motions	2024-12-15 18:34:05 +01:00
Dominik Demuth	adbf2af72b	fix wrong header path	2024-12-13 14:05:58 +01:00
Dominik Demuth	ee195c313c	add wobbling in a cone as model	2024-12-13 14:04:30 +01:00
Dominik Demuth	e6331749b2	add namespaces and cmakelists for sub-dirs	2024-12-13 14:03:33 +01:00
Dominik Demuth	e90c4c9543	split simulating and evaluating and stuff	2024-12-01 19:24:19 +01:00
Dominik Demuth	5b5aacff0b	change saving	2024-12-01 19:23:58 +01:00
Dominik Demuth	60ef1b0bcf	add 4 pulse ste	2024-11-30 16:15:38 +01:00