E×Bドリフト

fortranで電磁場中の荷電粒子(プロトン)の軌道を計算.

今回はZ方向の磁場 Bz と Y方向の電場 Ey 中にあるイオンの軌道を計算してみる.

もちろん、E×Bドリフトによって +x 方向にドリフトする。

方程式

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

設定

• z軸に 1T の磁場
• 初速度は vx=vy=vz=50m/s
• 時間幅 dt は結構小さく 1e-12 秒
• 時間・空間均一な Ey = 0.3 V/m があるとする

ソースコード

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]
real(8) :: a(3), t=0.d0
real(8) :: Bx=0.d0, By=0.d0, Bz=1.d0
real(8) :: Ex=0.d0, Ey=3.d-1, Ez=0.d0

do i = 1, maxstep
    if(i==1) write(101,'(7ES12.3)') t, r, v
    call leapflog
    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 :: q = 1.602d-19
    real(8), parameter :: m = 1.66054d-27
    real(8), intent(in)  :: v_in(3)
    real(8), intent(out) :: a_out(3)

    a_out(1) = q/m * (Ex + v_in(2)*Bz - v_in(3)*By)
    a_out(2) = q/m * (Ey + v_in(3)*Bx - v_in(1)*Bz)
    a_out(3) = q/m * (Ez + v_in(1)*By - v_in(2)*Bx)
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*1e6):($3*1e6) w l notitle,\
"fort.101" u ($2*1e6):($3*1e6) 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

いい感じに E×Bドリフトを計算できた。