add 4 pulse ste

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

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()

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

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

src/simulation/sims.cpp

@@ -87,9 +87,10 @@ void run_ste(
     ) {
     const int num_walker = static_cast<int>(parameter[std::string("num_walker")]);
 
-    const int num_mix_times = static_cast<int>(parameter[std::string("tmix_steps")]);
+    const int num_mix_times = static_cast<int>(parameter["tmix_steps"]);
     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 = 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<double, std::vector<double>> cc_dict;
@@ -103,7 +104,7 @@ void run_ste(
     std::vector<double> 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;