3元1次連立方程式を解く(Perl版)

2017-11-17 :  PCクリニック
Python、C言語、Perl、グルコサミン、Firefox
前(2017-11-15)の記事「3元1次連立方程式を解く」では、
3元1次連立方程式を解くGSL Shellのコードを作った。

これに従って、
(やりたかった)Perl版コードを作った。
まずは、単純な移植
##### Solve_by_Inv-Mat.pl
## 行列 A
( $a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33 ) =
( 0, 3, 6, 1, 4, 7, 2, 5, 9 );

$detA = $a11*$a22*$a33 + $a21*$a32*$a13 + $a31*$a12*$a23 -
$a11*$a32*$a23 - $a31*$a22*$a13 - $a21*$a12*$a33;

## 逆行列 Inv_A == G を、2段階で
( $g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33 ) =
($a22*$a33-$a23*$a32, $a13*$a32-$a12*$a33, $a12*$a23-$a13*$a22,
$a23*$a31-$a21*$a33, $a11*$a33-$a13*$a31, $a13*$a21-$a11*$a23,
$a21*$a32-$a22*$a31, $a12*$a31-$a11*$a32, $a11*$a22-$a12*$a21);
( $g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33 ) =
($g11/$detA, $g12/$detA, $g13/$detA,
$g21/$detA, $g22/$detA, $g23/$detA,
$g31/$detA, $g32/$detA, $g33/$detA );


## 一方、右辺 B
( $b1, $b2, $b3 ) = ( 1, 0, -2 );

## そして、逆行列により解を算出
$X = $g11*$b1+$g12*$b2+$g13*$b3;
$Y = $g21*$b1+$g22*$b2+$g23*$b3;
$Z = $g31*$b1+$g32*$b2+$g33*$b3;

## 確認
print( $X, ", ", $Y, ", ", $Z, ";\n" ); #> -2.3333333, 2.3333333, -1;
出来ました。

当たり前?


本日はここまで。


時には、Perl 学習もある?


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


171020
スポンサーサイト

3元1次連立方程式を解く

2017-11-15 :  PCクリニック
Python、C言語、Perl、グルコサミン、Firefox
昨年(2016-02-02)の記事「GSL Shell 学習:matrix.solve( A, b )」で、
3元1次連立方程式を解くには、matrix.solveで十分であった。

今回、
訳あって、自前コードで処理したくなった。


参考サイトとしては、下記の2つ。

連立方程式の解き方」から、

■連立1次方程式の解き方

  未知数3個、方程式3個
  a11*x + a12*y + a13*z = b1
  a21*x + a22*y + a23*z = b2
  a31*x + a32*y + a33*z = b3

は、行列を用いて、

  A * X = B

の様に書ける。

それで、
方法は2通りあるが、

(1) 逆行列を用いる方法

  X = Inv_A * B

とやれば良い。


それから、
逆行列」から、

3×3行列の逆行列の公式

  A = { { a11, a12, a13 }, { a21, a22, a23 }, { a31, a32, a33 } }

について、

  detA = a11*a22*a33 + a21*a32*a13 + a31*a12*a23 -
      a11*a32*a23 - a31*a22*a13 - a21*a12*a33

がゼロでないとき、

逆行列:
 Inv_A = 1 / detA *
   { { a22*a33 - a23*a32 , a13*a32 - a12*a33 , a12*a23 - a13*a22 },
    { a23*a31 - a21*a33 , a11*a33 - a13*a31 , a13*a21 - a11*a23 },
    { a21*a32 - a22*a31 , a12*a31 - a11*a32 , a11*a22 - a12*a21 } }


以上の公式に従って、
昨年の記事のデータで、
素直にコーディング
----- Solve_by_Inv-Mat.gsl
-- 行列 A
a11, a12, a13, a21, a22, a23, a31, a32, a33 =
0, 3, 6, 1, 4, 7, 2, 5, 9

detA = a11*a22*a33 + a21*a32*a13 + a31*a12*a23 -
a11*a32*a23 - a31*a22*a13 - a21*a12*a33

-- 逆行列 Inv_A == G を、2段階で
g11, g12, g13, g21, g22, g23, g31, g32, g33 =
a22*a33-a23*a32, a13*a32-a12*a33, a12*a23-a13*a22,
a23*a31-a21*a33, a11*a33-a13*a31, a13*a21-a11*a23,
a21*a32-a22*a31, a12*a31-a11*a32, a11*a22-a12*a21
g11, g12, g13, g21, g22, g23, g31, g32, g33 =
g11/detA, g12/detA, g13/detA,
g21/detA, g22/detA, g23/detA,
g31/detA, g32/detA, g33/detA

-- 一方、右辺 B
b1, b2, b3 = 1, 0, -2

-- そして、逆行列により解を算出
X = g11*b1+g12*b2+g13*b3
Y = g21*b1+g22*b2+g23*b3
Z = g31*b1+g32*b2+g33*b3

-- 確認
print( X, Y, Z ) --> -2.3333333, 2.3333333, -1
出来ました。

当たり前?


本日はここまで。


これで、(やりたかったコト) Perl 版を作ろう。


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


171020

DICOM ファイルを TIFF に変換(再改訂版)

