!> Particle mover with Newtonian/relativistic Boris scheme for Lorentz dynamics
!> By Jannis Teunissen, Bart Ripperda, Oliver Porth, and Fabio Bacchini (2016-2020)
module mod_particle_lorentz
  use mod_particle_base

  private

  public :: Lorentz_init
  public :: Lorentz_create_particles
  integer, parameter :: Boris=1, Vay=2, HC=3, LM=4

  ! Variables
  public :: bp, ep, vp, jp, rhop

contains

  subroutine Lorentz_init()
    use mod_global_parameters
    integer :: idir, nwx

    if(physics_type/='mhd') call mpistop(&
       "Lorentz particles need magnetic field!")
    if(ndir/=3) call mpistop("Lorentz particles need ndir=3!")

    ! The first 6 gridvars are always B and E
    ngridvars = ndir*2
    nwx=0
    ! density
    if(particles_etah>0) nwx = 1

    allocate(bp(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      bp(idir) = nwx
    end do
    allocate(ep(ndir)) ! electric field
    do idir = 1, ndir
      nwx = nwx + 1
      ep(idir) = nwx
    end do
    allocate(vp(ndir)) ! fluid velocity
    do idir = 1, ndir
      nwx = nwx + 1
      vp(idir) = nwx
    end do
    allocate(jp(ndir)) ! current
    do idir = 1, ndir
      nwx = nwx + 1
      jp(idir) = nwx
    end do 
    nwx = nwx + 1 ! density
    rhop = nwx

    ngridvars=nwx

    particles_fill_gridvars => lorentz_fill_gridvars

    if (associated(particles_define_additional_gridvars)) then
      call particles_define_additional_gridvars(ngridvars)
    end if

    select case(integrator_type_particles)
    case('Boris','boris')
      integrator = Boris
    case('Vay','vay')
      integrator = Vay
    case('HC','hc','higueracary')
      integrator = HC
    case('LM','lm','lapentamarkidis')
      integrator = LM
    case default
      integrator = Boris
    end select

    particles_integrate => Lorentz_integrate_particles

  end subroutine Lorentz_init

  subroutine Lorentz_create_particles()

    use mod_global_parameters
    use mod_usr_methods, only: usr_create_particles, usr_update_payload,&
        usr_check_particle

    integer          :: n, idir, igrid, ipe_particle, nparticles_local
    double precision :: lfac
    double precision :: x(3, num_particles)
    double precision :: v(3, num_particles)
    double precision :: q(num_particles)
    double precision :: m(num_particles)
    double precision :: rrd(num_particles,ndir)
    double precision :: defpayload(ndefpayload)
    double precision :: usrpayload(nusrpayload)
    logical          :: follow(num_particles), check

    follow = .false.
    x      = 0.0d0

    if (mype==0) then
      if (.not. associated(usr_create_particles)) then
        ! Randomly distributed
        do idir=1,ndir
          do n = 1, num_particles
            rrd(n,idir) = rng%unif_01()
          end do
        end do
        do n=1, num_particles
          x(1,n) = xprobmin1 + rrd(n+1,1) * (xprobmax1 - xprobmin1)
          x(2,n) = xprobmin2 + rrd(n+1,2) * (xprobmax2 - xprobmin2)
          x(3,n) = xprobmin3 + rrd(n+1,3) * (xprobmax3 - xprobmin3)
        end do
      else
        call usr_create_particles(num_particles, x, v, q, m, follow)
      end if
    end if

    call MPI_BCAST(x,3*num_particles,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)
    call MPI_BCAST(v,3*num_particles,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)
    call MPI_BCAST(q,num_particles,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)
    call MPI_BCAST(m,num_particles,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)
    call MPI_BCAST(follow,num_particles,MPI_LOGICAL,0,icomm,ierrmpi)

    nparticles_local = 0

    ! Find ipe and igrid responsible for particle
    do n = 1, num_particles
      call find_particle_ipe(x(:, n),igrid,ipe_particle)

      particle(n)%igrid = igrid
      particle(n)%ipe   = ipe_particle

      if (ipe_particle == mype) then
        check = .true.

        ! Check for user-defined modifications or rejection conditions
        if (associated(usr_check_particle)) call usr_check_particle(igrid, x(:,&
           n), v(:,n), q(n), m(n), follow(n), check)
        if (check) then
          call push_particle_into_particles_on_mype(n)
        else
          cycle
        end if

        nparticles_local = nparticles_local + 1

        call get_lfac_from_velocity(v(:,n), lfac)

        allocate(particle(n)%self)
        particle(n)%self%x      = x(:,n)
        particle(n)%self%u      = v(:,n) * lfac
        particle(n)%self%q      = q(n)
        particle(n)%self%m      = m(n)
        particle(n)%self%follow = follow(n)
        particle(n)%self%index  = n
        particle(n)%self%time   = global_time
        particle(n)%self%dt     = 0.0d0

        ! initialise payloads for Lorentz module
        allocate(particle(n)%payload(npayload))
        call Lorentz_update_payload(igrid,x(:,n),v(:,n)*lfac,q(n),m(n),&
           defpayload,ndefpayload,0.d0)
        particle(n)%payload(1:ndefpayload) = defpayload
        if (associated(usr_update_payload)) then
          call usr_update_payload(igrid,x(:,n),v(:,n)*lfac,q(n),m(n),&
             usrpayload,nusrpayload,0.d0)
          particle(n)%payload(ndefpayload+1:npayload) = usrpayload
        end if
      end if
    end do

    call MPI_ALLREDUCE(nparticles_local,nparticles,1,MPI_INTEGER,MPI_SUM,icomm,&
       ierrmpi)

  end subroutine Lorentz_create_particles

  subroutine lorentz_fill_gridvars
    use mod_global_parameters
    use mod_usr_methods, only: usr_particle_fields
    use mod_geometry

    integer                                   :: igrid, iigrid, idir, idim
    double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
       1:ndir) :: beta
    double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
       1:nw)   :: w,wold
    double precision                          :: current(ixGlo1:ixGhi1,&
       ixGlo2:ixGhi2,ixGlo3:ixGhi3,7-2*ndir:3)
    integer                                   :: idirmin
    double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
       1:ndir) :: vE, bhat
    double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
       ixGlo3:ixGhi3)        :: kappa, kappa_B, absB, tmp

    ! Fill electromagnetic quantities
    call fill_gridvars_default()

  end subroutine lorentz_fill_gridvars

  !> Relativistic particle integrator
  subroutine Lorentz_integrate_particles(end_time)
    use mod_global_parameters
    use mod_geometry
    use mod_usr_methods, only: usr_create_particles, usr_update_payload
    double precision, intent(in)      :: end_time
    double precision                  :: defpayload(ndefpayload)
    double precision                  :: usrpayload(nusrpayload)
    integer                           :: ipart, iipart, igrid
    double precision                  :: lfac, q, m, dt_p
    double precision                  :: xp(ndir), xpm(ndir), xpc(ndir),&
        xpcm(ndir)
    double precision                  :: up(ndir), upc(ndir), tp
    double precision, dimension(ndir) :: b, e, bc, ec, vfluid, current
    double precision                  :: rho, rhoold, td

    do iipart=1,nparticles_active_on_mype
      ipart = particles_active_on_mype(iipart);
      q     = particle(ipart)%self%q
      m     = particle(ipart)%self%m
      igrid = particle(ipart)%igrid
      xp    = particle(ipart)%self%x
      up    = particle(ipart)%self%u
      tp    = particle(ipart)%self%time   

      dt_p = Lorentz_get_particle_dt(particle(ipart), end_time)
      particle(ipart)%self%dt = dt_p

      ! Push particle over first half time step
      ! all pushes done in Cartesian coordinates
      call partcoord_to_cartesian(xp,xpc)
      call partvec_to_cartesian(xp,up,upc)
      call get_lfac(upc,lfac)
      xpcm = xpc + 0.5d0 * dt_p * upc/lfac
      call partcoord_from_cartesian(xpm,xpcm)
      ! Fix xp if the 0,2*pi boundary was crossed in cylindrical/spherical coords
      call fix_phi_crossing(xpm,igrid)
      
      ! Get E, B at n+1/2 position
      call get_bfield(igrid, xpm, tp+dt_p/2.d0, b)
      call get_efield(igrid, xpm, tp+dt_p/2.d0, e)

