クーロン力による散乱



原子核に向かう荷電粒子の運動方程式

電荷+Zeの非常に重い原子核に向かって、質量m,電荷+qeの荷電粒子が初速度V0を持って無限遠から衝突パラメーター（無限遠での粒子の運動方向を延長して得られる直線と標的（ここでは原子核）との距離）bで近付くとき、原子核を原点とする座標系での荷電粒子の運動方程式は

md2r/dt2 = kqZe2r/r3・・・(1)

（k:クーロン力の比例定数、r:荷電粒子の位置ベクトル）

と表せます。（粒子間にはクーロン力のみが働くものとし、また原子核の反跳は無視）

ここで、k,m,eの値をそのまま代入するのは困難なので、数値計算においては 規格化と呼ばれる作業をおこないます。
T=t/t0、R=r/r0として(1)式に代入すると

d2R/dT2 = q(kZe2t20/mr30) R/R3・・・(2)

となります。この計算では、t0=(kZe2/mr30)-1/2と なるようにしてみましょう。 荷電粒子の位置、および速度のx成分、y成分が満たす微分方程式を書き下せば

dVx/dT = qX/R3 ・・・(3)
dX / dT = VX ・・・(4)
dVy/dT = qY/R3 ・・・(5)
dY / dT = Vy ・・・(6)
（ただしR=(X2+Y2)1/2・・・(7))

となります。


初期座標(X0,Y0)、初速度V0、 電荷qを与えたとき、荷電粒子の軌道を計算するプログラムは以下のようになります。 ここでは保存量であるエネルギーと角運動量の誤差も計算しています。

E = V2/2 + q/R ・・・(8)
L=| R×V|・・・(9)


program main
  implicit none
  
  integer(kind=4),parameter:: neq=4              ←解く方程式の数

  integer(kind=4),parameter:: nmax=10000

  
  real(kind=8),parameter:: eps=1.0d-10
  
  external ex2
  
  common q                                        ←共通のqを用いる
  
  real(kind=8),dimension(neq):: y0,yn
  real(kind=8) x0,xn,dx                          
  real(kind=8) x_0,y_0,r0                         ←新たな変数の定義
  real(kind=8) q,v0,r,E0,E,L0,L                   ←新たな変数の定義
  
  real(kind=8),dimension(neq,12):: work
  
  integer(kind=4) ierr,i,n
  
  character(len=12) fname
  
  write(*,'(a,$)') ' q,v0 ? '
  read(*,*) q,v0                                   ←電荷、初速度を入力
  
  write(*,'(a,$)') ' x0,y0 ? '
  read(*,*) x_0,y_0                                ←初期座標を入力
  
  
  write(*,'(a,$)') 'file name ?'
  read(*,*) fname                            ←データを書き込むファイル名を入力
  
  open(10,file=fname)
  
  n=0
  x0=0.0d0
  dx=0.01d0
  
  y0(1) = v0                                       ←時刻0でのVx    
  y0(2) = x_0                                      ←時刻0でのX    
  y0(3) = 0.                                       ←時刻0でのVy    
  y0(4) = y_0                                      ←時刻0でのY    
  
  r0 = sqrt( x_0*x_0 +  y_0*y_0 )                  ←時刻0でのRを計算(7)式    
  E0 = 0.5*v0*v0 + q/r0                            ←時刻0でのEを計算(8)式    
  L0 = y_0*v0                                      ←時刻0でのLを計算(9)式    
  
  write(10,'(i5,1x,f7.3,6d23.15)') n,x0,y0(1),y0(2),y0(3),y0(4),0.0,0.0
  
  do n=1,nmax
     xn=x0+dx
     call runge(neq,ex2,x0,xn,y0,yn,eps,ierr,work)
     x0=xn
     
     do i=1,neq
        y0(i)=yn(i)
     end do
     
     r = sqrt( y0(2)*y0(2) +  y0(4)*y0(4) )         ←各々の時刻でのRを計算    
     E = 0.5*( y0(1)*y0(1) +  y0(3)*y0(3) ) + q/r   ←各々の時刻でのEを計算    
     L = abs(y0(2)*y0(3) -  y0(4)*y0(1))            ←各々の時刻でのLを計算    
     
     if (abs(r)>abs(r0))exit                        ←原子核から最初の距離よりも離れれば計算終了
     
     
     write(10,'(i5,1x,f7.3,6d23.15)') n,x0,y0(1),y0(2),y0(3),y0(4),abs((E-E0)/E0),abs((L-L0)/L0)                                  
                                                    ←E,Lは相対誤差を出力
  end do
  
  close(10)
  
end program main

real(kind=8) function ex2(i,x,y)
  implicit none
  
  integer(kind=4),parameter:: neq=4                  ←解く方程式の数
  
  integer(kind=4) i

  common q                                            ←共通のqを用いる

  real(kind=8) x
  real(kind=8),dimension(neq):: y
  real(kind=8) q                                     ←新たな変数の定義


  if (i==1) then
     ex2= q*y(2)/((y(2)*y(2)+y(4)*y(4))*sqrt(y(2)*y(2)+y(4)*y(4)))   ←(3)式
  else if (i==2) then
     ex2= y(1)                                                       ←(4)式
  else if (i==3) then
     ex2 = q*y(4)/((y(2)*y(2)+y(4)*y(4))*sqrt(y(2)*y(2)+y(4)*y(4)))  ←(5)式
  else if (i==4) then
     ex2 = y(3)                                                      ←(6)式
  else
     ex2=0.0d0
  end if
  
  return
end function ex2


プログラム中では、y(1)→Vx、y(2)→X、y(3)→Vy、y(4)→Y に対応しています。


計算例1：V0= 1.0、X0= -10.0でY0を変えた時の軌道と速度の変化

q= 1.0の時






q=-1.0の時






計算例2：V0= 1.0、X0= -10.0、Y0= 0.5でqを変えた時の軌道と速度の変化






計算例3：q= 1.0、V0= 1.0、X0= -10.0、Y0= 0.5のときのエネルギーと角運動量の誤差