2017-11-13 :  PCクリニック
Python、C言語、Perl、グルコサミン、Firefox
前(2017-11-10)の記事「DICOM ファイルを TIFF に変換(改訂版)」を
書いたばかりだが、
その後に、
更なる新たな .dcm ファイルを入手した。

だが、このファイル形式には
対応できなかった。

データ部分のサイズを指定する、
変数 g で表しているコード値が 0x7FE0 のモノが2つあった。

これの識別は難しいので、

g が 40 で、
e が 0x0010 と 0x0011 の時の val 値から、
データサイズを取得したら、
後は、
ファイルの最後から
データ本体部分を取得すること
にした。

プログラム名を“DICOM_to_TIFF_Rev2.gsl”とした:
----- DICOM_to_TIFF_Rev2.gsl -----
require'pl'; _=path.currentdir()
require'cv2_a'; SF = string.format
uint16 = |x| ffi.cast( 'uint16_t*', x )
uint32 = |x| ffi.cast( 'uint32_t*', x )

ReadUInt16 = || uint16( IN:read(2) )[0]
ReadUInt32 = || uint32( IN:read(4) )[0]
ReadChars = |n| IN:read(n)

Fname = 'x.dcm'; IN = io.open( Fname, 'rb' ) -- サンプルデータ --
FSZ = lfs.attributes( Fname, 'size' ) -- ファイルサイズ

-- 先頭 128 バイトは Preamble ( all ゼロ )
-- 次の4バイトは:'DICM'

zzz = IN:read(128); fid = IN:read(4)
if fid ~= 'DICM' then
print( 'Not a DICM' ); _=io.read(1); os.exit()
end

g = ReadUInt16()
------- g の判定は、以下のループの中に纏めた。
while g==2 or g==8 or g==16 or g==24 or g==32 or
g==40 or g==64 or g==0x7005 or g==0x7FE0 do
e = ReadUInt16();
vr = ReadChars(2)
if vr == 'AE' or vr == 'AS' or vr == 'AT' or vr == 'CS' or
vr == 'DA' or vr == 'DS' or vr == 'DT' or vr == 'FL' or
vr == 'FD' or vr == 'IS' or vr == 'LO' or vr == 'PN' or
vr == 'SH' or vr == 'SL' or vr == 'SS' or vr == 'ST' or
vr == 'TM' or vr == 'UI' or vr == 'UL' or vr == 'US'
then
length = ReadUInt16()
else
_tmp = ReadUInt16() -- Read the reserved byte
if string.sub(vr,1,1) == 'O' then
-- g==2 での'OB' と、0x7FE での'OW'
length = ReadUInt32()
else
length = ffi.cast('uint16_t*',vr)[0] + _tmp*256*256 -- *****
end
end
----- print( g, e, vr, length, SF('%X',IN:seek()) ) -- デバッグ用

if g==0x7FE0 then break end
val = ReadChars(length)
---------------------------
if g==40 then
if e==0x0010 then yy = uint16(val)[0] end
if e==0x0011 then xx = uint16(val)[0] end
if e==0x0100 then bpp= uint16(val)[0] end
end
---------------------------
g = ReadUInt16()
end

if bpp ~= 16 then print( 'Not 16bits' ); _=io.read(1); os.exit() end

----- 「g==0x7FE0」が2つあるので、全体から逆算する:
f_pos=IN:seek()
print( z, xx, yy, bpp, length + f_pos ) -- 進行確認用

nokori = FSZ - f_pos
if nokori > xx*yy*2 then _ = ReadChars( nokori- xx*yy*2 ) end

AT = uint16( IN:read('*a') ) -- Read '*a' == length == xx*yy*2
IN:close()

------------------------------
tif = cv2.cvCreateMat( yy, xx, CV_16UC1 )
BT = ffi.cast( 'uint16_t*', tif.Byte )

for n=0,xx*yy-1 do
if AT[0]>32767 then
BT[n] = AT[n] + 0x8000 --========
else
BT[n] = AT[n]
end
end

if string.find( Fname, 'dcm' ) ~= nil then
FN = string.gsub( Fname, '.dcm', '.tif' )
else
FN = Fname..'.tif'
end
cv2.cvSaveImage( 'o-'..FN , tif, 0 )
これで、?????

これで、
.dcmファイルを 16 bit TIFFファイルに変換できた。


前回のものより更なる汎用化ができた?


本日はここまで。


DICOM 学習は続く?


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


171021

DICOM ファイルを TIFF に変換(改訂版)

2017-11-10 :  PCクリニック
Python、C言語、Perl、グルコサミン、Firefox
先日(2017-10-16)の記事「DICOM ファイルを TIFF に変換」では、
.nii ファイルを経由せずに、
.dcm ファイルから 直接 TIFF ファイルに変換する
プログラム“DICOM_to_TIFF.gsl”を作った。

しかし、
新たな .dcm ファイルを入手したが、
対応できないモノだった。

更なる DICOM ファイルの仕様を解読し、
何とか読み込めるメドついた。

ポイントは、前回コードの
変数 g で表しているコード値が 2 以外の時でも
変数 vr で表しているタイプが指定されているモノでした。

これに対応した処理方式と、以前のデータに対する処理方式を
併せて、キモチ汎用版を作った。

