Comments

src/dynamics/base.h

@@ -6,23 +6,43 @@
 #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);
 };
 

src/dynamics/decoupled.h

@@ -7,6 +7,10 @@
 
 #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,

src/experiments/base.h

@@ -6,24 +6,41 @@
 #include <unordered_map>
 #include <vector>
 
+/// Time series of a single walker: NMR frequency, accumulated phase, and time.
 struct Trajectory {
-  std::vector<double> time;
-  std::vector<double> phase;
-  std::vector<double> omega;
+  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;
 };
 

src/experiments/spectrum.h

@@ -7,6 +7,14 @@
 #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,
@@ -19,10 +27,10 @@ public:
   void merge(const Experiment &other) override;
 
 private:
-  int m_num_acq{0};
-  std::vector<double> m_t_fid;
-  std::vector<double> m_echo_times;
-  std::map<double, std::vector<double>> m_fid_dict;
+  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};

src/experiments/ste.h

@@ -7,6 +7,16 @@
 #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,
@@ -20,12 +30,12 @@ public:
 
 private:
   int m_num_mix_times{0};
-  double m_tpulse4{0};
-  std::vector<double> m_evolution_times;
-  std::vector<double> m_mixing_times;
-  std::map<double, std::vector<double>> m_cc_dict;
-  std::map<double, std::vector<double>> m_ss_dict;
-  std::vector<double> m_f2;
+  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};

src/motions/base.cpp

@@ -10,6 +10,8 @@ BaseMotion::BaseMotion(std::string name, const double delta, const double 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;

src/motions/base.h

@@ -8,6 +8,17 @@
 #include <unordered_map>
 
 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 ~BaseMotion() = default;
@@ -15,12 +26,20 @@ public:
   BaseMotion(std::string, double, double);
   explicit BaseMotion(std::string);
 
+  /// Draw a uniform random orientation on the unit sphere.
   coordinates::SphericalPos draw_position(std::mt19937_64 &rng);
+
+  /// Compute quadrupolar frequency from orientation angles.
   [[nodiscard]] double omega_q(double, double) const;
   [[nodiscard]] double omega_q(const coordinates::SphericalPos &) const;
 
+  /// Set up initial state for a new walker (called once per walker).
   virtual void initialize(std::mt19937_64 &rng) = 0;
+
+  /// Perform one jump, returning the new omega_q.
   virtual double jump(std::mt19937_64 &rng) = 0;
+
+  /// 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> &);
@@ -35,14 +54,15 @@ public:
 
   [[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.};
+  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.};
+  double m_initial_omega{0.}; ///< Frequency at initial position
 };
 
 std::ostream &operator<<(std::ostream &os, const BaseMotion &m);

src/motions/conemotion.h

@@ -5,10 +5,15 @@
 #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;
@@ -16,8 +21,8 @@ public:
   getParameters() const override;
 
 protected:
-  double m_angle{0};
-  coordinates::SphericalPos m_axis{1, 0};
+  double m_angle{0};                       ///< Half-opening angle of the cone (radians)
+  coordinates::SphericalPos m_axis{1, 0};  ///< Symmetry axis orientation
 };
 } // namespace motions
 

src/motions/coordinates.h

@@ -2,7 +2,10 @@
 #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;
@@ -14,7 +17,10 @@ struct CartesianPos {
   double z;
 };
 
-SphericalPos rotate(const SphericalPos &, double, double);
+/// 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);
+
 CartesianPos spherical_to_xyz(const SphericalPos &);
 SphericalPos xyz_to_spherical(const CartesianPos &);
 } // namespace coordinates

src/motions/diffusivemotion.h

@@ -5,14 +5,19 @@
 #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.};
+  coordinates::SphericalPos m_prev_pos{0., 0.}; ///< Current orientation
 };
 } // namespace motions
 

src/sims.h

@@ -10,11 +10,18 @@
 #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);
 

src/times/base.h

@@ -6,6 +6,16 @@
 #include <unordered_map>
 
 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 ~BaseDistribution() = default;
@@ -21,19 +31,23 @@ public:
   [[nodiscard]] virtual std::unordered_map<std::string, double>
   getParameters() const;
 
+  /// Draw a walker-specific correlation time from the distribution.
   virtual void initialize(std::mt19937_64 &rng) = 0;
-  [[nodiscard]] double tau_wait(std::mt19937_64 &rng) const;
-  [[nodiscard]] virtual std::unique_ptr<BaseDistribution> clone() const = 0;
 
+  /// 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.};
+  double m_tau{1.};      ///< Mean correlation time (center of distribution)
+  double m_tau_jump{1.}; ///< Walker-specific correlation time (drawn in initialize)
 };
 } // namespace times
 

src/utils/functions.cpp

@@ -3,6 +3,9 @@
 #include <iostream>
 #include <vector>
 
+/// 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;
@@ -20,9 +23,7 @@ int nearest_index(const std::vector<double> &x_ref, const double x,
 
 double lerp(const std::vector<double> &x_ref, const std::vector<double> &y_ref,
             const double x, const int i) {
-  /*
-   * Linear interpolation between two
-   */
+  // 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];

src/utils/registry.h

@@ -8,6 +8,14 @@
 #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>()>;
@@ -21,6 +29,8 @@ public:
     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);
@@ -39,6 +49,7 @@ public:
   }
 };
 
+/// 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>(); });