##### プログラムはほとんど example-03-1.py と同じ #####
#### ここからしばらく一緒です ####
import numpy as np
import matplotlib.pyplot as plt
import xrayutilities as xu

# GaN, AlN の結晶構造をCIFから読む
GaN = xu.materials.Crystal.fromCIF("GaN.cif")
AlN = xu.materials.Crystal.fromCIF("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

#### ここまで一緒です。 ####

# さらに混晶も準備する
Alx = 0.2                                                # Al組成 x = 0.2
AlGaN = xu.materials.material.Alloy(GaN, AlN, Alx);      # AlN 0.2 + GaN (1-0.2)

## まずは 基板(厚さ無限のGaN Layer)、AlGaN層を準備
sub    = xu.simpack.Layer( GaN, float('inf') )            # 基板 (GaN)
lAlGaN = xu.simpack.Layer( AlGaN, 1500, relaxation=1.0 )  # 膜 (Al0.2GaN)
                                                          # relaxation=1.0 => 0.0 にすると回折ピークは広角側に動く
#
## 準備した Layer を積み上げて AlGaN層 / GaN 基板という積層構造を構築
Epi_Sub = xu.simpack.PseudomorphicStack001( "AlGaN/GaN", sub+lAlGaN )

#### この先はまた一緒になります。(コメントや画面表示用の文字列とかは別にして)

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

# 動的モデルの設定とシミュレーション
## リファレンスに基板(GaN)だけの計算もする
Sub_Model         = xu.simpack.DynamicalModel( sub,     energy=energy, resolution_width=0.0001 )
Int_Sub_Model     = Sub_Model.simulate( twotheta/2, hkl=(0, 0, 2) )
## こちらは Epi/Sub = AlGaN/GaN の計算
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.plot( twotheta/2, Int_Sub_Model )
plt.title( 'Dynamical Simulation of AlGaN/GaN and GaN Sub')
plt.yscale('log')
plt.xlabel(r'$\omega\ (deg)$')
plt.ylabel('Intensity (arb. u.)')
plt.show()