!      call get_vec(vp, igrid, xpm, tp+dt_p/2.d0, vfluid)
!      call get_vec(jp, igrid, xpm, tp+dt_p/2.d0, current)
!      select case (coordinate)
!      case (Cartesian,Cartesian_stretched,spherical)
!        e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
!        e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
!        e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)
!      case (cylindrical)
!        e(r_) = -vfluid(phi_)*b(z_)+vfluid(z_)*b(phi_) + particles_eta*current(r_)
!        e(phi_) = vfluid(r_)*b(z_)-vfluid(z_)*b(r_) + particles_eta*current(phi_)
!        e(z_) = -vfluid(r_)*b(phi_)+vfluid(phi_)*b(r_) + particles_eta*current(z_)
!      end select
!      if (particles_etah > zero) then
!        call interpolate_var(igrid,ixG^LL,ixM^LL,ps(igrid)%w(ixG^T,1),ps(igrid)%x,xpm,rho)
!        if (time_advance) then
!          td = (tp+dt_p/2.d0 - global_time) / dt
!          call interpolate_var(igrid,ixG^LL,ixM^LL,pso(igrid)%w(ixG^T,1),ps(igrid)%x,xpm,rhoold)
!          rho = rhoold * (1.d0-td) + rho * td
!        end if
!        select case (coordinate)
!        case (Cartesian,Cartesian_stretched,spherical)
!          e(1) = e(1) + particles_etah/rho * (current(2)*b(3) - current(3)*b(2))
!          e(2) = e(2) + particles_etah/rho * (-current(1)*b(3) + current(3)*b(1))
!          e(3) = e(3) + particles_etah/rho * (current(1)*b(2) - current(2)*b(1))
!        case (cylindrical)
!          e(r_) = e(r_) + particles_etah/rho * (current(phi_)*b(z_) - current(z_)*b(phi_))
!          e(phi_) = e(phi_) + particles_etah/rho * (-current(r_)*b(z_) + current(z_)*b(r_))
!          e(z_) = e(z_) + particles_etah/rho * (current(r_)*b(phi_) - current(phi_)*b(r_))
!        end select
!      end if

      ! Convert fields to Cartesian frame
      call partvec_to_cartesian(xpm,b,bc)
      call partvec_to_cartesian(xpm,e,ec)

      ! 'Kick' particle (update velocity) based on the chosen integrator
      call Lorentz_kick(upc,ec,bc,q,m,dt_p)

      ! Push particle over second half time step
      ! all pushes done in Cartesian coordinates
      call get_lfac(upc,lfac)
      xpc = xpcm + 0.5d0 * dt_p * upc/lfac
      call partcoord_from_cartesian(xp,xpc)
      ! Fix xp if the 0,2*pi boundary was crossed in cylindrical/spherical coords
      call fix_phi_crossing(xp,igrid)
      call partvec_from_cartesian(xp,up,upc)

      ! Store updated x,u
      particle(ipart)%self%x = xp
      particle(ipart)%self%u = up

      ! Time update
      tp = tp + dt_p
      particle(ipart)%self%time = tp

      ! Update payload
      call Lorentz_update_payload(igrid,xp,up,q,m,defpayload,ndefpayload,tp)
      particle(ipart)%payload(1:ndefpayload) = defpayload
      if (associated(usr_update_payload)) then
        call usr_update_payload(igrid,xp,up,q,m,usrpayload,nusrpayload,tp)
        particle(ipart)%payload(ndefpayload+1:npayload) = usrpayload
      end if

    end do ! ipart loop

  end subroutine Lorentz_integrate_particles

  !> Momentum update subroutine for full Lorentz dynamics
  subroutine Lorentz_kick(upart,e,b,q,m,dtp)
    use mod_global_parameters
    use mod_geometry
    double precision, intent(in)      :: e(ndir), b(ndir), q, m, dtp
    double precision, intent(inout)   :: upart(ndir)
    double precision                  :: lfac, cosphi, sinphi, phi1, phi2, r,&
        re, sigma, td
    double precision, dimension(ndir) :: emom, uprime, tau, s, tmp, uplus,&
        xcart1, xcart2, ucart2, radmom
    double precision, dimension(ndir) :: upartk, vbar, Fk, C1, C2, dupartk
    double precision                  :: abserr, tol, lfack, J11, J12, J13,&
        J21, J22, J23, J31, J32, J33
    double precision                  :: iJ11, iJ12, iJ13, iJ21, iJ22, iJ23,&
        iJ31, iJ32, iJ33, Det
    integer                           :: nk, nkmax

    ! Perform momentum update based on the chosen integrator
    select case(integrator)

    ! Boris integrator (works in Cartesian and cylindrical)
    case(Boris)
      ! Momentum update
      emom = q * e * dtp /(2.0d0 * m)
      uprime = upart + emom
      !!!!!! TODO: Adjust and document losses
