cpp/test.py

import pathlib
import subprocess

import numpy as np
import matplotlib.pyplot as plt



def run_sims(taus):
    for tau in taus:
        with pathlib.Path('./config.txt').open('a') as f:
            f.write(f'tau={tau}\n')
        subprocess.run(['./cmake-build-debug/rwsim', '-ste', './config.txt'])


def dampening(x: np.ndarray, apod: float) -> np.ndarray:
    return np.exp(-apod * x)


def pulse_attn(freq: np.ndarray, t_pulse: float):
    # cf. Schmitt-Rohr/Spieß eq. 2.126; omega_1 * t_p = pi/2
    pi_half_squared = np.pi**2 / 4
    omega = 2 * np.pi * freq

    numerator = np.sin(np.sqrt(pi_half_squared + omega**2 * t_pulse**2 / 2))
    denominator = np.sqrt(pi_half_squared + omega**2 * t_pulse**2 / 4)

    return np.pi * numerator / denominator / 2


def post_process_spectrum(taus, apod, tpulse):
    reduction_factor = np.zeros((taus.size, 5)) # hard-coded t_echo :(

    for i, tau in enumerate(taus):
        try:
            raw_data = np.loadtxt(f'fid_tau={tau:.6e}.dat')
        except OSError:
            continue

        t = raw_data[:, 0]
        timesignal = raw_data[:, 1:]

        timesignal *= dampening(t, apod)[:, None]
        timesignal[0, :] /= 2

        # 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]

        reduction_factor[i, :] = 2*timesignal[0, :]

        plt.plot(freq, spec)
        plt.show()

    plt.semilogx(taus, reduction_factor, '.')
    plt.show()


def post_process_ste(taus):
    for i, tau in enumerate(taus):

        try:
            raw_data_cc = np.loadtxt(f'coscos_tau={tau:.6e}.dat')
            raw_data_ss = np.loadtxt(f'sinsin_tau={tau:.6e}.dat')
        except OSError:
            continue

        t_mix = raw_data_cc[:, 0]

        plt.semilogx(t_mix, raw_data_cc[:, 1:])
        plt.show()

        plt.semilogx(t_mix, raw_data_ss[:, 1:])
        plt.show()

        plt.plot(raw_data_cc[0, 1:])
        plt.plot(raw_data_ss[0, 1:])
        plt.show()

        plt.plot(raw_data_cc[-1, 1:]/raw_data_cc[0, 1:])
        plt.plot(raw_data_ss[-1, 1:]/raw_data_ss[0, 1:])
        plt.show()


if __name__ == '__main__':
    tau_values = np.logspace(-3, -2, 1)
    lb = 2e3
    pulse_length = 2e-6

    run_sims(tau_values)
    post_process_ste(tau_values)
    # post_process_spectrum(tau_values, lb, pulse_length)