Merge pull request #203 from njoy/fix/mf8-mt45-discrete

Adding the possibility of using NT=8 in MF8 MT457 DiscreteSpectrum

.github/workflows/ContinuousIntegration.yml

@@ -16,8 +16,8 @@ jobs:
     runs-on: ${{matrix.os}}
     strategy:
       matrix:
-        os: [ ubuntu-20.04, macos-11 ]
-        cxx: [ clang++, g++-10 ]
+        os: [ ubuntu-22.04, macos-12 ]
+        cxx: [ clang++, g++-11 ]
         build_type: [ Debug, Release ]
 
     steps:

ReleaseNotes.md

@@ -4,7 +4,9 @@ Given here are some release notes for ENDFtk.
 ## [ENDFtk v1.1.1](https://github.com/njoy/ENDFtk/pull/xxx)
 This update removes the regions() and pairs() interface functions on the TAB1 record interface functions that are unused. The removal of these interface functions has no impact on the Python interface as these interface functions were not included on the Python side.
 
-In addition, a minor bug in the rectangular matrix covariance block was corrected. The values for the row and column energies are lifted out of a larger array using the std::ranges::take and std::ranges::drop function. For the column energies, we forgot to properly end the sequence. As a result, the end() iterator of the range did not point to the end of the column energies but to the end of the covariance values. This has been corrected.
+In addition, the following issues were corrected:
+  - A minor bug in the rectangular matrix covariance block was corrected. The values for the row and column energies are lifted out of a larger array using the std::ranges::take and std::ranges::drop function. For the column energies, we forgot to properly end the sequence. As a result, the end() iterator of the range did not point to the end of the column energies but to the end of the covariance values, which is now corrected.
+  - In MF8 MT457 DiscreteSpectrum, NT=8 (which can occur for electrons, i.e. STYP=8), was explicitly disallowed by ENDFtk. This was corrected and a new constructor reflecting this usage was added as well.
 
 A few changes were also made to remove some range-v3 code in MF1 MT451. These changes have no impact on functionality.
 

python/src/section/8/457/DiscreteSpectrum.python.cpp

@@ -60,6 +60,24 @@ void wrapDiscreteSpectrum( python::module& module, python::module& ) {
     "    type         the transition type\n"
     "    ris          internal pair formation coefficient or positron intensity"
   )
+  .def(
+
+    python::init< double, const std::array< double, 2 >&,
+                  const std::array< double, 2 >&, double,
+                  const std::array< double, 2 >&,
+                  const std::array< double, 2 >& >(),
+    python::arg( "chain" ), python::arg( "energy" ), python::arg( "intensity" ),
+    python::arg( "type" ), python::arg( "ris" ), python::arg( "ricc" ),
+    "Initialise the component\n\n"
+    "Arguments:\n"
+    "    self         the component\n"
+    "    chain        the ENDF decay chain\n"
+    "    energy       the energy and uncertainty\n"
+    "    intensity    the relative intensity and uncertainty\n"
+    "    type         the transition type\n"
+    "    ris          internal pair formation coefficient or positron intensity\n"
+    "    ricc         the total internal conversion coefficient"
+  )
   .def(
 
     python::init< double, const std::array< double, 2 >&,

python/test/MF8/Test_ENDFtk_MF8_MT457_DiscreteSpectrum.py

@@ -19,6 +19,10 @@ class Test_ENDFtk_MF8_MT457_DiscreteSpectrum( unittest.TestCase ) :
                      ' 3.000000+0 0.000000+0 1.420112-6 2.85306-10 1.000000+0 2.000000+03580 8457     \n'
                      ' 7.010000+5 1.106180-2 1.378999-1 2.000000-3 3.315735+5 5.608008+33580 8457     \n' )
 
+    chunk_styp_8 = ( ' 4.863000+4 5.000000+1          0          0          8          03580 8457     \n'
+                     ' 3.000000+0 1.000000+0 1.420112-6 2.85306-10 1.000000+0 2.000000+03580 8457     \n'
+                     ' 7.010000+5 1.106180-2                                            3580 8457     \n' )
+
     invalid = ( ' 4.863000+4 5.000000+1          0          0         10          03580 8457     \n'
                 ' 3.000000+0 0.000000+0 1.420112-6 2.85306-10 1.000000+0 2.000000+03580 8457     \n'
                 ' 7.010000+5 1.106180-2 1.378999-1 2.000000-3                      3580 8457     \n' )
@@ -139,6 +143,44 @@ def verify_chunk_styp_0( self, chunk ) :
             # verify string
             self.assertEqual( self.chunk_styp_0, chunk.to_string( 3580, 8, 457 ) )
 