プログラム名を“DICOM_to_TIFF_Rev.gsl”とした:
----- DICOM_to_TIFF_Rev.gsl -----
require'pl'; require'cv2_a'; SF = string.format
uint16 = |x| ffi.cast( 'uint16_t*', x )
uint32 = |x| ffi.cast( 'uint32_t*', x )

ReadUInt16 = || uint16( IN:read(2) )[0]
ReadUInt32 = || uint32( IN:read(4) )[0]
ReadChars = |n| IN:read(n)

Fname = 'x.dcm'; IN = io.open( Fname, 'rb' ) -- サンプルデータ --

-- 先頭 128 バイトは Preamble ( all ゼロ )
-- 次の4バイトは:'DICM'

zzz = IN:read(128); fid = IN:read(4)
if fid ~= 'DICM' then
print( 'Not a DICM' ); _=io.read(1); os.exit()
end

g = ReadUInt16()
------- g の判定は、以下のループの中に纏めた。
while g==2 or g==8 or g==16 or g==24 or g==32 or
g==40 or g==64 or g==0x7005 or g==0x7FE0 do
e = ReadUInt16();
vr = ReadChars(2)
if vr == 'AE' or vr == 'AS' or vr == 'AT' or vr == 'CS' or
vr == 'DA' or vr == 'DS' or vr == 'DT' or vr == 'FL' or
vr == 'FD' or vr == 'IS' or vr == 'LO' or vr == 'PN' or
vr == 'SH' or vr == 'SL' or vr == 'SS' or vr == 'ST' or
vr == 'TM' or vr == 'UI' or vr == 'UL' or vr == 'US'
then
length = ReadUInt16()
else
_tmp = ReadUInt16() -- Read the reserved byte
if string.sub(vr,1,1) == 'O' then -- g==2 の'OB' と、0x7FE の'OW'
length = ReadUInt32()
else
length = ffi.cast('uint16_t*',vr)[0] + _tmp*256*256 -- *****
end
end
----- print( g, e, vr, length, SF('%X',IN:seek()) ) -- デバッグ用

if g==0x7FE0 then break end
val = ReadChars(length)
---------------------------
if g==40 then
if e==0x0010 then yy = uint16(val)[0] end
if e==0x0011 then xx = uint16(val)[0] end
if e==0x0100 then bpp= uint16(val)[0] end
end
---------------------------
g = ReadUInt16()
end

if bpp ~= 16 then print( 'Not 16bits' ); _=io.read(1); os.exit() end

AT = uint16( IN:read('*a') ) -- Read '*a' == length == xx*yy*2
IN:close()

------------------------------
tif = cv2.cvCreateMat( yy, xx, CV_16UC1 )
BT = ffi.cast( 'uint16_t*', tif.Byte )

for n=0,xx*yy-1 do
if AT[0]>32767 then
BT[n] = AT[n] + 0x8000 --========
else
BT[n] = AT[n]
end
end

if string.find( Fname, 'dcm' ) ~= nil then
FN = string.gsub( Fname, '.dcm', '.tif' )
else
FN = Fname..'.tif'
end
cv2.cvSaveImage( 'o-'..FN , tif, 0 )
これで、?????

これで、
.dcmファイルを 16 bit TIFFファイルに変換できた。
前回のものよりキモチ汎用で。


本日はここまで。


DICOM 学習は続く?


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


171018

GSL Shell コード・テンプレート

2017-11-08 :  PCクリニック
Python、C言語、Perl、グルコサミン、Firefox
昨年暮れ(2016-12-26)に「viva “GSL Shell”」を書いた。

その後も、<紙>にとってのメイン言語 と云うことで、

新規に作るプログラムは、原則 GSL Shell で作り、
更に過去に作った、Perl や、R や、Python プログラムも
GSL Shell に移植してきた。


ここで、
今後のコーディングの見本として、

Code Template を作ってみた。

プログラム名は「Code_Template.gsl」:
----- Code_Template.gsl     <== プログラムの名前

--[[
   プログラムの仕様・コメント 等々

  プログラム・ファイルのコード:'Shift-JIS'

  コマンドライン引数:
   arg 引数リスト    #arg 引数の個数
arg[0] 当プログラム名

]]
--===========================
require'cv2_a' -- 'OpenCV' 使用: cv2.cvLoadImageM 等々
require'gsl_a' -- 'GSL' 使用: gsl.gsl_rng_get 等々
gp = require'gnuplot' -- 'gnuplot' 使用: gp.Plot 等々
require'pl' -- 'Penlight'使用: path.currentdir 等々

require'winman' -- 'afxLua/winman'使用: winman.send 等々
Iup = require'IUP' -- 'IUP' 使用: Iup.Open 等々

lfs = require'lfs' -- 'LFS' 使用: lfs.attributes 等々
--- 'Penlight' 経由で 'lfs' を使うなら
_ = path.currentdir()

--===========================
-- use'math' -- 'math.log' を 'log' と モジュール名の省略
-- use'bit' -- 'bit.band' を 'band'   ※ 使用は要注意

f = |x| x^2 - 1 --=== function f(x) return x^2-1 end

--===========================

