First EM pusher

First implementation of Electromagnetic pusher.

Some testing is still required.

Documentation needs to be upgraded to match the changes in this branch.

src/modules/mesh/0D/moduleMesh0D.f90

@@ -19,6 +19,7 @@ MODULE moduleMesh0D
       PROCEDURE, PASS:: randPos     => randPos0D
       PROCEDURE, NOPASS:: fPsi      => fPsi0D
       PROCEDURE, PASS:: gatherEF    => gatherEF0D
+      PROCEDURE, PASS:: gatherMF    => gatherMF0D
       PROCEDURE, PASS:: elemK       => elemK0D
       PROCEDURE, PASS:: elemF       => elemF0D
       PROCEDURE, PASS:: phy2log     => phy2log0D
@@ -123,6 +124,17 @@ MODULE moduleMesh0D
 
     END FUNCTION gatherEF0D
 
+    PURE FUNCTION gatherMF0D(self, xi) RESULT(MF)
+      IMPLICIT NONE
+
+      CLASS(meshVol0D), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8):: MF(1:3)
+
+      MF = 0.D0
+
+    END FUNCTION gatherMF0D
+
     PURE FUNCTION elemK0D(self) RESULT(localK)
       IMPLICIT NONE
 

src/modules/mesh/1DCart/moduleMesh1DCart.f90

@@ -77,6 +77,7 @@ MODULE moduleMesh1DCart
       PROCEDURE, PASS::   elemF       => elemFSegm
       PROCEDURE, NOPASS:: inside      => insideSegm
       PROCEDURE, PASS::   gatherEF    => gatherEFSegm
+      PROCEDURE, PASS::   gatherMF    => gatherMFSegm
       PROCEDURE, PASS::   getNodes    => getNodesSegm
       PROCEDURE, PASS::   phy2log     => phy2logSegm
       PROCEDURE, PASS::   nextElement => nextElementSegm
@@ -388,6 +389,28 @@ MODULE moduleMesh1DCart
 
     END FUNCTION gatherEFSegm
 
+    PURE FUNCTION gatherMFSegm(self, xi) RESULT(MF)
+      IMPLICIT NONE
+
+      CLASS(meshVol1DCartSegm), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8):: fPsi(1:2)
+      REAL(8):: MF_Nodes(1:2, 1:3)
+      REAL(8):: MF(1:3)
+      REAL(8):: invJ
+
+      MF_Nodes(1:2,1) = (/ self%n1%emData%B(1), &
+                           self%n2%emData%B(1) /)
+      MF_Nodes(1:2,2) = (/ self%n1%emData%B(2), &
+                           self%n2%emData%B(2) /)
+      MF_Nodes(1:2,3) = (/ self%n1%emData%B(3), &
+                           self%n2%emData%B(3) /)
+
+      fPsi  = self%fPsi(xi)
+      MF = MATMUL(fPsi, MF_Nodes)
+
+    END FUNCTION gatherMFSegm
+
     !Get nodes from 1D volume
     PURE FUNCTION getNodesSegm(self) RESULT(n)
       IMPLICIT NONE

src/modules/mesh/1DRad/moduleMesh1DRad.f90

@@ -78,6 +78,7 @@ MODULE moduleMesh1DRad
       PROCEDURE, PASS::   elemF       => elemFRad
       PROCEDURE, NOPASS:: inside      => insideRad
       PROCEDURE, PASS::   gatherEF    => gatherEFRad
+      PROCEDURE, PASS::   gatherMF    => gatherMFRad
       PROCEDURE, PASS::   getNodes    => getNodesRad
       PROCEDURE, PASS::   phy2log     => phy2logRad
       PROCEDURE, PASS::   nextElement => nextElementRad
@@ -400,6 +401,28 @@ MODULE moduleMesh1DRad
 
     END FUNCTION gatherEFRad
 
+    PURE FUNCTION gatherMFRad(self, xi) RESULT(MF)
+      IMPLICIT NONE
+
+      CLASS(meshVol1DRadSegm), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8):: fPsi(1:2)
+      REAL(8):: MF_Nodes(1:2, 1:3)
+      REAL(8):: MF(1:3)
+      REAL(8):: invJ
+
+      MF_Nodes(1:2,1) = (/ self%n1%emData%B(1), &
+                           self%n2%emData%B(1) /)
+      MF_Nodes(1:2,2) = (/ self%n1%emData%B(2), &
+                           self%n2%emData%B(2) /)
+      MF_Nodes(1:2,3) = (/ self%n1%emData%B(3), &
+                           self%n2%emData%B(3) /)
+
+      fPsi  = self%fPsi(xi)
+      MF = MATMUL(fPsi, MF_Nodes)
+
+    END FUNCTION gatherMFRad
+
     !Get nodes from 1D volume
     PURE FUNCTION getNodesRad(self) RESULT(n)
       IMPLICIT NONE

