Correction due to angle

Moving forward to a complete model for SEE. Now the yield has a
correction for the incident angle on the surface (still to check if this
 is correct)

Velocity distribution is still a dummy one.

src/modules/common/moduleTable.f90

@@ -30,6 +30,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 +57,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)

src/modules/mesh/moduleMeshBoundary.f90

@@ -212,19 +212,24 @@ 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):: energy !incident energy corrected by threshold and 
       INTEGER:: nNewElectrons
       REAL(8), ALLOCATABLE:: weight(:)
       INTEGER:: p
-      REAL(8), DIMENSION(1:3):: rElectron, XiElectron!Position of new electrons
       INTEGER:: cell
       TYPE(particle), POINTER:: secondaryElectron
 
@@ -235,64 +240,81 @@ MODULE moduleMeshBoundary
         !Convert to relative energy
         eRel = part%species%m*vRel**2*5.D-1
 
-        !Get number of secondary electrons macro-particles
+        !If energy is abound threshold calculate secondary electrons
-        yield = part%weight*bound%yield%get(eRel)
+        IF (eRel >= bound%energyThreshold) THEN
-        nNewElectrons = FLOOR(yield / bound%electron%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
+          !position of impacting ion
-          ALLOCATE(weight(1:nNewElectrons))
+          rElectron  = edge%intersection(part%r)
-          weight(1:nNewElectrons) = bound%electron%weight
+          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
+          energy = (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) !Check equation!
+          
+          !Convert the number to macro-particles
+          nNewElectrons = FLOOR(yield / bound%electron%weight)
+
+          !If the weight of new macro-particles is below the yield, correct adding a particle
+          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
+            secondaryElectron%v = 2.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
 
-        !position of impacting ion
-        rElectron  = edge%intersection(part%r)
-        XiElectron = mesh%cells(part%cell)%obj%phy2log(rElectron)
-
-        !New cell of origin
-        IF (ASSOCIATED(edge%e1)) THEN
-          cell = edge%e1%n
-
-        ELSEIF (ASSOCIATED(edge%e2)) THEN
-          cell = edge%e2%n
-
-        END IF
-
-        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
-          secondaryElectron%v = 2.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
-
         !Absorb impacting particle
         CALL absorption(edge, part)
 

src/modules/moduleBoundary.f90

@@ -51,6 +51,7 @@ MODULE moduleBoundary
     !Yield as a function of ion energy
     TYPE(table1D):: yield
     CLASS(speciesGeneric), POINTER:: electron !Electron species for secondary emission
+    REAL(8):: energyThreshold
     CONTAINS
 
   END TYPE boundarySEE
@@ -155,6 +156,7 @@ MODULE moduleBoundary
       CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary
       CHARACTER(:), ALLOCATABLE, INTENT(in):: tableFile
       INTEGER:: speciesID
+      INTEGER:: i
 
       ALLOCATE(boundarySEE:: boundary)
       SELECT TYPE(boundary)
@@ -162,6 +164,16 @@ MODULE moduleBoundary
         CALL boundary%yield%init(tableFile)
         CALL boundary%yield%convert(eV2J/(m_ref*v_ref**2), 1.D0)
         boundary%electron => species(speciesID)%obj
+        !Search for the threshold energy in the table
+        DO i = 1, SIZE(boundary%yield%f)
+          IF (boundary%yield%f(i) > 0.D0) THEN
+            boundary%energyThreshold = boundary%yield%x(i)
+
+            EXIT
+
+          END IF
+
+        END DO
 
       END SELECT