----- ディレクトリ内ファイル一覧取得:
require'pl'
f_list = dir.getfiles( Dir )


----- ファイル(不)存在チェック:
require'afx'
if not afx.fileexists( f ) then
-- ~
end


-- D&D でファイル名取得:
require'pl'
Dir, Fname = path.splitpath( arg[1] )
lfs.chdir( Dir ) -- これ必要???


----- 入力テキストファイルについて、
-- ファイル全体を一括で(文字列として)読み込む:
require'pl'
wstr = file.read( 'hoge.txt' )

-- ファイル全体を行毎の文字列配列で読み込む:
require'pl'
w = utils.readlines( 'hoge.txt' )
-- 更に、一行分を分割(空白区切り):
v = stringx.split(w[i])

-- CSVファイル全体をテーブルとして読み込む:
t = csv.read( 'hoge.csv' )
-- 更に、数値のみなら行列に変換:
m = matrix.def( t )
----- 非数値データの行(テーブル要素)を除くには:
table.remove( t, 1 ) -- とかで。

----- 入力バイナリファイルについて、
-- ファイル全体を4バイト整数配列として読み込む:
IN=io.open('hoge','rb'); buf=IN:read('*a'); IN:close()
pt=ffi.cast('int32_t*',buf)


----- 出力テキストファイルについて、
-- 'fprintf' で書き出し:
require'pl'; OT=io.open('f.txt','w')
utils.fprintf( OT, ' %g %d %s\n' , floatv, intv, str )
OT:close()

--======================
----- 文字列が数値を表すか の判定
if tonumber( str ) ~= nil then
-- ~
end

--======================

----- Python の 'range' 相当
require'pl'

for k in seq.range( 0, n-1 ) do -- for k in range( n ) :
-- ~ -- # ~
end -- # 無し

for k in seq.range( i, j-1 ) do -- for k in range( i, j ) :
-- ~ -- # ~
end -- # 無し


plc = require'pl.comprehension'.new()
T = plc('x for x=1,5')() -- T は { 1, 2, 3, 4, 5 }
--- plc 'x for x=1,5' () -- これでも可
V = matrix.vec( plc'x for x=1,5'() ) -- これで 'matrix'ベクトル


--======================
----- 三角関数 逆三角関数
-- sin、cos、tan、・・・・・ asin、acos、atan、atan2

----- 特殊関数
-- sf.fact(n) ・・・ 階乗(factorial) n!
-- sf.gamma(x)・・・ ガンマ関数(gamma function)
-- x は正数。 use:real Lanczos method ( ランチョス法 )


--======================
--[[ --- 過去記事/備忘録 ---

備忘録:(2017-04-08)記事
Lua 言語仕様:備忘録(忘れやすい事)4言語対比表
「H29-03-04_Lua言語仕様_備忘録.html」H29-03-22 /_GSL/

関数フィッティング:(2017-04-05)記事
GSL Shell 学習:Fitting サブ(改)
「SUB_GSL_Fit.lua」H28-02-04 /_GSL/H28-12-10_FitGSL/

非線形1変数方程式の求解:(2017-06-21)記事
GSL Shell:非線形1変数方程式の求解モジュール自作
「NL_1var_Rev.gsl」H28-05-07 /_GSL_tcc/H29-05-06_Solve_Newton/

非線形2変数方程式の求解:(2017-06-24)記事
GSL Shell:非線形2変数方程式の求解モジュール自作
「NL_2var_Rev2.gsl」H28-05-11 /_GSL_tcc/H29-05-06_Solve_Newton/

スプライン補間:(2016-09-19)記事
GSL Shell 学習:スプライン補間
「GDT_Spline.gsl」H28-08-12 /_GSL/H28-08-12_Spline/

'OpenCV' サンプルコード:(2017-07-30)記事
OpenCV:cvFindContours 正常に出来た
「cvFitEllipse2_7.gsl」H29-06-28 /_GSL_CV/H29-06-27_FindContours_FitEllipse2/

'GSL' サンプルコード:(2017-08-23)記事
GSL 学習:配列要素の並び替え
「gsl_sort.gsl」H29-07-17 /_GSL_gslib/H29-07-17_Sorting/

'gnuplot' 用モジュール定義:(2017-11-06)記事
Lua & gnuplot:gnuplot モジュール改訂2
「gnuplot.lua」H29-06-20 /_gnuplot/H29-06-20_gp_lua/

'IUP' サンプルコード:(2017-02-12)記事
IUP で GUI アプリ:コード改訂
「IUP_Exact_Dbl_Entou_Rev.gsl」H29-01-18 /_GSL/H29-01-18_IUP_3/

'afxLua' サンプルコード:'afx'、'cio'、'clipboard'、'winman'、'winreg'
検索「afxLua
]]
取り敢えずは、こんなところ?



本日はここまで。


Lua ( GSL Shell ) 学習は続く。


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


171009,13

Lua & gnuplot:gnuplot モジュール改訂2

2017-11-06 :  PCクリニック
Python、C言語、Perl、グルコサミン、Firefox
6月(2017-06-26)の記事「Lua & gnuplot:gnuplot モジュール改訂」で、
  ・・・・・
  ・・・・・
  今回インターフェイスは変えずに、若干の改訂を行った。

  1.生成する数値のフォーマットを変更

  ・・・・・
  ・・・・・

  2.フォント指定を変更

  ・・・・・
  ・・・・・

