from itertools import chain
from numbers import Integral, Real
import string

import numpy as np

import openmc.checkvalue as cv
import openmc.data

# Supported keywords for continuous-energy cross section plotting
PLOT_TYPES = ['total', 'scatter', 'elastic', 'inelastic', 'fission',
              'absorption', 'capture', 'nu-fission', 'nu-scatter', 'unity',
              'slowing-down power', 'damage']

# Supported keywords for multi-group cross section plotting
PLOT_TYPES_MGXS = ['total', 'absorption', 'scatter', 'fission',
                   'kappa-fission', 'nu-fission', 'prompt-nu-fission',
                   'deleyed-nu-fission', 'chi', 'chi-prompt', 'chi-delayed',
                   'inverse-velocity', 'beta', 'decay rate', 'unity']
# Create a dictionary which can be used to convert PLOT_TYPES_MGXS to the
# openmc.XSdata attribute name needed to access the data
_PLOT_MGXS_ATTR = {line: line.replace(' ', '_').replace('-', '_')
                   for line in PLOT_TYPES_MGXS}
_PLOT_MGXS_ATTR['scatter'] = 'scatter_matrix'

# Special MT values
UNITY_MT = -1
XI_MT = -2

# MTs to combine to generate associated plot_types
_INELASTIC = [mt for mt in openmc.data.SUM_RULES[3] if mt != 27]
PLOT_TYPES_MT = {
    'total': openmc.data.SUM_RULES[1],
    'scatter': [2] + _INELASTIC,
    'elastic': [2],
    'inelastic': _INELASTIC,
    'fission': [18],
    'absorption': [27],
    'capture': [101],
    'nu-fission': [18],
    'nu-scatter': [2] + _INELASTIC,
    'unity': [UNITY_MT],
    'slowing-down power': [2] + [XI_MT],
    'damage': [444]
}

# Types of plots to plot linearly in y
PLOT_TYPES_LINEAR = {'nu-fission / fission', 'nu-scatter / scatter',
                     'nu-fission / absorption', 'fission / absorption'}

# Minimum and maximum energies for plotting (units of eV)
_MIN_E = 1.e-5
_MAX_E = 20.e6


