diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90
index 9494517..492eaf5 100644
--- a/src/modules/init/moduleInput.f90
+++ b/src/modules/init/moduleInput.f90
@@ -826,6 +826,7 @@ MODULE moduleInput
       real(8):: eThreshold !Energy threshold
       integer:: speciesID, electronSecondaryID
       character(:), allocatable:: speciesName, crossSection, electronSecondary
+      integer:: s_incident
 
       ! Read models of particles
       object = 'boundaries.particles'
@@ -900,7 +901,12 @@ MODULE moduleInput
             END IF
 
           case('quasiNeutrality')
-            call initQuasiNeutrality(self)
+            call config%get(object // '.incident', speciesName, found)
+            if (.not. found) call criticalError("Incident species name not found for quasiNeutrality boundary model", 'readBoundary')
+
+            s_incident = speciesName2Index(speciesName)
+
+            call initQuasiNeutrality(self, s_incident)
 
           CASE DEFAULT
             CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary')
diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index 9133cca..25d7175 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -675,8 +675,9 @@ MODULE moduleMesh
 
     end subroutine initIonization
 
-    module subroutine initQuasiNeutrality(boundary)
+    module subroutine initQuasiNeutrality(boundary, s_incident)
       class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
+      integer, intent(in):: s_incident
 
     end subroutine initQuasiNeutrality
 
@@ -762,7 +763,7 @@ MODULE moduleMesh
 
   ! Ensures quasi-neutrality by changing the reflection coefficient
   type, public, extends(boundaryParticleGeneric):: boundaryQuasiNeutrality
-    real(8):: alpha ! Reflection parameter
+    real(8), allocatable:: alpha(:) ! Reflection parameter
     integer:: s_incident ! species index of the incident species
     type(meshEdgePointer), allocatable:: edges(:) !Array with edges
     contains
diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90
index 040f633..ba01430 100644
--- a/src/modules/mesh/moduleMesh@boundaryParticle.f90
+++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90
@@ -95,18 +95,23 @@ submodule(moduleMesh) boundaryParticle
 
     END SUBROUTINE initIonization
 
-    module subroutine initQuasiNeutrality(boundary)
+    module subroutine initQuasiNeutrality(boundary, s_incident)
       implicit none
 
       class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
+      integer, intent(in):: s_incident
 
       allocate(boundaryQuasiNeutrality:: boundary)
 
       select type(boundary)
       type is(boundaryQuasiNeutrality)
-        boundary%alpha = 0.d0
         allocate(boundary%edges(0))
 
+        boundary%s_incident = s_incident
+
+        allocate(boundary%alpha(1:mesh%numEdges)) ! TODO: Change this so only the edges associated to the boundary are here
+        boundary%alpha = 0.d0
+
         boundary%update => quasiNeutrality_update
 
       end select
@@ -333,7 +338,7 @@ submodule(moduleMesh) boundaryParticle
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 
-      if (random() <= self%alpha) then
+      if (random() <= self%alpha(edge%n)) then
         call genericReflection(edge, part)
 
       else
@@ -347,12 +352,52 @@ submodule(moduleMesh) boundaryParticle
       implicit none
 
       class(boundaryParticleGeneric), intent(inout):: self
+      integer:: e, n, s
+      integer, allocatable:: nodes(:)
+      class(meshEdge), pointer:: edge
+      real(8), allocatable:: density_nodes(:)
+      class(meshNode), pointer:: node
+      real(8):: density_incident, density_rest
+      real(8):: alpha
 
       select type(self)
       type is(boundaryQuasiNeutrality)
-        self%alpha = 0.1d0
+        do e = 1, size(self%edges)
+          edge => mesh%edges(e)%obj
 
-        print *, self%alpha
+          density_incident = 0.d0
+          density_rest     = 0.d0
+
+          nodes = edge%getNodes(edge%nNodes)
+          allocate(density_nodes(1:edge%nNodes))
+          do s = 1, nSpecies
+            do n = 1, edge%nNodes
+              node => mesh%nodes(n)%obj
+              density_nodes(n) = node%output(s)%den
+
+            end do
+
+            if (s == self%s_incident) then
+              density_incident = edge%gatherF((/0.d0,0.d0,0.d0/), edge%nNodes, density_nodes)
+
+            else
+              density_rest = density_rest + edge%gatherF((/0.d0,0.d0,0.d0/), edge%nNodes, density_nodes)
+
+            end if
+
+          end do
+
+          alpha = 1.d0 - density_incident/density_rest
+
+          ! Limit alpha between 0 and 1
+          alpha = min(alpha, 1.d0)
+          alpha = max(alpha, 0.d0)
+
+          self%alpha(edge%n) = alpha
+
+          deallocate(density_nodes)
+
+        end do
 
       end select