Merge pull request #451 from vsnever/enhancement/zero_rates_if_zero_p…

…arams

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

@@ -64,11 +64,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))
@@ -125,11 +122,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))
@@ -185,11 +179,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

@@ -96,14 +96,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)))
@@ -188,14 +182,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)))
@@ -281,14 +269,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

@@ -118,8 +118,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

@@ -70,11 +70,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))
@@ -136,11 +133,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

@@ -68,11 +68,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))
@@ -133,11 +130,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))
@@ -197,11 +191,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))