親知らずを抜いた話
※注意※
タイトルから想像される通り，痛い話です．
火曜日の午後．
親知らずの1本が横に生えてきてしかも虫歯だということで，大学病院で抜いてもらった．
初めに麻酔を打った．
何本か打った．
歯茎から唇までが奇妙な感触になった．
麻酔が効くまで少し待ってから手術に入った．
歯茎を切って歯を取り出すと言ってしまえば簡単に聞こえるが，実際は歯を割ってから取り出さなければならない．
そのため，削る音，ヒビを入れる音が交互に聞こえた．
目は開けていたが口が見えるわけないので，想像するしかなかった．
麻酔は効いていたはずなのに，何かを歯に押し込まれるときが一番痛かった．
我慢できると伝えたため，追加で打たれることはなかった．
歯がバキュームに吸い付く音がした．
4針縫って終わった．
正確には覚えていないが，作業自体は30分くらいであったと思う．
麻酔が効いていたのと，傷を刺激してはいけないと思ったので，その日は口が半開きだった．
アホっぽい顔をしていたと思う．
これまでも何度か歯を抜いたが，今回が一番痛かった．
今もまだ痛む．
やはり歯は大切にせねば．
Pythonで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 ** -9
ga = -1.7608 * 10 ** 11
M0 = 1
range_om = 100 # 掃引の点数
d_om = 0.01 # 掃引の間隔
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]
T1, T2 = 5 * 10 ** -8, 5 * 10 ** -8
H0, H1 = 3.3 * 10 ** -3, 10 ** -5
v0 = [0, 0, 1]
t = np.arange(0, time, time/N)
spc = [0] * range_om
plt.figure()
for i in range(range_om):
om = (9.5 + d_om * (i - np.floor(range_om/2))) * 10 ** 9
v = odeint(func, v0, t, args=(om, H0, H1, T1, T2))
sp = np.fft.fft(v[:,1])
freq = np.fft.fftfreq(v[:,1].shape[-1])
spc[i] = np.abs(sp)
# plt.plot(freq,np.abs(sp))
for i2 in range(N):
spc_tot = np.sum(spc,axis=0)
plt.plot(freq,spc_tot)
plt.show()
2回のforループでスペクトルの行列生成，および各周波数に対するスペクトルの総和を取っているが，もっと簡単な方法があるような気がする．
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()
以下にH1/H0 = 300および0.003のときの磁化ベクトルの軌跡を示す．
上の図は青点から始まり，赤点に至るまでベクトルがBeffの矢印を軸に回転していることがわかりやすい．
下の図では磁化ベクトルは動いていないように見えるが，拡大するとその回転運動を確認できる．
上はパルス測定，下はCW（連続波）測定に相当する．
今後はCW測定の場合について，周波数（om）の変化によるフーリエ変換したスペクトルの変化を追うことを考える．
PythonでBloch方程式を可視化する③〜回転座標系〜
前回までは空間に固定された座標系を使っていた．
今回は「回転座標系」と呼ばれる，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追記
パラメータを調整することで，磁化ベクトルの特異的な運動が明確になることが分かったので，可視化部分を復活させた．
PythonでBloch方程式を可視化する②〜odeint編〜
前回（https://kousukejapan.hatenablog.jp/entry/2019/03/16/222053）はEuler法でBloch方程式を解いたが，odeintを使うことを学んだのでそっちでもやってみた．
import numpy as np
# from math import sin, cos
from scipy.integrate import odeint
from matplotlib import cm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
M0 = 1
def func(v, t, ga, om, H0, H1, T1, T2):
return [ga * (v[1] * H0 - v[2] * H1 * np.sin(om * t)) - v[0] / T2,
ga * (v[2] * H1 *np.cos(om * t) - v[0] * H0) - v[1] / T2,
ga * H1 * (v[0] *np.sin(om * t) - v[1] *np.cos(om * t)) - (v[2] - M0) / T1]
ga, om, T1, T2 = 10 ** 11, 10 ** 9, 5 * 10 ** -10, 10 ** -10
H0, H1 = 10 ** -2, 5 * 10 ** -1
v0 = [0, 0, 1]
t = np.arange(0, 10 ** -8, 10 ** -13)
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='b')
ax.set_xlabel('x(t)')
ax.set_ylabel('y(t)')
ax.set_zlabel('z(t)')
forループがなくなった分スッキリした気がする．
精度はどうなんだろう？
2019-03-25追記
numpyに三角関数が入っているから，mathモジュールは要らなかった．
PythonでBloch方程式を可視化する
磁化ベクトルの微分方程式であるBloch方程式をPythonで解いてみた．
微分方程式の数値解は以下のpdfを参考に，Euler法で実装した．
http://www.isc.meiji.ac.jp/~mizutani/python/print/diffeq/differential_eq.pdf
import numpy as np
# from math import sin, cos
from matplotlib import cm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
ga, om, T1, T2 = 10 ** 11, 10 ** 9, 5 * 10 ** -10, 10 ** -10
H0, H1 = 10 ** -2, 5 * 10 ** -1
M0 = 1
def dx(t, x, y, z):
return(ga * (y*H0 - z*H1*np.sin(om*t)) - x/T2)
def dy(t, x, y, z):
return(ga * (z*H1*np.cos(om*t) - x*H0) - y/T2)
def dz(t, x, y, z):
return(ga * H1 * (x*np.sin(om*t) - y*np.cos(om*t)) - (z - M0)/T1)
def euler_orbit(x0, T, dt, vectorfield):
width = len(x0)
x = x0
t = 0
orbit = [] orbit.append(list(x))
while t <= T:
vx = list(map(lambda v: v(t, *x), vectorfield))
for i in range(width):
x[i] += dt * vx[i]
orbit.append(list(x)) t += dt
return(orbit)
vector = [dx, dy, dz]
x0 = [0, 0, 1]
dt = 10 ** -13
T = 10 ** -7
orbit = euler_orbit(x0, T, dt, vector)
nporbit = np.array(orbit)
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(nporbit[:,0], nporbit[:,1], nporbit[:,2])
ax.scatter3D(0, 0, 1,color='b')
ax.scatter3D(nporbit[-1,0], nporbit[-1,1], nporbit[-1,2],color='b')
ax.set_xlabel('x(t)')
ax.set_ylabel('y(t)')
ax.set_zlabel('z(t)')
plt.show()
自分なりの工夫としては，単位球を描かせて解ベクトルがその中を動いているかを確認できるようにしたことと，初期値（[0,0,1]）と終値を青点で表示したこと．前者については，MathematicaのSphereみたいに一発で書くことはできないのかな．
これは実験室座標系（固定されている座標系）で，本来は回転座標系（電磁波の周波数で回転する座標系）に移らないといけない．方程式をいじるだけなのだが，解を見るとどうも上手くいっていない．さて，どうしたものか．
2019-03-25追記
numpyに三角関数が含まれているから，mathモジュールは必要なかった．
