slxrd.py ダウンロード/コピー

slxrd.py をダウンロード

slxrd.py
slxrd.py
  1"""
  2X線回折シミュレーションスクリプト
  3
  4概要:
  5    多層膜構造からのX線回折強度を計算し、CSVファイルに出力し、プロットします。
  6    ブラッグの法則と複素ラウエ関数を用いて、異なるユニットセルが積層された構造の回折パターンをシミュレートします。
  7
  8関連リンク:
  9    slxrd_usage
 10"""
 11
 12import csv
 13from pprint import pprint
 14from math import sin, cos, asin, acos
 15from cmath import sin as csin, cos as ccos, sqrt as csqrt, exp as cexp
 16import numpy as np
 17from matplotlib import pyplot as plt
 18
 19
 20#=============================
 21# 定数
 22#=============================
 23pi          = 3.14159265358979323846
 24pi2         = pi + pi
 25pi2i        = pi2 * 1j
 26torad       = 0.01745329251944 # rad/deg";
 27todeg       = 57.29577951472   # deg/rad";
 28
 29
 30#=============================
 31# 大域変数
 32#=============================
 33wl    = 1.54 # angstrom
 34outcsv = 'slxrd.csv'
 35
 36#=============================
 37# 構造
 38#=============================
 39# 単位格子: c軸長、構造因子、繰り返し数
 40cells = [
 41            {'a': 4.0, 'F': 10.0, 'N': 5},
 42            {'a': 4.1, 'F': 20.0, 'N': 5},
 43        ]
 44# 層構造: cellsの番号で指定
 45#layers = [ 0 ]
 46layers = [ 0, 1 ]
 47#layers = [ 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1 ]
 48#layers = [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ]
 49
 50# 計算する2θ範囲、データ数
 51range2Q = [5.0, 150.0]
 52n2Q = 5001
 53step2Q = (range2Q[1] - range2Q[0]) / (n2Q - 1)
 54
 55
 56def cLaue(Q2, N, alatt, wl):
 57    """
 58    概要:
 59        複素ラウエ関数を計算します。
 60    詳細説明:
 61        N個のユニットセルからなる単一周期構造の回折強度を記述するのに用いられるラウエ関数を複素数で計算します。
 62        この関数は、周期構造からの回折の重ね合わせ効果を表現します。
 63    引数:
 64        :param Q2: 回折角度2θ（度）。
 65        :type Q2: float
 66        :param N: ユニットセルの繰り返し数。
 67        :type N: int
 68        :param alatt: ユニットセルの格子定数（Å）。
 69        :type alatt: float
 70        :param wl: X線の波長（Å）。
 71        :type wl: float
 72    戻り値:
 73        :returns: 複素ラウエ関数の値。
 74        :rtype: complex
 75    """
 76    delta = 2.0 * alatt * csin(Q2 / 2.0 * torad) / wl
 77    return (1.0 - cexp(pi2i * delta * N)) / (1.0 - cexp(pi2i * delta))
 78
 79def diffQ(d, wl):
 80    """
 81    概要:
 82        ブラッグの法則に基づき回折角度θを計算します。
 83    詳細説明:
 84        格子面間隔dとX線波長wlから、ブラッグの法則 (2d sinθ = nλ) を用いて、回折角度θを計算します。
 85        sinθの値が物理的に許容範囲外の場合 (abs(sinθ) > 1.0) はNoneを返します。
 86    引数:
 87        :param d: 格子面間隔（Å）。
 88        :type d: float
 89        :param wl: X線の波長（Å）。
 90        :type wl: float
 91    戻り値:
 92        :returns: 回折角度θ（度）。回折が生じない場合はNone。
 93        :rtype: float or None
 94    """
 95    sinQ = wl / 2.0 / d
 96    if abs(sinQ) > 1.0:
 97        return None
 98    Q = todeg * asin(sinQ)
 99    return Q
