PythonでBloch方程式を可視化する④〜ここまでのまとめ〜

明示していなかったかもしれないが，ここまでの目標は，Bloch方程式に基づく磁気共鳴スペクトルを出力することであった．

やりたいことの第一段階までは実装できたので，まとめておく．

import numpy as np
from scipy.integrate import odeint
from matplotlib import cm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

N = 1024
time = 10 ** -9
ga = -1.7608 * 10 ** 11 # 電子の磁気回転比/T * s ** -1
M0 = 1

# Bloch方程式（回転座標系）
def func(v, t, om, H0, H1, T1, T2):
return [ga*H0*v[1]-v[1]*om-v[1]/T2,
-ga*H0*v[0]+v[0]*om+ga*H1*v[2]-v[1]/T2,
-ga*H1*v[1]-(v[2]-M0)/T1]
# v[0]: 磁化ベクトルのx軸成分
# v[1]: 磁化ベクトルのy軸成分
# v[2]: 磁化ベクトルのz軸成分

om = 9.5 * 10 ** 9 # 測定周波数/Hz
T1, T2 = 5 * 10 ** -8, 5 * 10 ** -8 # 縦緩和時間，横緩和時間/s
H0, H1 = 3.3 * 10 ** -3, 10 ** 0 # 外部磁場および振動磁場の振幅/T
# H0, H1 = 3.3 * 10 ** -3, 10 ** -5
v0 = [0, 0, 1] # 初期磁化 (Mz = 1)
t = np.arange(0, time, time/N)

v = odeint(func, v0, t, args=(om, H0, H1, T1, T2))

# 単位球上に磁化ベクトルの軌跡を表示する
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
# （１）単位球を描く
U = np.linspace(0, 2 * np.pi, 100)
V = np.linspace(0, np.pi, 100)
X = np.outer(np.cos(U), np.sin(V))
Y = np.outer(np.sin(U), np.sin(V))
Z = np.outer(np.ones(np.size(U)), np.cos(V))
ax.plot_surface(X, Y, Z, cmap=cm.coolwarm, alpha=0.5)
# （２）有効磁場（磁化ベクトルの回転軸）を描く
norm = np.sqrt(H1 ** 2 + (H0-om/ga) ** 2)
ax.quiver(0, 0, 0, H1/norm, 0, (H0-om/ga)/norm, length=1, color='k', linewidth=2)
ax.text(H1/norm, 0, (H0-om/ga)/norm, "$\it{B_{eff}}$")
# （３）磁化ベクトルの運動を描く
ax.plot(v[:, 0], v[:, 1], v[:, 2])
ax.scatter3D(0, 0, 1, color='b')
ax.scatter3D(v[-1,0], v[-1,1], v[-1,2], color='r')
ax.set_xlabel('x(t)')
ax.set_ylabel('y(t)')
ax.set_zlabel('z(t)')

# 磁化ベクトルのy軸（回転座標系）成分を時間の関数として表示する
plt.figure()
plt.plot(t, v[:,1])

# 上記成分をフーリエ変換する
plt.figure()
sp = np.fft.fft(v[:,1])
freq = np.fft.fftfreq(v[:,1].shape[-1])
plt.plot(freq,np.abs(sp))

plt.show()