From 7f73b69dc2f42933344d75e887a63a76697e9c89 Mon Sep 17 00:00:00 2001
From: JGonzalez <jorge.gonzalez@upm.es>
Date: Sun, 27 Jul 2025 17:14:38 +0200
Subject: [PATCH] Fix injection

Half-Maxwellian distribution should inject particles correctly
---
 src/modules/common/moduleRandom.f90 | 17 +++++++++++++++++
 src/modules/moduleInject.f90        |  5 +----
 2 files changed, 18 insertions(+), 4 deletions(-)

diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90
index ae5c548..7d9cec6 100644
--- a/src/modules/common/moduleRandom.f90
+++ b/src/modules/common/moduleRandom.f90
@@ -66,6 +66,23 @@ MODULE moduleRandom
 
     END FUNCTION randomMaxwellian
 
+    !Returns a random number in a Maxwellian distribution of mean 0 and width 1
+    FUNCTION randomHalfMaxwellian() RESULT(rnd)
+      IMPLICIT NONE
+
+      REAL(8):: rnd
+      REAL(8):: x
+
+      rnd = 0.D0
+      x = 0.D0
+      DO WHILE (x == 0.D0)
+        CALL RANDOM_NUMBER(x)
+      END DO
+
+      rnd = DSQRT(-DLOG(x))
+
+    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 b4e0ed2..b6fc64b 100644
--- a/src/modules/moduleInject.f90
+++ b/src/modules/moduleInject.f90
@@ -298,10 +298,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