Remove near-zero thresholds in atomic rates and make them return zero if the parameters <= 0.

CHANGELOG.md

@@ -16,6 +16,7 @@ New:
 * **Beam dispersion calculation has changed from sigma(z) = sigma + z * tan(alpha) to sigma(z) = sqrt(sigma^2 + (z * tan(alpha))^2) for consistancy with the Gaussian beam model. Attention!!! The results of BES and CX spectroscopy are affected by this change. (#414)**
 * Improved beam direction calculation to allow for natural broadening of the BES line shape due to beam divergence. (#414)
 * Add kwargs to invert_regularised_nnls to pass them to scipy.optimize.nnls. (#438)
+* All interpolated atomic rates now return 0 if plasma parameters <= 0, which matches the behaviour of emission models. (#450)
 
 Bug fixes:
 * Fix deprecated transforms being cached in LaserMaterial after laser.transform update (#420)

cherab/openadas/rates/atomic.pyx

@@ -50,11 +50,8 @@ cdef class IonisationRate(CoreIonisationRate):
     cpdef double evaluate(self, double density, double temperature) except? -1e999:
 
         # need to handle zeros, also density and temperature can become negative due to cubic interpolation
-        if density < 1.e-300:
-            density = 1.e-300
-
-        if temperature < 1.e-300:
-            temperature = 1.e-300
+        if density <= 0 or temperature <= 0:
+            return 0
 
         # calculate rate and convert from log10 space to linear space
         return 10 ** self._rate.evaluate(log10(density), log10(temperature))
@@ -97,11 +94,8 @@ cdef class RecombinationRate(CoreRecombinationRate):
     cpdef double evaluate(self, double density, double temperature) except? -1e999:
 
         # need to handle zeros, also density and temperature can become negative due to cubic interpolation
-        if density < 1.e-300:
-            density = 1.e-300
-
-        if temperature < 1.e-300:
-            temperature = 1.e-300
+        if density <= 0 or temperature <= 0:
+            return 0
 
         # calculate rate and convert from log10 space to linear space
         return 10 ** self._rate.evaluate(log10(density), log10(temperature))
@@ -143,11 +137,8 @@ cdef class ThermalCXRate(CoreThermalCXRate):
     cpdef double evaluate(self, double density, double temperature) except? -1e999:
 
         # need to handle zeros, also density and temperature can become negative due to cubic interpolation
-        if density < 1.e-300:
-            density = 1.e-300
-
-        if temperature < 1.e-300:
-            temperature = 1.e-300
+        if density <= 0 or temperature <= 0:
+            return 0
 
         # calculate rate and convert from log10 space to linear space
         return 10 ** self._rate.evaluate(log10(density), log10(temperature))

cherab/openadas/rates/beam.pyx

@@ -78,14 +78,8 @@ cdef class BeamStoppingRate(CoreBeamStoppingRate):
         """
 
         # need to handle zeros, also density and temperature can become negative due to cubic interpolation
-        if energy < 1.e-300:
-            energy = 1.e-300
-
-        if density < 1.e-300:
-            density = 1.e-300
-
-        if temperature < 1.e-300:
-            temperature = 1.e-300
+        if energy <= 0 or density <= 0 or temperature <= 0:
+            return 0
 
         # calculate rate and convert from log10 space to linear space
         return 10 ** (self._npl_eb.evaluate(log10(energy), log10(density)) + self._tp.evaluate(log10(temperature)))
@@ -152,14 +146,8 @@ cdef class BeamPopulationRate(CoreBeamPopulationRate):
         """
 
         # need to handle zeros, also density and temperature can become negative due to cubic interpolation
-        if energy < 1.e-300:
-            energy = 1.e-300
-
-        if density < 1.e-300:
-            density = 1.e-300
-
-        if temperature < 1.e-300:
-            temperature = 1.e-300
+        if energy <= 0 or density <= 0 or temperature <= 0:
+            return 0
 
         # calculate rate and convert from log10 space to linear space
         return 10 ** (self._npl_eb.evaluate(log10(energy), log10(density)) + self._tp.evaluate(log10(temperature)))
@@ -228,14 +216,8 @@ cdef class BeamEmissionPEC(CoreBeamEmissionPEC):
         """
 
         # need to handle zeros, also density and temperature can become negative due to cubic interpolation
-        if energy < 1.e-300:
-            energy = 1.e-300
-
-        if density < 1.e-300:
-            density = 1.e-300
-
-        if temperature < 1.e-300:
-            temperature = 1.e-300
+        if energy <= 0 or density <= 0 or temperature <= 0:
+            return 0
 
         # calculate rate and convert from log10 space to linear space
         return 10 ** (self._npl_eb.evaluate(log10(energy), log10(density)) + self._tp.evaluate(log10(temperature)))

cherab/openadas/rates/cx.pyx

@@ -88,8 +88,8 @@ cdef class BeamCXPEC(CoreBeamCXPEC):
         cdef double rate
 
         # need to handle zeros for log-log interpolation
-        if energy < 1.e-300:
-            energy = 1.e-300
+        if energy <= 0:
+            return 0
 
         rate = 10 ** self._eb.evaluate(log10(energy))
 

cherab/openadas/rates/pec.pyx

@@ -56,11 +56,8 @@ cdef class ImpactExcitationPEC(CoreImpactExcitationPEC):
     cpdef double evaluate(self, double density, double temperature) except? -1e999:
 
         # need to handle zeros, also density and temperature can become negative due to cubic interpolation
-        if density < 1.e-300:
-            density = 1.e-300
-
-        if temperature < 1.e-300:
-            temperature = 1.e-300
+        if density <= 0 or temperature <= 0:
+            return 0
 
         # calculate rate and convert from log10 space to linear space
         return 10 ** self._rate.evaluate(log10(density), log10(temperature))
@@ -108,11 +105,8 @@ cdef class RecombinationPEC(CoreRecombinationPEC):
     cpdef double evaluate(self, double density, double temperature) except? -1e999:
 
         # need to handle zeros, also density and temperature can become negative due to cubic interpolation
-        if density < 1.e-300:
-            density = 1.e-300
-
-        if temperature < 1.e-300:
-            temperature = 1.e-300
+        if density <= 0 or temperature <= 0:
+            return 0
 
         # calculate rate and convert from log10 space to linear space
         return 10 ** self._rate.evaluate(log10(density), log10(temperature))

cherab/openadas/rates/radiated_power.pyx

@@ -49,11 +49,8 @@ cdef class LineRadiationPower(CoreLineRadiationPower):
     cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999:
 
         # need to handle zeros, also density and temperature can become negative due to cubic interpolation
-        if electron_density < 1.e-300:
-            electron_density = 1.e-300
-
-        if electron_temperature < 1.e-300:
-            electron_temperature = 1.e-300
+        if electron_density <= 0 or electron_temperature <= 0:
+            return 0
 
         # calculate rate and convert from log10 space to linear space
         return 10 ** self._rate.evaluate(log10(electron_density), log10(electron_temperature))
@@ -95,11 +92,8 @@ cdef class ContinuumPower(CoreContinuumPower):
     cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999:
 
         # need to handle zeros, also density and temperature can become negative due to cubic interpolation
-        if electron_density < 1.e-300:
-            electron_density = 1.e-300
-
-        if electron_temperature < 1.e-300:
-            electron_temperature = 1.e-300
+        if electron_density <= 0 or electron_temperature <= 0:
+            return 0
 
         # calculate rate and convert from log10 space to linear space
         return 10 ** self._rate.evaluate(log10(electron_density), log10(electron_temperature))
@@ -140,11 +134,8 @@ cdef class CXRadiationPower(CoreCXRadiationPower):
     cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999:
 
         # need to handle zeros, also density and temperature can become negative due to cubic interpolation
-        if electron_density < 1.e-300:
-            electron_density = 1.e-300
-
-        if electron_temperature < 1.e-300:
-            electron_temperature = 1.e-300
+        if electron_density <= 0 or electron_temperature <= 0:
+            return 0
 
         # calculate rate and convert from log10 space to linear space
         return 10 ** self._rate.evaluate(log10(electron_density), log10(electron_temperature))