100
101
102if __name__ == "__main__":
103# 構成単位格子の回折角度の計算
104    print("")
105    print("Diffraction angles from each unit cell")
106    for ic in range(len(cells)):
107        print("cell #{}: a={} F={} N={}"
108                .format(ic, cells[ic]['a'], cells[ic]['F'], cells[ic]['N']))
109        for l in range(1, 20):
110            d00l = cells[ic]['a'] / l
111            Q_00l = diffQ(d00l, wl)
112            if Q_00l is None:
113                break
114            Q2_00l = 2.0 * Q_00l
115            if range2Q[0] <= Q2_00l <= range2Q[1]:
116                sinQ = sin(torad * Q_00l)
117                LN = sin(pi * cells[ic]['N'] * cells[ic]['a'] * 2.0 * sinQ / wl)
118                L1 = sin(pi * 1 * cells[ic]['a'] * 2.0 * sinQ / wl)
119                LN2 = (LN * LN.conjugate()).real
120                L12 = (L1 * L1.conjugate()).real
121                I = LN2 / L12
122                print("  2Q(00%2d)=%10.6g" % (l, Q2_00l), 
123                      "  I = %8.4g = %8.4g / %8.4g" % (I, LN2, L12))
124
125# 繰り返し単位格子単位の回折角度の計算
126    print("")
127    print("Diffraction angles from each periodic unit cells unit")
128    for ic in range(len(cells)):
129        print("unit #{}: a={} F={} N={}"
130                .format(ic, cells[ic]['a'], cells[ic]['F'], cells[ic]['N']))
131        for n in range(1, 20):
132            sinQ = (n + 0.5) * wl / 2.0 / cells[ic]['N'] / cells[ic]['a']
133            if abs(sinQ) > 1.0:
134                break
135            Q2 = 2.0 * todeg * asin(sinQ)
136            LN = sin(pi * cells[ic]['N'] * cells[ic]['a'] * 2.0 * sinQ / wl)
137            L1 = sin(pi * 1 * cells[ic]['a'] * 2.0 * sinQ / wl)
138            LN2 = (LN * LN.conjugate()).real
139            L12 = (L1 * L1.conjugate()).real
140            I = LN2 / L12
141            if range2Q[0] <= Q2 <= range2Q[1]:
142                print("  2Q(%2d)=%10.6g" % (n, Q2), 
143                      "  I = %8.4g = %8.4g / %8.4g" % (I, LN2, L12))
144            
145
146# 回折強度の計算
147    x = []
148    y = []
149    for i in range(n2Q):
150        Q2 = range2Q[0] + step2Q * i
151    
152        t = 0.0
153        F = 0.0
154        for il in range(len(layers)):
155            cell = cells[layers[il]]
156            lf = cLaue(Q2, cell['N'], cell['a'], wl)
157            delta = 2.0 * t * sin(Q2 / 2.0 * torad) / wl
158            F += cells[ic]['F'] * lf * cexp(pi2i * delta)
159            t += cell['a'] * cell['N']
160
161        I = (F * F.conjugate()).real
162
163        x.append(Q2)
164        y.append(I)
165#    print(i, ": ", Q2, lf2)
166
167
168# CSVファイルへ書き出し
169    print("")
170    print("Write to [{}]".format(outcsv))
171    try:
172        f = open(outcsv, 'w')
173        fout = csv.writer(f, lineterminator='\n')
174        fout.writerow(('2Q', 'Intensity'))
175        for i in range(0, len(x)):
176            fout.writerow((x[i], y[i]))
177        f.close()
178    except IOError:
179        print("Error: Can not write to [{}]".format(outcsv))
180
181
182# グラフプロット
183    plt.plot(x, y, linestyle = 'solid')
184    plt.xlabel('2Q')
185    plt.ylabel('Intensity')
186
187    plt.pause(0.001)
188    print("Press ENTER to exit")
189    input()