PythonでBloch方程式を可視化する③〜回転座標系〜

プログラム Python

前回までは空間に固定された座標系を使っていた．

今回は「回転座標系」と呼ばれる，1つの軸（z軸とする）周りにマイクロ波と同じ周波数で回転する座標系でのBloch方程式を扱う．

下記コード中 v[:,1] の時間発展を追い，フーリエ変換を行った．

なお，連続波実験の条件に合わせてH1（電磁波の振幅）の値を小さくした．このとき磁化ベクトルはほとんど動かないため，Bloch球上での可視化を削除した．（追記参照）

 

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

N = 1024
time = 10 ** -8
M0 = 1
def func(v, t, ga, 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]
ga, om, T1, T2 = 10 ** 11, 9 * 10 ** 9, 5 * 10 ** -9, 10 ** -9
H0, H1 = 3 * 10 ** -2, 10 ** -3
v0 = [0, 0, 1]
t = np.arange(0, time, time/N)

v = odeint(func, v0, t, args=(ga, 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)

ax = fig.gca(projection='3d')
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)')
ax.view_init(75,180)

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.plot(freq,sp.real,freq,sp.imag)
plt.show()

 

前々回の最後に，座標変換による挙動の異変について触れたが，その原因はおそらく飽和の近似が破れていたためであると考えられる．

書いている途中で，窓関数のことを思い出したが，今後の課題としたい．

 

2019-03-27追記

パラメータを調整することで，磁化ベクトルの特異的な運動が明確になることが分かったので，可視化部分を復活させた．