src/modules/mesh/2DCart/moduleMesh2DCart.f90

@@ -85,6 +85,7 @@ MODULE moduleMesh2DCart
       PROCEDURE, PASS::   elemF       => elemFQuad
       PROCEDURE, NOPASS:: inside      => insideQuad
       PROCEDURE, PASS::   gatherEF    => gatherEFQuad
+      PROCEDURE, PASS::   gatherMF    => gatherMFQuad
       PROCEDURE, PASS::   getNodes    => getNodesQuad
       PROCEDURE, PASS::   phy2log     => phy2logQuad
       PROCEDURE, PASS::   nextElement => nextElementQuad
@@ -114,6 +115,7 @@ MODULE moduleMesh2DCart
       PROCEDURE, PASS::   elemF       => elemFTria
       PROCEDURE, NOPASS:: inside      => insideTria
       PROCEDURE, PASS::   gatherEF    => gatherEFTria
+      PROCEDURE, PASS::   gatherMF    => gatherMFTria
       PROCEDURE, PASS::   getNodes    => getNodesTria
       PROCEDURE, PASS::   phy2log     => phy2logTria
       PROCEDURE, PASS::   nextElement => nextElementTria
@@ -505,6 +507,33 @@ MODULE moduleMesh2DCart
 
     END FUNCTION gatherEFQuad
 
+    PURE FUNCTION gatherMFQuad(self,xi) RESULT(MF)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCartQuad), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8):: fPsi(1:4)
+      REAL(8):: MF_Nodes(1:4,1:3)
+      REAL(8):: MF(1:3)
+
+      MF_Nodes(1:4,1) = (/self%n1%emData%B(1), &
+                          self%n2%emData%B(1), &
+                          self%n3%emData%B(1), &
+                          self%n4%emData%B(1) /)
+      MF_Nodes(1:4,2) = (/self%n1%emData%B(2), &
+                          self%n2%emData%B(2), &
+                          self%n3%emData%B(2), &
+                          self%n4%emData%B(2) /)
+      MF_Nodes(1:4,3) = (/self%n1%emData%B(3), &
+                          self%n2%emData%B(3), &
+                          self%n3%emData%B(3), &
+                          self%n4%emData%B(3) /)
+
+      fPsi  =  self%fPsi(xi)
+      MF = MATMUL(fPsi(:), MF_Nodes)
+
+    END FUNCTION gatherMFQuad
+
     !Gets nodes from quadrilateral element
     PURE FUNCTION getNodesQuad(self) RESULT(n)
       IMPLICIT NONE
@@ -816,6 +845,30 @@ MODULE moduleMesh2DCart
 
     END FUNCTION gatherEFTria
 
+    PURE FUNCTION gatherMFTria(self,xi) RESULT(MF)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCartTria), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8):: fPsi(1:3)
+      REAL(8):: MF_Nodes(1:3,1:3)
+      REAL(8):: MF(1:3)
+
+      MF_Nodes(1:3,1) = (/self%n1%emData%B(1), &
+                          self%n2%emData%B(1), &
+                          self%n3%emData%B(1) /)
+      MF_Nodes(1:3,2) = (/self%n1%emData%B(2), &
+                          self%n2%emData%B(2), &
+                          self%n3%emData%B(2) /)
+      MF_Nodes(1:3,3) = (/self%n1%emData%B(3), &
+                          self%n2%emData%B(3), &
+                          self%n3%emData%B(3) /)
+
+      fPsi  =  self%fPsi(xi)
+      MF = MATMUL(fPsi, MF_Nodes)
+
+    END FUNCTION gatherMFTria
+
     !Gets node indexes from triangular element
     PURE FUNCTION getNodesTria(self) RESULT(n)
       IMPLICIT NONE

src/modules/mesh/2DCyl/moduleMesh2DCyl.f90

@@ -86,6 +86,7 @@ MODULE moduleMesh2DCyl
       PROCEDURE, PASS::   elemF       => elemFQuad
       PROCEDURE, NOPASS:: inside      => insideQuad
       PROCEDURE, PASS::   gatherEF    => gatherEFQuad
