様々な定数
来月から，一つの新しい時代が始まろうとしている。これまで，SI単位は，七つの基本単位が定義され，それらを元に様々な単位が定義されていた。 これが五月二十日から，新しい定義に移行するのである。新しい定義では，三つの基本単位については変更はないが，四つの基本単位は，プランク定数，ボルツマン定数，アボガドロ定数，電気素量から定義されることになった。これまでも光速からメートルが定義されていたが，これで五つの基本的な物理定数が真の意味で定数となり，ここから五つの基本単位ができることとなる。カンデラはなぜ必要なのかはよく分からないが，秒はこれまで通りセシウムを使って定義されている。基本的な物理定数で一定の周期をもつものがあれば，秒もこれを使って定義できるのかも知れないが，今のところそれは難しいようである。また，いずれは精度の高いストロンチウムに置き換えられるのでは無いかと予想している。 これらの物理定数をRを使って使うことも多くなるかも知れないので，以下に定義を書いておこう。 完全な定数になったわけなので、覚えないといけないかな。

nucs<-9192631770 #1/s hyperfine structure transition frezuency of 133Cs
cl<-299792458 #m/s speed of light
hp<-6.62607015e-34 #Js Plank constant
ec<-1.602176634e-19 #C elementary charge
kb<-1.380649e-23 #J/K Boltzmann constant
na<-6.02214076e23 #1/mol Avogadro constant
kcd<-683 #W/sr luminous efficacy of radiation of frequency 540e12Hz

diagの挙動
Rで角運動量演算子を行列表現するときに、diag(seq(j,-j,-1))というようにしていた。しかし、jが0のときには、エラーが出ることが分かった。調べてみると、diagは、成分が二個以上配列を渡されたときには、それを対角成分とする行列を返すが、成分が一個のときには、それをスカラーとみなして、そのスカラーの次元を持つ単位行列を返すことが分かった。そして、0に対しては、0次元の行列をつくろうとする。これがエラーの原因であることが分かった。臨む結果を得るために、if文で条件分岐することにしたが、良い書き方がないかな。

ruby/tkの互換性

昔ruby/tkで作ったプログラムを動かそうと思ったら、エラーが出て動かなかった。rubyのversionによる違いを修正しても動かないので、tk部分だと予想した。エラー箇所をよくよく見てみると、

t=TkcOval.new(@area,x-r,y-r,x+r,y+r,'fill'=>'black')

のところのfillが問題だということが分かった。これを

t=TkcOval.new(@area,x-r,y-r,x+r,y+r)
t.fill='black'

と直したら動くようになったが、一度変数を定義しないと属性を指定できないのは面倒に感じる。

imagemagickではダメだった

複数のpdfを結合して一つのpdfにしたいと思い調べてみた。convertを使って

convert c*.pdf out.pdf

とすれば良いらしいので、やってみたが、長い時間がかかってout.pdfを作ったが、エラーが出てうまくいかなかった。そこで、pdftkを使うことにした。インストールしようとしたら、すでに入っていた。コマンドラインで

pdftk c*.pdf cat output out.pdf

としたら、短時間でファイルが生成された。こちらはうまく読めた。少しコマンドが複雑だが、pdftkは便利そうだ。

固有値と変換行列

最近はPCの性能が上がり、便利なフリーソフトが増えてきたので、PCを使って簡単にいろいろなことができるようになってきている。 様々な学習にもPCを利用すると、理解が深まることがある。 私は式変形にはmaximaを、数値計算にはRを使って、理解の助けにしている。 例として、角運動量をRを使っていじってみよう。

角運動量演算子は、Jzの固有関数を使うと行列表示できる。 これをRで表すと、以下のようになる。

mjz<-function(j){return( diag(seq(j,-j,-1)) )}
mjp<-function(j){
 jp<-matrix(0,nrow=2*j+1,ncol=2*j+1,byrow=TRUE)
 for(i in 1:(j*2)){jp[i,i+1]<-sqrt(i*(j*2+1-i))}
 return(jp)
}
mjm<-function(j){return( t(mjp(j)) )}
mjx<-function(j){return( (mjp(j)+mjm(j))/2 )}
mjy<-function(j){return( (mjp(j)-mjm(j))/2i )}
me<-function(j){return( diag(2*j+1) )}

ここで、mjpとmjmはそれぞれJ+とJ-を、meは2J+1次元の単位行列を表す。 Jxを対角化するには、例えばeigen(mjx(3/2))とすると良い。 valuesに対角成分に対応する固有値が、vectorsには変換係数が出てくる。 このような振る舞いを見てみると、いろいろと理解が深まる。

2020/1/29追記 2019/3/21に指摘したバグも含めて、単純化したものも公開しておきます。そのプログラムがすぐに見つからなくて、公開した方が自分が便利なので。

mjz<-function(j){ diag(j:-j,2*j+1) }
mjp<-function(j){ matrix(sqrt(cumsum(j:-j)*2)*!1:(2*j+1)^2%%(2*j+2),2*j+1) }
mjm<-function(j){ t(mjp(j)) }
mjx<-function(j){ (mjp(j)+mjm(j))/2 }
mjy<-function(j){ (mjp(j)-mjm(j))/2i }
me<-function(j){ diag(2*j+1) }

2020/3/17追記 さらに、シンプルに書くことができたので、更新しておきます。