と、若干の改訂を行っているが、

今回、gnuplot.luaモジュールについて、
更なる若干の改訂を行った。

1.カラーパレット定義 関数を追加
  ・gp.Palette( )
  “シキノート”サイトの
  「gnuplotのカラーマップ」記事を参考にさせて頂いた。

   GNUPLOTのカラーマップ 2015年5月16日
   一番良いと思ったカラーマップは、matlabで使われている”jet”だと感じました。
  と云うことで、これを頂きました。


2.matrix データ定義 関数を追加
  ・gp.Mat( ~ )
  6月(2017-06-30)の記事「gnuplot 学習:matrix 形式データ
  ・・・・・
  ・・・・・
  このことには要注意だが、
  matrix形式データが定義できる様にしたい。
   先日(2017-06-26)の記事:
  「Lua & gnuplot:gnuplot モジュール改訂」で書いた、
  改訂“gnuplot”モジュールの追加改訂を行おう。
  ・・・・・
  ・・・・・
  に対応した追加作業。


結果は:
----- gnuplot.lua -----
local gp = {}

local SF = string.format
local SS = string.sub
local function split_path( p ) return string.match( p, "(.-)([^\\]-([^%.]+))$" ) end

function gp.Open( fn )
if fn == nil then
GP = io.popen( 'gnuplot -persistent', 'w' ) ----- ファイル名無し:画面出力
else
local Dir, FN, ext = split_path( fn )
if ext ~= 'plt' then FN = FN..'.plt' end ----- .plt 無し なら付加する
GP = io.open( FN, 'w' )
end
GP:write( 'set term wxt size 800, 600\n' ) ----- 窓サイズ
GP:write( 'set term wxt font "Arial,12"\n' ) ----- フォント
GP:write( 'set size 1.0, 0.98\n' ) ----- グラフサイズ
end

function gp.Pdf( fn )
local Dir, FN, ext = split_path( fn )
if ext ~= 'pdf' then FN = FN..'.pdf' end ----- .pdf 無し なら付加する
GP:write( 'set term pdf size 8, 6\n' ) ----- ページサイズ
GP:write( 'set term pdf font "Arial,12"\n' ) ----- フォント
GP:write( SF('set output "%s"\n', FN ) ) ----- PDFファイル名
end

function gp.Close() ----- ファイルclose
GP:close()
end

function gp.Com( c ) ----- 任意の コマンド用
GP:write( SF('%s\n', c ) ) ----- コマンド&パラメータ
end

function gp.Set( s ) ----- set コマンド用
GP:write( SF('set %s\n', s ) ) ----- パラメータ
end

function gp.UnSet( s ) ----- unset コマンド用
GP:write( SF('unset %s\n', s ) ) ----- パラメータ
end

function gp.Table( fn )
if fn == nil then
GP:write( SF('set table\n', fn ) ) ----- table コマンド用
elseif SS(fn,1,1) == '$' then
GP:write( SF('set table %s\n', fn ) )
else
GP:write( SF('set table "%s"\n', fn ) )
end
end

function gp.Title( t, s )
if s == nil then s = 36 end
GP:write( SF('set title "%s"\n', t ) ) ----- タイトル
GP:write( SF('set title font "Arial,%d"\n', s ) ) ----- フォント
end

function gp.Xlog()
GP:write( 'set logscale x\n' ) ----- 対数目盛
end

function gp.Xlabel( L, s )
if s == nil then s = 20 end
GP:write( SF('set xlabel "%s"\n', L ) ) ----- ラベル
GP:write( SF('set xlabel font "Arial,%d"\n', s ) ) ----- フォント
end

function gp.Xrange( s )
GP:write( SF('set xrange %s\n', s ) ) ----- 範囲
end

function gp.Ylog()
GP:write( 'set logscale y\n' ) ----- 対数目盛
end

function gp.Ylabel( L, s )
if s == nil then s = 20 end
GP:write( SF('set ylabel "%s"\n', L ) ) ----- ラベル
GP:write( SF('set ylabel font "Arial,%d"\n', s ) ) ----- フォント
end

function gp.Yrange( s )
GP:write( SF('set yrange %s\n', s ) ) ----- 範囲
end

function gp.KeyOff()
GP:write( 'set key off\n' ) ----- 凡例無し
end

function gp.Palette()
GP:write( 'set palette defined ( 0 "#000090", 1 "#000FFF", 2 "#0090FF", 3 "#0FFFEE", 4 "#90FF70", 5 "#FFEE00", 6 "#FF7000", 7 "#ee0000", 8 "#7F0000")\n' )
end ----- カラーパレット:jet

function gp.Plot( p ) ----- 引数そのママ
GP:write( SF('plot %s\n', p ) ) ----- Plot コマンド
end

function gp.RePlot( p )
GP:write( SF('replot %s\n', p ) ) ----- replot
end

function gp.Data( b, n, x, y ) ----- b:origin
local px, py = x, y
if b==0 then ----- 1 or 0
px, py = ffi.cast('double*', x ), ffi.cast('double*', y )
end
for k=b,n-(1-b) do
GP:write( SF('%g %g\n', px[k], py[k] ) )
end
GP:write( 'e\n' )
end

