一様な電磁場での荷電粒子の挙動のプログラミング
一様な電磁場での電子・陽子の動きをfortran77でプログラミングしようと思ったのですが、なかなかうまくいきません。数値設定も全て自分でやっているのであまり自信はないのですが・・・
何かおかしなところや直した方がいいところがあれば助言の方よろしくお願いします。。。
implicit none
real*8 me,re(3),ve(3),E(3),B(3),qe,Te,dte
real*8 mp,rp(3),vp(3),Tp,dtp,time,qp
c 初期値の設定
me=9.1d-31 ! 電子の質量
mp=1.67d-27 ! 陽子の質量
qe=-1.6d-19 ! 電子の電荷
qp=1.6d-19 ! 陽子の電荷
re(1)=0.0
re(2)=0.0 ! 電子の位置
re(3)=0.0
ve(1)=1.0d4
ve(2)=0.0 ! 電子の速度
ve(3)=0.0
rp(1)=0.0
rp(2)=0.0 ! 陽子の位置
rp(3)=0.0
vp(1)=1.0d4
vp(2)=0.0 ! 陽子の速度
vp(3)=0.0
E(1)=0.0
E(2)=1.0d-3 ! 電場
E(3)=0.0
B(1)=0.0
B(2)=0.0 ! 磁場
B(3)=1.0d-5
dte=-1.0d-5
dtp=1.0d-4
open(unit=10,file='densi.dat',status='unknown')
do time=0.0,0.1,dte
write(10,*) time,re,ve
call lorentz(re,ve,E,B,qe,me,dte,Te)
enddo
close(10)
open(unit=20,file='yousi.dat',status='unknown')
do time=0.0,0.1,dtp
write(20,*) time,rp,vp
call lorentz(rp,vp,E,B,qp,mp,dtp,Tp)
enddo
close(20)
stop
end
c ローレンツの計算
subroutine lorentz(r,v,E,B,q,m,dt,T)
real*8 r(3),v(3),E(3),B(3),q,m,dt,T
v(1)=v(1)+q*(E(1)+(v(2)*B(3)-v(3)*B(2)))*dt/m
v(2)=v(2)+q*(E(2)+(v(3)*B(1)-v(1)*B(3)))*dt/m
v(3)=v(3)+q*(E(3)+(v(1)*B(2)-v(2)*B(1)))*dt/m
r(1)=r(1)+v(1)*dt
r(2)=r(2)+v(2)*dt
r(3)=r(3)+v(3)*dt
T=(m*((v(1)**2)+(v(2)**2)+(v(3)**2)))/2
return
end
補足
ーがつくときはどのような座標系だと思いますか？