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で、これにはパッケージが無くて使えないという罠があった。

raspberry pi3を通じて

windowsで制御している装置のデータを、LANを通じて別のPCから見れるようにする際に、そのwindows PCはインターネットに接続したくなかったので、間にraspberry pi3を入れてみた。

raspberry pi3とそのwindows PCを有線LANで繋ぐ計画だったのだが、装置の制御に有線LANを使っていたので、まずはIPを調べる必要があった。その結果、PCが1で、装置が10と20のようだった。重複しないように、raspberry piは53にしてみた。なんとなく素数にしたくなってしまう。routerは未設定にした。これで、制御PCと装置とraspberry piのみのLANが構築できる。

次は、windows10でwindows共有をするのだが、これは非常に苦労した。ネットの情報を元にして、いろいろなところをいじったり、何度も再起動していたら、なんとか共有をすることができた。特に最近のwindowsはwindows共有をするのが非常に難しい。

あとはraspberry piにその共有フォルダをマウントするように/etc/fstabを書き換えた。このとき、単にautoとしてもダメで、下記のようにしたらうまく行くようだ。

//ip.1/Data /home/pi/public cifs username=user,password=pass,noauto,x-systemd.automount 0 0

そして、raspberry piを無線LANで繋ぎ、共有したフォルダをftpでanonymousから読み取りのみでアクセスできるようにしておいた。本当はsftpにしたいのだが、windowsからアクセスするのに面倒らしいので、ftpで我慢している。

今の所、装置との干渉も気にならないようだし、自動mountもうまく働いているので、このまま問題が起きないと良いな。あとはUSBの穴を塞いで。