module globals
	implicit none
	real(8) ,save :: x=0.d0,y=0.d0,u=0.d0,v=0.d0 !parameters of electron
	real(8) ,parameter :: e=1.602176d-19	!electron's charge(C)
	real(8) ,parameter :: me=9.109381d-31	!electron's mass(kg)
	real(8) ,parameter :: pi=2*acos(0.d0)
	real(8) ,parameter :: w1=29.d0	!長い底
	real(8) ,parameter :: w2=26.5d0	!短い底
	real(8) ,parameter :: h=9.1146859d0		!高さ
	real(8) ,parameter :: t=3.0d0		!磁極間の距離
 	real(8) ,parameter :: d=2.d-1	!width of detector
	real(8) ,parameter :: L1=80.d0	!distance between target and detector
	real(8) ,parameter :: L2=50.d0	!distance between magnet and detector
	real(8) ,parameter :: Ebeam=90.d0	!90MeV
	real(8) ,parameter :: Eelec=0.511d0	!0.511MeV
	real(8) ,parameter :: gamma=(Ebeam/Eelec)
	real(8) ,parameter :: LightV=2.99792458d10*sqrt(1.d0-1.d0/gamma**2)	!velocity of electron(cm/sec)
	!磁石の形を書くための変数
	real(8) ,parameter :: dw=w1-w2,theta=atan(dw/(2.0d0*h))&
		,sya=sqrt((dw/2.0d0)**2+h**2),endx1=w2/sya*h&
		,endy1=-w2/sya*dw/2.0d0,endx2=w1/sya*h&
		,endy2=-w1/sya*dw/2.0d0+sya
end module globals

module subprogs
	use globals
	implicit none
contains
	!磁石内の磁場(T=Weber/m^2=10^4Gauss)
	function B(x1,y1) result(z)	
		real(8) x1,y1,z
		if(y1>=tan(pi/2.d0-2*theta)*(x1-endx1)+endy1) then 
			z=0.63d0
		else
			z=0.d0
		endif
	endfunction
	function f(v1,b) result(z)
		real(8) v1,b,z
		z=e*v1*b/(me*gamma)
	endfunction
	function g(u1,b) result(z)
		real(8) u1,b,z
		z=-e*u1*b/(me*gamma)
	endfunction
	subroutine magnetform	!plot magnetform
		open(30,file='e_orbit_magnetform.dat')
		write(30,*) 0,0
		write(30,*) 0,sya
		write(30,*) endx2,endy2
		write(30,*) endx1,endy1
		write(30,*) 0,0
		write(30,*) ''
		write(30,*) 0,t/cos(theta)
		write(30,*) (t*dw/h+w2)/sya*h,&
				-(t*dw/h+w2)/sya*dw/2.0d0+t/cos(theta)
		write(30,*) ''
		write(30,*) 0,sya-t/cos(theta)
		write(30,*) ((h-t)*dw/h+w2)/sya*h,&
			-((h-t)*dw/h+w2)/sya*dw/2.0d0+sya-t/cos(theta)
		close(30)
	endsubroutine
	subroutine RK	!calculate electron's motion
	integer ::i,N=1000
	real(8) k1,k2,k3,k4,l1,l2,l3,l4,&
		m1,m2,m3,m4,n1,n2,n3,n4  !parameters for RK
	real(8) dt,tmax
	tmax=3.d-9	!3nsec
	dt=tmax/N
	write(10,*) x,y
	do i=1,N
		k1=dt*u
		l1=dt*f(v,B(x,y))
		m1=dt*v
		n1=dt*g(u,B(x,y))
		k2=dt*(u+l1/2.d0)
		l2=dt*f(v+n1/2.d0,B(x+k1/2.d0,y+m1/2.d0))
		m2=dt*(v+n1/2.d0)
		n2=dt*g(u+l1/2.d0,B(x+k1/2.d0,y+m1/2.d0))
		k3=dt*(u+l2/2.d0)
		l3=dt*f(v+n2/2.d0,B(x+k2/2.d0,y+m2/2.d0))
		m3=dt*(v+n2/2.d0)
		n3=dt*g(u+l2/2.d0,B(x+k2/2.d0,y+m2/2.d0))
		k4=dt*(u+l3)
		l4=dt*f(v+n3,B(x+k3,y+m3))
		m4=dt*(v+n3)
		n4=dt*g(u+l3,B(x+k3,y+m3))

		x=x+1/6.d0*(k1+2*k2+2*k3+k4)
		u=u+1/6.d0*(l1+2*l2+2*l3+l4)
		y=y+1/6.d0*(m1+2*m2+2*m3+m4)
		v=v+1/6.d0*(n1+2*n2+2*n3+n4)
		if(mod(i,10)==0) write(10,*) x,y
	enddo
	write(10,*) ''
	endsubroutine
end module subprogs

program e_orbit_rk
	use globals
	use subprogs
	implicit none
	integer i,j
	real(8) a,amin,amax
	real(8) in(0:6)
	open(10,file="e_orbit_rk.dat")
	call magnetform
	in(0)=t/cos(theta)	!入射位置
	do i=1,6
		in(i)=in(i-1)+d
	enddo

	amin=tan(atan(-2.d0*d/L1)+atan(0.13327))	!入射の傾き
	amax=tan(atan(2.d0*d/L1)+atan(0.13327))
	
	do i=0,0
		do j=0,1
			if(j==0) a=amin
			if(j==1) a=amax
			x=0.d0
			y=in(i)
			u=LightV*cos(atan(a))
			v=LightV*sin(atan(a))
			call RK
		enddo
		write(*,*) "i=",i
	enddo
	close(10)
endprogram