NMRRWSims/cpp
1
0
Fork 0
Code Issues 3 Pull Requests Actions Packages Projects Releases Wiki Activity

add more motions

Browse Source
This commit is contained in:
Dominik Demuth 2024-12-15 18:34:05 +01:00
parent adbf2af72b
commit bfb1cb314c
26 changed files with 243 additions and 71 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
angle_dependence.py
View File

@ -35,6 +35,18 @@ for fit_files in pathlib.Path('.').glob(f'IsotropicAngle/angle=*/Delta/tau=*/ste
tau_cc.append(fit_values[nearest_idx, 1])
tau_ss.append(fit_values[nearest_idx, 4])
angles = np.array(angles)
tau_cc = np.array(tau_cc)
tau_ss = np.array(tau_ss)
tau2 = np.array(tau2)
sortidx = np.argsort(angles)
angles = angles[sortidx]
tau_cc = tau_cc[sortidx]
tau_ss = tau_ss[sortidx]
tau2 = tau2[sortidx]
plt.semilogy(angles, tau_cc, '-.')
plt.semilogy(angles, tau_ss, '--')

41
process_spectrum.py Normal file
View File

@ -0,0 +1,41 @@
import pathlib
import matplotlib.pyplot as plt
from python.spectrum import *
tpulse = 3e-6
gb = 4e3
fig_spec, ax_spec = plt.subplots()
# for conf_file in pathlib.Path('.').glob(f'SixSiteOctahedral/angle=59/Delta/tau=*/timesignal_*_parameter.txt'):
for conf_file in pathlib.Path('.').glob(f'RandomJumpOnCone/angle=128.*/Delta/tau=*/timesignal_*_parameter.txt'):
vary_string = post_process_spectrum(conf_file, tpulse=tpulse, apod=gb)
print(conf_file)
# ax_spec
#
# ax_tau_cc.semilogy(tau_cc[:, 0], tau_cc[:, 1], label=vary_string)
# ax_tau_cc.axhline(tau_2[:, 1], color='k', linestyle='--')
# ax_beta_cc.plot(*beta_cc.T, label=vary_string)
# ax_finfty_cc.plot(*finfty_cc.T, label=vary_string)
# ax_tau_ss.semilogy(tau_ss[:, 0], tau_ss[:, 1], label=vary_string)
# ax_tau_ss.axhline(tau_2[:, 1], color='k', linestyle='--')
# ax_beta_ss.plot(*beta_ss.T, label=vary_string)
# ax_finfty_ss.plot(*finfty_ss.T, label=vary_string)
#
# np.savetxt(
# conf_file.with_name(f'ste_fit_{vary_string}.dat'),
# np.c_[
# tau_cc, beta_cc[:, 1], finfty_cc[:, 1],
# tau_ss[:, 1], beta_ss[:, 1], finfty_ss[:, 1],
# ],
# header=f'Fit STE {vary_string}\n'
# f'F2: tau={tau_2[0, 1]} beta={beta_2[0, 1]} finfty={finfty_2[0, 1]}\n'
# f'tevo\ttaucc\tbetacc\tfinftycc\ttauss\tbetass\tfinftyss',
# )
#
# for ax in [ax_tau_cc, ax_beta_cc, ax_finfty_cc, ax_tau_ss, ax_beta_ss, ax_finfty_ss]:
# ax.legend()
# plt.show()

4
process_sims.py → process_ste.py
View File

@ -1,7 +1,3 @@
import pathlib
import matplotlib.pyplot as plt
from python.ste import *
fig_tau_cc, ax_tau_cc = plt.subplots()

1
python/helpers.py
View File

@ -1,7 +1,6 @@
from __future__ import annotations
import pathlib
import re
import subprocess
from itertools import product

52
python/spectrum.py
View File

@ -1,8 +1,7 @@
import numpy
import numpy as np
from matplotlib import pyplot
import matplotlib.pyplot as plt
from python.helpers import read_parameter_file
# parameter for spectrum simulations
lb = 2e3
@ -41,30 +40,39 @@ def pulse_attn(freq: np.ndarray, t_pulse: float) -> np.ndarray:
return np.pi * numerator / denominator / 2
def post_process_spectrum(taus, apod, tpulse):
reduction_factor = np.zeros((taus.size, 5)) # hard-coded t_echo :(
def post_process_spectrum(parameter_file, apod, tpulse):
parameter = read_parameter_file(parameter_file)
for i, tau in enumerate(taus):
try:
raw_data = np.loadtxt(f'fid_tau={tau:.6e}.dat')
except OSError:
continue
# files have form ste_arg=0.000000e+01_parameter.txt, first remove ste part then parameter.txt to get variables
varied_string = str(parameter_file).partition('_')[-1].rpartition('_')[0]
t = raw_data[:, 0]
timesignal = raw_data[:, 1:]
# make evolution times
tevo = np.linspace(parameter['techo_start'], parameter['techo_stop'], num=int(parameter['techo_steps']))
timesignal *= dampening(t, apod)[:, None]
timesignal[0, :] /= 2
if varied_string:
raw_data = np.loadtxt(parameter_file.with_name(f'timesignal_{varied_string}.dat'))
else:
raw_data = np.loadtxt(parameter_file.with_name(f'timesignal.dat'))
# FT to spectrum
freq = np.fft.fftshift(np.fft.fftfreq(t.size, d=1e-6))
spec = np.fft.fftshift(np.fft.fft(timesignal, axis=0), axes=0).real
spec *= pulse_attn(freq, t_pulse=tpulse)[:, None]
t = raw_data[:, 0]
timesignal = raw_data[:, 1:]
reduction_factor[i, :] = 2*timesignal[0, :]
timesignal *= dampening(t, apod)[:, None]
timesignal[0, :] /= 2
plt.plot(freq, spec)
plt.show()
# FT to spectrum
freq = np.fft.fftshift(np.fft.fftfreq(t.size, d=parameter['dwell_time']))
spec = np.fft.fftshift(np.fft.fft(timesignal, axis=0), axes=0).real
spec *= pulse_attn(freq, t_pulse=tpulse)[:, None]
plt.semilogx(taus, reduction_factor, '.')
#
#
# reduction_factor[i, :] = 2*timesignal[0, :]
plt.plot(freq, spec)
plt.gca().set_title(varied_string)
plt.show()
#
# plt.semilogx(taus, reduction_factor, '.')
# plt.show()

6
python/ste.py
View File

@ -1,6 +1,5 @@
import pathlib
import numpy
import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit
@ -74,7 +73,10 @@ def fit_ste(
# make evolution times
tevo = np.linspace(parameter['tevo_start'], parameter['tevo_stop'], num=int(parameter['tevo_steps']))
raw_data = np.loadtxt(parameter_file.with_name(f'{prefix}_{varied_string}.dat'))
if varied_string:
raw_data = np.loadtxt(parameter_file.with_name(f'{prefix}_{varied_string}.dat'))
else:
raw_data = np.loadtxt(parameter_file.with_name(f'{prefix}.dat'))
t_mix = raw_data[:, 0]
decay = raw_data[:, 1:]

1
src/CMakeLists.txt
View File

@ -6,7 +6,6 @@ add_subdirectory(utils)
add_library(simulation STATIC sims.cpp sims.h)
target_link_libraries(simulation PRIVATE utils)
add_executable(
rwsim
main.cpp

2
src/motions/CMakeLists.txt
View File

@ -19,5 +19,7 @@ add_library(
bimodalangle.h
sixsitejump.cpp
sixsitejump.h
nsiteconejump.cpp
nsiteconejump.h
)

20
src/motions/base.cpp
View File

@ -1,14 +1,20 @@
#include "base.h"
#include <iostream>
#include "coordinates.h"
#include "bimodalangle.h"
#include "isosmallangle.h"
#include "random.h"
#include "foursitejump.h"
#include "nsiteconejump.h"
#include "sixsitejump.h"
#include "rjoac.h"
#include <stdexcept>
#include "sixsitejump.h"
namespace motions {
BaseMotion::BaseMotion(std::string name, const double delta, const double eta, std::mt19937_64& rng) : m_name(std::move(name)), m_delta(delta), m_eta(eta), m_rng(rng) {
@ -43,8 +49,8 @@ namespace motions {
if (input == "FourSiteTetrahedral")
return new FourSiteTetrahedron(rng);
if (input == "SixSiteOctahedral")
return new SixSiteOctahedron(rng);
if (input == "SixSiteOctahedralC3")
return new SixSiteOctahedronC3(rng);
if (input == "IsotropicAngle")
return new SmallAngle(rng);
@ -55,6 +61,12 @@ namespace motions {
if (input == "BimodalAngle")
return new BimodalAngle(rng);
if (input == "NSiteConeJump")
return new NSiteJumpOnCone(rng);
if (input == "RandomJumpOnCone")
return new RandomJumpOnCone(rng);
throw std::invalid_argument("Invalid input " + input);
}
@ -74,4 +86,4 @@ namespace motions {
os << m.getName();
return os;
}
}
}

2
src/motions/base.h
View File

@ -28,7 +28,7 @@ namespace motions {
[[nodiscard]] double getEta() const { return m_eta; }
void setEta(const double eta) { m_eta = eta; }
[[nodiscard]] std::string getName() const { return m_name; }
[[nodiscard]] double getInitOmega() const { return m_initial_omega; };
[[nodiscard]] double getInitOmega() const { return m_initial_omega; }
[[nodiscard]] virtual std::string toString() const = 0;

6
src/motions/bimodalangle.cpp
View File

@ -8,18 +8,18 @@ namespace motions {
BaseMotion(std::string("BimodalAngle"), delta, eta, rng),
m_angle1(angle1 * M_PI / 180.0),
m_angle2(angle2 * M_PI / 180.0),
m_prob(prob) {};
m_prob(prob) {}
BimodalAngle::BimodalAngle(std::mt19937_64 &rng) : BaseMotion(std::string("BimodalAngle"), rng) {}
void BimodalAngle::initialize() {
m_prev_pos = draw_position();
m_initial_omega = omega_q(m_prev_pos);
};
}
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 = coordinates::rotate(m_prev_pos, angle, gamma);
m_prev_pos = rotate(m_prev_pos, angle, gamma);
return omega_q(m_prev_pos);
}

2
src/motions/bimodalangle.h
View File

@ -6,7 +6,7 @@
#include "coordinates.h"
namespace motions {
class BimodalAngle : public BaseMotion {
class BimodalAngle final : public BaseMotion {
public:
BimodalAngle(double, double, double, double, double, std::mt19937_64& );
explicit BimodalAngle(std::mt19937_64&);

2
src/motions/conewobble.cpp
View File

@ -16,7 +16,7 @@ namespace motions {
double WobbleCone::jump() {
const double real_angle = m_uni_dist(m_rng) * m_angle;
const double phi = 2 * M_PI * m_uni_dist(m_rng);
return omega_q(coordinates::rotate(m_axis, real_angle, phi));
return omega_q(rotate(m_axis, real_angle, phi));
}
void WobbleCone::setParameters(const std::unordered_map<std::string, double> &parameters) {

2
src/motions/conewobble.h
View File

@ -7,7 +7,7 @@
#include <random>
namespace motions {
class WobbleCone: public BaseMotion {
class WobbleCone final: public BaseMotion {
public:
WobbleCone(double, double, double, std::mt19937_64&);
explicit WobbleCone(std::mt19937_64&);

6
src/motions/coordinates.cpp
View File

@ -22,9 +22,9 @@ namespace coordinates {
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
cos_alpha * x + sin_alpha * (-cos_beta * y - sin_beta * x * z) / norm,
cos_alpha * y + sin_alpha * (cos_beta * x - sin_beta * y * z) / norm,
cos_alpha * z + sin_alpha * sin_beta * norm,
};
return xyz_to_spherical(new_pos);

2
src/motions/foursitejump.cpp
View File

@ -16,7 +16,7 @@ namespace motions {
const double alpha = 2. * M_PI * m_uni_dist(m_rng);
for (int i = 1; i<4; i++) {
auto corner_pos = coordinates::rotate(pos, m_beta, alpha + (i-1) * 2*M_PI/3.);
auto corner_pos = rotate(pos, m_beta, alpha + (i-1) * 2*M_PI/3.);
m_corners[i] = omega_q(corner_pos);
}
m_initial_omega = FourSiteTetrahedron::jump();

6
src/motions/isosmallangle.cpp
View File

@ -6,17 +6,17 @@
namespace motions {
SmallAngle::SmallAngle(const double delta, const double eta, const double chi, std::mt19937_64 &rng) :
BaseMotion(std::string("IsotropicAngle"), delta, eta, rng), m_chi(chi * M_PI / 180.0) {};
BaseMotion(std::string("IsotropicAngle"), delta, eta, rng), m_chi(chi * M_PI / 180.0) {}
SmallAngle::SmallAngle(std::mt19937_64 &rng) : BaseMotion(std::string("IsotropicAngle"), rng) {}
void SmallAngle::initialize() {
m_prev_pos = draw_position();
m_initial_omega = omega_q(m_prev_pos);
};
}
double SmallAngle::jump() {
const double gamma{2 * M_PI * m_uni_dist(m_rng)};
m_prev_pos = coordinates::rotate(m_prev_pos, m_chi, gamma);
m_prev_pos = rotate(m_prev_pos, m_chi, gamma);
return omega_q(m_prev_pos);
}

59
src/motions/nsiteconejump.cpp Normal file
View File

@ -0,0 +1,59 @@
//
// Created by dominik on 12/14/24.
//
#include "nsiteconejump.h"
#include "coordinates.h"
#include <cmath>
#include <iostream>
#include <ostream>
#include <vector>
namespace motions {
NSiteJumpOnCone::NSiteJumpOnCone(const double delta, const double eta, const int num_sites, const double chi, std::mt19937_64 &rng) :
BaseMotion("NSiteJumpOnCone", delta, eta, rng),
m_num_sites(num_sites),
m_chi(chi*M_PI/180.) {}
NSiteJumpOnCone::NSiteJumpOnCone(std::mt19937_64 &rng) : BaseMotion("NSiteJumpOnCone", rng) { }
void NSiteJumpOnCone::initialize() {
m_sites = std::vector<double>(m_num_sites);
m_chooser = std::uniform_int_distribution<>{1, m_num_sites - 1};
m_axis = draw_position();
const double alpha = m_uni_dist(m_rng) * 2 * M_PI;
const double steps = 2*M_PI / m_num_sites;
for (int i = 0; i < m_num_sites; i++) {
m_sites[i] = omega_q(rotate(m_axis, m_chi, i * steps + alpha));
}
}
double NSiteJumpOnCone::jump() {
m_cone_idx += m_chooser(m_rng);
m_cone_idx %= m_num_sites;
return m_sites[m_cone_idx];
}
void NSiteJumpOnCone::setParameters(const std::unordered_map<std::string, double> &parameters) {
BaseMotion::setParameters(parameters);
m_num_sites = static_cast<int>(parameters.at("num_sites"));
m_chi = parameters.at("chi") * M_PI/180.;
}
std::unordered_map<std::string, double> NSiteJumpOnCone::getParameters() const {
auto parameter = BaseMotion::getParameters();
parameter["num_sites"] = m_num_sites;
parameter["chi"] = m_chi * 180. / M_PI;
return parameter;
}
std::string NSiteJumpOnCone::toString() const {
return std::to_string(m_num_sites) + "SiteJumpOnCone/angle=" + std::to_string(m_chi*180/M_PI);
}
} // motions

33
src/motions/nsiteconejump.h Normal file
View File

@ -0,0 +1,33 @@
#ifndef NSITEJUMPONCONE_H
#define NSITEJUMPONCONE_H
#include "base.h"
#include <random>
#include <vector>
namespace motions {
class NSiteJumpOnCone final : public BaseMotion {
public:
NSiteJumpOnCone(double, double, int, double, std::mt19937_64&);
explicit NSiteJumpOnCone(std::mt19937_64&);
void initialize() override;
double jump() override;
[[nodiscard]] std::string toString() const override;
void setParameters(const std::unordered_map<std::string, double> &) override;
[[nodiscard]] std::unordered_map<std::string, double> getParameters() const override;
private:
int m_num_sites{2};
std::vector<double> m_sites{};
double m_chi{1./2.};
int m_cone_idx = 0;
coordinates::SphericalPos m_axis{1, 0};
std::uniform_int_distribution<> m_chooser;
};
} // motions
#endif //NSITEJUMPONCONE_H

2
src/motions/rjoac.cpp
View File

@ -11,7 +11,7 @@ namespace motions {
double RandomJumpOnCone::jump() {
const double phi = 2 * M_PI * m_uni_dist(m_rng);
return omega_q(coordinates::rotate(m_axis, m_angle, phi));
return omega_q(rotate(m_axis, m_angle, phi));
}
void RandomJumpOnCone::setParameters(const std::unordered_map<std::string, double> &parameters) {

2
src/motions/rjoac.h
View File

@ -7,7 +7,7 @@
#include <random>
namespace motions {
class RandomJumpOnCone: public BaseMotion {
class RandomJumpOnCone final: public BaseMotion {
public:
RandomJumpOnCone(double, double, double, std::mt19937_64&);
explicit RandomJumpOnCone(std::mt19937_64&);

35
src/motions/sixsitejump.cpp
View File

@ -1,34 +1,41 @@
#include "sixsitejump.h"
#include <iomanip>
#include <iostream>
#include "coordinates.h"
namespace motions {
SixSiteOctahedron::SixSiteOctahedron(const double delta, const double eta, std::mt19937_64& rng) :
SixSiteOctahedronC3::SixSiteOctahedronC3(const double delta, const double eta, const double chi, std::mt19937_64& rng) :
m_chi{chi*M_PI/180.},
BaseMotion(std::string{"SixSiteOctahedral"}, delta, eta, rng) {}
SixSiteOctahedron::SixSiteOctahedron(std::mt19937_64& rng) : BaseMotion(std::string{"SixSiteOctahedral"}, rng) {}
void SixSiteOctahedron::initialize() {
const auto pos = draw_position();
m_corners[0] = omega_q(pos);
m_corners[1] = omega_q(rotate(pos, M_PI, 0));
SixSiteOctahedronC3::SixSiteOctahedronC3(std::mt19937_64& rng) : BaseMotion(std::string{"SixSiteOctahedralC3"}, rng) {}
void SixSiteOctahedronC3::initialize() {
const coordinates::SphericalPos c3_axis = draw_position();
const auto [x, y, z] = spherical_to_xyz(c3_axis);
const double alpha = 2. * M_PI * m_uni_dist(m_rng);
for (int i = 2; i<6; i++) {
auto corner_pos = coordinates::rotate(pos, M_PI_2, alpha + (i-2) * M_PI_2);
m_corners[i] = omega_q(corner_pos);
const double m_chi_opposite = M_PI - m_chi;
for (int i = 0; i<3; i++) {
m_corners[2*i] = omega_q(rotate(c3_axis, m_chi, alpha + i * 2./3.*M_PI));
m_corners[2*i+1] = omega_q(rotate(c3_axis, m_chi_opposite, alpha + i * 2./3.*M_PI + M_PI/3.));
}
m_initial_omega = SixSiteOctahedron::jump();
m_initial_omega = SixSiteOctahedronC3::jump();
}
double SixSiteOctahedron::jump() {
double SixSiteOctahedronC3::jump() {
m_corner_idx += m_chooser(m_rng);
m_corner_idx %= 6;
return m_corners[m_corner_idx];
}
std::string SixSiteOctahedron::toString() const {
return {"SixSiteOctahedral"};
std::string SixSiteOctahedronC3::toString() const {
return {"SixSiteOctahedral/angle=" + std::to_string(m_chi / M_PI * 180.)};
}
}

8
src/motions/sixsitejump.h
View File

@ -9,10 +9,10 @@
#include <array>
namespace motions {
class SixSiteOctahedron final : public BaseMotion {
class SixSiteOctahedronC3 final : public BaseMotion {
public:
SixSiteOctahedron(double, double, std::mt19937_64&);
explicit SixSiteOctahedron(std::mt19937_64&);
SixSiteOctahedronC3(double, double, double, std::mt19937_64&);
explicit SixSiteOctahedronC3(std::mt19937_64&);
void initialize() override;
double jump() override;
@ -20,7 +20,7 @@ namespace motions {
[[nodiscard]] std::string toString() const override;
private:
const double m_beta{std::acos(-1/3.)};
const double m_chi{0.95531661812450927816385710251575775424341469501000549095969812932191204590}; // 54.74 deg
std::array<double, 6> m_corners{};
int m_corner_idx{0};

2
src/times/base.h
View File

@ -14,7 +14,7 @@ namespace times {
[[nodiscard]] double getTau() const { return m_tau; }
void setTau(const double tau) { m_tau = tau; }
[[nodiscard]] std::string getName() const { return m_name; };
[[nodiscard]] std::string getName() const { return m_name; }
virtual void setParameters(const std::unordered_map<std::string, double>&);
[[nodiscard]] virtual std::unordered_map<std::string, double> getParameters() const;

4
src/utils/io.cpp
View File

@ -25,7 +25,7 @@ std::pair<std::string, double> get_optional_parameter(std::vector<std::string>::
stripped_arg = it->substr(1, it->size());
}
std::transform(stripped_arg.begin(), stripped_arg.end(), stripped_arg.begin(), [](unsigned char c) { return std::tolower(c); });
const auto stripped_value = std::stod(*(++it), nullptr);
const auto stripped_value = std::stod(*++it, nullptr);
return std::make_pair(stripped_arg, stripped_value);
}
@ -95,7 +95,7 @@ Arguments parse_args(const int argc, char* argv[]) {
std::unordered_map<std::string, double> read_parameter(const std::filesystem::path& infile) {
if (!std::filesystem::exists(infile)) {
if (!exists(infile)) {
std::cerr << "File " << infile << " does not exist" << std::endl;
exit(1);
}

2
start_sims.py
View File

@ -1,3 +1,5 @@
import numpy as np
from python.helpers import *
#Simulation parameter