電子の運動シミュレーション

1. 運動方程式

 電子の進行方向(z方向)での運動の様子についてのみ考えると、電子の運動方程式は、

である。ここでは、電場は空間成分と時間成分に分離できるものとする。さらに、独立変数である時間tをz座標で置き換える。例えば微分は、

となる。また、粒子の状態パラメーターも、運動量Pや時間tの代わりに、

  運動エネルギー (W)

  位相      (φ)

であらわす。相互の関係式は、

で与えられる。

 以上より、解くべき運動方程式は、zを独立変数とした以下の連立方程式である。

●(問1)この式(4)を導きなさい。

2. シミュレーション

 粒子シミュレーションでは、まず、粒子の1個について、状態を表すパラメータ(例えば、 W、φ)を考え、粒子がz方向に進行するにつれてそのパラメータが、どの様に変化するかを計算する。つまり、Wとφの初期値W0とφ0を決め、後は、式(4)を積分していけばよい。

単純にはこの積分は以下のように計算できる。ある位置zから微少距離z→(z+δz)進むする時、このWとφの変位量は以下の式で与えらるので、

この変化量を加速器の全長にわたって、繰り返し足していけば、最終的な値が得られる。

●(問2)E(z)は一定として、式(6)を計算するプログラムを作りなさい。

3. 粒子分布

 実験で測定することができるのは、1個の電子の運動ではなく、 電子ビームの分布であるので、これと比較するためには、 シミュレーションでも分布を計算する必要がある。 そこで、複数(N個)の粒子について、初期値を適当に変えて上と同様に計算し、 エネルギーや位相に関して統計を取る必要がある。

    ↓ ←運動方程式の積分

 

 このままでは、分布の様子がわかりにくいので、分布関数を計算する。 例えば、エネルギーの分布関数f(W)は、エネルギーがW±δWの範囲にある粒子の数を表す。また、分布関数は2次元でも考えることができ、f(W, φ) とすると、これはエネルギーがW±δW、位相がφ±δφの範囲にある粒子の数を表す。

●(問3) W-f平面で分布関数f(W)のグラフを書きなさい。


エネルギー分布計算プログラム例 FORTRAN, C, C with extra feature

	program rfgun
	parameter (clight=2.99792458E8,emass=511.)
	parameter(np=1000,nw=15,nph=18)
	real phi(np),w(np),charge(np),cur_w(0:nw),cur_ph(0:nph)
	integer nstep,i,j
	real pi,freq,omega,volt,gap,ez,deltaz,vz,gamma,dw
C Initialize
	pi=4*atan(1.0)
	freq=100.0E6
	omega=2.0*pi*freq
	volt=30.0
	nstep=1000
	gap=0.01
	ez=volt/gap
	deltaz=gap/nstep
	dw=volt/nw
C Initial distribution
	do 10 i=1,np
	  phi(i)=(pi/np)*(i-1)
	  w(i)=0.0
	  charge(i)=sin(phi(i))**1.5
 10	continue
C Integration
	do 20 i=1,nstep
	  do 30 j=1,np
	    w(j)=w(j)+ez*deltaz*sin(phi(j))
	    gamma=(w(j)+emass)/emass
	    if(gamma.LE.1.0) goto 30
	    vz=clight*sqrt(1.0-1.0/gamma/gamma)
	    phi(j)=phi(j)+omega*deltaz/vz
 30	  continue
 20	continue
C Energy bin
	do 40 i=0,nw
	  cur_w(i)=0.0
 40	continue
	do 50 j=1,np
	  ind=w(j)/dw
	  if(ind.ge.0 .and. ind.le.nw) then
	    cur_w(ind)=cur_w(ind)+charge(j)
	  endif
 50	continue
	write(*,'(''  Energy   Intensity''/(x,2f7.2))') 
     &	(i*dw,cur_w(i),i=1,nw)
C Phase bin
	do i=0,nph
	  cur_ph(i)=0.0
	enddo
	do j=1,np
	  ind=nint(mod(phi(j)/pi/2,1.0)*nph)
	  cur_ph(ind)=cur_ph(ind)+charge(j)
	enddo
	write(*,'(/''  phase   Intensity''/(x,2f7.2))') 
     &	(i*360./nph,cur_ph(i),i=1,nph)
	stop
	end