From 60ef1b0bcff16c762adc5bbd69802615575b3ef9 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Sat, 30 Nov 2024 16:15:38 +0100 Subject: [PATCH] add 4 pulse ste --- config.txt | 8 +++----- main.py | 24 +++++++++++++++++++----- python/helpers.py | 1 + python/ste.py | 9 ++------- src/simulation/sims.cpp | 16 +++++++++++----- 5 files changed, 36 insertions(+), 22 deletions(-) diff --git a/config.txt b/config.txt index c2a1337..b62c25f 100644 --- a/config.txt +++ b/config.txt @@ -5,13 +5,10 @@ delta=126e3 eta=0.0 # Distribution part tau=1e-3 -angle1=2 -angle2=30 -probability1=0 # Spectrum part dwell_time=1e-8 num_acq=4096 -techo_start=0e-6 +techo_start=1e-6 techo_stop=40e-6 techo_steps=5 # STE part @@ -20,4 +17,5 @@ tevo_stop=120e-6 tevo_steps=121 tmix_start=1e-5 tmix_stop=1e1 -tmix_steps=61 +tmix_steps=31 +tpulse4=10e-6 diff --git a/main.py b/main.py index 582ae99..4f2a07b 100644 --- a/main.py +++ b/main.py @@ -8,7 +8,7 @@ motion = 'IsotropicAngle' distribution = 'Delta' # parameter = {} parameter = { - "angle": [10, 109.47], + "angle": [10], } parameter = prepare_rw_parameter(parameter) @@ -33,13 +33,13 @@ ax_finfty_ss.set_title('f_infty_ss') for variation in parameter: print(f"\nRun RW for {motion}/{distribution} with arguments {variation}\n") - run_sims(motion, distribution, ste=True, spectrum=False, **variation) + # run_sims(motion, distribution, ste=True, spectrum=False, **variation) conf_file = find_config_file(motion, distribution, variation) - vary_string, tau_cc, beta_cc, finfty_cc = fit_and_save_ste(conf_file, 'coscos', plot_decays=False, verbose=False) - _, tau_ss, beta_ss, finfty_ss = fit_and_save_ste(conf_file, 'sinsin', plot_decays=False, verbose=False) - _, tau_2, beta_2, finfty_2 = fit_and_save_ste(conf_file, 'f2', plot_decays=True, verbose=True) + vary_string, tau_cc, beta_cc, finfty_cc = fit_ste(conf_file, f'coscos', plot_decays=False, verbose=False) + _, tau_ss, beta_ss, finfty_ss = fit_ste(conf_file, f'sinsin', plot_decays=False, verbose=False) + _, tau_2, beta_2, finfty_2 = fit_ste(conf_file, f'f2', plot_decays=True, verbose=True) ax_tau_cc.semilogy(tau_cc[:, 0], tau_cc[:, 1], label=vary_string) ax_tau_cc.axhline(tau_2[:, 1], color='k', linestyle='--') @@ -50,6 +50,20 @@ for variation in parameter: ax_beta_ss.plot(*beta_ss.T, label=vary_string) ax_finfty_ss.plot(*finfty_ss.T, label=vary_string) + + np.savetxt( + 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() \ No newline at end of file diff --git a/python/helpers.py b/python/helpers.py index 42c5786..a5505a3 100644 --- a/python/helpers.py +++ b/python/helpers.py @@ -5,6 +5,7 @@ import re import subprocess from itertools import product + def prepare_rw_parameter(parameter: dict) -> list: """ Converts a dictionary of iterables to list of dictionaries diff --git a/python/ste.py b/python/ste.py index ddf5a45..6d6c494 100644 --- a/python/ste.py +++ b/python/ste.py @@ -59,12 +59,12 @@ def fit_decay(x: np.ndarray, y: np.ndarray, tevo: np.ndarray, verbose: bool = Tr return tau_fit, beta_fit, finfty_fit -def fit_and_save_ste( +def fit_ste( parameter_file: pathlib.Path, prefix: str, plot_decays: bool = True, verbose: bool = True, -) -> tuple[str, np.ndarray, np.ndarray, np.ndarray]: + ) -> tuple[str, np.ndarray, np.ndarray, np.ndarray]: # read simulation parameters parameter = read_parameter_file(parameter_file) @@ -89,9 +89,4 @@ def fit_and_save_ste( print(f'Fit {prefix}') tau, beta, finfty = fit_decay(t_mix, decay, tevo, verbose=verbose) - - np.savetxt(f'tau_{prefix}_{varied_string}.dat', tau) - np.savetxt(f'beta_{prefix}_{varied_string}.dat', beta) - np.savetxt(f'finfty_{prefix}_{varied_string}.dat', finfty) - return varied_string, tau, beta, finfty diff --git a/src/simulation/sims.cpp b/src/simulation/sims.cpp index 5994c12..e6adf20 100644 --- a/src/simulation/sims.cpp +++ b/src/simulation/sims.cpp @@ -87,9 +87,10 @@ void run_ste( ) { const int num_walker = static_cast(parameter[std::string("num_walker")]); - const int num_mix_times = static_cast(parameter[std::string("tmix_steps")]); + const int num_mix_times = static_cast(parameter["tmix_steps"]); const std::vector evolution_times = linspace(parameter["tevo_start"], parameter["tevo_stop"], static_cast(parameter["tevo_steps"])); const std::vector mixing_times = logspace(parameter["tmix_start"], parameter["tmix_stop"], num_mix_times); + const double tpulse4= parameter["tpulse4"]; // make ste decay vectors and set them to zero std::map> cc_dict; @@ -103,7 +104,7 @@ void run_ste( std::vector f2(num_mix_times); // 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()); + const double tmax = *std::max_element(evolution_times.begin(), evolution_times.end()) * 2 + *std::max_element(mixing_times.begin(), mixing_times.end()) + 2*tpulse4; // set parameter in distribution and motion model dist.setParameters(parameter); @@ -127,6 +128,7 @@ void run_ste( f2[f2_idx] += traj_omega[f2_pos] * motion.getInitOmega() / num_walker; } + for (auto& [t_evo_j, _] : cc_dict) { auto& cc_j = cc_dict[t_evo_j]; auto& ss_j = ss_dict[t_evo_j]; @@ -138,16 +140,20 @@ void run_ste( const double ss_tevo = std::sin(dephased); for (int mix_idx = 0; mix_idx < num_mix_times; mix_idx++) { - // get phase at end of mixing time 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); + // get phase at position of 4th pulse + const double time_pulse4 = time_end_mix + tpulse4; + current_pos = nearest_index(traj_time, time_pulse4, current_pos); + const double phase_4pulse = lerp(traj_time, traj_phase, time_pulse4, current_pos); + // get phase at echo position - const double time_echo = mixing_times[mix_idx] + 2 * t_evo_j; + const double time_echo = time_pulse4 + tpulse4 + 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; + double rephased = lerp(traj_time, traj_phase, time_echo, current_pos) + phase_mix_end - 2*phase_4pulse; cc_j[mix_idx] += cc_tevo * std::cos(rephased) / num_walker; ss_j[mix_idx] += ss_tevo * std::sin(rephased) / num_walker;