diff --git a/data/see/constant.dat b/data/see/constant.dat
new file mode 100644
index 0000000..060c8ad
--- /dev/null
+++ b/data/see/constant.dat
@@ -0,0 +1,4 @@
+#Relative energy (eV) yield ()
+0.000 1.000E-01
+1.000 1.000E-01
+
diff --git a/data/see/constant_energy.dat b/data/see/constant_energy.dat
new file mode 100644
index 0000000..d9de3ba
--- /dev/null
+++ b/data/see/constant_energy.dat
@@ -0,0 +1,3 @@
+#Secondary energy (eV) Probability (a. u.)
+0.1 0.5
+1.0 0.5
diff --git a/doc/user-manual/.gitignore b/doc/user-manual/.gitignore
index a8475d7..db32e70 100644
--- a/doc/user-manual/.gitignore
+++ b/doc/user-manual/.gitignore
@@ -20,3 +20,4 @@ fpakc_UserManual-blx.bib
 *.gls
 *.ist
 fpakc_UserManual.run.xml
+*.bib.sav
diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90
index ae5c548..9134815 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
@@ -87,7 +105,6 @@ MODULE moduleRandom
 
         END IF
       END DO
-      ! rnd = MINLOC(DABS(rnd0b - cumWeight), 1)
 
     END FUNCTION randomWeighted
 
diff --git a/src/modules/common/moduleTable.f90 b/src/modules/common/moduleTable.f90
index 2dc1c75..a6b8c57 100644
--- a/src/modules/common/moduleTable.f90
+++ b/src/modules/common/moduleTable.f90
@@ -8,6 +8,8 @@ MODULE moduleTable
       PROCEDURE, PASS:: init    => initTable1D
       PROCEDURE, PASS:: get     => getValueTable1D
       PROCEDURE, PASS:: convert => convertUnits
+      PROCEDURE, PASS:: invertXF
+      PROCEDURE, PASS:: cumsumX
 
   END TYPE table1D
 
@@ -30,6 +32,8 @@ MODULE moduleTable
         READ(id, '(A)', iostat = stat) dummy
         !If EOF or error, exit file
         IF (stat /= 0) EXIT
+        !If empty line, exit file
+        IF (dummy == "") EXIT
         !Skip comment
         IF (INDEX(dummy,'#') /= 0) CYCLE
         !Add row
@@ -55,6 +59,7 @@ MODULE moduleTable
         !TODO: Make this a function
         IF (INDEX(dummy,'#') /= 0) CYCLE
         IF (stat /= 0) EXIT
+        IF (dummy == "") EXIT
         !Add data
         !TODO: substitute with extracting information from dummy
         BACKSPACE(id)
@@ -123,4 +128,50 @@ MODULE moduleTable
 
     END SUBROUTINE convertUnits
 
+    ! Invert the arrays of x and f
+    SUBROUTINE invertXF(self)
+      IMPLICIT NONE
+
+      CLASS(table1D), INTENT(inout):: self
+      REAL(8), ALLOCATABLE:: xTemp(:)
+      REAL(8):: xTempMin, xTempMax
+
+      ! Store temporal x data
+      xTemp     = self%x
+      xTempMin  = self%xMin
+      xTempMax  = self%xMax
+
+      ! Replace x data with f data
+      self%x    = self%f
+      self%xMin = self%fMin
+      self%xMax = self%fMax
+
+      ! Put temporal x data in f data
+      self%f    = xTemp
+      self%fMin = xTempMin
+      self%fMax = xTempMax
+
+    END SUBROUTINE invertXF
+
+    ! Makes the x axis its cumulative sum
+    SUBROUTINE cumSumX(self)
+      IMPLICIT NONE
+
+      CLASS(table1D), INTENT(inout):: self
+      INTEGER:: nTotal, n
+      REAL(8), ALLOCATABLE:: cumX(:)
+
+      nTotal = SIZE(self%x)
+      ALLOCATE(cumX(1:nTotal))
+      cumX(1) = self%x(1)
+      DO n = 2, nTotal
+        cumX(n) = self%x(n) + cumX(n-1)
+
+      END DO
+
+      self%x = cumX / cumX(nTotal)
+      DEALLOCATE(cumX)
+
+    END SUBROUTINE cumSumX
+
 END MODULE moduleTable
diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90
index b5f6881..31d9959 100644
--- a/src/modules/init/moduleInput.f90
+++ b/src/modules/init/moduleInput.f90
@@ -95,7 +95,10 @@ MODULE moduleInput
       TYPE(json_file), INTENT(inout):: config
       CHARACTER(LEN=*), INTENT(in):: step
 
-      IF (config%failed()) CALL criticalError("Error reading the JSON input file", TRIM(step))
+      IF (config%failed()) THEN
+        CALL criticalError("Error reading the JSON input file", TRIM(step))
+
+      END IF
 
     END SUBROUTINE checkStatus
 
