Lorentz_kick Subroutine

  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