GSL Shell 学習:結晶格子の座標変換

2017-11-03 :  理科部 部活
Python、C言語、Perl、グルコサミン、Firefox
7~8年前に、
結晶格子の座標変換」と云うことで、何度か記事を書いていた。

  ・・・・・
  ・・・・・
  結晶格子の単位格子(unit cell)は、
  良く利用(?)している CIFファイルなどに記述されている。

  3軸方向の長さ、a, b, c と、
  各軸間のなす角度、alpha, beta, gamma で表されている。

  そこで、結晶の各軸を表す単位ベクトルを、
  3次元(xyz)空間のベクトルとして表したい!

そこで、

条件として、(最初は)
  ベクトルc をz軸に重ね、
  ベクトルb をyz平面(x座標の値がゼロ)とする。
  これで、
  ベクトル a, b, c は { x, y, z } で表されるはず。

次に、(“結晶格子の座標変換2”)
  ベクトルa をx軸に重ね、
  ベクトルb をxy平面(z座標の値がゼロ)とする
  ことの方が、一般的(?)のような???

更に、(“結晶格子の座標変換3”)
  ま、なにはともあれ、
  「結晶格子の座標変換」も「結晶格子の座標変換2」も、
  意味があったのですね。

そして、(“結晶格子の座標変換4”)
  つまり、「Materials Studio」を使って、Import/Export することで、
  格子定数(.car)から、直交座標系ベクトル(.msi)に変換できる。
  ってことですか。
  と云うことは、
  <紙>のいい加減な Perl スクリプト は無駄だった?


ですが、
今回、(高速インタプリタ言語 GSL Shell に依る)コードを作りたくなった。

仕様としては:
-- INPUT 1
a, b, c = 5.4, 5.4, 8.8
alpha, beta, gamma = 90, 90, 120

-- INPUT 2
va = { 4.67654, -2.70000, 0.00000 }
vb = { 0.00000, 5.40000, 0.00000 }
vc = { 0.00000, 0.00000, 8.80000 }
の様な内容のファイルを読み込んで、
-- INPUT 1
a, b, c = 5.4, 5.4, 8.8
alpha, beta, gamma = 90, 90, 120

-- A=X, B in XY
va = { 5.40000, 0.00000, 0.00000 }
vb = { -2.70000, 4.67654, 0.00000 }
vc = { 0.00000, 0.00000, 8.80000 }

-- C=Z, B in YZ
va = { 4.67654, -2.70000, 0.00000 }
vb = { 0.00000, 5.40000, 0.00000 }
vc = { 0.00000, 0.00000, 8.80000 }

--=================================
-- INPUT 2
va = { 4.67654, -2.70000, 0.00000 }
vb = { 0.00000, 5.40000, 0.00000 }
vc = { 0.00000, 0.00000, 8.80000 }

-- output
a, b, c = 5.4, 5.4, 8.8
alpha, beta, gamma = 90, 90, 120
の様な内容のファイルを出力する。


で、以下の様な GSL Shell プログラムを作った。
----- Conv_Cell_Vec.gsl -----
inF = 'INPUT.txt'; otF = 'OUTPUT.txt' -- 入出力ファイル
--=================================

require'pl'; SF=string.format
sqrt=math.sqrt; sin=math.sin; cos=math.cos; acos=math.acos
degrad = |x| x*math.pi/180; raddeg = |x| x*180/math.pi
--============ ファイル入力 =======
w = utils.readlines( inF ); nn=#w -- 全体を行毎の文字列配列で読込
for n=1,nn do assert( load( w[n] ) )() end

ka, kb, kc = degrad(alpha), degrad(beta), degrad(gamma)
CZvc = { 0, 0, c }
CZvb = { 0, b*sin(ka), b*cos(ka) }
az = a*cos(kb); ay = a*(cos(kc)-cos(ka)*cos(kb)) / sin(ka)
CZva = { sqrt((a*sin(kb))^2 - ay^2), ay, az }

