diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90
index cd553a8..657ca28 100644
--- a/src/modules/common/moduleRandom.f90
+++ b/src/modules/common/moduleRandom.f90
@@ -66,6 +66,24 @@ MODULE moduleRandom
 
     END FUNCTION randomMaxwellian
 
+    FUNCTION randomHalfMaxwellian() RESULT(rnd)
+      USE moduleConstParam, ONLY: PI
+      IMPLICIT NONE
+
+      REAL(8):: rnd
+      REAL(8):: x, y
+
+      rnd = 0.D0
+      x = 0.D0
+      DO WHILE (x == 0.D0)
+        CALL RANDOM_NUMBER(x)
+      END DO
+      y = random(-PI, PI)
+
+      rnd = DSQRT(-2.D0*DLOG(x))*DCOS(y)
+
+    END FUNCTION randomHalfMaxwellian
+
     !Returns a random number weighted with the cumWeight array
     FUNCTION randomWeighted(cumWeight, sumWeight) RESULT(rnd)
       IMPLICIT NONE
diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90
index 4e57083..644cfa7 100644
--- a/src/modules/moduleInject.f90
+++ b/src/modules/moduleInject.f90
@@ -254,10 +254,7 @@ MODULE moduleInject
       REAL(8):: v
       v = 0.D0
 
-      DO WHILE (v <= 0.D0)
-        v = self%vTh*randomMaxwellian()
-
-      END DO
+      v = self%vTh*randomHalfMaxwellian()
 
     END FUNCTION randomVelHalfMaxwellian
 
@@ -334,17 +331,19 @@ MODULE moduleInject
 
         direction = self%n
 
-        partInj(n)%v = 0.D0
-
+        !Sample initial velocity
         partInj(n)%v = self%vMod*direction +  (/ self%v(1)%obj%randomVel(), &
                                                  self%v(2)%obj%randomVel(), &
                                                  self%v(3)%obj%randomVel() /)
 
-        !If velocity is not in the right direction, invert it
-        IF (DOT_PRODUCT(direction, partInj(n)%v) < 0.D0) THEN
-          partInj(n)%v = - partInj(n)%v
+        !For each direction, velocities have to agree with the direction of injection
+        DO i = 1, 3
+          DO WHILE (partInj(n)%v(i)*direction(i) < 0)
+            partInj(n)%v(i) = self%vMod*direction(i) + self%v(i)%obj%randomVel()
 
-        END IF
+          END DO
+
+        END DO
 
         !Obtain natural coordinates of particle in cell
         partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)