1from collections.abc import MutableSequence, Iterable
2import io
3
4import numpy as np
5from numpy.polynomial import Polynomial
6import pandas as pd
7
8import openmc.checkvalue as cv
9from .data import NEUTRON_MASS
10from .endf import get_head_record, get_cont_record, get_tab1_record, get_list_record
11try:
12    from .reconstruct import wave_number, penetration_shift, reconstruct_mlbw, \
13        reconstruct_slbw, reconstruct_rm
14    _reconstruct = True
15except ImportError:
16    _reconstruct = False
17
18
19class Resonances:
20    """Resolved and unresolved resonance data
21
22    Parameters
23    ----------
24    ranges : list of openmc.data.ResonanceRange
25        Distinct energy ranges for resonance data
26
27    Attributes
28    ----------
29    ranges : list of openmc.data.ResonanceRange
30        Distinct energy ranges for resonance data
31    resolved : openmc.data.ResonanceRange or None
32        Resolved resonance range
33    unresolved : openmc.data.Unresolved or None
34        Unresolved resonance range
35
36    """
37
38    def __init__(self, ranges):
39        self.ranges = ranges
40
41    def __iter__(self):
42        for r in self.ranges:
43            yield r
44
45    @property
46    def ranges(self):
47        return self._ranges
48
49    @property
50    def resolved(self):
51        resolved_ranges = [r for r in self.ranges
52                           if not isinstance(r, Unresolved)]
53        if len(resolved_ranges) > 1:
54            raise ValueError('More than one resolved range present')
55        elif len(resolved_ranges) == 0:
56            return None
57        else:
58            return resolved_ranges[0]
59
60    @property
61    def unresolved(self):
62        for r in self.ranges:
63            if isinstance(r, Unresolved):
64                return r
65        else:
66            return None
67
68    @ranges.setter
69    def ranges(self, ranges):
70        cv.check_type('resonance ranges', ranges, MutableSequence)
71        self._ranges = cv.CheckedList(ResonanceRange, 'resonance ranges',
72                                      ranges)
73
74    @classmethod
75    def from_endf(cls, ev):
76        """Generate resonance data from an ENDF evaluation.
77
78        Parameters
79        ----------
80        ev : openmc.data.endf.Evaluation
81            ENDF evaluation
82
83        Returns
84        -------
85        openmc.data.Resonances
86            Resonance data
87
88        """
89        file_obj = io.StringIO(ev.section[2, 151])
90
91        # Determine whether discrete or continuous representation
92        items = get_head_record(file_obj)
93        n_isotope = items[4]  # Number of isotopes
94
95        ranges = []
96        for iso in range(n_isotope):
97            items = get_cont_record(file_obj)
98            abundance = items[1]
99            fission_widths = (items[3] == 1)  # fission widths are given?
100            n_ranges = items[4]  # number of resonance energy ranges
101
102            for j in range(n_ranges):
103                items = get_cont_record(file_obj)
104                resonance_flag = items[2]  # flag for resolved (1)/unresolved (2)
105                formalism = items[3]  # resonance formalism
106
107                if resonance_flag in (0, 1):
108                    # resolved resonance region
109                    erange = _FORMALISMS[formalism].from_endf(ev, file_obj, items)
110
111                elif resonance_flag == 2:
112                    # unresolved resonance region
113                    erange = Unresolved.from_endf(file_obj, items, fission_widths)
114
115                # erange.material = self
116                ranges.append(erange)
117
118        return cls(ranges)
119
120
121class ResonanceRange:
122    """Resolved resonance range
123
124    Parameters
125    ----------
126    target_spin : float
127        Intrinsic spin, :math:`I`, of the target nuclide
128    energy_min : float
129        Minimum energy of the resolved resonance range in eV
130    energy_max : float
131        Maximum energy of the resolved resonance range in eV
132    channel : dict
133        Dictionary whose keys are l-values and values are channel radii as a
134        function of energy
135    scattering : dict
136        Dictionary whose keys are l-values and values are scattering radii as a
137        function of energy
138
139    Attributes
140    ----------
141    channel_radius : dict
142        Dictionary whose keys are l-values and values are channel radii as a
143        function of energy
144    energy_max : float
145        Maximum energy of the resolved resonance range in eV
146    energy_min : float
147        Minimum energy of the resolved resonance range in eV
148    scattering_radius : dict
149        Dictionary whose keys are l-values and values are scattering radii as a
150        function of energ
151    target_spin : float
152        Intrinsic spin, :math:`I`, of the target nuclide
153
154    """
155    def __init__(self, target_spin, energy_min, energy_max, channel, scattering):
156        self.target_spin = target_spin
157        self.energy_min = energy_min
158        self.energy_max = energy_max
159        self.channel_radius = channel
160        self.scattering_radius = scattering
161
162        self._prepared = False
163        self._parameter_matrix = {}
164
165    def __copy__(self):
166        cls = type(self)
167        new_copy = cls.__new__(cls)
168        new_copy.__dict__.update(self.__dict__)
169        new_copy._prepared = False
170        return new_copy
171
172    @classmethod
173    def from_endf(cls, ev, file_obj, items):
174        """Create resonance range from an ENDF evaluation.
175
176        This factory method is only used when LRU=0, indicating that only a
177        scattering radius appears in MF=2, MT=151. All subclasses of
178        ResonanceRange override this method with their own.
179
180        Parameters
181        ----------
182        ev : openmc.data.endf.Evaluation
183            ENDF evaluation
184        file_obj : file-like object
185            ENDF file positioned at the second record of a resonance range
186            subsection in MF=2, MT=151
187        items : list
188            Items from the CONT record at the start of the resonance range
189            subsection
190
191        Returns
192        -------
193        openmc.data.ResonanceRange
194            Resonance range data
195
196        """
197        energy_min, energy_max = items[0:2]
198
199        # For scattering radius-only, NRO must be zero
200        assert items[4] == 0
201
202        # Get energy-independent scattering radius
203        items = get_cont_record(file_obj)
204        target_spin = items[0]
205        ap = Polynomial((items[1],))
206
207        # Calculate channel radius from ENDF-102 equation D.14
208        a = Polynomial((0.123 * (NEUTRON_MASS*ev.target['mass'])**(1./3.) + 0.08,))
209
210        return cls(target_spin, energy_min, energy_max, {0: a}, {0: ap})
211
212    def reconstruct(self, energies):
213        """Evaluate cross section at specified energies.
214
215        Parameters
216        ----------
217        energies : float or Iterable of float
218            Energies at which the cross section should be evaluated
219
220        Returns
221        -------
222        3-tuple of float or numpy.ndarray
223            Elastic, capture, and fission cross sections at the specified
224            energies
225
226        """
227        if not _reconstruct:
228            raise RuntimeError("Resonance reconstruction not available.")
229
230        # Pre-calculate penetrations and shifts for resonances
231        if not self._prepared:
232            self._prepare_resonances()
233
234        if isinstance(energies, Iterable):
235            elastic = np.zeros_like(energies)
236            capture = np.zeros_like(energies)
237            fission = np.zeros_like(energies)
238
239            for i, E in enumerate(energies):
240                xse, xsg, xsf = self._reconstruct(self, E)
241                elastic[i] = xse
242                capture[i] = xsg
243                fission[i] = xsf
244        else:
245            elastic, capture, fission = self._reconstruct(self, energies)
246
247        return {2: elastic, 102: capture, 18: fission}
248
249
250class MultiLevelBreitWigner(ResonanceRange):
251    """Multi-level Breit-Wigner resolved resonance formalism data.
252
253    Multi-level Breit-Wigner resolved resonance data is identified by LRF=2 in
254    the ENDF-6 format.
255
256    Parameters
257    ----------
258    target_spin : float
259        Intrinsic spin, :math:`I`, of the target nuclide
260    energy_min : float
261        Minimum energy of the resolved resonance range in eV
262    energy_max : float
263        Maximum energy of the resolved resonance range in eV
264    channel : dict
265        Dictionary whose keys are l-values and values are channel radii as a
266        function of energy
267    scattering : dict
268        Dictionary whose keys are l-values and values are scattering radii as a
269        function of energy
270
271    Attributes
272    ----------
273    atomic_weight_ratio : float
274        Atomic weight ratio of the target nuclide given as a function of
275        l-value. Note that this may be different than the value for the
276        evaluation as a whole.
277    channel_radius : dict
278        Dictionary whose keys are l-values and values are channel radii as a
279        function of energy
280    energy_max : float
281        Maximum energy of the resolved resonance range in eV
282    energy_min : float
283        Minimum energy of the resolved resonance range in eV
284    parameters : pandas.DataFrame
285        Energies, spins, and resonances widths for each resonance
286    q_value : dict
287        Q-value to be added to incident particle's center-of-mass energy to
288        determine the channel energy for use in the penetrability factor. The
289        keys of the dictionary are l-values.
290    scattering_radius : dict
291        Dictionary whose keys are l-values and values are scattering radii as a
292        function of energy
293    target_spin : float
294        Intrinsic spin, :math:`I`, of the target nuclide
295
296    """
297
298    def __init__(self, target_spin, energy_min, energy_max, channel, scattering):
299        super().__init__(target_spin, energy_min, energy_max, channel,
300                         scattering)
301        self.parameters = None
302        self.q_value = {}
303        self.atomic_weight_ratio = None
304
305        # Set resonance reconstruction function
306        if _reconstruct:
307            self._reconstruct = reconstruct_mlbw
308        else:
309            self._reconstruct = None
310
311    @classmethod
312    def from_endf(cls, ev, file_obj, items):
313        """Create MLBW data from an ENDF evaluation.
314
315        Parameters
316        ----------
317        ev : openmc.data.endf.Evaluation
318            ENDF evaluation
319        file_obj : file-like object
320            ENDF file positioned at the second record of a resonance range
321            subsection in MF=2, MT=151
322        items : list
323            Items from the CONT record at the start of the resonance range
324            subsection
325
326        Returns
327        -------
328        openmc.data.MultiLevelBreitWigner
329            Multi-level Breit-Wigner resonance parameters
330
331        """
332
333        # Read energy-dependent scattering radius if present
334        energy_min, energy_max = items[0:2]
335        nro, naps = items[4:6]
336        if nro != 0:
337            params, ape = get_tab1_record(file_obj)
338
339        # Other scatter radius parameters
340        items = get_cont_record(file_obj)
341        target_spin = items[0]
342        ap = Polynomial((items[1],))  # energy-independent scattering-radius
343        NLS = items[4]  # number of l-values
344
345        # Read resonance widths, J values, etc
346        channel_radius = {}
347        scattering_radius = {}
348        q_value = {}
349        records = []
350        for l in range(NLS):
351            items, values = get_list_record(file_obj)
352            l_value = items[2]
353            awri = items[0]
354            q_value[l_value] = items[1]
355            competitive = items[3]
356
357            # Calculate channel radius from ENDF-102 equation D.14
358            a = Polynomial((0.123 * (NEUTRON_MASS*awri)**(1./3.) + 0.08,))
359
360            # Construct scattering and channel radius
361            if nro == 0:
362                scattering_radius[l_value] = ap
363                if naps == 0:
364                    channel_radius[l_value] = a
365                elif naps == 1:
366                    channel_radius[l_value] = ap
367            elif nro == 1:
368                scattering_radius[l_value] = ape
369                if naps == 0:
370                    channel_radius[l_value] = a
371                elif naps == 1:
372                    channel_radius[l_value] = ape
373                elif naps == 2:
374                    channel_radius[l_value] = ap
375
376            energy = values[0::6]
377            spin = values[1::6]
378            gt = np.asarray(values[2::6])
379            gn = np.asarray(values[3::6])
380            gg = np.asarray(values[4::6])
381            gf = np.asarray(values[5::6])
382            if competitive > 0:
383                gx = gt - (gn + gg + gf)
384            else:
385                gx = np.zeros_like(gt)
386
387            for i, E in enumerate(energy):
388                records.append([energy[i], l_value, spin[i], gt[i], gn[i],
389                                gg[i], gf[i], gx[i]])
390
391        columns = ['energy', 'L', 'J', 'totalWidth', 'neutronWidth',
392                   'captureWidth', 'fissionWidth', 'competitiveWidth']
393        parameters = pd.DataFrame.from_records(records, columns=columns)
394
395        # Create instance of class
396        mlbw = cls(target_spin, energy_min, energy_max,
397                   channel_radius, scattering_radius)
398        mlbw.q_value = q_value
399        mlbw.atomic_weight_ratio = awri
400        mlbw.parameters = parameters
401
402        return mlbw
403
404    def _prepare_resonances(self):
405        df = self.parameters.copy()
406
407        # Penetration and shift factors
408        p = np.zeros(len(df))
409        s = np.zeros(len(df))
410
411        # Penetration and shift factors for competitive reaction
412        px = np.zeros(len(df))
413        sx = np.zeros(len(df))
414
415        l_values = []
416        competitive = []
417
418        A = self.atomic_weight_ratio
419        for i, E, l, J, gt, gn, gg, gf, gx in df.itertuples():
420            if l not in l_values:
421                l_values.append(l)
422                competitive.append(gx > 0)
423
424            # Determine penetration and shift corresponding to resonance energy
425            k = wave_number(A, E)
426            rho = k*self.channel_radius[l](E)
427            rhohat = k*self.scattering_radius[l](E)
428            p[i], s[i] = penetration_shift(l, rho)
429
430            # Determine penetration at modified energy for competitive reaction
431            if gx > 0:
432                Ex = E + self.q_value[l]*(A + 1)/A
433                rho = k*self.channel_radius[l](Ex)
434                rhohat = k*self.scattering_radius[l](Ex)
435                px[i], sx[i] = penetration_shift(l, rho)
436            else:
437                px[i] = sx[i] = 0.0
438
439        df['p'] = p
440        df['s'] = s
441        df['px'] = px
442        df['sx'] = sx
443
444        self._l_values = np.array(l_values)
445        self._competitive = np.array(competitive)
446        for l in l_values:
447            self._parameter_matrix[l] = df[df.L == l].values
448
449        self._prepared = True
450
451
452class SingleLevelBreitWigner(MultiLevelBreitWigner):
453    """Single-level Breit-Wigner resolved resonance formalism data.
454
455    Single-level Breit-Wigner resolved resonance data is is identified by LRF=1
456    in the ENDF-6 format.
457
458    Parameters
459    ----------
460    target_spin : float
461        Intrinsic spin, :math:`I`, of the target nuclide
462    energy_min : float
463        Minimum energy of the resolved resonance range in eV
464    energy_max : float
465        Maximum energy of the resolved resonance range in eV
466    channel : dict
467        Dictionary whose keys are l-values and values are channel radii as a
468        function of energy
469    scattering : dict
470        Dictionary whose keys are l-values and values are scattering radii as a
471        function of energy
472
473    Attributes
474    ----------
475    atomic_weight_ratio : float
476        Atomic weight ratio of the target nuclide given as a function of
477        l-value. Note that this may be different than the value for the
478        evaluation as a whole.
479    channel_radius : dict
480        Dictionary whose keys are l-values and values are channel radii as a
481        function of energy
482    energy_max : float
483        Maximum energy of the resolved resonance range in eV
484    energy_min : float
485        Minimum energy of the resolved resonance range in eV
486    parameters : pandas.DataFrame
487        Energies, spins, and resonances widths for each resonance
488    q_value : dict
489        Q-value to be added to incident particle's center-of-mass energy to
490        determine the channel energy for use in the penetrability factor. The
491        keys of the dictionary are l-values.
492    scattering_radius : dict
493        Dictionary whose keys are l-values and values are scattering radii as a
494        function of energy
495    target_spin : float
496        Intrinsic spin, :math:`I`, of the target nuclide
497
498    """
499
500    def __init__(self, target_spin, energy_min, energy_max, channel, scattering):
501        super().__init__(target_spin, energy_min, energy_max, channel,
502                         scattering)
503
504        # Set resonance reconstruction function
505        if _reconstruct:
506            self._reconstruct = reconstruct_slbw
507        else:
508            self._reconstruct = None
509
510
511class ReichMoore(ResonanceRange):
512    """Reich-Moore resolved resonance formalism data.
513
514    Reich-Moore resolved resonance data is identified by LRF=3 in the ENDF-6
515    format.
516
517    Parameters
518    ----------
519    target_spin : float
520        Intrinsic spin, :math:`I`, of the target nuclide
521    energy_min : float
522        Minimum energy of the resolved resonance range in eV
523    energy_max : float
524        Maximum energy of the resolved resonance range in eV
525    channel : dict
526        Dictionary whose keys are l-values and values are channel radii as a
527        function of energy
528    scattering : dict
529        Dictionary whose keys are l-values and values are scattering radii as a
530        function of energy
531
532    Attributes
533    ----------
534    angle_distribution : bool
535        Indicate whether parameters can be used to compute angular distributions
536    atomic_weight_ratio : float
537        Atomic weight ratio of the target nuclide given as a function of
538        l-value. Note that this may be different than the value for the
539        evaluation as a whole.
540    channel_radius : dict
541        Dictionary whose keys are l-values and values are channel radii as a
542        function of energy
543    energy_max : float
544        Maximum energy of the resolved resonance range in eV
545    energy_min : float
546        Minimum energy of the resolved resonance range in eV
547    num_l_convergence : int
548        Number of l-values which must be used to converge the calculation
549    scattering_radius : dict
550        Dictionary whose keys are l-values and values are scattering radii as a
551        function of energy
552    parameters : pandas.DataFrame
553        Energies, spins, and resonances widths for each resonance
554    target_spin : float
555        Intrinsic spin, :math:`I`, of the target nuclide
556
557    """
558
559    def __init__(self, target_spin, energy_min, energy_max, channel, scattering):
560        super().__init__(target_spin, energy_min, energy_max, channel,
561                         scattering)
562        self.parameters = None
563        self.angle_distribution = False
564        self.num_l_convergence = 0
565
566        # Set resonance reconstruction function
567        if _reconstruct:
568            self._reconstruct = reconstruct_rm
569        else:
570            self._reconstruct = None
571
572    @classmethod
573    def from_endf(cls, ev, file_obj, items):
574        """Create Reich-Moore resonance data from an ENDF evaluation.
575
576        Parameters
577        ----------
578        ev : openmc.data.endf.Evaluation
579            ENDF evaluation
580        file_obj : file-like object
581            ENDF file positioned at the second record of a resonance range
582            subsection in MF=2, MT=151
583        items : list
584            Items from the CONT record at the start of the resonance range
585            subsection
586
587        Returns
588        -------
589        openmc.data.ReichMoore
590            Reich-Moore resonance parameters
591
592        """
593        # Read energy-dependent scattering radius if present
594        energy_min, energy_max = items[0:2]
595        nro, naps = items[4:6]
596        if nro != 0:
597            params, ape = get_tab1_record(file_obj)
598
599        # Other scatter radius parameters
600        items = get_cont_record(file_obj)
601        target_spin = items[0]
602        ap = Polynomial((items[1],))
603        angle_distribution = (items[3] == 1)  # Flag for angular distribution
604        NLS = items[4]  # Number of l-values
605        num_l_convergence = items[5]  # Number of l-values for convergence
606
607        # Read resonance widths, J values, etc
608        channel_radius = {}
609        scattering_radius = {}
610        records = []
611        for i in range(NLS):
612            items, values = get_list_record(file_obj)
613            apl = Polynomial((items[1],)) if items[1] != 0.0 else ap
614            l_value = items[2]
615            awri = items[0]
616
617            # Calculate channel radius from ENDF-102 equation D.14
618            a = Polynomial((0.123 * (NEUTRON_MASS*awri)**(1./3.) + 0.08,))
619
620            # Construct scattering and channel radius
621            if nro == 0:
622                scattering_radius[l_value] = apl
623                if naps == 0:
624                    channel_radius[l_value] = a
625                elif naps == 1:
626                    channel_radius[l_value] = apl
627            elif nro == 1:
628                if naps == 0:
629                    channel_radius[l_value] = a
630                    scattering_radius[l_value] = ape
631                elif naps == 1:
632                    channel_radius[l_value] = scattering_radius[l_value] = ape
633                elif naps == 2:
634                    channel_radius[l_value] = apl
635                    scattering_radius[l_value] = ape
636
637            energy = values[0::6]
638            spin = values[1::6]
639            gn = values[2::6]
640            gg = values[3::6]
641            gfa = values[4::6]
642            gfb = values[5::6]
643
644            for i, E in enumerate(energy):
645                records.append([energy[i], l_value, spin[i], gn[i], gg[i],
646                                gfa[i], gfb[i]])
647
648        # Create pandas DataFrame with resonance data
649        columns = ['energy', 'L', 'J', 'neutronWidth', 'captureWidth',
650                   'fissionWidthA', 'fissionWidthB']
651        parameters = pd.DataFrame.from_records(records, columns=columns)
652
653        # Create instance of ReichMoore
654        rm = cls(target_spin, energy_min, energy_max,
655                 channel_radius, scattering_radius)
656        rm.parameters = parameters
657        rm.angle_distribution = angle_distribution
658        rm.num_l_convergence = num_l_convergence
659        rm.atomic_weight_ratio = awri
660
661        return rm
662
663    def _prepare_resonances(self):
664        df = self.parameters.copy()
665
666        # Penetration and shift factors
667        p = np.zeros(len(df))
668        s = np.zeros(len(df))
669
670        l_values = []
671        lj_values = []
672
673        A = self.atomic_weight_ratio
674        for i, E, l, J, gn, gg, gfa, gfb in df.itertuples():
675            if l not in l_values:
676                l_values.append(l)
677            if (l, abs(J)) not in lj_values:
678                lj_values.append((l, abs(J)))
679
680            # Determine penetration and shift corresponding to resonance energy
681            k = wave_number(A, E)
682            rho = k*self.channel_radius[l](E)
683            rhohat = k*self.scattering_radius[l](E)
684            p[i], s[i] = penetration_shift(l, rho)
685
686        df['p'] = p
687        df['s'] = s
688
689        self._l_values = np.array(l_values)
690        for (l, J) in lj_values:
691            self._parameter_matrix[l, J] = df[(df.L == l) &
692                                              (abs(df.J) == J)].values
693
694        self._prepared = True
695
696
697class RMatrixLimited(ResonanceRange):
698    """R-matrix limited resolved resonance formalism data.
699
700    R-matrix limited resolved resonance data is identified by LRF=7 in the
701    ENDF-6 format.
702
703    Parameters
704    ----------
705    energy_min : float
706        Minimum energy of the resolved resonance range in eV
707    energy_max : float
708        Maximum energy of the resolved resonance range in eV
709    particle_pairs : list of dict
710        List of particle pairs. Each particle pair is represented by a
711        dictionary that contains the mass, atomic number, spin, and parity of
712        each particle as well as other characteristics.
713    spin_groups : list of dict
714        List of spin groups. Each spin group is characterized by channels,
715        resonance energies, and resonance widths.
716
717    Attributes
718    ----------
719    reduced_width : bool
720        Flag indicating whether channel widths in eV or reduced-width amplitudes
721        in eV^1/2 are given
722    formalism : int
723        Flag to specify which formulae for the R-matrix are to be used
724    particle_pairs : list of dict
725        List of particle pairs. Each particle pair is represented by a
726        dictionary that contains the mass, atomic number, spin, and parity of
727        each particle as well as other characteristics.
728    spin_groups : list of dict
729        List of spin groups. Each spin group is characterized by channels,
730        resonance energies, and resonance widths.
731
732    """
733
734    def __init__(self, energy_min, energy_max, particle_pairs, spin_groups):
735        super().__init__(0.0, energy_min, energy_max, None, None)
736        self.reduced_width = False
737        self.formalism = 3
738        self.particle_pairs = particle_pairs
739        self.spin_groups = spin_groups
740
741    @classmethod
742    def from_endf(cls, ev, file_obj, items):
743        """Read R-Matrix limited resonance data from an ENDF evaluation.
744
745        Parameters
746        ----------
747        ev : openmc.data.endf.Evaluation
748            ENDF evaluation
749        file_obj : file-like object
750            ENDF file positioned at the second record of a resonance range
751            subsection in MF=2, MT=151
752        items : list
753            Items from the CONT record at the start of the resonance range
754            subsection
755
756        Returns
757        -------
758        openmc.data.RMatrixLimited
759            R-matrix limited resonance parameters
760
761        """
762        energy_min, energy_max = items[0:2]
763
764        items = get_cont_record(file_obj)
765        reduced_width = (items[2] == 1)  # reduced width amplitude?
766        formalism = items[3]  # Specify which formulae are used
767        n_spin_groups = items[4]  # Number of Jpi values (NJS)
768
769        particle_pairs = []
770        spin_groups = []
771
772        items, values = get_list_record(file_obj)
773        n_pairs = items[5]//2  # Number of particle pairs (NPP)
774        for i in range(n_pairs):
775            first = {'mass': values[12*i],
776                     'z': int(values[12*i + 2]),
777                     'spin': values[12*i + 4],
778                     'parity': values[12*i + 10]}
779            second = {'mass': values[12*i + 1],
780                      'z': int(values[12*i + 3]),
781                      'spin': values[12*i + 5],
782                      'parity': values[12*i + 11]}
783
784            q_value = values[12*i + 6]
785            penetrability = values[12*i + 7]
786            shift = values[12*i + 8]
787            mt = int(values[12*i + 9])
788
789            particle_pairs.append(ParticlePair(
790                first, second, q_value, penetrability, shift, mt))
791
792        # loop over spin groups
793        for i in range(n_spin_groups):
794            items, values = get_list_record(file_obj)
795            J = items[0]
796            if J == 0.0:
797                parity = '+' if items[1] == 1.0 else '-'
798            else:
799                parity = '+' if J > 0. else '-'
800                J = abs(J)
801            kbk = items[2]
802            kps = items[3]
803            n_channels = items[5]
804            channels = []
805            for j in range(n_channels):
806                channel = {}
807                channel['particle_pair'] = particle_pairs[
808                    int(values[6*j]) - 1]
809                channel['l'] = values[6*j + 1]
810                channel['spin'] = values[6*j + 2]
811                channel['boundary'] = values[6*j + 3]
812                channel['effective_radius'] = values[6*j + 4]
813                channel['true_radius'] = values[6*j + 5]
814                channels.append(channel)
815
816            # Read resonance energies and widths
817            items, values = get_list_record(file_obj)
818            n_resonances = items[3]
819            records = []
820            m = n_channels//6 + 1
821            for j in range(n_resonances):
822                energy = values[6*m*j]
823                records.append([energy] + [values[6*m*j + k + 1]
824                                           for k in range(n_channels)])
825
826            # Determine column names
827            columns = ['energy']
828            for channel in channels:
829                mt = channel['particle_pair'].mt
830                if mt == 2:
831                    columns.append('neutronWidth')
832                elif mt == 18:
833                    columns.append('fissionWidth')
834                elif mt == 102:
835                    columns.append('captureWidth')
836                else:
837                    columns.append('width (MT={})'.format(mt))
838
839            # Create Pandas dataframe with resonance parameters
840            parameters = pd.DataFrame.from_records(records, columns=columns)
841
842            # Construct SpinGroup instance and add to list
843            sg = SpinGroup(J, parity, channels, parameters)
844            spin_groups.append(sg)
845
846            # Optional extension (Background R-Matrix)
847            if kbk > 0:
848                items, values = get_list_record(file_obj)
849                lbk = items[4]
850                if lbk == 1:
851                    params, rbr = get_tab1_record(file_obj)
852                    params, rbi = get_tab1_record(file_obj)
853
854            # Optional extension (Tabulated phase shifts)
855            if kps > 0:
856                items, values = get_list_record(file_obj)
857                lps = items[4]
858                if lps == 1:
859                    params, psr = get_tab1_record(file_obj)
860                    params, psi = get_tab1_record(file_obj)
861
862        rml = cls(energy_min, energy_max, particle_pairs, spin_groups)
863        rml.reduced_width = reduced_width
864        rml.formalism = formalism
865
866        return rml
867
868
869class ParticlePair:
870    def __init__(self, first, second, q_value, penetrability,
871                 shift, mt):
872        self.first = first
873        self.second = second
874        self.q_value = q_value
875        self.penetrability = penetrability
876        self.shift = shift
877        self.mt = mt
878
879
880class SpinGroup:
881    """Resonance spin group
882
883    Attributes
884    ----------
885    spin : float
886        Total angular momentum (nuclear spin)
887    parity : {'+', '-'}
888        Even (+) or odd(-) parity
889    channels : list of openmc.data.Channel
890        Available channels
891    parameters : pandas.DataFrame
892        Energies/widths for each resonance/channel
893
894    """
895
896    def __init__(self, spin, parity, channels, parameters):
897        self.spin = spin
898        self.parity = parity
899        self.channels = channels
900        self.parameters = parameters
901
902    def __repr__(self):
903        return '<SpinGroup: Jpi={}{}>'.format(self.spin, self.parity)
904
905
906class Unresolved(ResonanceRange):
907    """Unresolved resonance parameters as identified by LRU=2 in MF=2.
908
909    Parameters
910    ----------
911    target_spin : float
912        Intrinsic spin, :math:`I`, of the target nuclide
913    energy_min : float
914        Minimum energy of the unresolved resonance range in eV
915    energy_max : float
916        Maximum energy of the unresolved resonance range in eV
917    channel : openmc.data.Function1D
918        Channel radii as a function of energy
919    scattering : openmc.data.Function1D
920        Scattering radii as a function of energy
921
922    Attributes
923    ----------
924    add_to_background : bool
925        If True, file 3 contains partial cross sections to be added to the
926        average unresolved cross sections calculated from parameters.
927    atomic_weight_ratio : float
928        Atomic weight ratio of the target nuclide
929    channel_radius : openmc.data.Function1D
930        Channel radii as a function of energy
931    energies : Iterable of float
932        Energies at which parameters are tabulated
933    energy_max : float
934        Maximum energy of the unresolved resonance range in eV
935    energy_min : float
936        Minimum energy of the unresolved resonance range in eV
937    parameters : list of pandas.DataFrame
938        Average resonance parameters at each energy
939    scattering_radius : openmc.data.Function1D
940        Scattering radii as a function of energy
941    target_spin : float
942        Intrinsic spin, :math:`I`, of the target nuclide
943
944    """
945
946    def __init__(self, target_spin, energy_min, energy_max, channel, scattering):
947        super().__init__(target_spin, energy_min, energy_max, channel,
948                         scattering)
949        self.energies = None
950        self.parameters = None
951        self.add_to_background = False
952        self.atomic_weight_ratio = None
953
954    @classmethod
955    def from_endf(cls, file_obj, items, fission_widths):
956        """Read unresolved resonance data from an ENDF evaluation.
957
958        Parameters
959        ----------
960        file_obj : file-like object
961            ENDF file positioned at the second record of a resonance range
962            subsection in MF=2, MT=151
963        items : list
964            Items from the CONT record at the start of the resonance range
965            subsection
966        fission_widths : bool
967            Whether fission widths are given
968
969        Returns
970        -------
971        openmc.data.Unresolved
972            Unresolved resonance region parameters
973
974        """
975        # Read energy-dependent scattering radius if present
976        energy_min, energy_max = items[0:2]
977        nro, naps = items[4:6]
978        if nro != 0:
979            params, ape = get_tab1_record(file_obj)
980
981        # Get SPI, AP, and LSSF
982        formalism = items[3]
983        if not (fission_widths and formalism == 1):
984            items = get_cont_record(file_obj)
985            target_spin = items[0]
986            if nro == 0:
987                ap = Polynomial((items[1],))
988            add_to_background = (items[2] == 0)
989
990        if not fission_widths and formalism == 1:
991            # Case A -- fission widths not given, all parameters are
992            # energy-independent
993            NLS = items[4]
994            columns = ['L', 'J', 'd', 'amun', 'gn0', 'gg']
995            records = []
996            for ls in range(NLS):
997                items, values = get_list_record(file_obj)
998                awri = items[0]
999                l = items[2]
1000                NJS = items[5]
1001                for j in range(NJS):
1002                    d, j, amun, gn0, gg = values[6*j:6*j + 5]
1003                    records.append([l, j, d, amun, gn0, gg])
1004            parameters = pd.DataFrame.from_records(records, columns=columns)
1005            energies = None
1006
1007        elif fission_widths and formalism == 1:
1008            # Case B -- fission widths given, only fission widths are
1009            # energy-dependent
1010            items, energies = get_list_record(file_obj)
1011            target_spin = items[0]
1012            if nro == 0:
1013                ap = Polynomial((items[1],))
1014            add_to_background = (items[2] == 0)
1015            NE, NLS = items[4:6]
1016            records = []
1017            columns = ['L', 'J', 'E', 'd', 'amun', 'amuf', 'gn0', 'gg', 'gf']
1018            for ls in range(NLS):
1019                items = get_cont_record(file_obj)
1020                awri = items[0]
1021                l = items[2]
1022                NJS = items[4]
1023                for j in range(NJS):
1024                    items, values = get_list_record(file_obj)
1025                    muf = items[3]
1026                    d = values[0]
1027                    j = values[1]
1028                    amun = values[2]
1029                    gn0 = values[3]
1030                    gg = values[4]
1031                    gfs = values[6:]
1032                    for E, gf in zip(energies, gfs):
1033                        records.append([l, j, E, d, amun, muf, gn0, gg, gf])
1034            parameters = pd.DataFrame.from_records(records, columns=columns)
1035
1036        elif formalism == 2:
1037            # Case C -- all parameters are energy-dependent
1038            NLS = items[4]
1039            columns = ['L', 'J', 'E', 'd', 'amux', 'amun', 'amuf', 'gx', 'gn0',
1040                       'gg', 'gf']
1041            records = []
1042            for ls in range(NLS):
1043                items = get_cont_record(file_obj)
1044                awri = items[0]
1045                l = items[2]
1046                NJS = items[4]
1047                for j in range(NJS):
1048                    items, values = get_list_record(file_obj)
1049                    ne = items[5]
1050                    j = items[0]
1051                    amux = values[2]
1052                    amun = values[3]
1053                    amuf = values[5]
1054                    energies = []
1055                    for k in range(1, ne + 1):
1056                        E = values[6*k]
1057                        d = values[6*k + 1]
1058                        gx = values[6*k + 2]
1059                        gn0 = values[6*k + 3]
1060                        gg = values[6*k + 4]
1061                        gf = values[6*k + 5]
1062                        energies.append(E)
1063                        records.append([l, j, E, d, amux, amun, amuf, gx, gn0,
1064                                        gg, gf])
1065            parameters = pd.DataFrame.from_records(records, columns=columns)
1066
1067        # Calculate channel radius from ENDF-102 equation D.14
1068        a = Polynomial((0.123 * (NEUTRON_MASS*awri)**(1./3.) + 0.08,))
1069
1070        # Determine scattering and channel radius
1071        if nro == 0:
1072            scattering_radius = ap
1073            if naps == 0:
1074                channel_radius = a
1075            elif naps == 1:
1076                channel_radius = ap
1077        elif nro == 1:
1078            scattering_radius = ape
1079            if naps == 0:
1080                channel_radius = a
1081            elif naps == 1:
1082                channel_radius = ape
1083            elif naps == 2:
1084                channel_radius = ap
1085
1086        urr = cls(target_spin, energy_min, energy_max, channel_radius,
1087                  scattering_radius)
1088        urr.parameters = parameters
1089        urr.add_to_background = add_to_background
1090        urr.atomic_weight_ratio = awri
1091        urr.energies = energies
1092
1093        return urr
1094
1095
1096_FORMALISMS = {0: ResonanceRange,
1097               1: SingleLevelBreitWigner,
1098               2: MultiLevelBreitWigner,
1099               3: ReichMoore,
1100               7: RMatrixLimited}
1101
1102_RESOLVED = (SingleLevelBreitWigner, MultiLevelBreitWigner,
1103             ReichMoore, RMatrixLimited)
1104