!      if (losses) then
!        call get_lfac(particle(ipart)%self%u,lfac)
!        re = abs(q)**2 / (m * const_c**2)
!        call cross(upart,b,tmp)
!        radmom = - third * re**2 * lfac &
!             * ( sum((e(:)+tmp(:)/lfac)**2)  &
!             -  (sum(e(:)*upart(:))/lfac)**2 ) &
!             * particle(ipart)%self%u / m / const_c * dt_p
!        uprime = uprime + radmom
!      end if

      call get_lfac(uprime,lfac)
      tau = q * b * dtp / (2.0d0 * lfac * m)
      s = 2.0d0 * tau / (1.0d0+sum(tau(:)**2))

      call cross(uprime,tau,tmp)
      call cross(uprime+tmp,s,tmp)
      uplus = uprime + tmp

      upart = uplus + emom
      !!!!!! TODO: Adjust and document losses
!      if(losses) then
!        call cross(uplus,b,tmp)
!        radmom = - third * re**2 * lfac &
!             * ( sum((e(:)+tmp(:)/lfac)**2)  &
!             -  (sum(e(:)*uplus(:))/lfac)**2 ) &
!             * uplus / m / const_c * dt_p
!        upart = upart + radmom
!      end if

    ! Vay integrator
    case(Vay)
      call get_lfac(upart,lfac)
      emom = q * e * dtp /(2.0d0 * m)
      tau = q * b * dtp / (2.0d0 * m)

      call cross(upart,tau,tmp)
      uprime = upart + 2.d0*emom + tmp/lfac

      call get_lfac(uprime,lfac)
      sigma = lfac**2 - sum(tau(:)*tau(:))
      lfac = sqrt((sigma + sqrt(sigma**2 + 4.d0 * (sum(tau(:)*tau(:)) + &
         sum(uprime(:)*tau(:)/c_norm)**2))) / 2.d0)
      
      call cross(uprime,tau,tmp)
      upart = (uprime + sum(uprime(:)*tau(:))*tau/lfac**2 + tmp/lfac) / &
         (1.d0+sum(tau(:)*tau(:))/lfac**2)

    ! Higuera-Cary integrator
    case(HC)
      call get_lfac(upart,lfac)
      emom = q * e * dtp /(2.0d0 * m)
      tau = q * b * dtp / (2.0d0 * m)
      uprime = upart + emom

      call get_lfac(uprime,lfac)
      sigma = lfac**2 - sum(tau(:)*tau(:))
      lfac = sqrt((sigma + sqrt(sigma**2 + 4.d0 * (sum(tau(:)*tau(:)) + &
         sum(uprime(:)*tau(:)/c_norm)**2))) / 2.d0)
      
      call cross(uprime,tau,tmp)
      upart = (uprime + sum(uprime(:)*tau(:))*tau/lfac**2 + tmp/lfac) / &
         (1.d0+sum(tau(:)*tau(:))/lfac**2) + emom + tmp/lfac

    ! Lapenta-Markidis integrator
    case(LM)
      ! Initialise iteration quantities
      call get_lfac(upart,lfac)
      upartk = upart
      tau(:) = b(:)

      ! START OF THE NONLINEAR CYCLE
      abserr = 1.d0
      tol=1.d-14
      nkmax=10
      nk=0
      do while(abserr > tol .and. nk < nkmax)

        nk=nk+1

        call get_lfac(upartk,lfack)
        vbar = (upart + upartk) / (lfac + lfack)
        call cross(vbar,tau,tmp)

        ! Compute residual vector
        Fk = upartk - upart - q*dtp/m * (e + tmp)

        ! Compute auxiliary coefficients
        C1 = (lfack + lfac - upartk(1:ndim) / lfack / c_norm**2 * (upartk + &
           upart)) / (lfack + lfac)**2
        C2 = - upartk / lfack / c_norm**2 / (lfack + lfac)**2

        ! Compute Jacobian
        J11 = 1. - q*dtp/m * (C2(1) * (upartk(2) + upart(2)) * tau(3) - C2(1) &
           * (upartk(3) + upart(3)) * tau(2))
        J12 = - q*dtp/m * (C1(2) * tau(3) - C2(2) * (upartk(3) + upart(3)) * &
           tau(2))
        J13 = - q*dtp/m * (C2(3) * (upartk(2) + upart(2)) * tau(3) - C1(3) * &
           tau(2))
        J21 = - q*dtp/m * (- C1(1) * tau(3) + C2(1) * (upartk(3) + upart(3)) * &
           tau(1))
        J22 = 1. - q*dtp/m * (- C2(2) * (upartk(1) + upart(1)) * tau(3) + &
           C2(2) * (upartk(3) + upart(3)) * tau(1))
        J23 = - q*dtp/m * (- C2(3) * (upartk(1) + upart(1)) * tau(3) + C1(3) * &
           tau(1))
        J31 = - q*dtp/m * (C1(1) * tau(2) - C2(1) * (upartk(2) + upart(2)) * &
           tau(1))
        J32 = - q*dtp/m * (C2(2) * (upartk(1) + upart(1)) * tau(2) - C1(2) * &
           tau(1))
        J33 = 1. - q*dtp/m * (C2(3) * (upartk(1) + upart(1)) * tau(2) - C2(3) &
           * (upartk(2) + upart(2)) * tau(1))

        ! Compute inverse Jacobian
        Det = J11*J22*J33 + J21*J32*J13 + J31*J12*J23 - J11*J32*J23 - &
           J31*J22*J13 - J21*J12*J33
        iJ11 = (J22*J33 - J23*J32) / Det
        iJ12 = (J13*J32 - J12*J33) / Det
        iJ13 = (J12*J23 - J13*J22) / Det
        iJ21 = (J23*J31 - J21*J33) / Det
        iJ22 = (J11*J33 - J13*J31) / Det
        iJ23 = (J13*J21 - J11*J23) / Det
        iJ31 = (J21*J32 - J22*J31) / Det
        iJ32 = (J12*J31 - J11*J32) / Det
        iJ33 = (J11*J22 - J12*J21) / Det

        ! Compute new upartk = upartk - J^(-1) * F(upartk)
        dupartk(1) = - (iJ11 * Fk(1) + iJ12 * Fk(2) + iJ13 * Fk(3))
        dupartk(2) = - (iJ21 * Fk(1) + iJ22 * Fk(2) + iJ23 * Fk(3))
        dupartk(3) = - (iJ31 * Fk(1) + iJ32 * Fk(2) + iJ33 * Fk(3))

        ! Check convergence
        upartk=upartk+dupartk
        abserr=sqrt(sum(dupartk(:)*dupartk(:)))

      end do
      ! END OF THE NONLINEAR CYCLE

      ! Update velocity
      upart = upartk

    end select

  end subroutine Lorentz_kick

  !> Update payload subroutine
  subroutine Lorentz_update_payload(igrid,xpart,upart,qpart,mpart,mypayload,&
     mynpayload,particle_time)
    use mod_global_parameters
    use mod_geometry
    integer, intent(in)           :: igrid,mynpayload
    double precision, intent(in)  :: xpart(1:ndir),upart(1:ndir),qpart,mpart,&
       particle_time
    double precision, intent(out) :: mypayload(mynpayload)
    double precision              :: b(3), e(3), tmp(3), lfac, vfluid(3),&
        current(3), rho, rhoold, td

    call get_bfield(igrid, xpart, particle_time, b)
    call get_efield(igrid, xpart, particle_time, e)

