PythonでBloch方程式を可視化する⑤〜周波数掃引〜

プログラム Python

前回までで実装したフーリエ変換を周波数を変えて行い，実際の測定に近づけていく．

なお，まだ窓関数は使っていない．

 

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ループでスペクトルの行列生成，および各周波数に対するスペクトルの総和を取っているが，もっと簡単な方法があるような気がする．