AXva = { a, 0, 0 }
AXvb = { b*cos(kc), b*sin(kc), 0 }
cx = c*cos(kb); cy = c*(cos(ka)-cos(kc)*cos(kb)) / sin(kc)
AXvc = { cx, cy, sqrt((c*sin(kb))^2) - cy^2 }

Ax, Ay, Az = va[1], va[2], va[3]; A = sqrt(Ax^2+Ay^2+Az^2)
Bx, By, Bz = vb[1], vb[2], vb[3]; B = sqrt(Bx^2+By^2+Bz^2)
Cx, Cy, Cz = vc[1], vc[2], vc[3]; C = sqrt(Cx^2+Cy^2+Cz^2)

kA = raddeg(acos((Bx*Cx+By*Cy+Bz*Cz)/(B*C)))
kB = raddeg(acos((Cx*Ax+Cy*Ay+Cz*Az)/(C*A)))
kC = raddeg(acos((Ax*Bx+Ay*By+Az*Bz)/(A*B)))

OT=io.open( otF, 'wt' )
for i=1,3 do utils.fprintf( OT, SF( '%s\n', w[i] ) ) end
utils.fprintf( OT, SF( '\n-- A=X, B in XY\n' ) )
utils.fprintf( OT, SF( 'va = { %11.7f, %11.7f, %11.7f }\n', AXva[1], AXva[2], AXva[3] ) )
utils.fprintf( OT, SF( 'vb = { %11.7f, %11.7f, %11.7f }\n', AXvb[1], AXvb[2], AXvb[3] ) )
utils.fprintf( OT, SF( 'vc = { %11.7f, %11.7f, %11.7f }\n', AXvc[1], AXvc[2], AXvc[3] ) )
utils.fprintf( OT, SF( '\n-- C=Z, B in YZ\n' ) )
utils.fprintf( OT, SF( 'va = { %11.7f, %11.7f, %11.7f }\n', CZva[1], CZva[2], CZva[3] ) )
utils.fprintf( OT, SF( 'vb = { %11.7f, %11.7f, %11.7f }\n', CZvb[1], CZvb[2], CZvb[3] ) )
utils.fprintf( OT, SF( 'vc = { %11.7f, %11.7f, %11.7f }\n', CZvc[1], CZvc[2], CZvc[3] ) )
utils.fprintf( OT, SF( '\n--===================================\n' ) )
for i=5,8 do utils.fprintf( OT, SF( '%s\n', w[i] ) ) end
utils.fprintf( OT, SF( '\n-- output\n' ) )
utils.fprintf( OT, SF( 'a, b, c = %g, %g, %g\n', A, B, C ) )
utils.fprintf( OT, SF( 'alpha, beta, gamma = %g, %g, %g\n', kA, kB, kC ) )
OT:close()
でしょうか?



本日はここまで。


Lua を使った、理科部 部活は続く?


見ていただいた序でとは厚かましい限りですが、
お帰りに投票して頂けるとなお嬉しいです。 ⇒ blogram投票ボタン


170915
関連記事
スポンサーサイト

コメントの投稿

管理者にだけ表示を許可する

おきてがみ/blogram
blogram投票ボタン



おきてがみ

最新記事
カレンダー
10 | 2017/11 | 12
- - - 1 2 3 4
5 6 7 8 9 10 11
12 13 14 15 16 17 18
19 20 21 22 23 24 25
26 27 28 29 30 - -
月別アーカイブ
カテゴリ
最新コメント
検索フォーム
リンク
プロフィール

<紙>

Author:<紙>
ようこそ。
「パソコンヲタクの雑記帳」
もろもろなことを綴っています。
パソコン ヲタクってねくら?
画像は kami でなく kani です。

カウンター(fc2、i2i) /Google Analytics


i2i(from 2010-08-24)
Total =
Today  =  
Yesterday=
アンチエイジング

Google Analytics
ブックマーク