!    call get_vec(bp, igrid, xpart,particle_time,b)
!    call get_vec(vp, igrid, xpart,particle_time,vfluid)
!    call get_vec(jp, igrid, xpart,particle_time,current)
!    select case (coordinate)
!    case (Cartesian,Cartesian_stretched,spherical)
!      e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
!      e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
!      e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)
!    case (cylindrical)
!      e(r_) = -vfluid(phi_)*b(z_)+vfluid(z_)*b(phi_) + particles_eta*current(r_)
!      e(phi_) = vfluid(r_)*b(z_)-vfluid(z_)*b(r_) + particles_eta*current(phi_)
!      e(z_) = -vfluid(r_)*b(phi_)+vfluid(phi_)*b(r_) + particles_eta*current(z_)
!    end select
!    if (particles_etah > zero) then
!      call interpolate_var(igrid,ixG^LL,ixM^LL,ps(igrid)%w(ixG^T,1),ps(igrid)%x,xpart,rho)
!      if (time_advance) then
!        td = (particle_time - global_time) / dt
!        call interpolate_var(igrid,ixG^LL,ixM^LL,pso(igrid)%w(ixG^T,1),ps(igrid)%x,xpart,rhoold)
!        rho = rhoold * (1.d0-td) + rho * td
!      end if
!      select case (coordinate)
!      case (Cartesian,Cartesian_stretched,spherical)
!        e(1) = e(1) + particles_etah/rho * (current(2)*b(3) - current(3)*b(2))
!        e(2) = e(2) + particles_etah/rho * (-current(1)*b(3) + current(3)*b(1))
!        e(3) = e(3) + particles_etah/rho * (current(1)*b(2) - current(2)*b(1))
!      case (cylindrical)
!        e(r_) = e(r_) + particles_etah/rho * (current(phi_)*b(z_) - current(z_)*b(phi_))
!        e(phi_) = e(phi_) + particles_etah/rho * (-current(r_)*b(z_) + current(z_)*b(r_))
!        e(z_) = e(z_) + particles_etah/rho * (current(r_)*b(phi_) - current(phi_)*b(r_))
!      end select
!    end if

    ! Payload update
    ! Lorentz factor
    if (mynpayload > 0) then
      call get_lfac(upart,lfac)     
      mypayload(1) = lfac
    end if

    ! current gyroradius
    if (mynpayload > 1) then
      call cross(upart,b,tmp)
      tmp = tmp / sqrt(sum(b(:)**2))
      mypayload(2) = mpart/abs(qpart)*sqrt(sum(tmp(:)**2)) / &
         sqrt(sum(b(:)**2))
    end if

    ! magnetic moment
    if (mynpayload > 2) then
      mypayload(3) = mpart*sum(tmp(:)**2)/(2.0d0*sqrt(sum(b(:)**2)))
    end if

    ! e.b
    if (mynpayload > 3) then
      mypayload(4) = sum(e(:)*b(:))
    end if

  end subroutine Lorentz_update_payload

  function Lorentz_get_particle_dt(partp, end_time) result(dt_p)
    use mod_global_parameters
    use mod_geometry
    type(particle_ptr), intent(in)   :: partp
    double precision, intent(in)     :: end_time
    double precision                 :: dt_p
    integer                          :: ipart, iipart, nout
    double precision,dimension(ndir) :: b,v
    double precision                 :: lfac,absb,dt_cfl
    double precision                 :: tout
    double precision, parameter      :: cfl = 0.5d0

    if (const_dt_particles > 0) then
      dt_p = const_dt_particles
      return
    end if

    call get_vec(bp, partp%igrid,partp%self%x,partp%self%time,b)
    absb = sqrt(sum(b(:)**2))
    call get_lfac(partp%self%u,lfac)

    ! CFL timestep
    ! make sure we step only one cell at a time:
    v(:) = abs(partp%self%u(:) / lfac)

