MODULE modulePusher

  CONTAINS
    !Push neutral particles in cartesian coordinates
    PURE SUBROUTINE pushCartNeutral(part, tauIn)
      USE moduleSpecies
      IMPLICIT NONE

      TYPE(particle), INTENT(inout):: part
      REAL(8), INTENT(in):: tauIn

      part%r = part%r + part%v*tauIn

    END SUBROUTINE pushCartNeutral

    PURE SUBROUTINE pushCartElectrostatic(part, tauIn)
      USE moduleSPecies
      USE moduleMesh
      IMPLICIT NONE

      TYPE(particle), INTENT(inout):: part
      REAL(8), INTENT(in):: tauIn
      REAL(8):: qmEFt(1:3)

      !Get the electric field at particle position
      qmEFT = mesh%cells(part%cell)%obj%gatherElectricField(part%Xi)
      qmEFt = qmEFt*part%species%qm*tauIn

      !Update velocity
      part%v = part%v + qmEFt

      !Update position
      part%r = part%r + part%v*tauIn

    END SUBROUTINE pushCartElectrostatic

    PURE SUBROUTINE pushCartElectromagnetic(part, tauIn)
      USE moduleSpecies
      USE moduleMesh
      USE moduleMath
      IMPLICIT NONE

      TYPE(particle), INTENT(inout):: part
      REAL(8), INTENT(in):: tauIn
      REAL(8):: tauInHalf
      REAL(8):: qmEFt(1:3)
      REAL(8):: B(1:3), BNorm
      REAL(8):: fn
      REAL(8):: v_minus(1:3), v_prime(1:3), v_plus(1:3)

      tauInHalf = tauIn *0.5D0
      !Half of the force o f the electric field
      qmEFT = mesh%cells(part%cell)%obj%gatherElectricField(part%Xi)
      qmEFt = qmEFt*part%species%qm*tauInHalf

      !Half step for electrostatic
      v_minus = part%v + qmEFt

      !Full step rotation
      B = mesh%cells(part%cell)%obj%gatherMagneticField(part%Xi)
      BNorm = NORM2(B)
      IF (BNorm > 0.D0) THEN
        fn = DTAN(part%species%qm * tauInHalf*BNorm) / BNorm
        v_prime = v_minus + fn * crossProduct(v_minus, B)
        v_plus  = v_minus + 2.D0 * fn / (1.D0 + fn**2 * B**2)*crossProduct(v_prime, B)

      END IF

      !Half step for electrostatic
      part%v = v_plus + qmEFt

      !Update position
      part%r = part%r + part%v*tauIn

    END SUBROUTINE pushCartElectromagnetic

    !Push one particle. Boris pusher for 2D Cyl Neutral particle
    PURE SUBROUTINE push2DCylNeutral(part, tauIn)
      USE moduleSpecies
      IMPLICIT NONE

      TYPE(particle), INTENT(inout):: part
      REAL(8), INTENT(in):: tauIn
      TYPE(particle):: part_temp
      REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
      REAL(8):: v_p_oh_star(2:3)

      part_temp = part
      !z
      part_temp%v(1) = part%v(1)
      part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
      !r,theta
      v_p_oh_star(2) = part%v(2)
      x_new = part%r(2) + v_p_oh_star(2)*tauIn
      v_p_oh_star(3) = part%v(3)
      y_new =             v_p_oh_star(3)*tauIn
      r = DSQRT(x_new**2+y_new**2)
      part_temp%r(2) = r
      IF (r > 0.D0) THEN
        sin_alpha = y_new/r
        cos_alpha = x_new/r
      ELSE
        sin_alpha = 0.D0
        cos_alpha = 1.D0
      END IF
      part_temp%v(2) =  cos_alpha*v_p_oh_star(2)+sin_alpha*v_p_oh_star(3)
      part_temp%v(3) = -sin_alpha*v_p_oh_star(2)+cos_alpha*v_p_oh_star(3)

      !Copy temporal particle to particle
      part=part_temp

    END SUBROUTINE push2DCylNeutral

    !Push one particle. Boris pusher for 2D Cyl Electrostatic particle
    PURE SUBROUTINE push2DCylElectrostatic(part, tauIn)
      USE moduleSpecies
      USE moduleMesh
      IMPLICIT NONE

      TYPE(particle), INTENT(inout):: part
      REAL(8), INTENT(in):: tauIn
      REAL(8):: v_p_oh_star(2:3)
      TYPE(particle):: part_temp
      REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
      REAL(8):: qmEFt(1:3)!charge*tauIn*EF/mass

      part_temp = part
      !Get electric field at particle position
      qmEFT = mesh%cells(part_temp%cell)%obj%gatherElectricField(part_temp%Xi)
      qmEFt = qmEFt*part_temp%species%qm*tauIn
      !z
      part_temp%v(1) = part%v(1) + qmEFt(1)
      part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
      !r,theta
      v_p_oh_star(2) = part%v(2) + qmEFt(2)
      x_new = part%r(2) + v_p_oh_star(2)*tauIn
      v_p_oh_star(3) = part%v(3) + qmEFt(3)
      y_new =             v_p_oh_star(3)*tauIn
      r = DSQRT(x_new**2+y_new**2)
      part_temp%r(2) = r
      IF (r > 0.D0) THEN
        sin_alpha = y_new/r
        cos_alpha = x_new/r
      ELSE
        sin_alpha = 0.D0
        cos_alpha = 1.D0
      END IF
      part_temp%v(2) =  cos_alpha*v_p_oh_star(2)+sin_alpha*v_p_oh_star(3)
      part_temp%v(3) = -sin_alpha*v_p_oh_star(2)+cos_alpha*v_p_oh_star(3)

      !Copy temporal particle to particle
      part=part_temp

    END SUBROUTINE push2DCylElectrostatic

    !Push one particle. Boris pusher for 1D Radial Neutral particle
    PURE SUBROUTINE push1DRadNeutral(part, tauIn)
      USE moduleSpecies
      IMPLICIT NONE

      TYPE(particle), INTENT(inout):: part
      REAL(8), INTENT(in):: tauIn
      REAL(8):: v_p_oh_star(1:2)
      TYPE(particle):: part_temp
      REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha

      part_temp = part
      !r,theta
      v_p_oh_star(1) = part%v(1)
      x_new = part%r(1) + v_p_oh_star(1)*tauIn
      v_p_oh_star(2) = part%v(2)
      y_new =             v_p_oh_star(2)*tauIn
      r = DSQRT(x_new**2+y_new**2)
      part_temp%r(1) = r
      IF (r > 0.D0) THEN
        sin_alpha = y_new/r
        cos_alpha = x_new/r
      ELSE
        sin_alpha = 0.D0
        cos_alpha = 1.D0
      END IF
      part_temp%v(1) =  cos_alpha*v_p_oh_star(1)+sin_alpha*v_p_oh_star(2)
      part_temp%v(2) = -sin_alpha*v_p_oh_star(1)+cos_alpha*v_p_oh_star(2)

      !Copy temporal particle to particle
      part=part_temp

    END SUBROUTINE push1DRadNeutral

    !Push one particle. Boris pusher for 1D Radial Electrostatic particle
    PURE SUBROUTINE push1DRadElectrostatic(part, tauIn)
      USE moduleSpecies
      USE moduleMesh
      IMPLICIT NONE

      TYPE(particle), INTENT(inout):: part
      REAL(8), INTENT(in):: tauIn
      REAL(8):: v_p_oh_star(1:2)
      TYPE(particle):: part_temp
      REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
      REAL(8):: qmEFt(1:3)!charge*tauIn*EF/mass

      part_temp = part
      !Get electric field at particle position
      qmEFT = mesh%cells(part_temp%cell)%obj%gatherElectricField(part_temp%Xi)
      qmEFt = qmEFt*part_temp%species%qm*tauIn
      !r,theta
      v_p_oh_star(1) = part%v(1) + qmEFt(1)
      x_new = part%r(1) + v_p_oh_star(1)*tauIn
      v_p_oh_star(2) = part%v(2) + qmEFt(2)
      y_new =             v_p_oh_star(2)*tauIn
      r = DSQRT(x_new**2+y_new**2)
      part_temp%r(1) = r
      IF (r > 0.D0) THEN
        sin_alpha = y_new/r
        cos_alpha = x_new/r
      ELSE
        sin_alpha = 0.D0
        cos_alpha = 1.D0
      END IF
      part_temp%v(1) =  cos_alpha*v_p_oh_star(1)+sin_alpha*v_p_oh_star(2)
      part_temp%v(2) = -sin_alpha*v_p_oh_star(1)+cos_alpha*v_p_oh_star(2)

      !Copy temporal particle to particle
      part=part_temp

    END SUBROUTINE push1DRadElectrostatic

    !Dummy pusher for 0D geometry
    PURE SUBROUTINE push0D(part, tauIn)
      USE moduleSpecies
      IMPLICIT NONE

      TYPE(particle), INTENT(inout):: part
      REAL(8), INTENT(in):: tauIn

    END SUBROUTINE push0D

END MODULE modulePusher