function gp.DataBlk( BN, b, n, x, y, ... ) ----- BN:Block名
local t, pt = {...}, {} ----- 3列以上は、b==0 の場合のみ可
local px, py = x, y ----- b==1 では無視
if b==0 then
px, py = ffi.cast('double*', x ), ffi.cast('double*', y )
for i,p in ipairs(t) do pt[i]=ffi.cast('double*', t[i] ) end
end
GP:write( SF('%s << EOD\n', BN ) )
for k=b,n-(1-b) do
GP:write( SF('%g %g', px[k], py[k] ) )
for i,_ in ipairs(pt) do GP:write( SF(' %g', pt[i][k] ) ) end
GP:write( '\n' )
end
GP:write( 'EOD\n' )
end

function gp.Mat( BN, M, yy, xx ) ----- BN:Block名
local m = ffi.cast('double*', M ) ----- M : yy行 xx列
GP:write( SF('%s << EOD\n', BN ) )
for y=0,yy-1 do
for x=0,xx-1 do
GP:write( SF('%g ', m[y*xx+x]) )
end
GP:write( '\n' )
end
GP:write( 'EOD\n' )
end

return gp
となった。


本日はここまで。


Lua ( GSL Shell ) / gnuplot 学習は続く。


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


170620,30

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

image_viewer 自作

2017-11-01 :  PCクリニック
Python、C言語、Perl、グルコサミン、Firefox
7ヶ月前(2017-04-03)の記事:
OpenCV 学習:任意のタイプの画像表示+alpha」に載せた、
DnD_img_mouse_CV.gsl は、“グレースケール限定”だった。


今回、この制限を解除すべく、

1.機能拡張を行う。

序でに、

2.'cv2' モジュールの代わりに 'cv2_a' モジュールを使う。

3.'lfs' モジュールの代わりに 'Penlight' モジュールを使う。


それで、
出来上がったのは、・・・・・
----- DnD_image_viewer.gsl -----
require 'cio'; require'pl'; require'cv2_a'
function CHK_IF( pt )
local t, n = 5, 0; while pt[n]==0 do n=n+1 end -- 必ず終わるハズ!
if math.abs(pt[n])<2^29 then t=4 end
return t
end
-------------------------------------
Dir, Fname = path.splitpath( arg[1] ) -- D&D でファイル名取得
lfs.chdir( Dir ) -- これがないと、ダメ?
-- ==================================
local onmouse = ffi.cast( 'CvMouseCallback', function(event,x,y,flags)
if event == 0 then -- cv2.EVENT_MOUSEMOVE:
cio.gotoxy(0,3)
if typ>10 then pt=3*(y*xx+x)
b,g,r = IMG[pt],IMG[pt+1],IMG[pt+2]
print( string.format('Pos:(%d,%d)=%g,%g.%g ', x, y, r,g,b ) )
else
print( string.format('Pos:(%d,%d)=%g ', x, y, IMG[y*xx+x] ))
end
end
end )
-- ==================================
img = cv2.cvLoadImageM( Fname, -1 ) -- ファイル内のママ
yy, xx, typ = img.yy, img.xx, img.type%256
if typ == 0 then IMG=img.Byte
elseif typ==2 then IMG=img.Word
elseif typ==4 then IMG=img.Int ----- これ存在し得ない???
elseif typ==5 then IMG=img.Int
typ=CHK_IF(IMG)
if typ==5 then IMG=ffi.cast('float*',IMG) end
elseif typ==6 then IMG=img.Dbl
else typ=16; IMG=img.Byte ----- RGB 画像  ※ 今回機能追加
end

FT={ '8I', '', '16I', '', '32I', '32F', '64F' }; FT[16+1]='BGR'

if typ>10 then --- カラー画像 ------

dsp = cv2.cvCreateMat( yy, xx, 16 )
for n=0,3*xx*yy-1 do dsp.Byte[n]=IMG[n] end

cv2.cvNamedWindow( charA('img'), 0 )
cv2.cvSetMouseCallback( charA('img'), onmouse, NULL )
cio.gotoxy(0,0)
print( 'ESC key <=== Exit' )
print(string.format('File Type: %s. Size: %dx%d', FT[typ+1], xx, yy))
print( 'Value : R, G, B' )
while true do
cv2.cvShowImage( charA('img'), dsp )
k = bit.band( 0xFF, cv2.cvWaitKey(50) ) -- key bindings
if k == 27 then break end -- esc to exit
end

else ----- GrayScale画像 -------

Mx, Mn = -1e+100, 1e+100
for n=0,xx*yy-1 do
if IMG[n]>Mx then Mx=IMG[n] end
if IMG[n] end

dsp = cv2.cvCreateMat( yy, xx, 0 )
for n=0,xx*yy-1 do dsp.Byte[n]=(IMG[n]-Mn)/(Mx-Mn)*255+0.5 end

