2020年
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 31

setup diary

2007|12|
2008|01|02|03|04|05|06|07|08|09|10|11|12|
2009|01|02|03|04|05|06|07|08|09|10|11|12|
2010|01|02|03|04|05|06|07|08|09|10|11|12|
2011|01|02|03|04|05|06|07|08|09|10|11|12|
2012|01|02|03|04|05|06|07|08|10|11|12|
2013|01|02|03|04|05|06|07|08|09|10|11|12|
2014|01|02|03|04|06|08|11|
2015|01|02|03|04|05|06|07|08|10|11|12|
2016|01|02|03|04|05|06|07|08|09|10|11|12|
2017|01|02|03|04|05|06|07|08|09|10|11|12|
2018|01|02|03|04|05|06|07|08|09|10|11|12|
2019|01|02|03|04|05|06|07|08|09|10|11|12|
2020|01|02|03|04|05|06|07|08|09|10|11|12|
2021|01|02|03|04|05|06|07|08|09|10|11|12|
2022|01|02|03|04|05|06|07|08|09|10|11|12|
2023|01|02|03|04|05|06|07|08|09|10|11|

2020-12-06 原子形状因子

_ 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)}

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