fixed problem in spherical_to_xyz

functions.cpp

@@ -3,8 +3,11 @@
 //#
 
 #include <vector>
+#include <array>
 #include <chrono>
 #include <iostream>
+#include <cmath>
+#include <utility>
 
 
 int nearest_index(const std::vector<double> &x_ref, const double x, int start=0) {
@@ -28,6 +31,21 @@ 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(
     /*
      * Prints roughly every 10 seconds how many runs were done and gives a time estimation

functions.h

@@ -2,12 +2,19 @@
 #define RWSIM_FUNCTIONS_H
 
 #include <vector>
+#include <array>
+#include <utility>
+#include <chrono>
 
 
 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

io.cpp

@@ -15,7 +15,39 @@
 #include <filesystem>
 
 
-std::unordered_map<std::string, double> parse_arguments(const std::filesystem::path& infile) {
+
+Arguments parse_args(const int argc, char **argv) {
+    if (argc < 2) {
+        throw std::runtime_error("Not enough arguments: missing parameter file");
+    }
+
+    Arguments args;
+
+    for (int i=1; i<argc; i++) {
+        if (std::string arg = argv[i]; arg[0] == '-') {
+            if (arg == "-ste") {
+                args.ste = true;
+            } else if (arg == "-spectrum") {
+                args.spectrum = true;
+            } else {
+                throw std::runtime_error("Unrecognized option: " + arg);
+            }
+        } else if (args.parameter_file.empty()) {
+            args.parameter_file = arg;
+        } else {
+            throw std::runtime_error("Too many options for parameter file");
+        }
+    }
+
+    if (args.parameter_file.empty()) {
+        throw std::runtime_error("Missing parameter file");
+    }
+
+    return args;
+}
+
+
+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);
@@ -32,6 +64,8 @@ std::unordered_map<std::string, double> parse_arguments(const std::filesystem::p
     size_t delim_pos;
 
     while (std::getline(instream, line)) {
+        if (line[0] == '#' || line.length() == 1) continue;
+
         line.erase(std::remove(line.begin(), line.end(), ' '), line.end());
         delim_pos = line.find('=');
         key = line.substr(0, delim_pos);

io.h

@@ -8,7 +8,15 @@
 #include <filesystem>
 #include <vector>
 
-std::unordered_map<std::string, double> parse_arguments(const std::filesystem::path&);
+struct Arguments {
+    std::string parameter_file{};
+    bool ste = false;
+    bool spectrum = false;
+};
+
+Arguments parse_args(int argc, char **argv);
+
+std::unordered_map<std::string, double> read_parameter(const std::filesystem::path&);
 
 void fid_write_out(const std::string&, const std::vector<double>&, const std::vector<double>&, double, double);
 void fid_write_out(const std::string&, const std::vector<double>&, const std::map<double, std::vector<double>>&, double tau);

main.cpp

@@ -5,25 +5,31 @@
 #include "io.h"
 #include "sims.h"
 #include "motions/random.h"
+#include "motions/tetrahedral.h"
 #include "times/delta.h"
 
 
-int main (int argc, char *argv[]) {
+int main (const int argc, char *argv[]) {
-
+    Arguments args;
-    if (argc < 2) {
+    try {
-        std::cerr << "Usage: " << argv[0] << " PARAMETER_FILE" << std::endl;
+        args = parse_args(argc, argv);
+    } catch (std::runtime_error& error) {
+        std::cerr << error.what() << std::endl;
         return 1;
     }
-    std::cout << argv[0] << std::endl;
 
-    std::unordered_map parameter { parse_arguments(argv[1]) };
+    std::unordered_map parameter { read_parameter(args.parameter_file) };
 
     std::random_device rd;
     std::mt19937_64 rng(rd());
 
-    auto motion = RandomJump(rng);
+    auto motion = TetrahedralJump(rng);
     auto dist = DeltaDistribution(rng);
 
+    if (args.spectrum) {
        run_spectrum(parameter, motion, dist);
-
+    }
+    if (args.ste) {
+        run_ste(parameter, motion, dist);
+    }
 }

motions/base.cpp

@@ -9,6 +9,9 @@
 
 #include "base.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.);
@@ -22,7 +25,7 @@ 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;
 
-    return M_PI * m_delta * (3 * cos_theta_square - 1 - m_eta * sin_theta_square * std::cos(2.*phi));
+    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() {

motions/tetrahedral.cpp

@@ -4,51 +4,50 @@
 #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();
+    const auto [cos_theta, phi] = draw_position();
-    //
+    m_corners[0] = omega_q(cos_theta, phi);
-    // m_corners[0] = omega_q(cos_theta, phi);
+    const double alpha = 2. * M_PI * m_uni_dist(m_rng);
-    //
-    // const double alpha = 2. * M_PI * m_uni_dist(m_rng);
-    // std::cout << alpha << std::endl;
 
-    // v1[0] = sin(phi0)*sin(teta0);
+    const auto vec = spherical_to_xyz(cos_theta, phi);
-    // v2[0] = cos(phi0)*sin(teta0);
+    const double norm = std::sqrt(1 - cos_theta * cos_theta) + 1e-15;
-    // v3[0] = cos(teta0);
+
-    //
+    for (int i = 1; i<4; i++) {
-    // abs = v1[0]*v1[0]+v2[0]*v2[0]+v3[0]*v3[0];  // Normalization
+        const double cos_alpha = std::cos(alpha + (i-1) * 2*M_PI / 3.);
-    // v1[0] = v1[0]/abs;
+        const double sin_alpha = std::sin(alpha + (i-1) * 2*M_PI / 3.);
-    // v2[0] = v2[0]/abs;
+        std::array<double, 3> rotated_position{};
-    // v3[0] = v3[0]/abs;
+        if (cos_theta != 1 && cos_theta != -1) {
-    //
+            // std::cout << cos_theta << std::endl;
-    // /*Calculating the components of the other orientationts. Compare Master Thesis M. Sattig,corrected Version (Share/Abschlussarbeiten)*/
+            rotated_position = {
-    //     norm = sqrt(1 - v3[0]*v3[0])+1E-12;
+                0*m_cos_beta * vec[0] + m_sin_beta * (-vec[0] * vec[2] * sin_alpha - vec[1] * cos_alpha) / norm,
-    //     normI = 1.0/(norm);
+                0*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
-    //     cosgama = cos(delta0);
+            };
-    //     singama = sin(delta0);
+        } else {
-    //     v1[1] = v1[0]*cosBETA + sinBETA*normI*(-v1[0]*v3[0]*singama - v2[0]*cosgama);
+            rotated_position = {
-    //     v2[1] = v2[0]*cosBETA + sinBETA*normI*(-v2[0]*v3[0]*singama + v1[0]*cosgama);
+                m_sin_beta * cos_alpha,
-    //     v3[1] = v3[0]*cosBETA + sinBETA*norm*singama;
+                m_sin_beta * sin_alpha,
-    //
+                m_cos_beta * cos_theta
-    //     cosgama = cos(gama + delta0);
+            };
-    //     singama = sin(gama + delta0);
+        }
-    //     v1[2] = v1[0]*cosBETA + sinBETA*normI*(-v1[0]*v3[0]*singama - v2[0]*cosgama);
+        auto [new_cos_theta, new_phi] = xyz_to_spherical(rotated_position);
-    //     v2[2] = v2[0]*cosBETA + sinBETA*normI*(-v2[0]*v3[0]*singama + v1[0]*cosgama);
+        m_corners[i] = omega_q(new_cos_theta, new_phi);
-    //     v3[2] = v3[0]*cosBETA + sinBETA*norm*singama;
+    }
-    //
-    //     cosgama = cos(2.0*gama + delta0);
-    //     singama = sin(2.0*gama + delta0);
-    //     v1[3] = v1[0]*cosBETA + sinBETA*normI*(-v1[0]*v3[0]*singama - v2[0]*cosgama);
-    //     v2[3] = v2[0]*cosBETA + sinBETA*normI*(-v2[0]*v3[0]*singama + v1[0]*cosgama);
-    //     v3[3] = v3[0]*cosBETA + sinBETA*norm*singama;
 }
 
 double TetrahedralJump::jump() {
-    return 0.;
+    m_corner_idx += m_chooser(m_rng);
+    m_corner_idx %= 4;
+
+    return m_corners[m_corner_idx];
 }

