93 lines
2.4 KiB
Python
93 lines
2.4 KiB
Python
|
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(tau_values):
|
||
|
|
|
||
|
|
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)
|