cv2.cvNamedWindow( charA('img'), 0 )
cv2.cvSetMouseCallback( charA('img'), onmouse, NULL )
cio.gotoxy(0,0)
print( 'ESC key <=== Exit' )
print(string.format('File Type: %s. Size: %dx%d', FT[typ+1], xx, yy))
print(string.format('Value : Max=%g, Min=%g', Mx, Mn))
while true do
cv2.cvShowImage( charA('img'), dsp )
k = bit.band( 0xFF, cv2.cvWaitKey(50) ) -- key bindings
if k == 27 then break end -- esc → exit
end

end

cv2.cvDestroyAllWindows() -- 後処理
これで、・・・・・

対象とする画像の 種類は、
 ・ 8ビット整数
 ・ 16ビット整数
 ・ 32ビット整数
 ・ 32ビット浮動小数点数
 ・ 64ビット浮動小数点数
 ・ RGB(8ビット整数×3)
となる。


本日はここまで。


GSL Shell ( Lua ) / OpenCV 学習は続く。


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


170813

GSL Shell プログラム:xyz2mol.gsl

2017-10-30 :  PCクリニック
Python、C言語、Perl、グルコサミン、Firefox
前(2017-10-27)の記事「Python プログラム:xyz2mol.py」で、
既存の“xyz2mol”プログラム
即ち、

   座標データだけから成るファイルを元に、
   原子間の結合データを生成

する処理コード
を纏めた。


そこで、
このPython プログラムを、
Lua(GSL Shell)に移植した。


取り敢えずは単純な移植なのだが、
序でなので、
(1) 入力ファイルは、ドラッグ&ドロップで指定する
(2) 原子数(と結合数)は最大 999 までとする
ところまでは機能強化した。

--------------- DnD_xyz2mol.gsl
require'pl'
--------------------
Dir, Fname = path.splitpath( arg[1] ) -- D&D でファイル名取得
lfs.chdir( Dir ) -- これ必須?
file_nm = string.gsub( Fname, '.xyz', '' ) -- 入力「.xyz」ファイル ( 除 拡張子 )
-- 出力ファイル名は、(file_nm).mol ( 結合情報付き )

Dx=1.10 -- 結合判定用:共有結合半径の和の Dx 倍以内
--====================================
-- 1...92 元素名(左詰め)
tt = ' H HeLiBeB C N O F NeNaMgAlSiP S ClAr KCaScTiV CrMnFeCoNiCuZnGaGeAsSeBrKrRbSrY ZrNbMoTcRuRhPdAgCdInSnSbTeI XeCsBaLaCePrNdPmSmEuGdTbDyHoErTmYbLuHfTaW ReOsIrPtAuHgTlPbBiPoAtRnFrRaAcThPaU '

-- 1..92 共有結合半径 ( From Winmostar )
rr = '0.29 1.30 1.23 0.90 0.82 0.77 0.75 0.73 0.72 1.60 1.54 1.36 1.18 1.11 1.06 1.02 0.99 1.90 2.03 1.74 1.44 1.32 1.22 1.18 1.17 1.17 1.16 1.15 1.17 1.25 1.26 1.22 1.20 1.16 1.14 2.00 2.16 1.91 1.62 1.45 1.34 1.30 1.27 1.25 1.25 1.28 1.34 1.48 1.44 1.41 1.40 1.36 1.33 2.20 2.35 1.98 1.69 1.70 1.70 1.70 1.70 1.70 1.70 1.70 1.70 1.70 1.70 1.70 1.70 1.70 1.70 1.44 1.34 1.30 1.28 1.26 1.27 1.30 1.34 1.49 1.48 1.47 1.45 1.46 1.45 2.71 2.40 2.10 1.70 1.70 1.70 1.70'
rTc = stringx.split( rr ) -- 文字列配列(空白区切)

An={}; Ac=matrix.new( 999, 1 )
Ax=matrix.new( 999, 1 ); Ay=matrix.new( 999, 1 ); Az=matrix.new( 999, 1 )
B1=matrix.new( 999, 1 ); B2=matrix.new( 999, 1 ); rT=matrix.new( 92, 1 )
for i=1, #rTc do rT[i]=rTc[i]+0 end

w = utils.readlines( file_nm .. '.xyz' ) -- 入力「.xyz」ファイル取り込み
aa=w[1]+0 -- Atom 数
wl2=w[2] -- セルベクトル情報 ( 文字列のママ )
for a=1,aa do v = stringx.split(w[a+2])
Ax[a]=v[2]+0; Ay[a]=v[3]+0; Az[a]=v[4]+0
an = string.sub(v[1]..' ',1,2); k = string.find( tt, an )
An[a]=an; Ac[a]=k/2
end

bb=0 -- ===== 結合計算 =======
for k=1,aa do
for j=1,k-1 do
dd2=(Ax[k]-Ax[j])^2+(Ay[k]-Ay[j])^2+(Az[k]-Az[j])^2
dr2=((rT[Ac[k]]+rT[Ac[j]])*Dx)^2
if dd2 end
end