def plot_xs(this, types, divisor_types=None, temperature=294., data_type=None,
            axis=None, sab_name=None, ce_cross_sections=None,
            mg_cross_sections=None, enrichment=None, plot_CE=True, orders=None,
            divisor_orders=None, **kwargs):
    """Creates a figure of continuous-energy cross sections for this item.

    Parameters
    ----------
    this : str or openmc.Material
        Object to source data from
    types : Iterable of values of PLOT_TYPES
        The type of cross sections to include in the plot.
    divisor_types : Iterable of values of PLOT_TYPES, optional
        Cross section types which will divide those produced by types
        before plotting. A type of 'unity' can be used to effectively not
        divide some types.
    temperature : float, optional
        Temperature in Kelvin to plot. If not specified, a default
        temperature of 294K will be plotted. Note that the nearest
        temperature in the library for each nuclide will be used as opposed
        to using any interpolation.
    data_type : {'nuclide', 'element', 'material', 'macroscopic'}, optional
        Type of object to plot. If not specified, a guess is made based on the
        `this` argument.
    axis : matplotlib.axes, optional
        A previously generated axis to use for plotting. If not specified,
        a new axis and figure will be generated.
    sab_name : str, optional
        Name of S(a,b) library to apply to MT=2 data when applicable; only used
        for items which are instances of openmc.Element or openmc.Nuclide
    ce_cross_sections : str, optional
        Location of cross_sections.xml file. Default is None.
    mg_cross_sections : str, optional
        Location of MGXS HDF5 Library file. Default is None.
    enrichment : float, optional
        Enrichment for U235 in weight percent. For example, input 4.95 for
        4.95 weight percent enriched U. Default is None. This is only used for
        items which are instances of openmc.Element
    plot_CE : bool, optional
        Denotes whether or not continuous-energy will be plotted. Defaults to
        plotting the continuous-energy data.
    orders : Iterable of Integral, optional
        The scattering order or delayed group index to use for the
        corresponding entry in types. Defaults to the 0th order for scattering
        and the total delayed neutron data. This only applies to plots of
        multi-group data.
    divisor_orders : Iterable of Integral, optional
        Same as orders, but for divisor_types
    **kwargs
        All keyword arguments are passed to
        :func:`matplotlib.pyplot.figure`.

    Returns
    -------
    fig : matplotlib.figure.Figure
        If axis is None, then a Matplotlib Figure of the generated
        cross section will be returned. Otherwise, a value of
        None will be returned as the figure and axes have already been
        generated.

    """
    import matplotlib.pyplot as plt

    cv.check_type("plot_CE", plot_CE, bool)

    if data_type is None:
        if isinstance(this, openmc.Nuclide):
            data_type = 'nuclide'
        elif isinstance(this, openmc.Element):
            data_type = 'element'
        elif isinstance(this, openmc.Material):
            data_type = 'material'
        elif isinstance(this, openmc.Macroscopic):
            data_type = 'macroscopic'
        elif isinstance(this, str):
            if this[-1] in string.digits:
                data_type = 'nuclide'
            else:
                data_type = 'element'
        else:
            raise TypeError("Invalid type for plotting")

    if plot_CE:
        # Calculate for the CE cross sections
        E, data = calculate_cexs(this, data_type, types, temperature, sab_name,
                                 ce_cross_sections, enrichment)
        if divisor_types:
            cv.check_length('divisor types', divisor_types, len(types))
            Ediv, data_div = calculate_cexs(this, divisor_types, temperature,
                                            sab_name, ce_cross_sections,
                                            enrichment)

            # Create a new union grid, interpolate data and data_div on to that
            # grid, and then do the actual division
            Enum = E[:]
            E = np.union1d(Enum, Ediv)
            data_new = np.zeros((len(types), len(E)))

            for line in range(len(types)):
                data_new[line, :] = \
                    np.divide(np.interp(E, Enum, data[line, :]),
                              np.interp(E, Ediv, data_div[line, :]))
                if divisor_types[line] != 'unity':
                    types[line] = types[line] + ' / ' + divisor_types[line]
            data = data_new
    else:
        # Calculate for MG cross sections
        E, data = calculate_mgxs(this, data_type, types, orders, temperature,
                                 mg_cross_sections, ce_cross_sections,
                                 enrichment)

        if divisor_types:
            cv.check_length('divisor types', divisor_types, len(types))
            Ediv, data_div = calculate_mgxs(this, data_type, divisor_types,
                                            divisor_orders, temperature,
                                            mg_cross_sections,
                                            ce_cross_sections, enrichment)

            # Perform the division
            for line in range(len(types)):
                data[line, :] /= data_div[line, :]
                if divisor_types[line] != 'unity':
                    types[line] += ' / ' + divisor_types[line]

    # Generate the plot
    if axis is None:
        fig, ax = plt.subplots()
    else:
        fig = None
        ax = axis
    # Set to loglog or semilogx depending on if we are plotting a data
    # type which we expect to vary linearly
    if set(types).issubset(PLOT_TYPES_LINEAR):
        plot_func = ax.semilogx
    else:
        plot_func = ax.loglog

    # Plot the data
    for i in range(len(data)):
        data[i, :] = np.nan_to_num(data[i, :])
        if np.sum(data[i, :]) > 0.:
            plot_func(E, data[i, :], label=types[i])

    ax.set_xlabel('Energy [eV]')
    if plot_CE:
        ax.set_xlim(_MIN_E, _MAX_E)
    else:
        ax.set_xlim(E[-1], E[0])
    if divisor_types:
        if data_type == 'nuclide':
            ylabel = 'Nuclidic Microscopic Data'
        elif data_type == 'element':
            ylabel = 'Elemental Microscopic Data'
        elif data_type == 'material' or data_type == 'macroscopic':
            ylabel = 'Macroscopic Data'
    else:
        if data_type == 'nuclide':
            ylabel = 'Microscopic Cross Section [b]'
        elif data_type == 'element':
            ylabel = 'Elemental Cross Section [b]'
        elif data_type == 'material' or data_type == 'macroscopic':
            ylabel = 'Macroscopic Cross Section [1/cm]'
    ax.set_ylabel(ylabel)
    ax.legend(loc='best')
    name = this.name if data_type == 'material' else this
    if len(types) > 1:
        ax.set_title('Cross Sections for ' + name)
    else:
        ax.set_title('Cross Section for ' + name)

    return fig