mjp<-function(j){ diag(sqrt(diffinv(j:-j)*2))[-1,-2*j-2] }

androidとPCの間でのデータ転送
最近、タブレットを使う人が増えてきているように感じる。しかし、ipadとPCの間で、直接データの転送をするのは、なかなか面倒らしい。androidの場合には、どうなのか調べてみた。 linuxにUSBで接続すると、dmesgでは認識しているけど、マウントされない。調べてみると、androidはMTPプロトコルを使っているかららしい。そこで、linuxにjmtpfsをインストールして、android側でMTPを有効にすると、自動的にマウントされるようになり、データの転送も自由にできる。アンマウントが あまりうまく行かないけど。

データフレームの扱い

Rを使って結晶内の特定の距離の範囲内にある原子を数え上げるプログラムを作っていたのだが、実行が異様に遅い。いろいろと試してみると、データフレームを結合するために使っているrbindのところで、ほとんどの時間を使っていることが分かった。webで検索してみると、dplyrというlibraryのbind_rowsを使うと速くなるというので、試してみたら、なんとか許容範囲の実行時間になって来た。

しかし、もっと速くできるのでは無いかといろいろと試行錯誤していたら、かなり速くすることに成功した。expand.gridを使うと、指定したすべての組み合わせのデータフレームを作ってくれるので、それを利用した。極端な例でこれを説明しよう。

球の中にある格子点の数を数えるプログラムを作るとしよう。座標をふって球の中にある場合に、データフレームに加えていくということをする場合には、以下のようなプログラムになる。

d<-data.frame(x=c(),y=c(),z=c())
for(x in -10:10){
for(y in -10:10){
for(z in -10:10){
 if(x^2+y^2+z^2<100){
  d<-rbind(d,data.frame(x=x,y=y,z=z))
 }
}}}

このrbindが時間をかなり使っているようだ。これをexpand.gridで格子点を作って、それを一度に条件判定をするプログラムは、以下のようになる。

d<-expand.grid(x=-10:10,y=-10:10,z=-10:10)
d<-subset(d,x^2+y^2+z^2<100)

このときの実行時間は、数百倍も違う。こんなに違うのかと、少し驚いたが、今後Rを使うときには、このようなことを考慮しようと思う。

久々のlubuntuのインストール
SSDをつんだノートにubuntu18.04を入れたら、以外に苦労した。まず、USBメモリにisoを入れて、そこからブートしてインストールを開始した。しかし、swapパーティションを無くすのを忘れたと思ってと中で止めたのがいけなかったか、それ以降苦労することになった。 partitionerが途中で止まって動かないようになってしまったので、レスキューモードで立ち上げたら、最後までインストール出来た。そのためにユーザーの設定がされない状態になってしまった。そこで、同じくレスキューモードでユーザーを作成した。 後で調べたら、今のubuntuはswapパーティションが無くなって、swapファイルになっていたので、最初のインストールを止める必要が無かったことが分かった。 さらに、aptでupdateをすると、「ハッシュサムが適合しません」と出てしまう。そのようなバグも報告されていたようである。/var/lib/apt/lists/* を消去することによって、この問題も解決した。そして、upgradeをしたら、それなりに使えるようになってきた。

Rのrglを使った原子軌道の角度部分の絵

原子軌道の角度部分の立体的な絵を、ソフトを使って書いてみた。webを調べるとpovrayを使って書いている人が多かったが、特殊関数などの定義が面倒だし、Rのrglを使って書いた方が良い気がして、その方針でやってみた。ソースは以下のとおりです。

library(misc3d)
library(gsl)
theta<-function(x,y,z)acos(z/sqrt(x^2+y^2+z^2))
phi<-function(x,y,z)atan( (sqrt(x^2+y^2)-x)/y )*2
th<- 0:100/100*pi
ph<- 0:200/200*2*pi
fc <- function(x,y,z)ifelse(fr(theta(x,y,z),phi(x,y,z))<0,'blue','red')
fx <- function(t,p)abs(fr(t,p))*sin(t)*cos(p)
fy <- function(t,p)abs(fr(t,p))*sin(t)*sin(p)
fz <- function(t,p)abs(fr(t,p))*cos(t)
for(l in 0:5){
 for(m in 0:l){
  fr <- function(t,p)legendre_Plm(l,m,cos(t))*cos(m*p)
  parametric3d(fx, fy, fz, th, ph,color=fc)
  play3d(spin3d(c(0, 0, 1), 10),5)
}}

ライブラリとしては、特殊関数を使うためにrglを使って、プロットにはmisc3dのparametric3dを使った。極座標のphiは、tan(2phi)から計算するようにして、条件分岐を減らすことができた。とりあえず、h軌道ぐらいまで書いて、五秒間10rpmで回すようにしてみたが、分かりやすくなったかな。

ルジャンドル陪多項式

Rを使って、球面調和関数展開をしようと思ったら、ルジャンドル陪多項式が見当たらない。調べてみたらgslにあった。r-cran-gslをaptで入れたら、

library(gsl)
legendre_Plm(l,abs(m),cos(theta))*exp(-1i*m*phi)

という感じで使うことができた。これで、必要なら軌道の絵なども書けるようになった。しかし、普段使っているノートは、少し古いubuntuで、これにはパッケージが無くて使えないという罠があった。