|
setup diary | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
使っているdebianで何気なく
aptitude search xray
としてみたら、
libxray-absorption-perl libxray-scattering-perl libxray-spacegroup-perl python-aws-xray-sdk python3-aws-xray-sdk
というパッケージが出てきた。二番目と三番目に興味を惹かれたが、まずは二番目をインストールしてみた。すると、/usr/share/perl5/Xray/Scattering/に原子形状を求めるperlスクリプトが入ったので、これをRに移植することにした。そこにあるcromann.dbから係数を読み込んでいるようだが、これをRで読めるようにするために、yaml形式に変更した。libyaml-perlとr-cran-yamlをインストールして、
use strict; use Storable; my $r_cromann = retrieve( "/usr/share/perl5/Xray/Scattering/cromann.db" ); use YAML (); print YAML::Dump($r_cromann);というスクリプトを実行して、その出力をcromann.yamlに入れた。後で分かったことだが、窒素NとイットリウムYはそのままではうまく動かないので、"n"と"y"に変える必要がある。そして、Rで
library("yaml")
crom<-yaml.load_file("cromann.yaml")
#s = sin(theta)/lambda, s = q/4/pi
f0 <- function(sym,s){ #atomic form factor
cef<-crom[[tolower(sym)]]
sapply(s^2,function(ss) sum( cef[1:4*2-1]*exp(-1*ss*cef[1:4*2]) ) + cef[9] ) }
と関数を定義した。ここで、sapplyを使っているのは、sとしてベクトルを与えても動くようにするためにである。これで、
for(ele in c("H","He","Li","Be","B","C","N","O","F","Ne")){
plot(0:100/4,f0(ele,0:100/4/4/pi),type='l',ylim=c(0,10))
par(new=TRUE)}
などとすると、原子形状因子をプロットできるようになった。
Rでデータ処理をしていたのだが、化学式の式量を計算する必要がある。以前はstringrというライブラリを使って計算するプログラムを書いたが、ライブラリのインストールが面倒だということで、ライブラリを使わずに式量を計算するプログラムを書いてみた。
fw<-function(fml){
aw<-read.table("formula.dat",row.names=1)
en<-gregexpr("([A-Z][a-z]?)([0-9.]*)",fml,perl=TRUE)[[1]]
cs<-attr(en,"capture.start")
cl<-attr(en,"capture.length")
sum(apply( matrix(substring(fml,cs,cs+cl-1),,2),1,
function(x) aw[x[1],]*ifelse(x[2]=="",1,as.numeric(x[2])) ))
}
以前書いたものとあまり変わらないが、attrを使ってキャプチャした文字列の場所を変数に入れて、substringで文字列を取り出している。substrを使うと、fmlを複数指定しないといけないので、一括で文字列を取り出す時には、substringの方が便利である。あとは、元素と数を行列の形にして、applyで掛け算をして、sumで足し合わせている。もう少し短く書けそうな気はするけど。