implicit real*8 (a-h,o-z)
real*8 y(6),work(226),m
integer iwork(5)
common qm,r_maa
external func
namelist /input/ q,m,r_maa,nt,dt,neqn,relerr,abserr,y

open(10,file='lorentz_prooton.in')
open(20,file='lorentz.dat')
read(10,input)

t=0.
qm=q/m
tout=0.
iflag=1

do i=1,nt
tout=t+dt

call ode(func,neqn,y,t,tout,relerr,abserr,iflag,work,iwork)
if(iflag /= 2)then
print *,iflag
iflag=2
endif

write(20,*) y(1:3)
enddo

stop
end

subroutine func(t,y,dy)
implicit real*8 (a-h,o-z)
real*8 y(6),dy(6),b(3),re(3)
common qm,r_maa
re(1)=y(1)/r_maa
re(2)=y(2)/r_maa
re(3)=y(3)/r_maa
call emag(re,b)
dy(1)=y(4)
dy(2)=y(5)
dy(3)=y(6)
dy(4)=qm*(y(5)*b(3)-y(6)*b(2))
dy(5)=qm*(y(6)*b(1)-y(4)*b(3))
dy(6)=qm*(y(4)*b(2)-y(5)*b(1))
return
end

subroutine emag(r,b)
implicit real*8 (a-h,o-z)
real*8 r(3), b(3)
b0=3.07e-05
x=r(1)
y=r(2)
z=r(3)
rr=sqrt(x**2+y**2+z**2)
b(1)=-3.*b0*z*x/(rr**5)
b(2)=-3.*b0*z*y/(rr**5)
b(3)=-b0*(2*z**2-x**2-y**2)/(rr**5)
return
end