def calculate_cexs(this, data_type, types, temperature=294., sab_name=None,
                   cross_sections=None, enrichment=None):
    """Calculates continuous-energy cross sections of a requested type.

    Parameters
    ----------
    this : {str, openmc.Nuclide, openmc.Element, openmc.Material}
        Object to source data from
    data_type : {'nuclide', 'element', 'material'}
        Type of object to plot
    types : Iterable of values of PLOT_TYPES
        The type of cross sections to calculate
    temperature : float, optional
        Temperature in Kelvin to plot. If not specified, a default
        temperature of 294K will be plotted. Note that the nearest
        temperature in the library for each nuclide will be used as opposed
        to using any interpolation.
    sab_name : str, optional
        Name of S(a,b) library to apply to MT=2 data when applicable.
    cross_sections : str, optional
        Location of cross_sections.xml file. Default is None.
    enrichment : float, optional
        Enrichment for U235 in weight percent. For example, input 4.95 for
        4.95 weight percent enriched U. Default is None
        (natural composition).

    Returns
    -------
    energy_grid : numpy.ndarray
        Energies at which cross sections are calculated, in units of eV
    data : numpy.ndarray
        Cross sections calculated at the energy grid described by energy_grid

    """

    # Check types
    cv.check_type('temperature', temperature, Real)
    if sab_name:
        cv.check_type('sab_name', sab_name, str)
    if enrichment:
        cv.check_type('enrichment', enrichment, Real)

    if data_type == 'nuclide':
        if isinstance(this, str):
            nuc = openmc.Nuclide(this)
        else:
            nuc = this
        energy_grid, xs = _calculate_cexs_nuclide(nuc, types, temperature,
                                                  sab_name, cross_sections)
        # Convert xs (Iterable of Callable) to a grid of cross section values
        # calculated on the points in energy_grid for consistency with the
        # element and material functions.
        data = np.zeros((len(types), len(energy_grid)))
        for line in range(len(types)):
            data[line, :] = xs[line](energy_grid)
    elif data_type == 'element':
        if isinstance(this, str):
            elem = openmc.Element(this)
        else:
            elem = this
        energy_grid, data = _calculate_cexs_elem_mat(elem, types, temperature,
                                                     cross_sections, sab_name,
                                                     enrichment)
    elif data_type == 'material':
        cv.check_type('this', this, openmc.Material)
        energy_grid, data = _calculate_cexs_elem_mat(this, types, temperature,
                                                     cross_sections)
    else:
        raise TypeError("Invalid type")

    return energy_grid, data


