単純ラーマー運動の数値計算

fortranで電磁場中の荷電粒子(プロトン)の軌道を計算. 今回は単純なz方向磁場中のラーマー運動を計算してみる.

方程式

dv/dt = (q/m) (E +v×B)

設定

• z軸に 1T の磁場
• 初速度は vx=vy=vz=50m/s
• 時間幅 dt は結構小さく 1e-12 秒

ソースコード

Leap-Frog法で実装

program orbit
implicit none

real(8), parameter :: dt = 1d-12
integer, parameter :: maxstep = 3000000

integer :: i, j, k

real(8) :: r(3) = [0.d0, 0.d0, 0.d0]
real(8) :: v(3) = [5.d1, 5.d1, 5.d1], a(3), t=0.d0
real(8) :: Bx=0.d0, By=0.d0, Bz=1.d0
real(8) :: Ex=0.d0, Ey=0.d0, Ez=0.d0

! calculate trajectry of proton
do i = 1, maxstep
    if(i==1) write(101,'(7ES12.3)') t, r, v
    call leapflog   ! update r & v
    if( mod(i,500)==0 ) write(101,'(7ES12.3)') t, r, v
end do

contains
! ================================================== !
subroutine leapflog
    call accl(v,a)       ! update of accel
    v = v + a * dt*0.5d0 ! half step update of velocity
    r = r + v * dt       ! full step update of position
    call accl(v,a)       ! update of accel
    v = v + a * dt*0.5d0 ! half step update of velocity
    t = t + dt           ! time update
end subroutine leapflog
! ================================================== !
subroutine accl(v_in, a_out)
    implicit none
    real(8), parameter :: k = 1.602d-19 / 1.66054d-27
    real(8), intent(in)  :: v_in(3)
    real(8), intent(out) :: a_out(3)

    a_out(1) = (Ex + v_in(2)*Bz - v_in(3)*By)
    a_out(2) = (Ey + v_in(3)*Bx - v_in(1)*Bz)
    a_out(3) = (Ez + v_in(1)*By - v_in(2)*Bx)
    a_out = a_out * k
end subroutine accl
! ================================================== !
end program orbit

コンパイル・実行

gfortran orbit.f90 -O2
./a.out
gnuplot plot_orbit.gp # plot

軌道プロット

set term pdfcairo color enh font "Computer Modern,14"
set output "iorbit.pdf"
set nogrid
set size ratio 1
set colorsequence classic
set autoscale

set xlabel "time"
set ylabel "kinetic energy"
plot "fort.101" u 1:(sqrt($5**2+$6**2+$7**2)) w l notitle

set ylabel "position"
plot \
"fort.101" u 1:2 w l t "x",\
"fort.101" u 1:3 w l t "y",\
"fort.101" u 1:4 axis x1y2 w l t "z"

set ylabel "position"
plot "fort.101" u 1:3 w l t "z"

set size square
set xlabel "x"
set ylabel "y"
plot \
"fort.101" u 2:3 w l notitle,\
"fort.101" u 2:3 every ::0::0 w p pt 7 ps 1 notitle

set zlabel "z"
splot \
"fort.101" u 2:3:4 w l notitle,\
"fort.101" u 2:3:4 every ::0::0 w p pt 7 ps 1 notitle

きれいな円軌道を描く. 半径は r = mv_perp / qB で大体同じくらいの値が得られた.