diff --git a/CMakeLists.txt b/CMakeLists.txt index 1ea903a..b556b3f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -26,6 +26,10 @@ add_executable(rwsim main.cpp motions/isosmallangle.h motions/coordinates.cpp motions/coordinates.h + motions/bimodalangle.cpp + motions/bimodalangle.h + times/lognormal.cpp + times/lognormal.h ) target_compile_options(rwsim PUBLIC -Werror -Wall -Wextra -Wconversion -O2) diff --git a/config.txt b/config.txt index 5c60c0d..5136da4 100644 --- a/config.txt +++ b/config.txt @@ -12,9 +12,26 @@ techo_start=0e-6 techo_stop=40e-6 techo_steps=5 # STE part -tevo_start=0e-6 +tevo_start=1e-6 tevo_stop=60e-6 tevo_steps=121 -tmix_start=1e-6 -tmix_stop=1e0 -tmix_steps=51 +tmix_start=1e-5 +tmix_stop=1e1 +tmix_steps=61 +tau=0.001 +tau=0.001 +tau=0.001 +tau=0.001 +tau=0.001 +tau=0.001 +tau=0.01 +tau=0.0031622776601683794 +tau=0.001 +tau=0.00031622776601683794 +tau=0.0001 +tau=3.1622776601683795e-05 +tau=1e-05 +tau=3.162277660168379e-06 +tau=1e-06 +tau=3.162277660168379e-07 +tau=1e-07 diff --git a/main.cpp b/main.cpp index 60152e6..20385b0 100644 --- a/main.cpp +++ b/main.cpp @@ -4,10 +4,12 @@ #include "io.h" #include "sims.h" +#include "motions/bimodalangle.h" #include "motions/isosmallangle.h" #include "motions/random.h" #include "motions/tetrahedral.h" #include "times/delta.h" +#include "times/lognormal.h" int main (const int argc, char *argv[]) { @@ -19,15 +21,22 @@ int main (const int argc, char *argv[]) { return 1; } + std::unordered_map parameter { read_parameter(args.parameter_file) }; + + std::random_device rd; std::mt19937_64 rng(rd()); - auto motion = TetrahedralJump(rng); + // auto motion = BimodalAngle(1, 1, 2, 30, 0.98, rng); + // auto motion = TetrahedralJump(rng); // auto motion = RandomJump(rng); - // auto motion = SmallAngle(1, 1, 123, rng); - auto dist = DeltaDistribution(rng); + auto motion = SmallAngle(1, 1, 1, rng); + + // auto dist = DeltaDistribution(rng); + auto dist = LogNormalDistribution(1, 2, rng); + if (args.spectrum) { run_spectrum(parameter, motion, dist); diff --git a/motions/base.cpp b/motions/base.cpp index d926deb..9b63f49 100644 --- a/motions/base.cpp +++ b/motions/base.cpp @@ -9,7 +9,7 @@ #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.); @@ -38,4 +38,4 @@ SphericalPos Motion::draw_position() { const double phi = 2.0 * M_PI * m_uni_dist(m_rng); return {cos_theta, phi}; -} \ No newline at end of file +} diff --git a/motions/base.h b/motions/base.h index 54f106b..a2aa782 100644 --- a/motions/base.h +++ b/motions/base.h @@ -7,7 +7,7 @@ #include "coordinates.h" #include - +#include class Motion { @@ -36,4 +36,5 @@ protected: std::uniform_real_distribution<> m_uni_dist; }; + #endif //RWSIM_MOTIONBASE_H diff --git a/motions/bimodalangle.cpp b/motions/bimodalangle.cpp new file mode 100644 index 0000000..f2db206 --- /dev/null +++ b/motions/bimodalangle.cpp @@ -0,0 +1,25 @@ +// +// Created by dominik on 8/23/24. +// + +#include "bimodalangle.h" +#include "base.h" + +BimodalAngle::BimodalAngle(const double delta, const double eta, const double angle1, const double angle2, const double prob, std::mt19937_64 &rng) : + Motion(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(rng) {} + +void BimodalAngle::initialize() { + m_prev_pos = draw_position(); +}; + +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); + + return omega_q(m_prev_pos); +} diff --git a/motions/bimodalangle.h b/motions/bimodalangle.h new file mode 100644 index 0000000..09211f7 --- /dev/null +++ b/motions/bimodalangle.h @@ -0,0 +1,25 @@ +// +// Created by dominik on 8/23/24. +// + +#ifndef BIMODALANGLE_H +#define BIMODALANGLE_H + +#include "base.h" + +class BimodalAngle : public Motion { +public: + BimodalAngle(double, double, double, double, double, std::mt19937_64& ); + explicit BimodalAngle(std::mt19937_64&); + + void initialize() override; + double jump() override; + +protected: + double m_angle1{0}; + double m_angle2{0}; + double m_prob{0}; + SphericalPos m_prev_pos{0., 0.}; +}; + +#endif //BIMODALANGLE_H diff --git a/motions/isosmallangle.cpp b/motions/isosmallangle.cpp index c1ac818..8c58ab8 100644 --- a/motions/isosmallangle.cpp +++ b/motions/isosmallangle.cpp @@ -3,12 +3,11 @@ // #include "isosmallangle.h" +#include "coordinates.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) { -} +SmallAngle::SmallAngle(const double delta, const double eta, const double chi, std::mt19937_64 &rng) : Motion(delta, eta, rng), m_chi(chi * M_PI / 180.0) {}; +SmallAngle::SmallAngle(std::mt19937_64 &rng) : Motion(rng) {} void SmallAngle::initialize() { m_prev_pos = draw_position(); @@ -16,7 +15,7 @@ void SmallAngle::initialize() { double SmallAngle::jump() { const double gamma{2 * M_PI * m_uni_dist(m_rng)}; - m_prev_pos = rotate(m_prev_pos, gamma, m_chi); + m_prev_pos = rotate(m_prev_pos, m_chi, gamma); return omega_q(m_prev_pos); } diff --git a/motions/tetrahedral.cpp b/motions/tetrahedral.cpp index 027bc07..3b65852 100644 --- a/motions/tetrahedral.cpp +++ b/motions/tetrahedral.cpp @@ -9,12 +9,10 @@ TetrahedralJump::TetrahedralJump(const double delta, const double eta, std::mt19 TetrahedralJump::TetrahedralJump(std::mt19937_64& rng) : Motion(rng) {} void TetrahedralJump::initialize() { - // const auto pos = SphericalPos{0., 0}; const auto pos = draw_position(); m_corners[0] = omega_q(pos); const double alpha = 2. * M_PI * m_uni_dist(m_rng); - // const double alpha = 0.; for (int i = 1; i<4; i++) { auto corner_pos = rotate(pos, m_beta, alpha + (i-1) * 2*M_PI/3.); diff --git a/sims.cpp b/sims.cpp index efbffde..892f2bb 100644 --- a/sims.cpp +++ b/sims.cpp @@ -112,6 +112,7 @@ void run_ste(std::unordered_map& parameter, Motion& motion, motion.setEta(parameter["eta"]); const auto start = std::chrono::system_clock::now(); + auto last_print_out = std::chrono::system_clock::now(); const time_t start_time = std::chrono::system_clock::to_time_t(start); std::cout << "Start tau = " << tau << "s : " << ctime(&start_time); @@ -148,6 +149,7 @@ void run_ste(std::unordered_map& parameter, Motion& motion, 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 @@ -164,12 +166,13 @@ void run_ste(std::unordered_map& parameter, Motion& motion, } -void make_trajectory(Motion& motion, const Distribution& dist, const double t_max, std::vector& out_time, std::vector& out_phase) { +void make_trajectory(Motion& motion, Distribution& dist, const double t_max, std::vector& out_time, std::vector& out_phase) { // Starting position double t_passed = 0; double phase = 0; motion.initialize(); + dist.initialize(); out_time.emplace_back(t_passed); out_phase.emplace_back(0); diff --git a/sims.h b/sims.h index ee51500..67b4d6d 100644 --- a/sims.h +++ b/sims.h @@ -13,6 +13,6 @@ void run_spectrum(std::unordered_map& parameter, Motion& motion, Distribution& dist); void run_ste(std::unordered_map& parameter, Motion& motion, Distribution& dist); -void make_trajectory(Motion&, const Distribution&, double, std::vector&, std::vector&); +void make_trajectory(Motion&, Distribution&, double, std::vector&, std::vector&); #endif //RWSIM_SIMS_H diff --git a/test.py b/test.py index 37a062d..4bfb0b0 100644 --- a/test.py +++ b/test.py @@ -58,7 +58,7 @@ def post_process_spectrum(taus, apod, tpulse): def post_process_ste(taus): - for i, tau in enumerate(tau_values): + for i, tau in enumerate(taus): try: raw_data_cc = np.loadtxt(f'coscos_tau={tau:.6e}.dat') diff --git a/times/base.cpp b/times/base.cpp index 5a51e94..b14e967 100644 --- a/times/base.cpp +++ b/times/base.cpp @@ -3,11 +3,11 @@ // #include "base.h" -Distribution::Distribution(const double tau, std::mt19937_64 &rng) : m_tau(tau), m_rng(rng) {} +Distribution::Distribution(const double tau, std::mt19937_64 &rng) : m_tau(tau), m_tau_jump(tau), m_rng(rng) {} Distribution::Distribution(std::mt19937_64 &rng) : m_rng(rng) {} double Distribution::tau_wait() const { - return std::exponential_distribution(1./m_tau)(m_rng); + return std::exponential_distribution(1./m_tau_jump)(m_rng); } diff --git a/times/base.h b/times/base.h index 91eab4f..bdd1e0b 100644 --- a/times/base.h +++ b/times/base.h @@ -17,11 +17,13 @@ public: [[nodiscard]] double getTau() const { return m_tau; } void setTau(const double tau) { m_tau = tau;} + virtual void initialize() = 0; virtual void draw_tau() = 0; [[nodiscard]] double tau_wait() const; -private: +protected: double m_tau{1.}; + double m_tau_jump{1.}; std::mt19937_64& m_rng; }; diff --git a/times/delta.cpp b/times/delta.cpp index d229d99..cec764d 100644 --- a/times/delta.cpp +++ b/times/delta.cpp @@ -6,5 +6,9 @@ DeltaDistribution::DeltaDistribution(const double tau, std::mt19937_64& rng) : Distribution(tau, rng) {} DeltaDistribution::DeltaDistribution(std::mt19937_64& rng) : Distribution(rng) {} +void DeltaDistribution::initialize() { + m_tau_jump = m_tau; +} + void DeltaDistribution::draw_tau() {} diff --git a/times/delta.h b/times/delta.h index 62e920f..4367db6 100644 --- a/times/delta.h +++ b/times/delta.h @@ -12,6 +12,7 @@ public: DeltaDistribution(double, std::mt19937_64&); explicit DeltaDistribution(std::mt19937_64 &rng); + void initialize() override; void draw_tau() override; }; diff --git a/times/lognormal.cpp b/times/lognormal.cpp new file mode 100644 index 0000000..8e2ca6f --- /dev/null +++ b/times/lognormal.cpp @@ -0,0 +1,18 @@ +// +// Created by dominik on 8/24/24. +// + +#include "lognormal.h" +#include + +LogNormalDistribution::LogNormalDistribution(const double tau, const double sigma, std::mt19937_64& rng) : Distribution(tau, rng), m_sigma(sigma), m_distribution(std::log(tau), sigma) {} +LogNormalDistribution::LogNormalDistribution(std::mt19937_64& rng) : Distribution(rng) {} + +void LogNormalDistribution::initialize() { + m_distribution = std::lognormal_distribution(std::log(m_tau), m_sigma); + m_tau_jump = m_distribution(m_rng); +} + +void LogNormalDistribution::draw_tau() { + m_tau_jump = m_distribution(m_rng); +} diff --git a/times/lognormal.h b/times/lognormal.h new file mode 100644 index 0000000..818c35c --- /dev/null +++ b/times/lognormal.h @@ -0,0 +1,21 @@ + +#ifndef LOGNORMAL_H +#define LOGNORMAL_H + +#include "base.h" +#include + +class LogNormalDistribution final : public Distribution { +public: + LogNormalDistribution(double, double, std::mt19937_64&); + explicit LogNormalDistribution(std::mt19937_64 &rng); + + void initialize() override; + void draw_tau() override; + +private: + double m_sigma{1}; + std::lognormal_distribution<> m_distribution; +}; + +#endif //LOGNORMAL_H