[gensan]ラザフォード散乱の古典軌道を求めるプログラム

学生実験「原子核散乱」(の二つあるテーマのうちの片方)では26 MeVのα粒子の金標的での散乱の微分散乱断面積を測定します。測定した散乱断面積を、古典的な散乱断面積をルンゲクッタ法を用いた軌道計算から求めて比較してみる、というのをレポート課題の一つとしています(必須ではないですが、なるべくたくさんの人に挑戦してもらいたい)。最近、物理学教室のカリキュラム変更があって計算機の授業が始まったこともあり、プログラミングは頑張ってくれる学生も多いのですが、プログラムを作ることで満足してしまって肝心の物理の議論まで辿り着かないレポートも散見されます。もう少しデータを理解する、という点に時間が取れるようサンプルプログラムを用意しました。

ラザフォード散乱の古典軌道計算のプログラムは、早野龍五・高橋忠幸著『計算物理』(共立出版)の5.5節に丁寧な解析とプログラムコードがあります。ただし、少し古い本でほとんどのプログラムがFortranで書かれていますので、C++で書き直したものを載せます。

rutherford.cc

使い方

プログラムはC++で書かれています。計算機実験の授業でCの使い方は習っているはずですが、C++もほとんど同じなので読めば理解出来ると想います。このプログラムでは、特にC++特有のことはほとんどしていないので、標準ライブラリの名前と標準出力の書き方が多少違う程度です。コンパイルはgccではなく、g++を使います。

$ g++ rutherford.cc

で、a.outというファイルができます。実行するときに、入力として衝突係数を入れます。衝突係数を10.0(単位はfm)のときの散乱角度を計算したければ、以下のようにします。

$ ./a.out 10.0

  Alpha beam energy = 26 MeV
  Impact paramter   = 10 fm

  Scatteing Angle   = 47.2239 deg

これで、散乱角度が47.2度と計算されます。同時に、output.gpというファイルができます。gnuplotの実行ファイルなので、Gnuplotで開けば散乱の軌道が表示されるはずです。

$ gnuplot output.gp

サンプルプログラムでは軌道を表示して散乱角度を計算するところまで載せてあります。計算された軌道から微分散乱断面積にするためには、「原子核散乱」の教科書にある(5)式を使ってください。

サンプルプログラムでは、ポテンシャルとして一番簡単な点電荷によるクーロンポテンシャルのみが書かれていますが、これをもとに実験データを再現するようなポテンシャルを探してください。

プログラムの解説

まず、全体の空間を定義します。サンプルプログラムでは、3000 fm (=3 pm)としています。クーロン力が十分小さく軌道が(必要な計算精度の範囲内で)変化しない程度であれば構いません。

単位系は、教科書に従ってガウス系を使っています。計算に必要な定数や質量などの値は教科書にも載せてあります。これらは以下のページで調べることができます。

http://physics.nist.gov/cuu/index.html

1ステップの大きさは、10.0 fm/cとしています。単純に軌道を求めるだけならばこんなに細かくする必要はありませんが、最近接距離近傍で吸収の効果をとりいれる際にもっと細かくしたくなるかもしれません。何れにしても、最近のパソコンならば計算時間が問題になるようなことはないでしょう。

関数solve()がルンゲクッタ法の式が書いてあります。各ステップごとの計算はさらにsolve内で呼ばれているsolve_1step()で定義してあります。つまりここに、微分方程式が書かれています。力は、関数potential()で定義したポテンシャルの微分です。微分計算には、Richardson's extrapolation法というのを使っています。

計算の詳細は、早野龍五・高橋忠幸著『計算物理』(共立出版)の5.5節を読んでください。おそらく絶版(?)で購入はできないと思うので貸し出します。