From 15917a0d75e556f9a49d4aa852f618e5a949bb67 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Fri, 23 Aug 2024 18:33:38 +0200 Subject: [PATCH] use struct for coordinates --- CMakeLists.txt | 4 ++++ functions.cpp | 14 ------------- functions.h | 6 ------ main.cpp | 3 +++ motions/base.cpp | 18 +++++++++------- motions/base.h | 8 +++++--- motions/coordinates.cpp | 43 +++++++++++++++++++++++++++++++++++++++ motions/coordinates.h | 23 +++++++++++++++++++++ motions/isosmallangle.cpp | 22 ++++++++++++++++++++ motions/isosmallangle.h | 24 ++++++++++++++++++++++ motions/random.cpp | 3 +-- motions/random.h | 1 + motions/tetrahedral.cpp | 37 +++++++-------------------------- motions/tetrahedral.h | 2 -- 14 files changed, 144 insertions(+), 64 deletions(-) create mode 100644 motions/coordinates.cpp create mode 100644 motions/coordinates.h create mode 100644 motions/isosmallangle.cpp create mode 100644 motions/isosmallangle.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 7949d5f..1ea903a 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/functions.cpp b/functions.cpp index 5e82827..95b4689 100644 --- a/functions.cpp +++ b/functions.cpp @@ -31,20 +31,6 @@ double lerp(const std::vector& x_ref, const std::vector& y_ref, return y_left + dydx * (x - x_left); } -std::array spherical_to_xyz(const double cos_theta, const double phi) { - const double sin_theta = std::sin(std::acos(cos_theta)); - const std::array xyz{sin_theta * std::cos(phi), sin_theta * std::sin(phi), cos_theta}; - - return xyz; -} - -std::pair xyz_to_spherical(const std::array& xyz) { - // std::cout << "length " < printSteps( /* diff --git a/functions.h b/functions.h index 8be3e94..8bd3e8d 100644 --- a/functions.h +++ b/functions.h @@ -2,8 +2,6 @@ #define RWSIM_FUNCTIONS_H #include -#include -#include #include @@ -11,10 +9,6 @@ int nearest_index(const std::vector&, double, int); double lerp(const std::vector&, const std::vector&, double, int); -std::array spherical_to_xyz(double, double); - -std::pair xyz_to_spherical(const std::array& xyz); - std::chrono::time_point printSteps(std::chrono::time_point, std::chrono::time_point, int, int); #endif diff --git a/main.cpp b/main.cpp index 33bc3e3..60152e6 100644 --- a/main.cpp +++ b/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) { diff --git a/motions/base.cpp b/motions/base.cpp index 5a86d61..d926deb 100644 --- a/motions/base.cpp +++ b/motions/base.cpp @@ -5,13 +5,11 @@ // #include #include #include -#include #include "base.h" +#include "coordinates.h" #include -#include - 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 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}; +} \ No newline at end of file diff --git a/motions/base.h b/motions/base.h index 6cceb51..54f106b 100644 --- a/motions/base.h +++ b/motions/base.h @@ -5,9 +5,10 @@ #ifndef RWSIM_MOTIONBASE_H #define RWSIM_MOTIONBASE_H - +#include "coordinates.h" #include -#include + + class Motion { public: @@ -16,8 +17,9 @@ public: Motion(double, double, std::mt19937_64&); explicit Motion(std::mt19937_64&); - std::pair 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; diff --git a/motions/coordinates.cpp b/motions/coordinates.cpp new file mode 100644 index 0000000..db6629d --- /dev/null +++ b/motions/coordinates.cpp @@ -0,0 +1,43 @@ +// +// Created by dominik on 8/22/24. +// + +#include "coordinates.h" +#include +#include +#include + +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)}; +} diff --git a/motions/coordinates.h b/motions/coordinates.h new file mode 100644 index 0000000..18d9649 --- /dev/null +++ b/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 diff --git a/motions/isosmallangle.cpp b/motions/isosmallangle.cpp new file mode 100644 index 0000000..c1ac818 --- /dev/null +++ b/motions/isosmallangle.cpp @@ -0,0 +1,22 @@ +// +// Created by dominik on 8/21/24. +// + +#include "isosmallangle.h" + +#include + +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); +} diff --git a/motions/isosmallangle.h b/motions/isosmallangle.h new file mode 100644 index 0000000..8973669 --- /dev/null +++ b/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 diff --git a/motions/random.cpp b/motions/random.cpp index fe68628..fd319bd 100644 --- a/motions/random.cpp +++ b/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()); } diff --git a/motions/random.h b/motions/random.h index ed4f9e2..17e003e 100644 --- a/motions/random.h +++ b/motions/random.h @@ -8,6 +8,7 @@ #include "base.h" #include + class RandomJump final : public Motion { public: RandomJump(double, double, std::mt19937_64&); diff --git a/motions/tetrahedral.cpp b/motions/tetrahedral.cpp index a23ea72..027bc07 100644 --- a/motions/tetrahedral.cpp +++ b/motions/tetrahedral.cpp @@ -4,44 +4,21 @@ #include #include "tetrahedral.h" -#include -#include - -#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 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); } } diff --git a/motions/tetrahedral.h b/motions/tetrahedral.h index cae0578..ec564cb 100644 --- a/motions/tetrahedral.h +++ b/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 m_corners{}; int m_corner_idx{0};