c Lorentz Force //  E = Exi + Ey j    B= B0 k  (units SI )
      real mass
      data c,q,mass/3.E8,1.602E-19,9.1e-31/
      data vx0,vy0,vz0/0.,-1.E7,0./
      beta(vx,vy,vz)=sqrt(vx**2+vy**2+vz**2)
      gamma(vx,vy,vz)=1./sqrt(1.-(beta(vx,vy,vz)/c)**2)
      dgdt(vx,vy,vz,vx0,vy0,vz0)=(gamma(vx,vy,vz)-gamma(vx0,vy0,vz0))/dt
      ex(x,y,z)=-1.e2
      ey(x,y,z)=1.e-12
      ez(x,y,z)=0.
      bx(x,y,z)=0.
      by(x,y,z)=0.
      bz(x,y,z)=1.e-4
      Eabs(x,y,z)=sqrt( ex(x,y,z)**2 + ey(x,y,z)**2 +ez(x,y,z)**2)
      ax(x,y,z,vy,vz)= (q/(mass*gamma(vx,vy,vz)))*
     $ (ex(x,y,z)+vy*Bz(x,y,z)
     $ -vz*by(x,y,z))-vx* dgdt(vx,vy,vz,vx0,vy0,vz0)/gamma(vx,vy,vz)
      ay(x,y,z,vx,vz)= (q/(mass*gamma(vx,vy,vz)))*
     $ (ey(x,y,z)-vx*bz(x,y,z)
     $ +vz*bx(x,y,z)) -vy* dgdt(vx,vy,vz,vx0,vy0,vz0)/gamma(vx,vy,vz)
      az(x,y,z,vx,vy)= (q/(mass*gamma(vx,vy,vz)))*
     $ (ez(x,y,z)-vy*bx(x,y,z)
     $ +vx*by(x,y,z)) -vz* dgdt(vx,vy,vz,vx0,vy0,vz0)/gamma(vx,vy,vz)
      tscale1=mass/(q*Bz(0.,0.,0.))
      tscale2=c*mass/(q*Eabs(0.,0.,0.))
      tscale=amin1(tscale1,tscale2)
      print*,'tscale=',tscale
      nstep=4000
c      dt=6.*tscale/float(nstep)
      dt=5.e-7/float(nstep)
      x0=1.
      y0=0.
      z0=0.
      x1=vx0*dt+x0
      y1=vy0*dt+y0
      z1=vz0*dt+z0
      kprint=int(float(nstep)/60.)
      kount=kprint
      do 10 i=2,nstep
      t=dt*float(i)
      vx=(x1-x0)/dt
      vy=(y1-y0)/dt
      vz=(z1-z0)/dt
c      if(gamma(vx,vy,vy).ge.1.1)goto  200
      x2=2.*x1-x0+dt**2*ax(x,y,z,vy,vz)
      y2=2.*y1-y0+dt**2*ay(x,y,z,vx,vz)
      z2=2.*z1-z0+dt**2*az(x,y,z,vx,vy)
      v=sqrt( ((x2-x1)/dt)**2 +  ((y2-y1)/dt)**2 +vz**2 )
      if(i.eq.kount)then
      kount=kount+kprint
      print 100 ,t,x2,y2,z2
      endif
      x0=x1
      x1=x2
      y0=y1
      y1=y2
      z0=z1
      z1=z2
      vx0=vx
      vy0=vy
      vz0=vz
10    continue
100   format(2x,'t,x,y,v=',4(3x,e11.3))
200   stop
      end