+      PROCEDURE, PASS::   gatherMF    => gatherMFQuad
       PROCEDURE, PASS::   getNodes    => getNodesQuad
       PROCEDURE, PASS::   phy2log     => phy2logQuad
       PROCEDURE, PASS::   nextElement => nextElementQuad
@@ -115,6 +116,7 @@ MODULE moduleMesh2DCyl
       PROCEDURE, PASS::   elemF       => elemFTria
       PROCEDURE, NOPASS:: inside      => insideTria
       PROCEDURE, PASS::   gatherEF    => gatherEFTria
+      PROCEDURE, PASS::   gatherMF    => gatherMFTria
       PROCEDURE, PASS::   getNodes    => getNodesTria
       PROCEDURE, PASS::   phy2log     => phy2logTria
       PROCEDURE, PASS::   nextElement => nextElementTria
@@ -526,6 +528,33 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION gatherEFQuad
 
+    PURE FUNCTION gatherMFQuad(self,xi) RESULT(MF)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCylQuad), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8):: fPsi(1:4)
+      REAL(8):: MF_Nodes(1:4,1:3)
+      REAL(8):: MF(1:3)
+
+      MF_Nodes(1:4,1) = (/self%n1%emData%B(1), &
+                          self%n2%emData%B(1), &
+                          self%n3%emData%B(1), &
+                          self%n4%emData%B(1) /)
+      MF_Nodes(1:4,2) = (/self%n1%emData%B(2), &
+                          self%n2%emData%B(2), &
+                          self%n3%emData%B(2), &
+                          self%n4%emData%B(2) /)
+      MF_Nodes(1:4,3) = (/self%n1%emData%B(3), &
+                          self%n2%emData%B(3), &
+                          self%n3%emData%B(3), &
+                          self%n4%emData%B(3) /)
+
+      fPsi  =  self%fPsi(xi)
+      MF = MATMUL(fPsi(:), MF_Nodes)
+
+    END FUNCTION gatherMFQuad
+
     !Gets nodes from quadrilateral element
     PURE FUNCTION getNodesQuad(self) RESULT(n)
       IMPLICIT NONE
@@ -845,6 +874,30 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION gatherEFTria
 
+    PURE FUNCTION gatherMFTria(self,xi) RESULT(MF)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCylTria), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8):: fPsi(1:3)
+      REAL(8):: MF_Nodes(1:3,1:3)
+      REAL(8):: MF(1:3)
+
+      MF_Nodes(1:3,1) = (/self%n1%emData%B(1), &
+                          self%n2%emData%B(1), &
+                          self%n3%emData%B(1) /)
+      MF_Nodes(1:3,2) = (/self%n1%emData%B(2), &
+                          self%n2%emData%B(2), &
+                          self%n3%emData%B(2) /)
+      MF_Nodes(1:3,3) = (/self%n1%emData%B(3), &
+                          self%n2%emData%B(3), &
+                          self%n3%emData%B(3) /)
+
+      fPsi  =  self%fPsi(xi)
+      MF = MATMUL(fPsi, MF_Nodes)
+
+    END FUNCTION gatherMFTria
+
     !Gets node indexes from triangular element
     PURE FUNCTION getNodesTria(self) RESULT(n)
       IMPLICIT NONE

src/modules/mesh/3DCart/moduleMesh3DCart.f90

@@ -78,6 +78,7 @@ MODULE moduleMesh3DCart
       PROCEDURE, PASS::   elemF       => elemFTetra
       PROCEDURE, NOPASS:: inside      => insideTetra
       PROCEDURE, PASS::   gatherEF    => gatherEFTetra
+      PROCEDURE, PASS::   gatherMF    => gatherMFTetra
       PROCEDURE, PASS::   getNodes    => getNodesTetra
       PROCEDURE, PASS::   phy2log     => phy2logTetra
       PROCEDURE, PASS::   nextElement => nextElementTetra
@@ -498,6 +499,33 @@ MODULE moduleMesh3DCart
 
     END FUNCTION gatherEFTetra
 