def _calculate_cexs_nuclide(this, types, temperature=294., sab_name=None,
                            cross_sections=None):
    """Calculates continuous-energy cross sections of a requested type.

    Parameters
    ----------
    this : openmc.Nuclide
        Nuclide object to source data from
    types : Iterable of str or Integral
        The type of cross sections to calculate; values can either be those
        in openmc.PLOT_TYPES or integers which correspond to reaction
        channel (MT) numbers.
    temperature : float, optional
        Temperature in Kelvin to plot. If not specified, a default
        temperature of 294K will be plotted. Note that the nearest
        temperature in the library for each nuclide will be used as opposed
        to using any interpolation.
    sab_name : str, optional
        Name of S(a,b) library to apply to MT=2 data when applicable.
    cross_sections : str, optional
        Location of cross_sections.xml file. Default is None.

    Returns
    -------
    energy_grid : numpy.ndarray
        Energies at which cross sections are calculated, in units of eV
    data : Iterable of Callable
        Requested cross section functions

    """

    # Load the library
    library = openmc.data.DataLibrary.from_xml(cross_sections)

    # Convert temperature to format needed for access in the library
    strT = "{}K".format(int(round(temperature)))
    T = temperature

    # Now we can create the data sets to be plotted
    energy_grid = []
    xs = []
    lib = library.get_by_material(this)
    if lib is not None:
        nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path'])
        # Obtain the nearest temperature
        if strT in nuc.temperatures:
            nucT = strT
        else:
            delta_T = np.array(nuc.kTs) - T * openmc.data.K_BOLTZMANN
            closest_index = np.argmin(np.abs(delta_T))
            nucT = nuc.temperatures[closest_index]

        # Prep S(a,b) data if needed
        if sab_name:
            sab = openmc.data.ThermalScattering.from_hdf5(sab_name)
            # Obtain the nearest temperature
            if strT in sab.temperatures:
                sabT = strT
            else:
                delta_T = np.array(sab.kTs) - T * openmc.data.K_BOLTZMANN
                closest_index = np.argmin(np.abs(delta_T))
                sabT = sab.temperatures[closest_index]

            # Create an energy grid composed the S(a,b) and the nuclide's grid
            grid = nuc.energy[nucT]
            sab_Emax = 0.
            sab_funcs = []
            if sab.elastic is not None:
                elastic = sab.elastic.xs[sabT]
                if isinstance(elastic, openmc.data.CoherentElastic):
                    grid = np.union1d(grid, elastic.bragg_edges)
                    if elastic.bragg_edges[-1] > sab_Emax:
                        sab_Emax = elastic.bragg_edges[-1]
                elif isinstance(elastic, openmc.data.Tabulated1D):
                    grid = np.union1d(grid, elastic.x)
                    if elastic.x[-1] > sab_Emax:
                        sab_Emax = elastic.x[-1]
                sab_funcs.append(elastic)
            if sab.inelastic is not None:
                inelastic = sab.inelastic.xs[sabT]
                grid = np.union1d(grid, inelastic.x)
                if inelastic.x[-1] > sab_Emax:
                        sab_Emax = inelastic.x[-1]
                sab_funcs.append(inelastic)
            energy_grid = grid
        else:
            energy_grid = nuc.energy[nucT]

        # Parse the types
        mts = []
        ops = []
        yields = []
        for line in types:
            if line in PLOT_TYPES:
                tmp_mts = [mtj for mti in PLOT_TYPES_MT[line] for mtj in
                           nuc.get_reaction_components(mti)]
                mts.append(tmp_mts)
                if line.startswith('nu'):
                    yields.append(True)
                else:
                    yields.append(False)
                if XI_MT in tmp_mts:
                    ops.append((np.add,) * (len(tmp_mts) - 2) + (np.multiply,))
                else:
                    ops.append((np.add,) * (len(tmp_mts) - 1))
            else:
                # Not a built-in type, we have to parse it ourselves
                cv.check_type('MT in types', line, Integral)
                cv.check_greater_than('MT in types', line, 0)
                tmp_mts = nuc.get_reaction_components(line)
                mts.append(tmp_mts)
                ops.append((np.add,) * (len(tmp_mts) - 1))
                yields.append(False)

        for i, mt_set in enumerate(mts):
            # Get the reaction xs data from the nuclide
            funcs = []
            op = ops[i]
            for mt in mt_set:
                if mt == 2:
                    if sab_name:
                        # Then we need to do a piece-wise function of
                        # The S(a,b) and non-thermal data
                        sab_sum = openmc.data.Sum(sab_funcs)
                        pw_funcs = openmc.data.Regions1D(
                            [sab_sum, nuc[mt].xs[nucT]],
                            [sab_Emax])
                        funcs.append(pw_funcs)
                    else:
                        funcs.append(nuc[mt].xs[nucT])
                elif mt in nuc:
                    if yields[i]:
                        # Get the total yield first if available. This will be
                        # used primarily for fission.
                        for prod in chain(nuc[mt].products,
                                          nuc[mt].derived_products):
                            if prod.particle == 'neutron' and \
                                prod.emission_mode == 'total':
                                func = openmc.data.Combination(
                                    [nuc[mt].xs[nucT], prod.yield_],
                                    [np.multiply])
                                funcs.append(func)
                                break
                        else:
                            # Total doesn't exist so we have to create from
                            # prompt and delayed. This is used for scatter
                            # multiplication.
                            func = None
                            for prod in chain(nuc[mt].products,
                                              nuc[mt].derived_products):
                                if prod.particle == 'neutron' and \
                                    prod.emission_mode != 'total':
                                    if func:
                                        func = openmc.data.Combination(
                                            [prod.yield_, func], [np.add])
                                    else:
                                        func = prod.yield_
                            if func:
                                funcs.append(openmc.data.Combination(
                                    [func, nuc[mt].xs[nucT]], [np.multiply]))
                            else:
                                # If func is still None, then there were no
                                # products. In that case, assume the yield is
                                # one as its not provided for some summed
                                # reactions like MT=4
                                funcs.append(nuc[mt].xs[nucT])
                    else:
                        funcs.append(nuc[mt].xs[nucT])
                elif mt == UNITY_MT:
                    funcs.append(lambda x: 1.)
                elif mt == XI_MT:
                    awr = nuc.atomic_weight_ratio
                    alpha = ((awr - 1.) / (awr + 1.))**2
                    xi = 1. + alpha * np.log(alpha) / (1. - alpha)
                    funcs.append(lambda x: xi)
                else:
                    funcs.append(lambda x: 0.)
            funcs = funcs if funcs else [lambda x: 0.]
            xs.append(openmc.data.Combination(funcs, op))
    else:
        raise ValueError(this + " not in library")

    return energy_grid, xs


