use struct for coordinates

CMakeLists.txt

@@ -22,6 +22,10 @@ add_executable(rwsim main.cpp
         ranges.h
         motions/tetrahedral.cpp
         motions/tetrahedral.h
+        motions/isosmallangle.cpp
+        motions/isosmallangle.h
+        motions/coordinates.cpp
+        motions/coordinates.h
 )
 
 target_compile_options(rwsim PUBLIC -Werror -Wall -Wextra -Wconversion -O2)

functions.cpp

@@ -31,20 +31,6 @@ double lerp(const std::vector<double>& x_ref, const std::vector<double>& y_ref,
     return y_left + dydx * (x - x_left);
 }
 
-std::array<double, 3> spherical_to_xyz(const double cos_theta, const double phi) {
-    const double sin_theta = std::sin(std::acos(cos_theta));
-    const std::array<double, 3> xyz{sin_theta * std::cos(phi), sin_theta * std::sin(phi), cos_theta};
-
-    return xyz;
-}
-
-std::pair<double, double> xyz_to_spherical(const std::array<double, 3>& xyz) {
-    // std::cout << "length " <<xyz[0]*xyz[0] + xyz[1]*xyz[1] + xyz[2]*xyz[2] << std::endl;
-    double cos_theta = xyz[2];
-    double phi = std::atan2(xyz[1], xyz[0]);
-
-    return std::make_pair(cos_theta, phi);
-}
 
 std::chrono::time_point<std::chrono::system_clock> printSteps(
     /*

functions.h

@@ -2,8 +2,6 @@
 #define RWSIM_FUNCTIONS_H
 
 #include <vector>
-#include <array>
-#include <utility>
 #include <chrono>
 
 
@@ -11,10 +9,6 @@ int nearest_index(const std::vector<double>&, double, int);
 
 double lerp(const std::vector<double>&, const std::vector<double>&, double, int);
 
-std::array<double, 3> spherical_to_xyz(double, double);
-
-std::pair<double, double> xyz_to_spherical(const std::array<double, 3>& xyz);
-
 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

main.cpp

@@ -4,6 +4,7 @@
 
 #include "io.h"
 #include "sims.h"
+#include "motions/isosmallangle.h"
 #include "motions/random.h"
 #include "motions/tetrahedral.h"
 #include "times/delta.h"
@@ -24,6 +25,8 @@ int main (const int argc, char *argv[]) {
     std::mt19937_64 rng(rd());
 
     auto motion = TetrahedralJump(rng);
+    // auto motion = RandomJump(rng);
+    // auto motion = SmallAngle(1, 1, 123, rng);
     auto dist = DeltaDistribution(rng);
 
     if (args.spectrum) {

motions/base.cpp

@@ -5,13 +5,11 @@
 // #include <cmath>
 #include <random>
 #include <cmath>
-#include <utility>
 
 #include "base.h"
+#include "coordinates.h"
 
 #include <iostream>
-#include <ostream>
-
 
 Motion::Motion(const double delta, const double eta, std::mt19937_64& rng) : m_delta(delta), m_eta(eta), m_rng(rng) {
     m_uni_dist = std::uniform_real_distribution(0., 1.);
@@ -28,10 +26,16 @@ double Motion::omega_q(const double cos_theta, const double phi) const {
     return M_PI * m_delta * (3. * cos_theta_square - 1. - m_eta * sin_theta_square * std::cos(2.*phi));
 }
 
-std::pair<double, double> Motion::draw_position() {
+double Motion::omega_q(const SphericalPos& pos) const {
+    // std::cout << pos.cos_theta << " " << pos.phi << std::endl;
+    auto [cos_theta, phi] = pos;
+
+    return omega_q(cos_theta, phi);
+}
+
+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 std::make_pair(cos_theta, phi);
+    return {cos_theta, phi};
 }
-

motions/base.h

@@ -5,9 +5,10 @@
 #ifndef RWSIM_MOTIONBASE_H
 #define RWSIM_MOTIONBASE_H
 
-
+#include "coordinates.h"
 #include <random>
-#include <utility>
+
+
 
 class Motion {
 public:
@@ -16,8 +17,9 @@ public:
     Motion(double, double, std::mt19937_64&);
     explicit Motion(std::mt19937_64&);
 
-    std::pair<double, double> draw_position();
+    SphericalPos draw_position();
     [[nodiscard]] double omega_q(double, double) const;
+    [[nodiscard]] double omega_q(const SphericalPos&) const;
 
     virtual void initialize() = 0;
     virtual double jump() = 0;

motions/coordinates.cpp

@@ -0,0 +1,43 @@
+//
+// Created by dominik on 8/22/24.
+//
+
+#include "coordinates.h"
+#include <cmath>
+#include <iostream>
+#include <ostream>
+
+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 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});
+    }
+
+    const double norm{std::sqrt(1 - cos_theta * cos_theta) + 1e-15};
+
+    auto [x, y , z] = spherical_to_xyz(old_pos);
+
+    const auto new_pos = CartesianPos{
+        cos_alpha * x + sin_alpha * (-x * z * sin_beta - y * cos_beta) / norm,
+        cos_alpha * y + sin_alpha * (-y * z * sin_beta + x * cos_beta) / norm,
+        cos_alpha * z + sin_alpha * norm * sin_beta
+    };
+
+    return xyz_to_spherical(new_pos);
+}
+
+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)};
+}

motions/coordinates.h

@@ -0,0 +1,23 @@
+//
+// Created by dominik on 8/22/24.
+//
+
+#ifndef COORDINATES_H
+#define COORDINATES_H
+
+struct SphericalPos {
+    double cos_theta;
+    double phi;
+};
+
+struct CartesianPos {
+    double x;
+    double y;
+    double z;
+};
+
+SphericalPos rotate(const SphericalPos&, double, double);
+CartesianPos spherical_to_xyz(const SphericalPos&);
+SphericalPos xyz_to_spherical(const CartesianPos&);
+
+#endif //COORDINATES_H

motions/isosmallangle.cpp

@@ -0,0 +1,22 @@
+//
+// Created by dominik on 8/21/24.
+//
+
+#include "isosmallangle.h"
+
+#include <algorithm>
+
+SmallAngle::SmallAngle(const double delta, const double eta, const double chi, std::mt19937_64 &rng) : Motion(delta, eta, rng), m_chi(chi) {};
+SmallAngle::SmallAngle(std::mt19937_64 &rng) : Motion(rng) {
+}
+
+void SmallAngle::initialize() {
+    m_prev_pos = draw_position();
+};
+
+double SmallAngle::jump() {
+    const double gamma{2 * M_PI * m_uni_dist(m_rng)};
+    m_prev_pos = rotate(m_prev_pos, gamma, m_chi);
+
+    return omega_q(m_prev_pos);
+}

motions/isosmallangle.h

@@ -0,0 +1,24 @@
+//
+// Created by dominik on 8/21/24.
+//
+
+#ifndef RWSIM_MOTIONISOSMALLANGLE_H
+#define RWSIM_MOTIONISOSMALLANGLE_H
+
+#include "base.h"
+#include "coordinates.h"
+
+class SmallAngle final : public Motion {
+public:
+    SmallAngle(double, double, double, std::mt19937_64& );
+    explicit SmallAngle(std::mt19937_64&);
+
+    void initialize() override;
+    double jump() override;
+
+private:
+    double m_chi{0};
+    SphericalPos m_prev_pos{0., 0.};
+};
+
+#endif //RWSIM_MOTIONISOSMALLANGLE_H

motions/random.cpp

@@ -12,6 +12,5 @@ RandomJump::RandomJump(std::mt19937_64 &rng) : Motion(rng) {}
 void RandomJump::initialize() {}
 
 double RandomJump::jump() {
-    const auto [cos_theta, phi] = draw_position();
-    return omega_q(cos_theta, phi);
+    return omega_q(draw_position());
 }

motions/random.h

@@ -8,6 +8,7 @@
 
 #include "base.h"
 #include <random>
+
 class RandomJump final : public Motion {
 public:
     RandomJump(double, double, std::mt19937_64&);

motions/tetrahedral.cpp

@@ -4,44 +4,21 @@
 #include <random>
 #include "tetrahedral.h"
 
-#include <iostream>
-#include <ostream>
-
-#include "../functions.h"
-
-
 TetrahedralJump::TetrahedralJump(const double delta, const double eta, std::mt19937_64& rng) : Motion(delta, eta, rng) {}
 
 TetrahedralJump::TetrahedralJump(std::mt19937_64& rng) : Motion(rng) {}
 
 void TetrahedralJump::initialize() {
-    const auto [cos_theta, phi] = draw_position();
-    m_corners[0] = omega_q(cos_theta, phi);
-    const double alpha = 2. * M_PI * m_uni_dist(m_rng);
+    // const auto pos = SphericalPos{0., 0};
+    const auto pos = draw_position();
+    m_corners[0] = omega_q(pos);
 
-    const auto vec = spherical_to_xyz(cos_theta, phi);
-    const double norm = std::sqrt(1 - cos_theta * cos_theta) + 1e-15;
+    const double alpha = 2. * M_PI * m_uni_dist(m_rng);
+    // const double alpha = 0.;
 
     for (int i = 1; i<4; i++) {
-        const double cos_alpha = std::cos(alpha + (i-1) * 2*M_PI / 3.);
-        const double sin_alpha = std::sin(alpha + (i-1) * 2*M_PI / 3.);
-        std::array<double, 3> rotated_position{};
-        if (cos_theta != 1 && cos_theta != -1) {
-            // std::cout << cos_theta << std::endl;
-            rotated_position = {
-                m_cos_beta * vec[0] + m_sin_beta * (-vec[0] * vec[2] * sin_alpha - vec[1] * cos_alpha) / norm,
-                m_cos_beta * vec[1] + m_sin_beta * (-vec[1] * vec[2] * sin_alpha + vec[0] * cos_alpha) / norm,
-                m_cos_beta * vec[2] + m_sin_beta * norm * sin_alpha
-            };
-        } else {
-            rotated_position = {
-                m_sin_beta * cos_alpha,
-                m_sin_beta * sin_alpha,
-                m_cos_beta * cos_theta
-            };
-        }
-        auto [new_cos_theta, new_phi] = xyz_to_spherical(rotated_position);
-        m_corners[i] = omega_q(new_cos_theta, new_phi);
+        auto corner_pos = rotate(pos, m_beta, alpha + (i-1) * 2*M_PI/3.);
+        m_corners[i] = omega_q(corner_pos);
     }
 }
 

motions/tetrahedral.h

@@ -20,8 +20,6 @@ public:
 
 private:
     const double m_beta{std::acos(-1/3.)};
-    const double m_cos_beta{-1./3.};
-    const double m_sin_beta{ 2. * std::sqrt(2.)/3. };
 
     std::array<double, 4> m_corners{};
     int m_corner_idx{0};