+    PURE FUNCTION gatherMFTetra(self, xi) RESULT(MF)
+      IMPLICIT NONE
+
+      CLASS(meshVol3DCartTetra), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8):: fPsi(1:4)
+      REAL(8):: MF_Nodes(1:4,1:3)
+      REAL(8):: MF(1:3)
+
+      MF_Nodes(1:4,1) = (/self%n1%emData%B(1), &
+                          self%n2%emData%B(1), &
+                          self%n3%emData%B(1), &
+                          self%n4%emData%B(1) /)
+      MF_Nodes(1:4,2) = (/self%n1%emData%B(2), &
+                          self%n2%emData%B(2), &
+                          self%n3%emData%B(2), &
+                          self%n4%emData%B(2) /)
+      MF_Nodes(1:4,3) = (/self%n1%emData%B(3), &
+                          self%n2%emData%B(3), &
+                          self%n3%emData%B(3), &
+                          self%n4%emData%B(3) /)
+
+      fPsi  = self%fPsi(xi)
+      MF = MATMUL(fPsi, MF_Nodes)
+
+    END FUNCTION gatherMFTetra
+
     PURE FUNCTION getNodesTetra(self) RESULT(n)
       IMPLICIT NONE
 

src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90

@@ -200,6 +200,21 @@ MODULE moduleMeshOutputGmsh2
           WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF(xi)*EF_ref
         END DO
         WRITE(20, "(A)") '$EndElementData'
+
+        WRITE(20, "(A)") '$NodeData'
+        WRITE(20, "(A)") '1'
+        WRITE(20, "(A)") '"Magnetic Field (T)"'
+        WRITE(20, *) 1
+        WRITE(20, *) time
+        WRITE(20, *) 3
+        WRITE(20, *) t
+        WRITE(20, *) 3
+        WRITE(20, *) self%numNodes
+        DO n=1, self%numNodes
+          WRITE(20, *) n, self%nodes(n)%obj%emData%B * B_ref
+        END DO
+        WRITE(20, "(A)") '$EndNodeData'
+
         CLOSE(20)
 
       END IF

src/modules/mesh/moduleMesh.f90

@@ -166,6 +166,7 @@ MODULE moduleMesh
       PROCEDURE(fPsi_interface),        DEFERRED, NOPASS:: fPsi
       PROCEDURE,                                  PASS::   scatter
       PROCEDURE(gatherEF_interface),    DEFERRED, PASS::   gatherEF
+      PROCEDURE(gatherMF_interface),    DEFERRED, PASS::   gatherMF
       PROCEDURE(elemK_interface),       DEFERRED, PASS::   elemK
       PROCEDURE(elemF_interface),       DEFERRED, PASS::   elemF
       PROCEDURE,                                  PASS::   findCell
@@ -194,6 +195,14 @@ MODULE moduleMesh
 
     END FUNCTION gatherEF_interface
 
+    PURE FUNCTION gatherMF_interface(self, xi) RESULT(MF)
+      IMPORT:: meshVol
+      CLASS(meshVol), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8):: MF(1:3)
+
+    END FUNCTION gatherMF_interface
+
     PURE FUNCTION getNodesVol_interface(self) RESULT(n)
       IMPORT:: meshVol
       CLASS(meshVol), INTENT(in):: self

src/modules/moduleEM.f90

@@ -62,6 +62,23 @@ MODULE moduleEM
 
     END FUNCTION gatherElecField
 
+    PURE FUNCTION gatherMagnField(part) RESULT(BField)
+      USE moduleSpecies
+      USE moduleMesh
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(in):: part
+      REAl(8):: xi(1:3) !Logical coordinates of particle in element
+      REAL(8):: BField(1:3)
+
+      BField = 0.D0
+
+      xi  = part%xi
+
+      BField = mesh%vols(part%vol)%obj%gatherMF(xi)
+
+    END FUNCTION gatherMagnField
+
     SUBROUTINE assembleSourceVector(vectorF)
       USE moduleMesh
       USE moduleRefParam

src/modules/moduleInput.f90

@@ -155,6 +155,7 @@ MODULE moduleInput
       Vol_ref   = L_ref**3                       !reference volume
       EF_ref    = qe*n_ref*L_ref/eps_0           !reference electric field
       Volt_ref  = EF_ref*L_ref                   !reference voltage
+      B_ref     = m_ref / (ti_ref * qe)          !reference magnetic field
 
     END SUBROUTINE readReference
 
@@ -177,8 +178,9 @@ MODULE moduleInput
       !simulation final and initial times in [t]
       REAL(8):: finalTime, initialTime
       CHARACTER(:), ALLOCATABLE:: pusherType, WSType, EMType
+      REAL(8):: B(1:3)
       INTEGER:: nTau, nSolver
