Allowing XS plotting module to handle excitation reactions

tools/ALARAJOYWrapper/README.md

@@ -18,7 +18,11 @@ This preprocessor uses [NJOY 2016](https://github.com/njoy/NJOY2016) Nuclear Dat
   * [Warnings](https://docs.python.org/3/library/warnings.html)
 - Generic Python packages
   * [Matplotlib.pyplot](https://matplotlib.org/3.5.3/api/_as_gen/matplotlib.pyplot.html)
-  * [NumPy](https://numpy.org/install/)
+  * [NumPy](https://numpy.org/install/) ([numpy.f2py](https://numpy.org/doc/stable/f2py/) is used for Fortran interfacing, which depends on the following):
+    * C-compiler (i.e. [gcc](https://hprc.tamu.edu/kb/Software/GNU-Compiler-Collection/#gcc-versions))
+    * Fortran-compiler (i.e. [gfortran](https://fortran-lang.org/learn/os_setup/install_gfortran/))
+    * [Meson](https://mesonbuild.com/)
+    * [Ninja](https://github.com/ninja-build/ninja.git)
   * [Pandas](https://pandas.pydata.org/docs/getting_started/install.html)
   * [PyYAML](https://pyyaml.org/wiki/PyYAMLDocumentation) (only needed if running `xs_plotting.py` as an independent script.)
 - Domain-specific packages

tools/ALARAJOYWrapper/njoy_endf_wrapper.f90

@@ -0,0 +1,59 @@
+! ===========================
+! This lightweight NJOY wrapper accesses the NJOY ENDF module's utility
+! function endf.terpa() to make use of NJOY's built-in interpolation 
+! functionality to apply the appropriate interpolation schemes. The
+! designation for these schemes are encoded within the TAB1 record header 
+! variable INT, representing one of the following interpolation schemes, as
+! laid out in the ENDF-6 manual
+! (https://www.nndc.bnl.gov/endfdocs/ENDF-102-2023.pdf), under section 0.5.2
+! ("Interpolation Laws"):
+!
+!     INT   | Interpolation Scheme
+!     1     | y is constant in x (constant, histogram)
+!     2     | y is linear in x (linear-linear)
+!     3     | y is linear in ln(x) (linear-log)
+!     4     | ln(y) is linear in x (log-linear)
+!     5     | ln(y) is linear in ln(x) (log-log)
+!     6     | special one-dimensional interpolation law, used for charged-
+!           | particle cross sections only
+!     11-25 | method of corresponding points (follow interpolation laws 1-5)
+!     21-25 | unit base interpolation (follow interpolation laws of 1-5)
+!
+! This module's incorporation and usage within ALARAJOYWrapper is managed by
+! njoy_tools.import_njoy_endf_wrapper(), which conditionally compiles this
+! Fortran file to an executable using numpy.f2py (if such an executable has
+! not already been created) and importing njoy_endf_wrapper as a Python
+! package. This allows the subroutine interpolate_tab1() to be callable within
+! ALARAJOYWrapper, which is necessary for the construction of pathway-specific
+! reaction cross-section from MF9 multiplicities multiplied by MF3 cumulative
+! cross-sections (see xs_plotting.extract_continuous_data() for specific use-
+! case implementation).
+! ===========================
+
+module njoy_endf_wrapper
+
+   use endf
+   implicit none
+
+contains
+
+   subroutine interpolate_tab1(tab1, x, y)
+
+      real(kind=8), intent(in) :: tab1(:)
+      real(kind=8), intent(in) :: x(:)
+      real(kind=8), intent(out) :: y(size(x))
+
+      integer :: i
+      integer :: ip, ir, idis
+      real(kind=8) :: xnext
+
+      ip = 2
+      ir = 1
+
+      do i = 1, size(x)
+         call terpa(y(i), x(i), xnext, idis, tab1, ip, ir)
+      end do
+
+   end subroutine interpolate_tab1
+
+end module njoy_endf_wrapper

tools/ALARAJOYWrapper/njoy_tools.py

@@ -4,6 +4,8 @@
 from pathlib import Path
 import re
 import numpy as np
+from sys import executable
+from shutil import which
 
 def set_directory():
     '''
@@ -699,4 +701,49 @@ def cleanup_njoy_files(element, A):
 
     output_dir = dir / 'njoy_outputs'
     output_dir.mkdir(exist_ok=True)
-    Path('output').rename(output_dir / f'njoy_output_{element}{A}.out')
+    Path('output').rename(output_dir / f'njoy_output_{element}{A}.out')
+
+def import_njoy_endf_wrapper():
+    """
+    Import the njoy_endf_wrapper module defined by the targeted NJOY-wrapper
+        njoy_endf_wrapper.f90. This wrapper utilizes the NJOY function
+        `endf.terpa()`, which interpolates TAB1 data according to the encoded
+        interpolation scheme(s). If the module cannot be accessed, compile
+        njoy_endf_wrapper.f90 to a CPython executable using NumPy.f2py and
+        subsequently import the module.
+
+    Arguments:
+        None
+
+    Returns:
+        njoy_endf_wrapper.njoy_endf_wrapper (fortran object): Python module of
+            the compiled njoy_endf_wrapper.f90 NJOY wrapper containing the
+            subroutine `interpolate_tab1()`, which can be used to apply
+            `endf.terpa()` to interpolate TAB1 according to the encoded
+            interpolation scheme(s).
+    """
+
+    try:
+        import njoy_endf_wrapper
+
+    except (ModuleNotFoundError, ImportError):
+        njoy_dir = Path(which('njoy')).parent
+        subprocess.run(
+            [
+                executable,
+                '-m',
+                'numpy.f2py',
+                '-c',
+                '-m',
+                'njoy_endf_wrapper',
+                str(Path(__file__).parent / 'njoy_endf_wrapper.f90'),
+                f'-I{njoy_dir / "fortran_modules"}',
+                f'-L{njoy_dir}',
+                '-lnjoy'
+            ],
+            check=True,
+            cwd=Path.cwd()
+        )
+        import njoy_endf_wrapper
+
+    return njoy_endf_wrapper.njoy_endf_wrapper

tools/ALARAJOYWrapper/preprocess_fendl3.py

@@ -11,6 +11,7 @@
 from pathlib import Path
 from collections import defaultdict
 from subprocess import TimeoutExpired
+from endf_parserpy import EndfParserPy
 
 def make_argparser():
     parser = argparse.ArgumentParser()
@@ -482,6 +483,9 @@ def store_results(
         dsv.write(f'{nGroups} {group_name}\n')
         for parent in sorted(all_rxns):
             element, A = tp.interpret_KZA(parent)
+            endf_dict = EndfParserPy().parsefile(
+                tendl_dir / f'{element}{A}.tendl'
+            )
             for daughter in all_rxns[parent]:
                 if parent != daughter:
                     for MT, rxn in all_rxns[parent][daughter].items():
@@ -501,8 +505,7 @@ def store_results(
                                 )
 
                                 continuous_dict = xp.extract_continuous_data(
-                                    tendl_dir / f'{element}{A}.tendl',
-                                    xp.flagged_num_to_int(MT)
+                                    endf_dict, MT
                                 )
 
                                 energies = njt.load_external_group_struct(
@@ -521,7 +524,7 @@ def store_results(
                                 plot_path = xp.set_plot_save_path(
                                     element, A, emitted, tendl_dir, group_name
                                 )
-                                
+
                                 plt.savefig(plot_path)
 
         if plotting:

tools/ALARAJOYWrapper/tendl_processing.py

@@ -23,6 +23,7 @@
     for val in arr
 }
 ISOMERIC_STATES = 'mnopqrstuvwxyz'
+PATH_SPECIFIC_MFS = (9,10)
 
 def parse_endf_file_level_data(endf_path, MF=3, endf_format='endf6-ext'):
     """
@@ -187,10 +188,9 @@ def determine_all_excitations(endf_path, MTs):
 
     isomer_dict = defaultdict(lambda: defaultdict(list))
 
-    path_specific_MFs = (9,10)
     mf_dict = {
         MF: parse_endf_file_level_data(endf_path, MF)[0]
-        for MF in path_specific_MFs
+        for MF in PATH_SPECIFIC_MFS
     }
 
     for MT in MTs:
@@ -199,7 +199,7 @@ def determine_all_excitations(endf_path, MTs):
             # Isomer pathways contained either in MF 9 ("Multiplicities for
             # Production of Radioactive Nuclides") and MF 10 ("Cross Sections
             # for Production of Radioactive Nuclides").
-            for MF in path_specific_MFs:
+            for MF in PATH_SPECIFIC_MFS:
                 pathways = mf_dict[MF].get(MT, {}).get('subsection', {})
                 for pathway_data in pathways.values():
                     isomer_dict[MT][MF].append(pathway_data['LFS'])