def _calculate_cexs_elem_mat(this, types, temperature=294.,
                             cross_sections=None, sab_name=None,
                             enrichment=None):
    """Calculates continuous-energy cross sections of a requested type.

    Parameters
    ----------
    this : openmc.Material or openmc.Element
        Object to source data from
    types : Iterable of values of PLOT_TYPES
        The type of cross sections to calculate
    temperature : float, optional
        Temperature in Kelvin to plot. If not specified, a default
        temperature of 294K will be plotted. Note that the nearest
        temperature in the library for each nuclide will be used as opposed
        to using any interpolation.
    cross_sections : str, optional
        Location of cross_sections.xml file. Default is None.
    sab_name : str, optional
        Name of S(a,b) library to apply to MT=2 data when applicable.
    enrichment : float, optional
        Enrichment for U235 in weight percent. For example, input 4.95 for
        4.95 weight percent enriched U. Default is None
        (natural composition).

    Returns
    -------
    energy_grid : numpy.ndarray
        Energies at which cross sections are calculated, in units of eV
    data : numpy.ndarray
        Cross sections calculated at the energy grid described by energy_grid

    """

    if isinstance(this, openmc.Material):
        if this.temperature is not None:
            T = this.temperature
        else:
            T = temperature
    else:
        T = temperature

    # Load the library
    library = openmc.data.DataLibrary.from_xml(cross_sections)

    if isinstance(this, openmc.Material):
        # Expand elements in to nuclides with atomic densities
        nuclides = this.get_nuclide_atom_densities()
        # For ease of processing split out the nuclide and its fraction
        nuc_fractions = {nuclide[1][0]: nuclide[1][1]
                         for nuclide in nuclides.items()}
        # Create a dict of [nuclide name] = nuclide object to carry forward
        # with a common nuclides format between openmc.Material and
        # openmc.Element objects
        nuclides = {nuclide[1][0]: nuclide[1][0]
                    for nuclide in nuclides.items()}
    else:
        # Expand elements in to nuclides with atomic densities
        nuclides = this.expand(1., 'ao', enrichment=enrichment,
                               cross_sections=cross_sections)
        # For ease of processing split out the nuclide and its fraction
        nuc_fractions = {nuclide[0]: nuclide[1] for nuclide in nuclides}
        # Create a dict of [nuclide name] = nuclide object to carry forward
        # with a common nuclides format between openmc.Material and
        # openmc.Element objects
        nuclides = {nuclide[0]: nuclide[0] for nuclide in nuclides}

    # Identify the nuclides which have S(a,b) data
    sabs = {}
    for nuclide in nuclides.items():
        sabs[nuclide[0]] = None
    if isinstance(this, openmc.Material):
        for sab_name in this._sab:
            sab = openmc.data.ThermalScattering.from_hdf5(
                library.get_by_material(sab_name, data_type='thermal')['path'])
            for nuc in sab.nuclides:
                sabs[nuc] = library.get_by_material(sab_name,
                        data_type='thermal')['path']
    else:
        if sab_name:
            sab = openmc.data.ThermalScattering.from_hdf5(sab_name)
            for nuc in sab.nuclides:
                sabs[nuc] = library.get_by_material(sab_name,
                        data_type='thermal')['path']

    # Now we can create the data sets to be plotted
    xs = {}
    E = []
    for nuclide in nuclides.items():
        name = nuclide[0]
        nuc = nuclide[1]
        sab_tab = sabs[name]
        temp_E, temp_xs = calculate_cexs(nuc, 'nuclide', types, T, sab_tab,
                                         cross_sections)
        E.append(temp_E)
        # Since the energy grids are different, store the cross sections as
        # a tabulated function so they can be calculated on any grid needed.
        xs[name] = [openmc.data.Tabulated1D(temp_E, temp_xs[line])
                    for line in range(len(types))]

    # Condense the data for every nuclide
    # First create a union energy grid
    energy_grid = E[0]
    for grid in E[1:]:
        energy_grid = np.union1d(energy_grid, grid)

    # Now we can combine all the nuclidic data
    data = np.zeros((len(types), len(energy_grid)))
    for line in range(len(types)):
        if types[line] == 'unity':
            data[line, :] = 1.
        else:
            for nuclide in nuclides.items():
                name = nuclide[0]
                data[line, :] += (nuc_fractions[name] *
                                  xs[name][line](energy_grid))

    return energy_grid, data


