Fix self-attenuating RelAct calculations as AD goes to zero.

Previous niave implementation would result in a inf for AD equal to zero, and probably bad values of derivative as it approached; replaced with more well-behaved approach if mu*ad is below a threshold.
Also fixup some debug checks.

Author: wcjohns

InterSpec/PeakDists_imp.hpp

@@ -374,7 +374,7 @@ void gauss_exp_integral( const T peak_mean,
       const double simple_answer = peak_amplitude * PeakDists::gauss_exp_integral( peak_mean, peak_sigma, skew, lower_energy, upper_energy );
       const double diff = fabs(val - simple_answer);
       const double frac_diff = diff / std::max(val, simple_answer);
-      assert( (frac_diff < 1.0E-5) || (diff < 1.0E-8) );
+      assert( (frac_diff < 1.0E-5) || (diff < 1.0E-7) || (fabs(peak_mean - 0.5*(lower_energy + upper_energy)) > (5.0*peak_sigma)) );
     }
 #endif
 

InterSpec/RelActCalc_imp.hpp

@@ -253,9 +253,26 @@ T eval_physical_model_eqn_imp( const double energy,
 
     assert( areal_density >= 0.0 );
 
-    if( (mu >= 0.0) && (areal_density >= 0.0) )
-      answer *= (1.0 - exp(-mu * areal_density)) / (mu * self_atten->areal_density);
-    
+    //The simple computation to handel self-attenuation would be:
+    //  answer *= (1.0 - exp(-mu * areal_density)) / (mu * areal_density);
+    //But if -mu*areal_density is zero, we get NaN, and even as it goes towards zero, we will greatly lose precision
+
+    const T epsilon = T(1.0e-8);
+    const T x = mu * areal_density;
+
+    if( x < epsilon )
+    {
+      // Use the Taylor series expansion for small x.
+      // f(x) = 1 - x/2 + x^2/6
+      answer *= (T(1.0) - 0.5*x + ((x*x) / 6.0)); // + -(x*x*x)/24.0 + ...
+    }else
+    {
+      // For larger x, use the standard formula, but with expm1
+      // to avoid catastrophic cancellation near zero.
+      // 1 - exp(-x) = -expm1(-x)
+      answer *= -expm1(-x) / x;
+    }//if( (mu*areal_density) < epsilon ) / else
+
     assert( !isnan(answer) && !isinf(answer) );
   }//if( self_atten->has_value() )
 

src/PeakFit.cpp

@@ -6412,17 +6412,16 @@ double fit_amp_and_offset( const float *x, const float *data, const size_t nbin,
 {
   //20250410: I believe `PeakFit::fit_amp_and_offset_imp(...)` should be the better function to use, as it uses SVD
   //          but currently leaving the below check in to make sure the re-implementation of this function is correct
-  //          TODO: remove this check after a while
+#define COMPARE_TO_OLD_FIT_WAY 0
   double * const dummy_channel_counts = nullptr;
-  std::vector<double> dummy_amplitudes, dummy_continuum_coeffs, dummy_amplitudes_uncerts, dummy_continuum_coeffs_uncerts;
-#if( !PERFORM_DEVELOPER_CHECKS )
+#if( !COMPARE_TO_OLD_FIT_WAY )
   return PeakFit::fit_amp_and_offset_imp( x, data, nbin, num_polynomial_terms, step_continuum,
                   ref_energy, means, sigmas, fixedAmpPeaks, skew_type, skew_parameters,
                                          amplitudes, continuum_coeffs,
                                          amplitudes_uncerts, continuum_coeffs_uncerts,
                                          dummy_channel_counts );
 #else
-
+  std::vector<double> dummy_amplitudes, dummy_continuum_coeffs, dummy_amplitudes_uncerts, dummy_continuum_coeffs_uncerts;
   const double dummy_chi2 = PeakFit::fit_amp_and_offset_imp( x, data, nbin, num_polynomial_terms, step_continuum,
                           ref_energy, means, sigmas, fixedAmpPeaks, skew_type, skew_parameters,
                                          dummy_amplitudes, dummy_continuum_coeffs,
@@ -6852,7 +6851,7 @@ double fit_amp_and_offset( const float *x, const float *data, const size_t nbin,
   }
 
   return chi2;
-#endif // !PERFORM_DEVELOPER_CHECKS )
+#endif // if( COMPARE_TO_OLD_FIT_WAY ) / else
 }//double fit_amp_and_offset(...)
 
 

src/RelActAutoGuiNuclide.cpp

@@ -1324,6 +1324,9 @@ void RelActAutoGuiNuclide::handleColorChange()
 std::variant<std::monostate, const SandiaDecay::Nuclide *, const SandiaDecay::Element *, const ReactionGamma::Reaction *> RelActAutoGuiNuclide::source() const
 {
   const string nucstr = m_nuclide_edit->text().toUTF8();
+  if( nucstr.empty() )
+    return std::monostate{};
+  
   const SandiaDecay::SandiaDecayDataBase * const db = DecayDataBaseServer::database();
   assert( db );
   if( !db )