"""
X線回折シミュレーションスクリプト
概要:
多層膜構造からのX線回折強度を計算し、CSVファイルに出力し、プロットします。
ブラッグの法則と複素ラウエ関数を用いて、異なるユニットセルが積層された構造の回折パターンをシミュレートします。
関連リンク:
:doc:`slxrd_usage`
"""
import csv
from pprint import pprint
from math import sin, cos, asin, acos
from cmath import sin as csin, cos as ccos, sqrt as csqrt, exp as cexp
import numpy as np
from matplotlib import pyplot as plt
#=============================
# 定数
#=============================
pi = 3.14159265358979323846
pi2 = pi + pi
pi2i = pi2 * 1j
torad = 0.01745329251944 # rad/deg";
todeg = 57.29577951472 # deg/rad";
#=============================
# 大域変数
#=============================
wl = 1.54 # angstrom
outcsv = 'slxrd.csv'
#=============================
# 構造
#=============================
# 単位格子: c軸長、構造因子、繰り返し数
cells = [
{'a': 4.0, 'F': 10.0, 'N': 5},
{'a': 4.1, 'F': 20.0, 'N': 5},
]
# 層構造: cellsの番号で指定
#layers = [ 0 ]
layers = [ 0, 1 ]
#layers = [ 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1 ]
#layers = [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ]
# 計算する2θ範囲、データ数
range2Q = [5.0, 150.0]
n2Q = 5001
step2Q = (range2Q[1] - range2Q[0]) / (n2Q - 1)
[ドキュメント]
def cLaue(Q2, N, alatt, wl):
"""
複素ラウエ関数を計算します。
詳細説明:
N個のユニットセルからなる単一周期構造の回折強度を記述するのに用いられるラウエ関数を複素数で計算します。
この関数は、周期構造からの回折の重ね合わせ効果を表現します。
:param Q2: 回折角度2θ（度）。
:type Q2: float
:param N: ユニットセルの繰り返し数。
:type N: int
:param alatt: ユニットセルの格子定数（Å）。
:type alatt: float
:param wl: X線の波長（Å）。
:type wl: float
:returns: 複素ラウエ関数の値。
:rtype: complex
"""
delta = 2.0 * alatt * csin(Q2 / 2.0 * torad) / wl
return (1.0 - cexp(pi2i * delta * N)) / (1.0 - cexp(pi2i * delta))
[ドキュメント]
def diffQ(d, wl):
"""
ブラッグの法則に基づき回折角度θを計算します。
詳細説明:
格子面間隔dとX線波長wlから、ブラッグの法則 (2d sinθ = nλ) を用いて、回折角度θを計算します。
sinθの値が物理的に許容範囲外の場合 (abs(sinθ) > 1.0) はNoneを返します。
:param d: 格子面間隔（Å）。
:type d: float
:param wl: X線の波長（Å）。
:type wl: float
:returns: 回折角度θ（度）。回折が生じない場合はNone。
:rtype: float or None
"""
sinQ = wl / 2.0 / d
if abs(sinQ) > 1.0:
return None
Q = todeg * asin(sinQ)
return Q
if __name__ == "__main__":
# 構成単位格子の回折角度の計算
print("")
print("Diffraction angles from each unit cell")
for ic in range(len(cells)):
print("cell #{}: a={} F={} N={}"
.format(ic, cells[ic]['a'], cells[ic]['F'], cells[ic]['N']))
for l in range(1, 20):
d00l = cells[ic]['a'] / l
Q_00l = diffQ(d00l, wl)
if Q_00l is None:
break
Q2_00l = 2.0 * Q_00l
if range2Q[0] <= Q2_00l <= range2Q[1]:
sinQ = sin(torad * Q_00l)
LN = sin(pi * cells[ic]['N'] * cells[ic]['a'] * 2.0 * sinQ / wl)
L1 = sin(pi * 1 * cells[ic]['a'] * 2.0 * sinQ / wl)
LN2 = (LN * LN.conjugate()).real
L12 = (L1 * L1.conjugate()).real
I = LN2 / L12
print(" 2Q(00%2d)=%10.6g" % (l, Q2_00l),
" I = %8.4g = %8.4g / %8.4g" % (I, LN2, L12))
# 繰り返し単位格子単位の回折角度の計算
print("")
print("Diffraction angles from each periodic unit cells unit")
for ic in range(len(cells)):
print("unit #{}: a={} F={} N={}"
.format(ic, cells[ic]['a'], cells[ic]['F'], cells[ic]['N']))
for n in range(1, 20):
sinQ = (n + 0.5) * wl / 2.0 / cells[ic]['N'] / cells[ic]['a']
if abs(sinQ) > 1.0:
break
Q2 = 2.0 * todeg * asin(sinQ)
LN = sin(pi * cells[ic]['N'] * cells[ic]['a'] * 2.0 * sinQ / wl)
L1 = sin(pi * 1 * cells[ic]['a'] * 2.0 * sinQ / wl)
LN2 = (LN * LN.conjugate()).real
L12 = (L1 * L1.conjugate()).real
I = LN2 / L12
if range2Q[0] <= Q2 <= range2Q[1]:
print(" 2Q(%2d)=%10.6g" % (n, Q2),
" I = %8.4g = %8.4g / %8.4g" % (I, LN2, L12))
# 回折強度の計算
x = []
y = []
for i in range(n2Q):
Q2 = range2Q[0] + step2Q * i
t = 0.0
F = 0.0
for il in range(len(layers)):
cell = cells[layers[il]]
lf = cLaue(Q2, cell['N'], cell['a'], wl)
delta = 2.0 * t * sin(Q2 / 2.0 * torad) / wl
F += cells[ic]['F'] * lf * cexp(pi2i * delta)
t += cell['a'] * cell['N']
I = (F * F.conjugate()).real
x.append(Q2)
y.append(I)
# print(i, ": ", Q2, lf2)
# CSVファイルへ書き出し
print("")
print("Write to [{}]".format(outcsv))
try:
f = open(outcsv, 'w')
fout = csv.writer(f, lineterminator='\n')
fout.writerow(('2Q', 'Intensity'))
for i in range(0, len(x)):
fout.writerow((x[i], y[i]))
f.close()
except IOError:
print("Error: Can not write to [{}]".format(outfile))
# グラフプロット
plt.plot(x, y, linestyle = 'solid')
plt.xlabel('2Q')
plt.ylabel('Intensity')
plt.pause(0.001)
print("Press ENTER to exit")
input()