222 lines
7.1 KiB
Python
222 lines
7.1 KiB
Python
import pathlib
|
|
import subprocess
|
|
|
|
import numpy as np
|
|
from scipy.optimize import curve_fit
|
|
import matplotlib.pyplot as plt
|
|
|
|
|
|
|
|
def run_sims(taus, ste: bool = True, spectrum: bool = False, exec_file: str = './build/rwsim', config_file: str = './config.txt'):
|
|
for tau in taus:
|
|
arguments = [exec_file, config_file]
|
|
if ste:
|
|
arguments += ['--ste']
|
|
if spectrum:
|
|
arguments += ['--spectrum']
|
|
|
|
arguments += ['BimodalAngle']
|
|
|
|
# arguments += ['-TAU', '1']
|
|
|
|
with pathlib.Path(config_file).open('a') as f:
|
|
f.write(f'tau={tau}\n')
|
|
|
|
subprocess.run(arguments)
|
|
|
|
|
|
|
|
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):
|
|
tevo = np.linspace(1e-6, 120e-6, num=8)
|
|
|
|
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]
|
|
|
|
fig_cc_raw, ax_cc_raw = plt.subplots()
|
|
fig_ss_raw, ax_ss_raw = plt.subplots()
|
|
|
|
# ax_cc_raw.semilogx(t_mix, raw_data_cc[:, 1:], '.')
|
|
ax_ss_raw.semilogx(t_mix, raw_data_ss[:, 1:], '.')
|
|
|
|
scaled_cc = (raw_data_cc[:, 1:]-raw_data_cc[-1, 1:])/(raw_data_cc[0, 1:]-raw_data_cc[-1, 1:])
|
|
scaled_ss = (raw_data_ss[:, 1:]-raw_data_ss[-1, 1:])/(raw_data_ss[0, 1:]-raw_data_ss[-1, 1:])
|
|
|
|
fig_tau, ax_tau = plt.subplots()
|
|
fig_beta, ax_beta= plt.subplots()
|
|
fig_finfty, ax_finfty = plt.subplots()
|
|
|
|
tau_cc_fit = []
|
|
beta_cc_fit = []
|
|
finfty_cc_fit = []
|
|
|
|
tau_ss_fit = []
|
|
beta_ss_fit = []
|
|
finfty_ss_fit = []
|
|
|
|
tau_plus_fit = []
|
|
beta_plus_fit = []
|
|
finfty_plus_fit = []
|
|
|
|
for j in range(1, raw_data_cc.shape[1]-1):
|
|
p0 = [
|
|
raw_data_cc[0, 1],
|
|
t_mix[np.argmin(np.abs(scaled_cc[:, j]-np.exp(-1)))],
|
|
0.5,
|
|
0.1
|
|
]
|
|
|
|
res = curve_fit(ste, t_mix, raw_data_cc[:, j+1], p0, bounds=[(0, 0, 0., 0), (np.inf, np.inf, 1, 1)])
|
|
m0, tauc, beta, finfty = res[0]
|
|
# print(f'Cos-Cos-Fit for {tevo[j]}: tau_c = {tauc:.6e}, beta={beta:.4e}, amplitude = {m0: .4e}, f_infty={finfty:.4e}')
|
|
l = ax_cc_raw.semilogx(t_mix, raw_data_cc[:, j+1], 'x', label=f'{tevo[j]}')
|
|
ax_cc_raw.semilogx(t_mix, ste(t_mix, *res[0]), linestyle='--', color = l[0].get_color())
|
|
|
|
tau_cc_fit.append(res[0][1])
|
|
beta_cc_fit.append(res[0][2])
|
|
finfty_cc_fit.append(res[0][3])
|
|
|
|
p0 = [
|
|
raw_data_cc[0, 1],
|
|
t_mix[np.argmin(np.abs(scaled_ss[:, j]-np.exp(-1)))],
|
|
0.5,
|
|
0.1
|
|
]
|
|
|
|
res = curve_fit(ste, t_mix, raw_data_ss[:, j+1], p0, bounds=[(0, 0, 0, 0), (np.inf, np.inf, 1, 1)])
|
|
m0, tauc, beta, finfty = res[0]
|
|
# print(f'Sin-Sin-Fit for {tevo[j]}: tau_c = {tauc:.6e}, beta={beta:.4e}, amplitude = {m0: .4e}, f_infty={finfty:.4e}')
|
|
|
|
tau_ss_fit.append(res[0][1])
|
|
beta_ss_fit.append(res[0][2])
|
|
finfty_ss_fit.append(res[0][3])
|
|
|
|
p0 = [
|
|
raw_data_cc[0, 1],
|
|
t_mix[np.argmin(np.abs(scaled_ss[:, j]-np.exp(-1)))],
|
|
0.5,
|
|
0.1
|
|
]
|
|
|
|
res = curve_fit(ste, t_mix, raw_data_ss[:, j+1] + raw_data_cc[:, j+1], p0, bounds=[(0, 0, 0, 0), (np.inf, np.inf, 1, 1)])
|
|
m0, tauc, beta, finfty = res[0]
|
|
# print(f'Sin-Sin-Fit for {tevo[j]}: tau_c = {tauc:.6e}, beta={beta:.4e}, amplitude = {m0: .4e}, f_infty={finfty:.4e}')
|
|
|
|
tau_plus_fit.append(res[0][1])
|
|
beta_plus_fit.append(res[0][2])
|
|
finfty_plus_fit.append(res[0][3])
|
|
|
|
# ax_tau.axhline(tau)
|
|
|
|
ax_tau.semilogy(tevo[1:], np.array(tau_cc_fit)/tau_cc_fit[0], 'C0-')
|
|
ax_beta.plot(tevo[1:], beta_cc_fit, 'C0-')
|
|
ax_finfty.plot(tevo[1:], finfty_cc_fit, 'C0-')
|
|
|
|
ax_tau.semilogy(tevo[1:], np.array(tau_ss_fit)/tau_ss_fit[0], 'C1-')
|
|
ax_beta.plot(tevo[1:], beta_ss_fit, 'C1-')
|
|
ax_finfty.plot(tevo[1:], finfty_ss_fit, 'C1-')
|
|
|
|
ax_tau.semilogy(tevo[1:], np.array(tau_plus_fit)/tau_plus_fit[0], 'C2-')
|
|
ax_beta.plot(tevo[1:], beta_plus_fit, 'C2-')
|
|
ax_finfty.plot(tevo[1:], finfty_plus_fit, 'C2-')
|
|
|
|
ax_cc_raw.legend()
|
|
plt.show()
|
|
|
|
# np.savetxt('cc_tauc.dat', list(zip(tevo[1:], tau_cc_fit)))
|
|
# np.savetxt('cc_beta.dat', list(zip(tevo[1:], beta_cc_fit)))
|
|
# np.savetxt('cc_finfty.dat', list(zip(tevo[1:], finfty_cc_fit)))
|
|
# np.savetxt('ss_tauc.dat', list(zip(tevo[1:], tau_ss_fit)))
|
|
# np.savetxt('ss_beta.dat', list(zip(tevo[1:], beta_ss_fit)))
|
|
# np.savetxt('ss_finfty.dat', list(zip(tevo[1:], finfty_ss_fit)))
|
|
|
|
|
|
|
|
# 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()
|
|
|
|
|
|
def ste(x, m0, t, beta, finfty):
|
|
return m0 * ((1-finfty) * np.exp(-(x/t)**beta) + finfty)
|
|
|
|
|
|
if __name__ == '__main__':
|
|
tau_values = np.geomspace(1e-2, 1e-6, num=1) # if num=1, tau is first value
|
|
|
|
# parameter for spectrum simulations
|
|
lb = 2e3
|
|
pulse_length = 2e-6
|
|
|
|
run_sims(tau_values)
|
|
post_process_ste(tau_values)
|
|
# post_process_spectrum(tau_values, lb, pulse_length) |