!    ! convert to angular velocity:
!    if(coordinate ==cylindrical.and.phi_>0) then
!      v(phi_) = abs(v(phi_)/partp%self%x(r_))
!    end if

    dt_cfl = min(bigdouble, rnode(rpdx1_,partp%igrid)/max(v(1), smalldouble),&
       rnode(rpdx2_,partp%igrid)/max(v(2), smalldouble),rnode(rpdx3_,&
       partp%igrid)/max(v(3), smalldouble))

    if(coordinate ==cylindrical.and.phi_>0) then
      ! phi-momentum leads to radial velocity:
      if(phi_ .gt. ndim) dt_cfl = min(dt_cfl, sqrt(rnode(rpdx1_,&
         partp%igrid)/partp%self%x(r_)) / v(phi_))
      ! limit the delta phi of the orbit (just for aesthetic reasons):
      dt_cfl = min(dt_cfl,0.1d0/max(v(phi_), smalldouble))
      ! take some care at the axis:
      dt_cfl = min(dt_cfl,(partp%self%x(r_)+smalldouble)/max(v(r_),&
          smalldouble))
    end if

    dt_cfl = dt_cfl * cfl

    ! bound by gyro-rotation:
    dt_p = abs( dtheta * partp%self%m * lfac / (partp%self%q * absb) )

    dt_p = min(dt_p, dt_cfl)

    ! Make sure we don't advance beyond end_time
    call limit_dt_endtime(end_time - partp%self%time, dt_p)

  end function Lorentz_get_particle_dt

end module mod_particle_lorentz