matplotlib for solid state physics

固体物理学にあらわれるグラフを、 matplotlib で描く

グラフの形を見たいだけなので、定数や係数は適当です。

初歩から学ぶ固体物理学 矢口裕之 (5.15)

図5.4 さまざまな温度におけるフェルミ分布関数

本で、T=100K から 1000K とあるのは、T=10K から 100K の誤りと思います。

図5.5 さまざまな温度におけるボース分布関数

# coding: utf-8

# 初歩から学ぶ固体物理学 矢口裕之 (5.15)

# 図5.4 さまざまな温度におけるフェルミ分布関数

import numpy as np

import matplotlib.pyplot as plt

E=np.arange(0, 600)

for T in 1, 10, 50, 100:

    f = 1 / (np.exp((E - 200)/T) + 1)

    plt.plot(E, f)

plt.show()

# 図5.5 さまざまな温度におけるボース分布関数

E=np.arange(0.001, 10, 0.1)

plt.ylim(0, 1000)

for T in 1, 100, 500, 1000:

    f = 1 / (np.exp((E - 0)/T) - 1)

    plt.plot(E, f)

plt.show()

初歩から学ぶ固体物理学 矢口裕之 (7.10)

図7.7 ２種類の原子からなる１次元の格子振動における波数と角振動数の関係

# coding: utf-8

# 初歩から学ぶ固体物理学 矢口裕之 (7.10)

# 図7.7 ２種類の原子からなる１次元の格子振動における波数と角振動数の関係

import numpy as np

import matplotlib.pyplot as plt

# a=K=1

Ma = 2

Mb = 1

k=np.arange(-np.pi, np.pi, 0.1)

delta = np.sqrt((1/Ma + 1/Mb)**2 - (4* np.sin(k/2)**2)/(Ma * Mb))

omega0 = np.sqrt((1/Ma + 1/Mb) + delta)

omega1 = np.sqrt((1/Ma + 1/Mb) - delta)

plt.plot(k, omega0)

plt.plot(k, omega1)

plt.show()

初歩から学ぶ固体物理学 矢口裕之 (10.26)

図 10.5 ほとんど自由な電子モデルから求められた電子のエネルギー

# coding: utf-8

# 初歩から学ぶ固体物理学 矢口裕之 (10.26)

# 図 10.5 ほとんど自由な電子モデルから求められた電子のエネルギー

import numpy as np

import matplotlib.pyplot as plt

# 逆格子ベクトル Gm

G=np.int(10)

# 周期ポテンシャルの k=Gm のフーリエ係数

V=np.int(4)

# 波数ベクトル

#k=np.arange(-2*G, G, 0.1)

k=np.arange(-1.2* G, 0.2* G, 0.1)

# h, me, V0 は無視しています。

delta=np.sqrt((k**2 - (k+G)**2)**2 + 4*V**2)

E0=k**2 + (k+G)**2 + delta

E1=k**2 + (k+G)**2 - delta

plt.ylim(0, 150)

plt.plot(k, E0)

plt.plot(k, E1)

plt.show()

トポロジカル絶縁体・超伝導体 野村健太郎 (3.41)

図 3.3 (a) ２バンド模型 (3.39) のエネルギー分散

以下を実行すると、ぐりぐり動くよ。

# coding: utf-8

# トポロジカル絶縁体・超伝導体 野村健太郎 (3.41)

# 図 3.3 (a) ２バンド模型 (3.39) のエネルギー分散

import matplotlib

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

m=np.int(-1)

kx, ky = np.mgrid[-np.pi:np.pi:0.1, -np.pi:np.pi:0.1]

# t=r=a=1

R=np.sqrt(np.sin(kx)** 2 + np.sin(ky)**2 +

          (m + (1-np.cos(kx)) + (1-np.cos(ky)))**2)

# ε(k)

e = 10

fig = plt.figure()

ax = Axes3D(fig)

ax.plot_surface(kx, ky, e + R)

ax.plot_surface(kx, ky, e - R)

plt.show()

以上