# いつものおまじない。numpy, matplotlib, xrayutilities を使用可能にする
import numpy as np
import matplotlib.pyplot as plt
import xrayutilities as xu

#### ここからが本体 ####

# GaN, AlN の結晶構造をCIFから読む
GaN = xu.materials.Crystal.fromCIF("GaN.cif")  # 所持している GaN の cif ファイルの名前と場所 "場所/名前.cif"
AlN = xu.materials.Crystal.fromCIF("AlN.cif")  # 無ければここにのせている GaN, AlN の cif をダウンロードしてください

# 今作った結晶(GaN, AlN)は、結晶構造の情報しか持ってないので
# 積層構造にした時の歪の計算ができるように、弾性に関わる情報を追加する
#
## どちらも Hexagonal 
GaN.symm_class_name = 'Hexagonal'
AlN.symm_class_name = 'Hexagonal'
#
## GaN, AlN の弾性定数 Cij の値を、名前のついた値として準備
C11_GaN, C12_GaN, C13_GaN, C33_GaN, C44_GaN = 390, 145, 106, 398, 105
C11_AlN, C12_AlN, C13_AlN, C33_AlN, C44_AlN = 396, 137, 108, 373, 116
#
## xu にサポートしてもらって 2階テンソルの形に組み上げる
cij_array_GaN = xu.materials.material.HexagonalElasticTensor( C11_GaN, C12_GaN, C13_GaN, C33_GaN, C44_GaN )
cij_array_AlN = xu.materials.material.HexagonalElasticTensor( C11_AlN, C12_AlN, C13_AlN, C33_AlN, C44_AlN )
#
## さらに4階テンソルの形にもしておく ij to ijkl (プログラムでは '2' を 'to' の意味でよく使う)
cijkl_tensor_GaN = xu.materials.material.Cij2Cijkl(cij_array_GaN)
cijkl_tensor_AlN = xu.materials.material.Cij2Cijkl(cij_array_AlN)
#
## Crystal オブジェクト(GaN, AlN)に4階テンソルを属性として設定
setattr(GaN, 'cijkl', cijkl_tensor_GaN)
setattr(AlN, 'cijkl', cijkl_tensor_AlN)
#
## 2階テンソルも elastic 属性に代入
GaN.elastic = cij_array_GaN
AlN.elastic = cij_array_AlN

# ここまでで、弾性定数の情報も持った GaN, AlN という結晶が準備できた。
# 次は、この物質のを積み上げていくためにその部品として Layer を作り、
# それを積み上げて積層構造にする
#
## まずは 基板(厚さ無限のGaN Layer)、AlN層を準備
sub = xu.simpack.Layer( GaN, float('inf') )            # 基板 (GaN)
lAlN = xu.simpack.Layer( AlN, 1500, relaxation=1.0 )   # 膜 (AlN)
#
## 準備した Layer を積み上げて AlN層 / GaN 基板という積層構造を構築
Epi_Sub = xu.simpack.PseudomorphicStack001( "AlN/GaN", sub+lAlN )

# ω-2θ 強度計算
## 波長、エネルギー、計算する角度範囲指定
wavelength = xu.wavelength('CuKa1')           # 波長は Cu Ka1 の波長(1.5405...A)
energy     = 12390 / wavelength               # エネルギーは波長から計算
twotheta   = np.linspace(32,38,200)           # 横軸になる角度範囲のベクトルを作っておく

# 動的モデルの設定とシミュレーション
Epi_Sub_Model     = xu.simpack.DynamicalModel( Epi_Sub, energy=energy, resolution_width=0.0001 )
Int_Epi_Sub_Model = Epi_Sub_Model.simulate( twotheta/2, hkl=(0, 0, 2) )

# 結果のプロット
plt.figure()
plt.plot( twotheta/2, Int_Epi_Sub_Model )
plt.title( 'Dynamical Simulation of AlN(1500nm)/GaN Sub')
plt.yscale('log')
plt.xlabel(r'$\omega\ (deg)$')
plt.ylabel('Intensity (arb. u.)')
plt.show()