def calculate_mgxs(this, data_type, types, orders=None, temperature=294.,
                   cross_sections=None, ce_cross_sections=None,
                   enrichment=None):
    """Calculates multi-group cross sections of a requested type.

    If the data for the nuclide or macroscopic object in the library is
    represented as angle-dependent data then this method will return the
    geometric average cross section over all angles.

    Parameters
    ----------
    this : str or openmc.Material
        Object to source data from
    data_type : {'nuclide', 'element', 'material', 'macroscopic'}
        Type of object to plot
    types : Iterable of values of PLOT_TYPES_MGXS
        The type of cross sections to calculate
    orders : Iterable of Integral, optional
        The scattering order or delayed group index to use for the
        corresponding entry in types. Defaults to the 0th order for scattering
        and the total delayed neutron data.
    temperature : float, optional
        Temperature in Kelvin to plot. If not specified, a default
        temperature of 294K will be plotted. Note that the nearest
        temperature in the library for each nuclide will be used as opposed
        to using any interpolation.
    cross_sections : str, optional
        Location of MGXS HDF5 Library file. Default is None.
    ce_cross_sections : str, optional
        Location of continuous-energy cross_sections.xml file. Default is None.
        This is used only for expanding an openmc.Element object passed as this
    enrichment : float, optional
        Enrichment for U235 in weight percent. For example, input 4.95 for
        4.95 weight percent enriched U. Default is None
        (natural composition).

    Returns
    -------
    energy_grid : numpy.ndarray
        Energies at which cross sections are calculated, in units of eV
    data : numpy.ndarray
        Cross sections calculated at the energy grid described by energy_grid

    """

    # Check types
    cv.check_type('temperature', temperature, Real)
    if enrichment:
        cv.check_type('enrichment', enrichment, Real)
    cv.check_iterable_type('types', types, str)

    cv.check_type("cross_sections", cross_sections, str)
    library = openmc.MGXSLibrary.from_hdf5(cross_sections)

    if data_type in ('nuclide', 'macroscopic'):
        mgxs = _calculate_mgxs_nuc_macro(this, types, library, orders,
                                         temperature)
    elif data_type in ('element', 'material'):
        mgxs = _calculate_mgxs_elem_mat(this, types, library, orders,
                                        temperature, ce_cross_sections,
                                        enrichment)
    else:
        raise TypeError("Invalid type")

    # Convert the data to the format needed
    data = np.zeros((len(types), 2 * library.energy_groups.num_groups))
    energy_grid = np.zeros(2 * library.energy_groups.num_groups)
    for g in range(library.energy_groups.num_groups):
        energy_grid[g * 2: g * 2 + 2] = \
            library.energy_groups.group_edges[g: g + 2]
    # Ensure the energy will show on a log-axis by replacing 0s with a
    # sufficiently small number
    energy_grid[0] = max(energy_grid[0], _MIN_E)

    for line in range(len(types)):
        for g in range(library.energy_groups.num_groups):
            data[line, g * 2: g * 2 + 2] = mgxs[line, g]

    return energy_grid[::-1], data