OT = io.open( file_nm .. '.mol', 'w' ) --=== 「.mol」ファイル出力
utils.fprintf( OT, '--- from %s.xyz --> %s.mol\n' , file_nm,file_nm )
utils.fprintf( OT, '%s\n' , wl2 )
utils.fprintf( OT, "By 'xyz2mol.gsl'\n" )
utils.fprintf( OT, '%3d%3d 0 0 0 0 1 V2000\n' , aa, bb )
for a=1,aa do
utils.fprintf( OT, ' %9.5f %9.5f %9.5f %-2s 0 0 0 0 0 0\n' ,
Ax[a], Ay[a], Az[a], An[a] )
end
for b=1,bb do
utils.fprintf( OT, '%3d%3d 1 0 0 0\n' , B1[b], B2[b] )
end
utils.fprintf( OT , 'M END\n' )
OT:close()
となった。


本日はここまで。


この GSL Shell プログラムを流用して、目的のモノを作ろう!


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


171009

Python プログラム:xyz2mol.py

2017-10-27 :  PCクリニック
Python、C言語、Perl、グルコサミン、Firefox
3年前(2014-08-08)になる古い記事ですが、
プログラム設計:msi2mol」正しくは、xyz2molでしたが。

それで、
この記事にある設計に従って、
Python プログラム xyz2mol.py を作っていた。

なぜ今頃?

実は同様なことを、最近やりたくなって、
当然 GSL Shell で作る予定だが、・・・

探していたら、過去に作っていたので、
これをベースに移植を考えた。


なので、当時のコードをここに整理した。

# ----------------------- xyz2mol.py
from math import *
from string import *

file_nm = 'test' # 入力「.xyz」ファイル名 ( 除 拡張子 )

# 出力ファイル名は、(file_nm).mol ( 結合情報付き )

Dx=1.10 # 結合判定:共有結合半径の和の Dx 倍以内

tt = " HHeLiBe B C N O FNeNaMgAlSi P SClAr KCaScTi VCrMnFeCoNiCuZnGaGeAsSeBrKrRbSr YZrNbMoTcRuRhPdAgCdInSnSbTe IXeCsBaLaCePrNdPmSmEuGdTbDyHoErTmYbLuHfTa WReOsIrPtAuHgTlPbBiPoAtRnFrRaAcThPa U"
# 1...92 元素名(右詰め)

rr = " 0.00 0.29 1.30 1.23 0.90 0.82 0.77 0.75 0.73 0.72 1.60 1.54 1.36 1.18 1.11 1.06 1.02 0.99 1.90 2.03 1.74 1.44 1.32 1.22 1.18 1.17 1.17 1.16 1.15 1.17 1.25 1.26 1.22 1.20 1.16 1.14 2.00 2.16 1.91 1.62 1.45 1.34 1.30 1.27 1.25 1.25 1.28 1.34 1.48 1.44 1.41 1.40 1.36 1.33 2.20 2.35 1.98 1.69 1.70 1.70 1.70 1.70 1.70 1.70 1.70 1.70 1.70 1.70 1.70 1.70 1.70 1.70 1.44 1.34 1.30 1.28 1.26 1.27 1.30 1.34 1.49 1.48 1.47 1.45 1.46 1.45 2.71 2.40 2.10 1.70 1.70 1.70 1.70" # 1..92 共有結合半径 ( From Winmostar )
rTc=rr.split() # 文字列配列(空白区切)

An=['']*300; Ax=[0]*300; Ay=[0]*300; Az=[0]*300;
Ac=[0]*300; B1=[0]*300; B2=[0]*300; rT=[0]*100; # 配列(リスト)初期化
for i in range(len(rTc)): rT[i]=float(rTc[i])

IN = open( file_nm + '.xyz' ) # === 入力「.xyz」ファイル取り込み ===
w=IN.readline(); aa=int(w) # Atom 数
wl2=IN.readline()
for a in range(aa):
w=IN.readline(); w=w.rstrip('\r\n'); v=w.split()
an=v[0]; Ax[a]=float(v[1]); Ay[a]=float(v[2]); Az[a]=float(v[3])
if( len(an)==1 ): an=" "+an
k=tt.find(an)
An[a]=v[0]; Ac[a]=k/2
IN.close()

bb=0 # ===== 結合関係計算 =====
for k in range(1,aa):
kb=k-1
for j in range(kb):
if( Ac[k]==Ac[j] ): continue
dd2=(Ax[k]-Ax[j])**2+(Ay[k]-Ay[j])**2+(Az[k]-Az[j])**2
dr2=((rT[Ac[k]]+rT[Ac[j]])*Dx)**2
if( dd2
OT = open( file_nm + '.mol', 'w' ) # === 「.mol」ファイル出力 ===

OT.write( '--- from ' + file_nm + '.xyz --> ' + file_nm + '.mol\n' )
OT.write( wl2 )
OT.write( "By 'xyz2mol.py'\n" )
OT.write( str(aa).rjust(3) + str(bb).rjust(3) + ' 0 0 0 0 1 V2000\n' ) # sprintf 形式
for a in range(aa):
OT.write( " %9.5f %9.5f %9.5f %-2s 0 0 0 0 0 0\n" % ( Ax[a], Ay[a], Az[a], An[a] ) ) # printf 形式
for b in range(bb):
OT.write( "%3d%3d 1 0 0 0\n" % ( B1[b], B2[b] ) )
OT.write( "M END\n" )
OT.close();
と云うこと?


本日はここまで。


このプログラムを参考にして GSL Shell で作ろう!


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


140215,171006
おきてがみ/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
ブックマーク