+        def verify_chunk_styp_8( self, chunk ) :
+
+            # verify content
+            self.assertEqual( 8, chunk.NT )
+            self.assertAlmostEqual( 3., chunk.RTYP )
+            self.assertAlmostEqual( 3., chunk.decay_chain )
+            self.assertAlmostEqual( 4.863000e+4, chunk.ER[0] )
+            self.assertAlmostEqual( 5.000000e+1, chunk.ER[1] )
+            self.assertAlmostEqual( 4.863000e+4, chunk.discrete_energy[0] )
+            self.assertAlmostEqual( 5.000000e+1, chunk.discrete_energy[1] )
+            self.assertAlmostEqual( 1.420112e-6, chunk.RI[0] )
+            self.assertAlmostEqual( 2.85306e-10, chunk.RI[1] )
+            self.assertAlmostEqual( 1.420112e-6, chunk.relative_intensity[0] )
+            self.assertAlmostEqual( 2.85306e-10, chunk.relative_intensity[1] )
+            self.assertAlmostEqual( 1., chunk.TYPE )
+            self.assertAlmostEqual( 1., chunk.type )
+            self.assertAlmostEqual( 1.0, chunk.RIS[0] )
+            self.assertAlmostEqual( 2.0, chunk.RIS[1] )
+            self.assertAlmostEqual( 1.0, chunk.internal_pair_formation_coefficient[0] )
+            self.assertAlmostEqual( 2.0, chunk.internal_pair_formation_coefficient[1] )
+            self.assertAlmostEqual( 7.010000e+5, chunk.RICC[0] )
+            self.assertAlmostEqual( 1.106180e-2, chunk.RICC[1] )
+            self.assertAlmostEqual( 7.010000e+5, chunk.total_internal_conversion_coefficient[0] )
+            self.assertAlmostEqual( 1.106180e-2, chunk.total_internal_conversion_coefficient[1] )
+            self.assertAlmostEqual( 0.0, chunk.RICK[0] )
+            self.assertAlmostEqual( 0.0, chunk.RICK[1] )
+            self.assertAlmostEqual( 0.0, chunk.internal_conversion_coefficient_k_shell[0] )
+            self.assertAlmostEqual( 0.0, chunk.internal_conversion_coefficient_k_shell[1] )
+            self.assertAlmostEqual( 0.0, chunk.RICL[0] )
+            self.assertAlmostEqual( 0.0, chunk.RICL[1] )
+            self.assertAlmostEqual( 0.0, chunk.internal_conversion_coefficient_l_shell[0] )
+            self.assertAlmostEqual( 0.0, chunk.internal_conversion_coefficient_l_shell[1] )
+
+            self.assertEqual( 3, chunk.NC )
+
+            # verify string
+            self.assertEqual( self.chunk_styp_8, chunk.to_string( 3580, 8, 457 ) )
+
         # the data is given explicitly
         chunk = DiscreteSpectrum( energy = [ 4.863000e+4, 5.000000e+1 ],
                                   chain = 3., type = 0.,
@@ -195,6 +237,25 @@ def verify_chunk_styp_0( self, chunk ) :
 
         verify_chunk_styp_0( self, copy )
 
+        # the data is given explicitly
+        chunk = DiscreteSpectrum( energy = [ 4.863000e+4, 5.000000e+1 ],
+                                  chain = 3., type = 1.,
+                                  intensity = [ 1.420112e-6, 2.85306e-10 ],
+                                  ris = [ 1.0, 2.0 ],
+                                  ricc = [ 7.010000e+5, 1.106180e-2 ] )
+
+        verify_chunk_styp_8( self, chunk )
+
+        # the data is read from a string
+        chunk = DiscreteSpectrum.from_string( self.chunk_styp_8, 3580, 8, 457 )
+
+        verify_chunk_styp_8( self, chunk )
+
+        # the data is copied
+        copy = DiscreteSpectrum( chunk )
+
+        verify_chunk_styp_8( self, copy )
+
     def test_failures( self ) :
 
         print( '\n' )

src/ENDFtk/section/1/451.hpp

@@ -4,8 +4,6 @@
 // system includes
 
 // other includes
-#include "range/v3/iterator/operations.hpp"
-#include "range/v3/range/conversion.hpp"
 #include "range/v3/view/all.hpp"
 #include "range/v3/view/split.hpp"
 #include "ENDFtk/macros.hpp"

src/ENDFtk/section/14.hpp

@@ -5,7 +5,6 @@
 #include <variant>
 
 // other includes
-#include "range/v3/algorithm/count_if.hpp"
 #include "range/v3/view/transform.hpp"
 #include "ENDFtk/macros.hpp"
 #include "ENDFtk/ControlRecord.hpp"
@@ -82,7 +81,7 @@ namespace section {
      */
     int NI() const {
 
-      return ranges::count_if( this->photons_,
+      return std::count_if( this->photons_.begin(), this->photons_.end(),
                                [] ( const auto& entry )
                                   { return entry.index() == 0; } );
     }

src/ENDFtk/section/2/151.hpp

@@ -7,7 +7,6 @@
 #include <optional>
 
 // other includes
-#include "range/v3/algorithm/count.hpp"
 #include "range/v3/algorithm/find.hpp"
 #include "range/v3/range/conversion.hpp"
 #include "range/v3/view/all.hpp"

src/ENDFtk/section/2/151/RMatrixLimited/BackgroundChannels.hpp

@@ -50,11 +50,9 @@ class ENDFTK_PYTHON_EXPORT BackgroundChannels {
    */
   unsigned int KBK() const {
 
-    return ranges::cpp20::count( this->backgroundRMatrices()
-                                     | ranges::cpp20::views::transform(
-                                           [] ( const auto& value )
-                                              { return bool( value ); } ),
-                                 true );
+    auto matrices = this->backgroundRMatrices();
+    return std::count_if( matrices.begin(), matrices.end(),
+                          [] ( auto&& matrix ) { return matrix.has_value(); } );
   }
 
   /**

src/ENDFtk/section/2/151/RMatrixLimited/TabulatedBackgroundRMatrix/src/verify.hpp

@@ -1,10 +1,7 @@
 template < typename Range >
 static bool compare( const Range& left, const Range& right ) {
 
-  return ranges::cpp20::count(
-               ranges::views::zip_with(
-                   [] ( double left, double right ) { return left == right; },
-                   left, right ), true ) == ranges::cpp20::distance( left );
+  return std::equal( left.begin(), left.end(), right.begin(), right.end() );
 }
 
 static void

src/ENDFtk/section/26.hpp

@@ -4,7 +4,6 @@
 // system includes
 
 // other includes
-#include "range/v3/algorithm/find_if_not.hpp"
 #include "range/v3/range/conversion.hpp"
 #include "range/v3/view/all.hpp"
 #include "range/v3/view/chunk.hpp"

src/ENDFtk/section/26/ContinuumEnergyAngle/src/verifyLANG.hpp

@@ -5,14 +5,13 @@ verifyLANG( int LANG, const Array& sequence ) {
   auto compare = [LANG] ( const auto& entry )
                         { return entry.LANG() == LANG; };
 
-  auto iter = ranges::cpp20::find_if_not( sequence, compare );
+  auto iter = std::find_if_not( sequence.begin(), sequence.end(), compare );
 
-  if ( iter != ranges::cpp20::end( sequence ) ) {
+  if ( iter != sequence.end() ) {
 
     Log::error( "All subsections must use the same LANG format option" );
     Log::info( "Expected LANG={} for the subsection with index={}", LANG,
-               ranges::cpp20::distance(
-                   ranges::cpp20::begin( sequence ), iter ) );
+               ranges::cpp20::distance( sequence.begin(), iter ) );
     throw std::exception();
   }
 }

src/ENDFtk/section/6.hpp

@@ -5,7 +5,6 @@
 #include <complex>
 
 // other includes
-#include "range/v3/algorithm/find_if_not.hpp"
 #include "range/v3/range/conversion.hpp"
 #include "range/v3/view/all.hpp"
 #include "range/v3/view/chunk.hpp"

src/ENDFtk/section/6/ContinuumEnergyAngle/src/verifyLANG.hpp

@@ -7,14 +7,13 @@ verifyLANG( int LANG, const Array& sequence ) {
                                                 { return entry.LANG(); },
                                              variant ) == LANG; };
 
-  auto iter = ranges::cpp20::find_if_not( sequence, compare );
+  auto iter = std::find_if_not( sequence.begin(), sequence.end(), compare );
 
-  if ( iter != ranges::cpp20::end( sequence ) ) {
+  if ( iter != sequence.end() ) {
 
     Log::error( "All subsections must use the same LANG format option" );
     Log::info( "Expected LANG={} for the subsection with index={}", LANG,
-               ranges::cpp20::distance(
-                   ranges::cpp20::begin( sequence ), iter ) );
+               ranges::cpp20::distance( sequence.begin(), iter ) );
     throw std::exception();
   }
 }

src/ENDFtk/section/7/2.hpp

@@ -6,7 +6,6 @@
 #include <vector>
 
 // other includes
-#include "range/v3/algorithm/find_if_not.hpp"
 #include "range/v3/view/concat.hpp"
 #include "range/v3/view/single.hpp"
 #include "range/v3/view/transform.hpp"
@@ -22,7 +21,7 @@ namespace ENDFtk {
 namespace section{
 
   template<>
-  class ENDFTK_PYTHON_EXPORT Type< 7, 2 > : 
+  class ENDFTK_PYTHON_EXPORT Type< 7, 2 > :
     protected BaseWithoutMT< Type< 7, 2 > > {
 
     friend BaseWithoutMT< Type< 7, 2 > >;

src/ENDFtk/section/7/2/CoherentElastic/src/verifyTemperatures.hpp

@@ -5,14 +5,13 @@ verifyTemperatures( int NP, const Array& temperatures ) {
   auto compare = [NP] ( const auto& entry )
                       { return entry.NPL() == NP; };
 
-  auto iter = ranges::cpp20::find_if_not( temperatures, compare );
+  auto iter = std::find_if_not( temperatures.begin(), temperatures.end(), compare );
 
-  if ( iter != ranges::cpp20::end( temperatures ) ) {
+  if ( iter != temperatures.end() ) {
 
     Log::error( "All S(E,T) functions must have the same size" );
     Log::info( "Expected NP={} for the temperature with index={}", NP,
-               ranges::cpp20::distance(
-                   ranges::cpp20::begin( temperatures ), iter ) );
+               ranges::cpp20::distance( temperatures.begin(), iter ) );
     throw std::exception();
   }
 }

src/ENDFtk/section/8/457/DiscreteSpectrum/src/ctor.hpp

@@ -53,6 +53,30 @@ DiscreteSpectrum( double chain,
                                   intensity,
                                   ris ) ) ) {}
 
+/**
+ *  @brief Constructor
+ *
+ *  @param[in] chain        the ENDF decay chain
+ *  @param[in] energy       the energy and uncertainty
+ *  @param[in] intensity    the relative intensity and uncertainty
+ *  @param[in] type         the transition type
+ *  @param[in] ris          internal pair formation coefficient or
+ *                          positron intensity
+ *  @param[in] ricc         the total internal conversion coefficient
+ */
+DiscreteSpectrum( double chain,
+                  const std::array< double, 2 >& energy,
+                  const std::array< double, 2 >& intensity,
+                  double type,
+                  const std::array< double, 2 >& ris,
+                  const std::array< double, 2 >& ricc ) :
+  // this can never fail, try-catch would be unreachable
+  DiscreteSpectrum(
+        ListRecord( energy[0], energy[1], 0, 0, 0,
+                    generateList( chain, type,
+                                  intensity,
+                                  ris, ricc ) ) ) {}
+
 /**
  *  @brief Constructor
  *

src/ENDFtk/section/8/457/DiscreteSpectrum/src/generateList.hpp

@@ -13,6 +13,16 @@ generateList( double rtyp, double type,
   return { rtyp, type, ri[0], ri[1], ris[0], ris[1] };
 }
 
+static std::vector< double >
+generateList( double rtyp, double type,
+              const std::array< double, 2 >& ri,
+              const std::array< double, 2 >& ris,
+              const std::array< double, 2 >& ricc ) {
+
+  return { rtyp, type, ri[0], ri[1], ris[0], ris[1],
+           ricc[0], ricc[1] };
+}
+
 static std::vector< double >
 generateList( double rtyp, double type,
               const std::array< double, 2 >& ri,

src/ENDFtk/section/8/457/DiscreteSpectrum/src/returnArray.hpp

@@ -1,9 +1,12 @@
 std::array< double, 2 >
 returnArray( int index ) const {
+
   if ( this->NT() > index ) {
+
     return {{ ListRecord::list()[index], ListRecord::list()[index+1] }};
   }
   else {
+
     return {{ 0.0, 0.0 }};
   }
-}  
+}

src/ENDFtk/section/8/457/DiscreteSpectrum/src/verifyNT.hpp

@@ -1,10 +1,10 @@
 static void
 verifyNT( int NT ) {
 
-  if ( ( NT != 4 ) and ( NT != 6 ) and ( NT != 12 ) ) {
+  if ( ( NT != 4 ) && ( NT != 6 ) && ( NT != 8 ) && ( NT != 12 ) ) {
 
     Log::error( "Encountered illegal NT value" );
-    Log::info( "NT can only be equal to 4, 6 or 12" );
+    Log::info( "NT can only be equal to 4, 6, 8 or 12" );
     Log::info( "NT value: {}", NT );
     throw std::exception();
   }