def _calculate_mgxs_nuc_macro(this, types, library, orders=None,
                              temperature=294.):
    """Determines the multi-group cross sections of a nuclide or macroscopic
    object.

    If the data for the nuclide or macroscopic object in the library is
    represented as angle-dependent data then this method will return the
    geometric average cross section over all angles.

    Parameters
    ----------
    this : openmc.Nuclide or openmc.Macroscopic
        Object to source data from
    types : Iterable of str
        The type of cross sections to calculate; values can either be those
        in openmc.PLOT_TYPES_MGXS
    library : openmc.MGXSLibrary
        MGXS Library containing the data of interest
    orders : Iterable of Integral, optional
        The scattering order or delayed group index to use for the
        corresponding entry in types. Defaults to the 0th order for scattering
        and the total delayed neutron data.
    temperature : float, optional
        Temperature in Kelvin to plot. If not specified, a default
        temperature of 294K will be plotted. Note that the nearest
        temperature in the library for each nuclide will be used as opposed
        to using any interpolation.

    Returns
    -------
    data : numpy.ndarray
        Cross sections calculated at the energy grid described by energy_grid

    """

    # Check the parameters and grab order/delayed groups
    if orders:
        cv.check_iterable_type('orders', orders, Integral,
                               min_depth=len(types), max_depth=len(types))
    else:
        orders = [None] * len(types)
    for i, line in enumerate(types):
        cv.check_type("line", line, str)
        cv.check_value("line", line, PLOT_TYPES_MGXS)
        if orders[i]:
            cv.check_greater_than("order value", orders[i], 0, equality=True)

    xsdata = library.get_by_name(this)

    if xsdata is not None:
        # Obtain the nearest temperature
        t = np.abs(xsdata.temperatures - temperature).argmin()

        # Get the data
        data = np.zeros((len(types), library.energy_groups.num_groups))
        for i, line in enumerate(types):
            if 'fission' in line and not xsdata.fissionable:
                continue
            elif line == 'unity':
                data[i, :] = 1.
            else:
                # Now we have to get the cross section data and properly
                # treat it depending on the requested type.
                # First get the data in a generic fashion
                temp_data = getattr(xsdata, _PLOT_MGXS_ATTR[line])[t]
                shape = temp_data.shape[:]
                # If we have angular data, then want the geometric
                # average over all provided angles.  Since the angles are
                # equi-distant, un-weighted averaging will suffice
                if xsdata.representation == 'angle':
                    temp_data = np.mean(temp_data, axis=(0, 1))

                # Now we can look at the shape of the data to identify how
                # it should be modified to produce an array of values
                # with groups.
                if shape in (xsdata.xs_shapes["[G']"],
                             xsdata.xs_shapes["[G]"]):
                    # Then the data is already an array vs groups so copy
                    # and move along
                    data[i, :] = temp_data
                elif shape == xsdata.xs_shapes["[G][G']"]:
                    # Sum the data over outgoing groups to create our array vs
                    # groups
                    data[i, :] = np.sum(temp_data, axis=1)
                elif shape == xsdata.xs_shapes["[DG]"]:
                    # Then we have a constant vs groups with a value for each
                    # delayed group. The user-provided value of orders tells us
                    # which delayed group we want. If none are provided, then
                    # we sum all the delayed groups together.
                    if orders[i]:
                        if orders[i] < len(shape[0]):
                            data[i, :] = temp_data[orders[i]]
                    else:
                        data[i, :] = np.sum(temp_data[:])
                elif shape in (xsdata.xs_shapes["[DG][G']"],
                               xsdata.xs_shapes["[DG][G]"]):
                    # Then we have an array vs groups with values for each
                    # delayed group. The user-provided value of orders tells us
                    # which delayed group we want. If none are provided, then
                    # we sum all the delayed groups together.
                    if orders[i]:
                        if orders[i] < len(shape[0]):
                            data[i, :] = temp_data[orders[i], :]
                    else:
                        data[i, :] = np.sum(temp_data[:, :], axis=0)
                elif shape == xsdata.xs_shapes["[DG][G][G']"]:
                    # Then we have a delayed group matrix. We will first
                    # remove the outgoing group dependency
                    temp_data = np.sum(temp_data, axis=-1)
                    # And then proceed in exactly the same manner as the
                    # "[DG][G']" or "[DG][G]" shapes in the previous block.
                    if orders[i]:
                        if orders[i] < len(shape[0]):
                            data[i, :] = temp_data[orders[i], :]
                    else:
                        data[i, :] = np.sum(temp_data[:, :], axis=0)
                elif shape == xsdata.xs_shapes["[G][G'][Order]"]:
                    # This is a scattering matrix with angular data
                    # First remove the outgoing group dependence
                    temp_data = np.sum(temp_data, axis=1)
                    # The user either provided a specific order or we resort
                    # to the default 0th order
                    if orders[i]:
                        order = orders[i]
                    else:
                        order = 0
                    # If the order is available, store the data for that order
                    # if it is not available, then the expansion coefficient
                    # is zero and thus we already have the correct value.
                    if order < shape[1]:
                        data[i, :] = temp_data[:, order]
    else:
        raise ValueError("{} not present in provided MGXS "
                         "library".format(this))

    return data


