setup diary

_ atomic form factor using R

使っている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)}

などとすると、原子形状因子をプロットできるようになった。