motions/tetrahedral.h

@@ -20,7 +20,14 @@ 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};
+
+    std::uniform_int_distribution<> m_chooser{1, 3};
+
 };
 
 

sims.cpp

@@ -10,50 +10,44 @@
 #include <cmath>
 #include <chrono>
 
-#include "functions.h"
 #include "motions/base.h"
 #include "times/base.h"
+#include "functions.h"
 #include "ranges.h"
 #include "sims.h"
-
-#include <bits/fs_fwd.h>
-
 #include "io.h"
 
 
 void run_spectrum(std::unordered_map<std::string, double>& parameter, Motion& motion, Distribution& dist) {
-    const int num_acq = static_cast<int>(parameter["num_acq"]);
     const int num_walker = static_cast<int>(parameter["num_walker"]);
 
-    const std::vector<double> correlation_times = logspace(parameter["tau_start"], parameter["tau_stop"], static_cast<int>(parameter["tau_steps"]));
-    motion.setDelta(parameter["delta"]);
-    motion.setEta(parameter["eta"]);
-
     // time axis for all time signals
-    const auto t_fid = arange(num_acq, parameter["dwell_time"]);
+    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_evo_i: echo_times) {
         fid_dict[t_evo_i] = std::vector<double>(num_acq);
         std::fill(fid_dict[t_evo_i].begin(), fid_dict[t_evo_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();
 
-    for (const auto tau_i: correlation_times) {
+    // set parameter in distribution and motion model
-        auto start = std::chrono::system_clock::now();
+    const double tau = parameter["tau"];
+    dist.setTau(tau);
+    motion.setDelta(parameter["delta"]);
+    motion.setEta(parameter["eta"]);
+
+    const auto start = std::chrono::system_clock::now();
     auto last_print_out = std::chrono::system_clock::now();
-        time_t start_time = std::chrono::system_clock::to_time_t(start);
+    const time_t start_time = std::chrono::system_clock::to_time_t(start);
-        std::cout << "Start tau = " <<  tau_i << "s : " << ctime(&start_time);
+    std::cout << "Start tau = " <<  tau << "s : " << ctime(&start_time);
 
-        dist.setTau(tau_i);
+    // let the walker walk
-
-        for (auto& [_, fid_j]: fid_dict) {
-            std::fill(fid_j.begin(), fid_j.end(), 0.);
-        }
-
-        // reset array for each correlation time
     for (int mol_i = 0; mol_i < num_walker; mol_i++){
         std::vector<double> traj_time{};
         std::vector<double> traj_phase{};
@@ -63,9 +57,8 @@ void run_spectrum(std::unordered_map<std::string, double>& parameter, Motion& mo
         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_tevo = lerp(traj_time, traj_phase, t_echo_j, current_pos);
+            const double phase_techo = lerp(traj_time, traj_phase, t_echo_j, current_pos);
 
-                // time axis by echo delay to get time in trajectory
 
             for (int acq_idx = 0; acq_idx < num_acq; acq_idx++) {
                 const double real_time = t_fid[acq_idx] + 2 * t_echo_j;
@@ -73,80 +66,93 @@ void run_spectrum(std::unordered_map<std::string, double>& parameter, Motion& mo
                 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_tevo) / num_walker;
+                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
-        fid_write_out("fid", t_fid, fid_dict, tau_i);
+    fid_write_out("fid", t_fid, fid_dict, tau);
 
-        auto end = std::chrono::system_clock::now();
+    const auto end = std::chrono::system_clock::now();
 
     std::chrono::duration<float> duration = end - start;
-        time_t end_time = std::chrono::system_clock::to_time_t(end);
+    const time_t end_time = std::chrono::system_clock::to_time_t(end);
-        std::cout << "End   tau = " << tau_i << "s : " << ctime(&end_time);
+    std::cout << "End   tau = " << tau << "s : " << ctime(&end_time);
     std::cout << "Duration: " << duration.count() << "s\n" << std::endl;
-    }
+
 }
 
+
 void run_ste(std::unordered_map<std::string, double>& parameter, Motion& motion, Distribution& dist) {
-    const int num_acq = static_cast<int>(parameter[std::string("num_acq")]);
     const int num_walker = static_cast<int>(parameter[std::string("num_walker")]);
 
-    const double tau = parameter["tau"];
+    const int num_mix_times = static_cast<int>(parameter[std::string("tmix_steps")]);
-    dist.setTau(tau);
-
-    motion.setDelta(parameter["delta"]);
-    motion.setEta(parameter["eta"]);
-
     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 = 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);
 
-    std::map<double, std::vector<double>> fid_dict;
+
+    // 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) {
-        fid_dict[t_evo_i] = std::vector<double>(num_acq);
+        cc_dict[t_evo_i] = std::vector<double>(num_mix_times);
-        std::fill(fid_dict[t_evo_i].begin(), fid_dict[t_evo_i].end(), 0.);
+        ss_dict[t_evo_i] = std::vector<double>(num_mix_times);
+        std::fill(ss_dict[t_evo_i].begin(), ss_dict[t_evo_i].end(), 0.);
     }
 
     // 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
+    const double tau = parameter["tau"];
+    dist.setTau(tau);
+    motion.setDelta(parameter["delta"]);
+    motion.setEta(parameter["eta"]);
+
     const auto start = 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);
 
-    for (auto& [_, fid_j]: fid_dict) {
+    // let the walker walk
-        std::fill(fid_j.begin(), fid_j.end(), 0.);
-    }
-
-    // reset array for each correlation time
     for (int mol_i = 0; mol_i < num_walker; mol_i++){
         std::vector<double> traj_time{};
         std::vector<double> traj_phase{};
 
         make_trajectory(motion, dist, tmax, traj_time, traj_phase);
 
-        for (auto& [t_evo_j, fid_j] : fid_dict) {
+        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 phase_tevo = lerp(traj_time, traj_phase, t_evo_j, current_pos);
+            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);
 
-            // time axis by echo delay to get time in trajectory
+            for (int mix_idx = 0; mix_idx < num_mix_times; mix_idx++) {
 
-            for (int acq_idx = 0; acq_idx < num_acq; acq_idx++) {
+                // get phase at end of mixing time
-                const double real_time = mixing_times[acq_idx] + t_evo_j;
+                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);
 
-                current_pos = nearest_index(traj_time, real_time, current_pos);
+                // get phase at echo position
-                const double phase_acq = lerp(traj_time, traj_phase, real_time, current_pos);
+                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;
 
-                fid_j[acq_idx] += cos(phase_acq - 2 * phase_tevo);
+                cc_j[mix_idx] += cc_tevo * std::cos(rephased) / num_walker;
+                ss_j[mix_idx] += ss_tevo * std::sin(rephased) / num_walker;
             }
         }
     }
 
-    // write fid to files
+    // write to files
-    fid_write_out("ste", mixing_times, fid_dict, tau);
+    fid_write_out("coscos", mixing_times, cc_dict, tau);
+    fid_write_out("sinsin", mixing_times, ss_dict, tau);
 
     const auto end = std::chrono::system_clock::now();
 
@@ -171,7 +177,6 @@ void make_trajectory(Motion& motion, const Distribution& dist, const double t_ma
     while (t_passed < t_max) {
         const double t = dist.tau_wait();
         t_passed += t;
-
         phase += motion.jump() * t;
 
         out_time.emplace_back(t_passed);

sims.h

@@ -12,6 +12,7 @@
 #include "times/base.h"
 
 void run_spectrum(std::unordered_map<std::string, double>& parameter, Motion& motion, Distribution& dist);
+void run_ste(std::unordered_map<std::string, double>& parameter, Motion& motion, Distribution& dist);
 void make_trajectory(Motion&, const Distribution&, double, std::vector<double>&, std::vector<double>&);
 
 #endif //RWSIM_SIMS_H