-      INTEGER:: i
+      INTEGER:: i, n
       CHARACTER(2):: iString
       CHARACTER(1):: tString
 
@@ -267,6 +269,31 @@ MODULE moduleInput
           !Read BC
           CALL readEMBoundary(config)
 
+        CASE("ConstantB")
+          !Read BC
+          CALL readEMBoundary(config)
+          !Read constant magnetic field
+          DO i = 1, 3
+            WRITE(istring, '(i2)') i
+            CALL config%get(object // '.B(' // istring // ')', B(i), found)
+            IF (.NOT. found) THEN
+              CALL criticalError('Constant magnetic field not provided in direction ' // iString, 'readSolver')
+
+            END IF
+
+          END DO
+          
+          !Non-dimensional magnetic field
+          B = B / B_ref
+
+          !Assign it to the nodes
+          DO n =1, mesh%numNodes
+            mesh%nodes(n)%obj%emData%B(1) = B(1)
+            mesh%nodes(n)%obj%emData%B(2) = B(2)
+            mesh%nodes(n)%obj%emData%B(3) = B(3)
+
+          END DO
+
         CASE DEFAULT
           CALL criticalError('EM Solver ' // EMType // ' not found', 'readSolver')
 

src/modules/moduleOutput.f90

@@ -11,6 +11,7 @@ MODULE moduleOutput
   TYPE emNode
     CHARACTER(:), ALLOCATABLE:: type
     REAL(8):: phi
+    REAL(8):: B(1:3)
 
   END TYPE emNode
 

src/modules/moduleRefParam.f90

@@ -3,7 +3,7 @@ MODULE moduleRefParam
   !Parameters that define the problem (inputs)
   REAL(8):: n_ref, m_ref, T_ref, r_ref, debye_ref, sigma_ref
   !Reference parameters for non-dimensional problem
-  REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref
+  REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref
 
 END MODULE moduleRefParam
 

src/modules/moduleSolver.f90

@@ -27,7 +27,7 @@ MODULE moduleSolver
 
   INTERFACE
     !Push a particle
-    PURE SUBROUTINE push_interafece(part, tauIn)
+    SUBROUTINE push_interafece(part, tauIn)
       USE moduleSpecies
 
       TYPE(particle), INTENT(inout):: part
@@ -80,7 +80,10 @@ MODULE moduleSolver
             self%pushParticle => pushCartNeutral
 
           CASE('Electrostatic')
-            self%pushParticle => pushCartCharged
+            self%pushParticle => pushCartElectrostatic
+
+          CASE('Electromagnetic')
+            self%pushParticle => pushCartElectromagnetic
 
           CASE DEFAULT
             CALL criticalError('Pusher ' // pusherType // ' not found for Cart','initPusher')
@@ -88,22 +91,17 @@ MODULE moduleSolver
           END SELECT
 
         CASE('Cyl')
-          IF (self%dimen == 2) THEN
+          SELECT CASE(pusherType)
-            SELECT CASE(pusherType)
+          CASE('Neutral')
-            CASE('Neutral')
+            self%pushParticle => push2DCylNeutral
-              self%pushParticle => push2DCylNeutral
 
-            CASE('Electrostatic')
+          CASE('Electrostatic')
-              self%pushParticle => push2DCylCharged
+            self%pushParticle => push2DCylElectrostatic
 
-            CASE DEFAULT
+          CASE DEFAULT
-              CALL criticalError('Pusher ' // pusherType // ' not found for Cyl','initPusher')
+            CALL criticalError('Pusher ' // pusherType // ' not found for Cyl','initPusher')
 
-            END SELECT
+          END SELECT
-
-          END IF
-
-        END SELECT
 
         CASE('Rad')
           SELECT CASE(pusherType)
@@ -111,7 +109,7 @@ MODULE moduleSolver
             self%pushParticle => push1DRadNeutral
 
           CASE('Electrostatic')
-            self%pushParticle => push1DRadCharged
+            self%pushParticle => push1DRadElectrostatic
 
           CASE DEFAULT
             CALL criticalError('Pusher ' // pusherType // ' not found for Rad','initPusher')
@@ -134,7 +132,7 @@ MODULE moduleSolver
       CHARACTER(:), ALLOCATABLE:: EMType
 
       SELECT CASE(EMType)
-      CASE('Electrostatic')
+      CASE('Electrostatic','ConstantB')
         self%solveEM => solveElecField
 
       END SELECT
@@ -191,56 +189,71 @@ MODULE moduleSolver
 
       TYPE(particle), INTENT(inout):: part
       REAL(8), INTENT(in):: tauIn
-      TYPE(particle):: part_temp
 
-      part_temp = part
+      part%r = part%r + part%v*tauIn
-
-      !x
-      part_temp%v(1) = part%v(1)
-      part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
-
-      !y
-      part_temp%v(2) = part%v(2)
-      part_temp%r(2) = part%r(2) + part_temp%v(2)*tauIn
-
-      !z
-      part_temp%v(3) = part%v(3)
-      part_temp%r(3) = part%r(3) + part_temp%v(3)*tauIn
-
-      part = part_temp
 
     END SUBROUTINE pushCartNeutral
 
-    !Push charged particles in 3D cartesian coordinates
+    PURE SUBROUTINE pushCartElectrostatic(part, tauIn)
-    PURE SUBROUTINE pushCartCharged(part, tauIn)
       USE moduleSPecies
       USE moduleEM
       IMPLICIT NONE
 
       TYPE(particle), INTENT(inout):: part
       REAL(8), INTENT(in):: tauIn
-      TYPE(particle):: part_temp
       REAL(8):: qmEFt(1:3)
 
-      part_temp = part
       !Get the electric field at particle position
-      qmEFt = part_temp%qm*gatherElecField(part_temp)*tauIn
+      qmEFt = part%qm*gatherElecField(part)*tauIn
 
-      !x
+      !Update velocity
-      part_temp%v(1) = part%v(1) + qmEFt(1)
+      part%v = part%v + qmEFt
-      part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
 
-      !y
+      !Update position
-      part_temp%v(2) = part%v(2) + qmEFt(2)
+      part%r = part%r + part%v*tauIn
-      part_temp%r(2) = part%r(2) + part_temp%v(2)*tauIn
 
-      !z
+    END SUBROUTINE pushCartElectrostatic
-      part_temp%v(3) = part%v(3) + qmEFt(3)
-      part_temp%r(3) = part%r(3) + part_temp%v(3)*tauIn
 
-      part = part_temp
+    SUBROUTINE pushCartElectromagnetic(part, tauIn)
+      USE moduleSPecies
+      USE moduleEM
+      USE moduleMath
+      IMPLICIT NONE
 
-    END SUBROUTINE pushCartCharged
+      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 = part%qm*gatherElecField(part)*tauInHalf
+
+      !Half step for electrostatic
+      v_minus = part%v + qmEFt
+
+      !Full step rotation
+      B = gatherMagnField(part)
+      BNorm = NORM2(B)
+      IF (BNorm > 0.D0) THEN
+        fn = DTAN(part%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)
+
+        PRINT *, v_minus, v_plus
+
+      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)
@@ -279,8 +292,8 @@ MODULE moduleSolver
 
     END SUBROUTINE push2DCylNeutral
 
-    !Push one particle. Boris pusher for 2D Cyl Charged particle
+    !Push one particle. Boris pusher for 2D Cyl Electrostatic particle
-    PURE SUBROUTINE push2DCylCharged(part, tauIn)
+    PURE SUBROUTINE push2DCylElectrostatic(part, tauIn)
       USE moduleSpecies
       USE moduleEM
       IMPLICIT NONE
@@ -318,7 +331,7 @@ MODULE moduleSolver
       !Copy temporal particle to particle
       part=part_temp
 
-    END SUBROUTINE push2DCylCharged
+    END SUBROUTINE push2DCylElectrostatic
 
     !Push one particle. Boris pusher for 1D Radial Neutral particle
     PURE SUBROUTINE push1DRadNeutral(part, tauIn)
@@ -355,8 +368,8 @@ MODULE moduleSolver
 
     END SUBROUTINE push1DRadNeutral
 
-    !Push one particle. Boris pusher for 1D Radial Charged particle
+    !Push one particle. Boris pusher for 1D Radial Electrostatic particle
-    PURE SUBROUTINE push1DRadCharged(part, tauIn)
+    PURE SUBROUTINE push1DRadElectrostatic(part, tauIn)
       USE moduleSpecies
       USE moduleEM
       IMPLICIT NONE
@@ -391,7 +404,7 @@ MODULE moduleSolver
       !Copy temporal particle to particle
       part=part_temp
 
-    END SUBROUTINE push1DRadCharged
+    END SUBROUTINE push1DRadElectrostatic
 
     !Dummy pusher for 0D geometry
     PURE SUBROUTINE push0D(part, tauIn)