Making some progress

Things are better organized.

data/see/constant_energy.dat

@@ -0,0 +1,3 @@
+#Secondary energy (eV) Probability (a. u.)
+0.1 0.5
+1.0 0.5

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, yield, electronSecondary
+      CHARACTER(:), ALLOCATABLE:: speciesName, crossSection, yield, energySpectrum, electronSecondary
       LOGICAL:: found
       INTEGER:: nTypes
 
@@ -886,11 +889,13 @@ MODULE moduleInput
           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, speciesID)
+            CALL initSEE(boundary(i)%bTypes(s)%obj, yield, energySpectrum, speciesID)
 
           CASE('axis')
             ALLOCATE(boundaryAxis:: boundary(i)%bTypes(s)%obj)

src/modules/mesh/moduleMeshBoundary.f90

@@ -232,7 +232,7 @@ MODULE moduleMeshBoundary
       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 
+      REAL(8):: energyIncident !incident energy corrected by threshold
       INTEGER:: nNewElectrons
       REAL(8), ALLOCATABLE:: weight(:)
       INTEGER:: p
@@ -258,7 +258,7 @@ MODULE moduleMeshBoundary
           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
+          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)

src/modules/moduleBoundary.f90

@@ -49,8 +49,8 @@ MODULE moduleBoundary
 
   !Secondary electron emission (by ion impact)
   TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundarySEE
-    !Yield as a function of ion energy
+    TYPE(table1D):: yield                     !Yield as a function of ion energy
-    TYPE(table1D):: yield
+    TYPE(table1D):: energySpectrum            !Spectrum of secondary particle emitted
     CLASS(speciesGeneric), POINTER:: electron !Electron species for secondary emission
     REAL(8):: energyThreshold
     CONTAINS
@@ -165,33 +165,30 @@ MODULE moduleBoundary
 
     END SUBROUTINE initIonization
 
-    SUBROUTINE initSEE(boundary, tableFile, speciesID)
+    SUBROUTINE initSEE(boundary, yieldFile, energySpectrumFile, speciesID)
       USE moduleRefParam
       USE moduleConstParam
       USE moduleSpecies
       IMPLICIT NONE
 
       CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary
-      CHARACTER(:), ALLOCATABLE, INTENT(in):: tableFile
+      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(tableFile)
+        CALL boundary%yield%init(yieldFile)
         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
+        ! Use first value from table as threshold
-        DO i = 1, SIZE(boundary%yield%f)
+        boundary%energyThreshold = boundary%yield%xMin
-          IF (boundary%yield%f(i) > 0.D0) THEN
+        CALL boundary%energySpectrum%init(energySpectrumFile)
-            boundary%energyThreshold = boundary%yield%x(i)
+        CALL boundary%energySpectrum%convert(eV2J/(m_ref*v_ref**2), 1.D0)
-
+        CALL boundary%energySpectrum%invertXF()
-            EXIT
+        CALL boundary%energySpectrum%cumSumX()
-
-          END IF
-
-        END DO
 
       END SELECT