とあるスクールで、講師が事前に用意したOriginのマクロを使ってXMCDのSum ruleから磁気モーメントを求めるという講習があったが、Originを使いたいと全く思わなかったため、一人でpythonスクリプトを黙々と作っていた。
ローレンツ関数やら誤差関数やらで適当に作った架空のXASとXMCD。
これを処理してSum ruleを適用する。
手順としては
- pre edgeを差し引く。
- 引いたものをpost edgeで割る。
- XMCDを(再度)作る。
- ピークの位置を求める。
- step関数を作る。
- step関数をXASから引いて積分する。
- との境目を決めてXMCDを積分する。
という感じだろう。
pre edge、post edgeを以下の関数を使って、求める。
import numpy as np def pre_edge( x, y, start, end ): pre_x = [] pre_y = [] for ix, xx in enumerate( x ): if start < xx: pre_x.append( xx ) pre_y.append( y[ ix ] ) if end < xx: break pre_x = np.array( pre_x ) pre_y = np.array( pre_y ) A = np.c_[ pre_x, np.ones( len( pre_x ) ) ] a, b = np.linalg.lstsq( A, pre_y, rcond=None )[0] return a * x + b def post_edge( x, y, start, end ): post_x = [] post_y = [] for ix, xx in enumerate( x ): if start < xx: post_x.append( xx ) post_y.append( y[ ix ] ) if end < xx: break post_x = np.array( post_x ) post_y = np.array( post_y ) A = np.c_[ post_x, np.ones( len( post_x ) ) ] a, b = np.linalg.lstsq( A, post_y, rcond=None )[0] return a * x + b
要は、各エネルギー領域を指定してあげて、その中のデータ点を用いて最小二乗法で直線を定めるというもの。
今の場合だと、
- pre edge: 680 ~ 700 eV
- pre edge: 760 ~ 800 eV
くらいで設定した。
とりあえずプラスhelicityのXASにのみ適用してみる。
# energy: array for energy points # spec_p: array for XAS spectrum with plus helicity pre = pre_edge( energy, spec_p, 680, 700 ) norm0 = spec_p - pre post = post_edge( energy, norm0, 760, 800 ) norm_p = norm0 / post fig = plt.figure( figsize=( 8, 2.5 ) ) ax1 = fig.add_subplot( 1, 3, 1 ) ax2 = fig.add_subplot( 1, 3, 2 ) ax3 = fig.add_subplot( 1, 3, 3 ) ax1.plot( energy, spec_p, "r", label="spec_p" ) ax1.plot( energy, pre, "green" ) ax1.set_xlabel( "Energy [eV]") ax2.plot( energy, norm0, "r", label="norm0" ) ax2.plot( energy, post, "purple" ) ax2.set_xlabel( "Energy [eV]") ax3.plot( energy, norm_p, "r", label="norm_p" ) ax3.set_xlabel( "Energy [eV]") plt.tight_layout() plt.show()
うまく引けている感じがする。
規格化後のXASとXMCD。
xas = ( norm_p + norm_m ) / 2.
xmcd = norm_p - norm_m
ピークの位置を、次の関数で求める。
最初のピークは710 ~ 730 eV、次のピークは730 ~ 750 eVの範囲にあるので、そのように設定。
def white_line_position( x, y, x0, x1 ): p = max([ yy for xx, yy in zip( x, y ) if x0 < xx < x1 ]) for ix, xx in enumerate( x ): if x0 < xx < x1: if abs( y[ ix ] - p ) < 1e-7: break return ix, xx _, e1 = white_line_position( energy, xas, 710, 730 ) _, e2 = white_line_position( energy, xas, 730, 750 ) print( e1, e2 ) # 720.36..., 740.54...
d状態以外の寄与を差し引くためのstep関数を次のように用意。
from scipy.special import erf def step_func( x, x0, d ): return ( erf( ( x - x0 ) / d ) + 1. ) / 2. def subtract_func( x, x01, x02, d, c ): return c * step_func( x, x01, d ) + ( c / 2. ) * step_func( x, x02, d)
元々の誤差関数はの値を持つため、step_funcでは]になるように調整した。
強度比は軌道モーメントがなければであるため、大きさを決める係数は一つで良い。(この関数は、(遷移金属原子における)d状態以外の寄与のためのモデルであり、それらは一般に軌道モーメントは小さい。)
求めたピークの位置を使って、step関数を作る。
step関数のfitting領域は760 eV以降のpost edgeに設定。
def index_energy( energy, e0 ): for i, e in enumerate( energy ): if e >= e0: if abs( energy[ i ] - e0 ) <= abs( energy[ i - 1 ] - e0 ): return i else: return i - 1 i_post = index_energy( energy, 760 ) from scipy.optimize import curve_fit c, _ = curve_fit( lambda e, c: subtract_func( e, e1, e2, de, c ), energy[ i_post : ], xas[ i_post : ] ) subtract = subtract_func( energy, e1, e2, de, c ) plt.plot( energy, subtract ) plt.plot( energy, xas ) plt.xlabel( "Energy [eV]" ) plt.ylabel( "Intensity [arb.unit]" ) plt.show()
それらしく階段関数が置けている。
curve_fitの使い方は以下のサイトを参照。
pythonでフィッティングをする - おっぱいそん!
階段関数を差し引いたXAS、そしてXMCDを次のように積分する。
ここではSimpson法を用いる。
from scipy.integrate import simps def integrate( x, y ): result = [ 0. ] for ix in range( len( x ) )[ 1: ]: result.append( simps( y[ :ix ], x[ :ix ] ) ) return np.array( result ) int_xas = integrate( energy, xas - subtract ) int_xmcd = integrate( energy, xmcd ) plt.plot( energy, int_xas ) plt.plot( energy, int_xmcd ) plt.xlabel( "Energy [eV]" ) plt.show() print( int_xas[-1] ) # 18.914...
ちょっと後ろの方でフニャッとしているが、こんなもんだろう。
との境目は、XASの一回微分がゼロになるところを持って来た。
def get_index( xs, x0 ): for i, x in enumerate( xs ): if x >= x0: if abs( xs[ i ] - x0 ) <= abs( xs[ i - 1 ] - x0 ): return i else: return i - 1 def index_zero( x, y, x_start, x_end ): iy_start = get_index( x, x_start ) iy_end = get_index( x, x_end ) y_zero = min( abs( y[ iy_start : iy_end ] ) ) for iy, yy in enumerate( y ): if not iy_start <= iy <= iy_end: continue if abs( yy ) - y_zero <= 1e-7: return iy i_zero = index_zero( energy[ 1 : ], np.diff( xas, n=1 ), 730, 735 ) plt.plot( energy, xas ) plt.plot( energy[ im ], xas[ im ], "o" ) plt.show() print( int_xmcd[ i_zero ] ) # -4.5744...
差分による近似的に微分した配列が得られるnp.diffの使い方は、以下のサイトを参照。
要素の差分、足し合わせを計算するNumPyのdiff関数とcumsum関数の使い方 - DeepAge
ひとまずはこれでSum ruleに必要なXAS及びXMCDの積分値を求めることが出来た。
ちなみに、適当に作ったスペクトルは以下のスクリプトを使用。
def lorentzian_norm( x, x0, gamma ): return gamma**2 / ( ( x - x0 )**2 + gamma**2 ) energy = np.linspace( 680, 800, 1000 ) de = energy[1] - energy[0] peak1 = 720 peak2 = 740 back = ( 0.01 * ( energy - 900 ) )**2 base = subtract_func( energy, peak1, peak2, 10 * de, 1 ) peak_p = 2 * lorentzian_norm( energy, peak1, 2 ) peak_p += 1 * lorentzian_norm( energy, peak2, 2 ) spec_p = base + peak_p + back peak_m = 3.0 * lorentzian_norm( energy, peak1, 2 ) peak_m += 0.5 * lorentzian_norm( energy, peak2, 2 ) spec_m = base + peak_m + back
ローレンツ関数の大きさ等は以下を参照。
ローレンツ関数とは?