def _calculate_mgxs_elem_mat(this, types, library, orders=None,
                             temperature=294., ce_cross_sections=None,
                             enrichment=None):
    """Determines the multi-group cross sections of an element or material
    object.

    If the data for the nuclide or macroscopic object in the library is
    represented as angle-dependent data then this method will return the
    geometric average cross section over all angles.

    Parameters
    ----------
    this : openmc.Element or openmc.Material
        Object to source data from
    types : Iterable of str
        The type of cross sections to calculate; values can either be those
        in openmc.PLOT_TYPES_MGXS
    library : openmc.MGXSLibrary
        MGXS Library containing the data of interest
    orders : Iterable of Integral, optional
        The scattering order or delayed group index to use for the
        corresponding entry in types. Defaults to the 0th order for scattering
        and the total delayed neutron data.
    temperature : float, optional
        Temperature in Kelvin to plot. If not specified, a default
        temperature of 294K will be plotted. Note that the nearest
        temperature in the library for each nuclide will be used as opposed
        to using any interpolation.
    ce_cross_sections : str, optional
        Location of continuous-energy cross_sections.xml file. Default is None.
        This is used only for expanding the elements
    enrichment : float, optional
        Enrichment for U235 in weight percent. For example, input 4.95 for
        4.95 weight percent enriched U. Default is None
        (natural composition).

    Returns
    -------
    data : numpy.ndarray
        Cross sections calculated at the energy grid described by energy_grid

    """

    if isinstance(this, openmc.Material):
        if this.temperature is not None:
            T = this.temperature
        else:
            T = temperature

        # Check to see if we have nuclides/elements or a macrocopic object
        if this._macroscopic is not None:
            # We have macroscopics
            nuclides = {this._macroscopic: (this._macroscopic, this.density)}
        else:
            # Expand elements in to nuclides with atomic densities
            nuclides = this.get_nuclide_atom_densities()

        # For ease of processing split out nuc and nuc_density
        nuc_fraction = [nuclide[1][1] for nuclide in nuclides.items()]
    else:
        T = temperature
        # Expand elements in to nuclides with atomic densities
        nuclides = this.expand(100., 'ao', enrichment=enrichment,
                               cross_sections=ce_cross_sections)

        # For ease of processing split out nuc and nuc_fractions
        nuc_fraction = [nuclide[1] for nuclide in nuclides]

    nuc_data = []
    for nuclide in nuclides.items():
        nuc_data.append(_calculate_mgxs_nuc_macro(nuclide[0], types, library,
                                                  orders, T))

    # Combine across the nuclides
    data = np.zeros((len(types), library.energy_groups.num_groups))
    for line in range(len(types)):
        if types[line] == 'unity':
            data[line, :] = 1.
        else:
            for n in range(len(nuclides)):
                data[line, :] += nuc_fraction[n] * nuc_data[n][line, :]

    return data