@@ -808,7 +811,7 @@ MODULE moduleInput
       REAL(8):: effTime
       REAL(8):: eThreshold !Energy threshold
       INTEGER:: speciesID, electronSecondaryID
-      CHARACTER(:), ALLOCATABLE:: speciesName, crossSection, electronSecondary
+      CHARACTER(:), ALLOCATABLE:: speciesName, crossSection, yield, energySpectrum, electronSecondary
       LOGICAL:: found
       INTEGER:: nTypes
 
@@ -883,6 +886,17 @@ MODULE moduleInput
 
             CALL initWallTemperature(boundary(i)%bTypes(s)%obj, Tw, cw)
 
+          CASE('secondaryEmission')
+            CALL config%get(object // '.yield', yield, found)
+            IF (.NOT. found) CALL criticalError("missing parameter 'yield' for secondary emission", 'readBoundary')
+            CALL config%get(object // '.energySpectrum', energySpectrum, found)
+            IF (.NOT. found) CALL criticalError("missing parameter 'energySpectrum' for secondary emission", 'readBoundary')
+            CALL config%get(object // '.electron', speciesName, found)
+            IF (.NOT. found) CALL criticalError("missing parameter 'electron' for secondary emission", 'readBoundary')
+            speciesID = speciesName2Index(speciesName)
+
+            CALL initSEE(boundary(i)%bTypes(s)%obj, yield, energySpectrum, speciesID)
+
           CASE('axis')
             ALLOCATE(boundaryAxis:: boundary(i)%bTypes(s)%obj)
 
diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90
index 091e52e..618c61e 100644
--- a/src/modules/mesh/moduleMeshBoundary.f90
+++ b/src/modules/mesh/moduleMeshBoundary.f90
@@ -199,7 +199,7 @@ MODULE moduleMeshBoundary
 
       END SELECT
 
-      !Removes ionizing electron regardless the number of pair created
+      !Removes ionizing electron regardless the number of pairs created
       part%n_in = .FALSE.
 
     END SUBROUTINE ionization
@@ -218,6 +218,118 @@ MODULE moduleMeshBoundary
 
     END SUBROUTINE symmetryAxis
 
+    !Secondary emission of electrons by impacting particle.
+    SUBROUTINE secondaryEmission(edge, part)
+      USE moduleSpecies
+      USE moduleRandom
+      USE moduleConstParam
+      IMPLICIT NONE
+
+      CLASS(meshEdge), INTENT(inout):: edge
+      CLASS(particle), INTENT(inout):: part
+      REAL(8):: vRel, eRel
+      REAL(8), DIMENSION(1:3):: rElectron, XiElectron!Position of new electrons (impacting particle)
+      REAL(8), DIMENSION(1:3):: rIncident !Vector from imapcting particle position to particle position
+      REAL(8):: theta !incident angle
+      REAL(8):: yield
+      REAL(8):: energyIncident !incident energy corrected by threshold
+      INTEGER:: nNewElectrons
+      REAL(8), ALLOCATABLE:: weight(:)
+      INTEGER:: p
+      INTEGER:: cell
+      TYPE(particle), POINTER:: secondaryElectron
+
+      SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj)
+      TYPE IS(boundarySEE)
+        !Get relative velocity
+        vRel = NORM2(part%v)
+        !Convert to relative energy
+        eRel = part%species%m*vRel**2*5.D-1
+
+        !If energy is abound threshold calculate secondary electrons
+        IF (eRel >= bound%energyThreshold) THEN
+
+          !position of impacting ion
+          rElectron  = edge%intersection(part%r)
+          XiElectron = mesh%cells(part%cell)%obj%phy2log(rElectron)
+
+          !Calculate incident angle
+          rIncident = part%r - rElectron
+          theta = ACOS(DOT_PRODUCT(rIncident, edge%normal) / (NORM2(rIncident) * NORM2(edge%normal)))
+
+          !Calculate incident energy
+          energyIncident = (eRel - bound%energyThreshold) * PI2 / (PI2 + theta**2) + bound%energyThreshold
+
+          !Get number of secondary electrons particles
+          yield = part%weight*bound%yield%get(eRel) * (1.D0 * theta**2 / PI2)
+          
+          !Convert the number to macro-particles
+          nNewElectrons = FLOOR(yield / bound%electron%weight)
+
+          !If the weight of new macro-particles is below the yield, add an additional particle with less weight
+          IF (REAL(nNewElectrons) * bound%electron%weight < yield) THEN
+            nNewElectrons = nNewElectrons + 1
+            ALLOCATE(weight(1:nNewElectrons))
+            weight(1:nNewElectrons-1) = bound%electron%weight
+            weight(nNewElectrons)     = yield - SUM(weight(1:nNewElectrons-1))
+
+          ELSE
+            ALLOCATE(weight(1:nNewElectrons))
+            weight(1:nNewElectrons) = bound%electron%weight
+          
+          END IF
+
+          !New cell of origin
+          IF (ASSOCIATED(edge%e1)) THEN
+            cell = edge%e1%n
+
+          ELSEIF (ASSOCIATED(edge%e2)) THEN
+            cell = edge%e2%n
+
+          END IF
+
+          !Create the new electron macro-particles
+          DO p = 1, nNewElectrons
+            !Create new macro-particle
+            ALLOCATE(secondaryElectron)
+
+            !Assign species to electron
+            secondaryElectron%species => bound%electron
+
+            !Assign position to particle
+            secondaryElectron%r  = rElectron
+            secondaryElectron%Xi = XiElectron
+
+            !Assign cell to electron
+            secondaryElectron%cell = cell
+
+            !Assign weight
+            secondaryElectron%weight = weight(p)
+
+            !Assume particle is inside the numerical domain
+            secondaryElectron%n_in = .TRUE.
+
+            !Assign velocity
+            !TODO: Calculate energy and velocity
+            !TODO: Better way to inject particles from surface. For this and for general injection
+            secondaryElectron%v = 1.D0 * edge%normal + 1.D-1 * (/ randomMaxwellian(), randomMaxwellian(), randomMaxwellian() /)
+
+            !Add particle to list
+            CALL partSurfaces%setLock()
+            CALL partSurfaces%add(secondaryElectron)
+            CALL partSurfaces%unsetLock()
+
+          END DO
+
+        END IF
+
+        !Absorb impacting particle
+        CALL absorption(edge, part)
+
+      END SELECT
+
+    END SUBROUTINE secondaryEmission
+
     !Points the boundary function to specific type
     SUBROUTINE pointBoundaryFunction(edge, s)
       USE moduleErrors
@@ -245,6 +357,9 @@ MODULE moduleMeshBoundary
       TYPE IS(boundaryAxis)
         edge%fBoundary(s)%apply => symmetryAxis
 
+      TYPE IS(boundarySEE)
+        edge%fBoundary(s)%apply => secondaryEmission
+
       CLASS DEFAULT
         CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction')
 
diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90
index 0b76105..091cccf 100644
--- a/src/modules/moduleBoundary.f90
+++ b/src/modules/moduleBoundary.f90
@@ -47,6 +47,16 @@ MODULE moduleBoundary
 
   END TYPE boundaryIonization
 
+  !Secondary electron emission (by ion impact)
+  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundarySEE
+    TYPE(table1D):: yield                     !Yield as a function of ion energy
+    TYPE(table1D):: energySpectrum            !Spectrum of secondary particle emitted
+    CLASS(speciesGeneric), POINTER:: electron !Electron species for secondary emission
+    REAL(8):: energyThreshold
+    CONTAINS
+
+  END TYPE boundarySEE
+
   !Symmetry axis
   TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis
     CONTAINS
@@ -155,4 +165,33 @@ MODULE moduleBoundary
 
     END SUBROUTINE initIonization
 
+    SUBROUTINE initSEE(boundary, yieldFile, energySpectrumFile, speciesID)
+      USE moduleRefParam
+      USE moduleConstParam
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: yieldFile
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: energySpectrumFile
+      INTEGER:: speciesID
+      INTEGER:: i
+
+      ALLOCATE(boundarySEE:: boundary)
+      SELECT TYPE(boundary)
+      TYPE IS(boundarySEE)
+        CALL boundary%yield%init(yieldFile)
+        CALL boundary%yield%convert(eV2J/(m_ref*v_ref**2), 1.D0)
+        boundary%electron => species(speciesID)%obj
+        ! Use first value from table as threshold
+        boundary%energyThreshold = boundary%yield%xMin
+        CALL boundary%energySpectrum%init(energySpectrumFile)
+        CALL boundary%energySpectrum%convert(eV2J/(m_ref*v_ref**2), 1.D0)
+        CALL boundary%energySpectrum%invertXF()
+        CALL boundary%energySpectrum%cumSumX()
+
+      END SELECT
+
+    END SUBROUTINE initSEE
+
 END MODULE moduleBoundary
diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90
index 5b96be2..41f7626 100644
--- a/src/modules/moduleInject.f90
+++ b/src/modules/moduleInject.f90
@@ -283,10 +283,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
 
@@ -312,7 +309,7 @@ MODULE moduleInject
 
       CLASS(injectGeneric), INTENT(in):: self
       INTEGER, SAVE:: nMin
-      INTEGER:: i, e
+      INTEGER:: i, j, e
       INTEGER:: n, sp
       CLASS(meshEdge), POINTER:: randomEdge
       REAL(8):: direction(1:3)
@@ -365,6 +362,7 @@ MODULE moduleInject
 
           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() /)