def write_metadata(self, nucs, libs, dirname):
     track_actinides = [n for n in nucs if nucname.znum(n) in nucname.act]
     with open(os.path.join(dirname, "manifest.txt"), "w") as f:
         f.write("\n".join([str(nucname.zzaaam(act)) for act in track_actinides]))
         f.write("\n")
     with open(os.path.join(dirname, "params.txt"), "w") as f:
         if self.rc.get("enrichment") is None:
             enrichment = self.rc.initial_heavy_metal.get(922350)
         else:
             enrichment = self.rc.enrichment
         if enrichment is not None:
             f.write("ENRICHMENT {}\n".format(enrichment))
         if self.rc.get("batches") is not None:
             f.write("BATCHES {}\n".format(self.rc.batches))
         if self.rc.get("pnl") is not None:
             f.write("PNL {}\n".format(self.rc.pnl))
         f.write("BURNUP {}\n".format(sum(libs["fuel"]["BUd"])))
         f.write("FLUX {:.0E}\n".format(np.mean(libs["fuel"]["phi_tot"][1:])))
     with open(os.path.join(dirname, "structural.txt"), "w") as f:
         clad_linear_density = pi * self.rc.clad_density * \
             (self.rc.clad_cell_radius ** 2 - self.rc.void_cell_radius ** 2)
         fuel_linear_density = pi * self.rc.fuel_density * \
             self.rc.fuel_cell_radius ** 2
         clad_frac = float(clad_linear_density / fuel_linear_density)
         cladrows = ["{} {:.8f}".format(nucname.zzaaam(n), f*clad_frac)
                     for n, f in self.rc.clad_material.comp.items()]
         f.write("\n".join(cladrows))
         f.write("\n")
     shutil.copyfile("TAPE9.INP", os.path.join(dirname, "TAPE9.INP"))

def parse_nucs(s):
    """Parses a string into a set of nuclides."""
    nset = set()
    nucs = s.split(',')

    for nuc in nucs:
        if len(nuc) == 0:
            continue
        elif '-' in nuc:
            nsplit = nuc.split()
            nlower = nucname.zzaaam(nsplit[0])
            nupper = nucname.zzaaam(nsplit[1])
            if 0 == nupper%10000:
                nupper += 10000
            else:
                nupper += 1
            tmpset = set(range(nlower, nupper))
        else:
            n = nucname.zzaaam(nuc)
            if 0 == n%10000:
                nrange = range(n, n + 10000)
            else:
                nrange = [n]
            tmpset = set(nrange)

        # Add the union 
        nset = (nset | tmpset)

    return nset

def parse_scattering_lengths(build_dir):
    """Converts to scattering lenth data to a numpy array."""
    build_filename = os.path.join(build_dir, "scattering_lengths.html")
    
    # Read in cinder data file
    with open(build_filename, 'r') as f:
        raw_data = f.read()

    sl_data = []

    # Iterate over all isotopes in the table
    for m in re.finditer(scat_len_pattern, raw_data):
        md = m.groupdict()

        slrow = (nucname.zzaaam(md['iso']),
                 nist_num(md['b_coherent']) * (1E-13),
                 nist_num(md['b_incoherent']) * (1E-13),
                 nist_num(md['xs_coherent']),
                 nist_num(md['xs_incoherent']),
                 nist_num(md['xs']))

        sl_data.append(slrow)

    sl_array = np.array(sl_data, dtype=sl_dtype)
    return sl_array

    def update(self, libs, dirname, fname):
        rownames = ["TIME", "phi_tot", "NEUT_PROD", "NEUT_DEST", "BUd"]
        trans_matrix = {}
        for mat, matlib in libs.items():
            if isinstance(mat, int):
                fname = str(nucname.zzaaam(mat))
            elif mat == 'fuel':
                fname = mat
            else:
                continue
            f = open(os.path.join(dirname, fname + ".txt"), "r")         
            lines = f.readlines()
            time_steps = len(lines[0].split())
            i = 5
            while i < range(len(lines)):
                nuc_array = lines[i].split()
                nuc_name = nuc_array[0]
                nuc_values = nuc_array[1:]
                if len(nuc_array) == time_steps:
                    trans_matrix[nuc_name] = nuc_values
                    trans_matrix[nuc_name].append(matlib['material'][-1].comp[temp_nuc]*1000)
                else:
                    if matlib['material'][-1].comp[nuc_name] > self.rc.track_nuc_threshold:
                        zero_array = [0.]*(time_steps-1)
                        trans_matrix[nuc_name] = zero_array
                        trans_matrix[nuc_name].append(matlib['material'][-1].comp[temp_nuc]*1000)
                i+=1;
            if time_steps == len(matlib['TIME']):
                return libs, trans_matrix
            for j in range(len(rownames)):
                store = matlib[rownames[j]][-1]
                matlib[rownames[j]] = lines[j].split()[1:]
                matlib[rownames[j]].append(store)

        return libs, trans_matrix    

def make_1g_xs_graphs(nuc, sig):
    global burn_times
    nuc_zz = nucname.zzaaam(nuc)
    plt.clf()

    reactions = ['sigma_t', 'sigma_s', 'sigma_a', 'sigma_f']
    markers = ['-', 'o', 's', 'x']
    ncol = 2

    if nuc_zz < 860000:
        reactions = reactions[:-1]
        markers = markers[:-1]
        ncol = 3

    for reaction, marker in zip(reactions, markers):
        r, s, diff = sig[reaction, nuc]
        plt.plot(burn_times, s, color='k', marker=marker, linestyle='-', label="Serpent $\\{0}$".format(reaction))
        #plt.errorbar(burn_times, r, diff*r, color='r', marker=marker, linestyle='-', label="RMG $\\{0}$".format(reaction))
        plt.plot(burn_times, r, color='r', marker=marker, linestyle='-', label="RMG $\\{0}$".format(reaction))

    plt.xlabel("burn time [days]")
    plt.ylabel(nuc + " One-Group Cross Sections [barns]")
    plt.ticklabel_format(scilimits=(-5, 5))
    plt.legend(loc=0, ncol=ncol)
    plt.savefig(nuc + '_1g_xs.png')
    plt.savefig(nuc + '_1g_xs.eps')
    plt.clf()

def sigma_f(nuc, temp=300.0, group_struct=None, phi_g=None, xs_cache=None):
    """Calculates the neutron fission cross section for a nuclide. 

    Parameters
    ----------
    nuc : int, str, Material, or dict-like 
        A nuclide or nuclide-atom fraction mapping.
    temp : float, optional
        Temperature [K] of material, defaults to 300.0.
    group_struct : array-like of floats, optional
        Energy group structure E_g [MeV] from highest-to-lowest energy, length G+1,
        defaults to xs_cache['E_g'].
    phi_g : array-like of floats, optional
        Group fluxes [n/cm^2/s] matching group_struct, length G, defaults to 
        xs_cache['phi_g'].
    xs_cache : XSCache, optional
        Cross section cache to use, defaults to pyne.xs.cache.xs_cache. 

    Returns
    -------
    sigma_f_g : ndarray 
        The fission cross-section, length G.

    Notes
    -----
    This always pulls the fission cross section out of the cache.

    """
    xs_cache = cache.xs_cache if xs_cache is None else xs_cache
    _prep_cache(xs_cache, group_struct, phi_g)
    if isinstance(nuc, collections.Iterable) and not isinstance(nuc, basestring):
        return _atom_weight_channel(sigma_f, nuc, temp=temp, xs_cache=xs_cache)
    nuc = nucname.zzaaam(nuc)
    key = (nuc, "f", temp)
    return xs_cache[key]

def test_FuelFabrication_7():
    # Reactor to use
    rp = ReactorParameters()
    r1g = Reactor1G(rp=rp, n="r1g")

    # Mass streams to use
    u235 = Material({922350: 1.0}, 1.0, name="U-235")
    u238 = Material({922380: 1.0}, 1.0, name="U-238")
    mats = {"U235": u235, "U238": u238}

    # Mass weights to use
    mws = {"U235": -1.0, "U238": -1.0}

    # Fuel Fabrication Facility
    ff = FuelFabrication(mats=mats, mws_in=mws, r=r1g, paramtrack=set(["Mass"]), n="ff")

    keys = ["U235", "U238"]
    assert_equal(set(ff.materials.keys()), set(keys))

    for iso in keys:
        assert_equal(ff.materials[iso].mass, 1.0)
        assert_equal(ff.materials[iso].comp[nucname.zzaaam(iso)], 1.0)

    assert_equal(ff.mass_weights_in, mws)

    assert_equal(ff.track_params, set(["Mass", "Weight_U235", "deltaR_U235", "Weight_U238", "deltaR_U238"]))

    assert_equal(ff.name, "ff")

    assert_equal(ff.reactor.name, "r1g")
    r1g.name = "r1g name"
    ff.initialize(mats, mws, r1g)
    assert_equal(ff.reactor.name, "r1g name")

def metastable_ratio(nuc, rx, temp=300.0, group_struct=None, phi_g=None, xs_cache=None):
    """Calculates the ratio between a reaction that leaves the nuclide in a 
    metastable state and the equivalent reaction that leaves the nuclide in 
    the ground state.  This allows the calculation of metastable cross sections 
    via sigma_ms = ratio * sigma_ground. 

    Parameters
    ----------
    nuc : int, str, Material, or dict-like 
        A nuclide or nuclide-atom fraction mapping.
    rx : str
        Reaction key. ('gamma', 'alpha', 'p', etc.)
    temp : float, optional
        Temperature [K] of material, defaults to 300.0.
    group_struct : array-like of floats, optional
        Energy group structure E_g [MeV] from highest-to-lowest energy, length G+1,
        defaults to xs_cache['E_g'].
    phi_g : array-like of floats, optional
        Group fluxes [n/cm^2/s] matching group_struct, length G, defaults to 
        xs_cache['phi_g'].
    xs_cache : XSCache, optional
        Cross section cache to use, defaults to pyne.xs.cache.xs_cache. 

    Returns
    -------
    ratio_rx_g : ndarray
        An array of the ratio of the metastable cross section for a reaction 
        to the ground state reaction.

    Notes
    -----
    This always pulls the absorption reaction cross section out of the cache.

    See Also
    --------
    pyne.xs.data_source.RX_TYPES
    pyne.xs.data_source.RX_TYPES_MAP 

    """
    if isinstance(nuc, int) or isinstance(nuc, basestring):
        xs_cache = cache.xs_cache if xs_cache is None else xs_cache
        _prep_cache(xs_cache, group_struct, phi_g)
        nuc = nucname.zzaaam(nuc)
        key = (nuc, rx + "_x_ratio", temp)
        if key in xs_cache:
            return xs_cache[key]

    # Get the cross-sections
    sigma_rx = sigma_a_reaction(nuc, rx, temp, group_struct, phi_g, xs_cache)
    sigma_rx_x = sigma_a_reaction(nuc, rx + "_x", temp, group_struct, phi_g, xs_cache)

    # Get the ratio
    ratio_rx_g = sigma_rx_x / sigma_rx
    ratio_rx_g[ratio_rx_g < 0.0] = 0.0
    ratio_rx_g[ratio_rx_g == np.inf] = 0.0
    ratio_rx_g[np.isnan(ratio_rx_g)] = 0.0

    if isinstance(nuc, int):
        xs_cache[key] = ratio_rx_g
    return ratio_rx_g

def make_atomic_weight_table(nuc_data, build_dir=""):
    """Makes an atomic weight table in the nuc_data library.

    Parameters
    ----------
    nuc_data : str
        Path to nuclide data file.
    build_dir : str
        Directory to place html files in.
    """
    # Grab raw data
    atomic_abund = parse_atomic_abund(build_dir)
    atomic_masses = parse_atmoic_mass_adjustment(build_dir)

    A = {}

    # Add normal isotopes to A
    for nuc_zz, mass, error in atomic_masses:
        try:
            nuc_name = nucname.name(nuc_zz)
        except RuntimeError:
            continue

        if nuc_zz in atomic_abund:
            A[nuc_zz] = nuc_name, nuc_zz, mass, error, atomic_abund[nuc_zz]
        else:
            A[nuc_zz] = nuc_name, nuc_zz, mass, error, 0.0

    # Add naturally occuring elements
    for element in nucname.name_zz:
        nuc_zz = nucname.zzaaam(element)
        A[nuc_zz] = element, nuc_zz, 0.0, 0.0, 0.0

    for nuc, abund in atomic_abund.items():
        zz = nuc / 10000
        element_zz = zz * 10000
        element = nucname.zz_name[zz]

        nuc_name, nuc_zz, nuc_mass, _error, _abund = A[nuc]
        elem_name, elem_zz, elem_mass, _error, _abund = A[element_zz]

        new_elem_mass = elem_mass + (nuc_mass * abund)
        A[element_zz] = element, element_zz, new_elem_mass, 0.0, 0.0

    A = sorted(A.values(), key=lambda x: x[1])
    # A = np.array(A, dtype=atomic_weight_dtype)

    # Open the HDF5 File
    kdb = tb.openFile(nuc_data, "a", filters=BASIC_FILTERS)

    # Make a new the table
    Atable = kdb.createTable("/", "atomic_weight", atomic_weight_desc, "Atomic Weight Data [amu]", expectedrows=len(A))
    Atable.append(A)

    # Ensure that data was written to table
    Atable.flush()

    # Close the hdf5 file
    kdb.close()

def sigma_a_reaction(nuc, rx, E_g=None, E_n=None, phi_n=None):
    """Calculates the neutron absorption reaction cross-section for a nuclide for a 
    new, lower resolution group structure using a higher fidelity flux.  Note that 
    g indexes G, n indexes N, and G < N.

    Parameters
    ----------
    nuc : int, str, Material, or dict-like 
        A nuclide or nuclide-atom fraction mapping for which to calculate the 
        absorption reaction cross-section.
    rx : str
        Reaction key. ('gamma', 'alpha', 'p', etc.)
    E_g : array-like of floats, optional
        New, lower fidelity energy group structure [MeV] that is of length G+1. 
    E_n : array-like of floats, optional
        Higher resolution energy group structure [MeV] that is of length N+1. 
    phi_n : array-like of floats, optional
        The high-fidelity flux [n/cm^2/s] to collapse the fission cross-section over.  
        Length N.  

    Returns
    -------
    sigma_rx_g : ndarray 
        An array of the collapsed absorption reaction cross section.

    Notes
    -----
    This always pulls the absorption reaction cross-section out of the nuc_data.    

    See Also
    --------
    pyne.xs.cache.ABSORPTION_RX
    pyne.xs.cache.ABSORPTION_RX_MAP 
    """
    _prep_cache(E_g, E_n, phi_n)

    if isinstance(nuc, collections.Iterable) and not isinstance(nuc, basestring):
        return _atom_weight_channel(sigma_a_reaction, nuc, rx)

    # Get the absorption XS
    nuc_zz = nucname.zzaaam(nuc)
    key_n = ('sigma_rx_n', nuc_zz, rx)
    key_g = ('sigma_rx_g', nuc_zz, rx)

    # Don't recalculate anything if you don't have to
    if key_g in xs_cache:
        return xs_cache[key_g]
    else:
        sigma_rx_n = xs_cache[key_n]

    # Perform the group collapse, knowing that the right data is in the cache
    sigma_rx_g = group_collapse(sigma_rx_n, xs_cache['phi_n'], phi_g=xs_cache['phi_g'], 
                                partial_energies=xs_cache['partial_energy_matrix'])

    # Put this value back into the cache, with the appropriate label
    xs_cache[key_g] = sigma_rx_g

    return sigma_rx_g

def load_nuc_file(path):
    """Takes a file that contains whitespace separated nuclide names and
    returns the zzaaam representation as a sorted list."""
    with open(path, 'r') as f:
        s = f.read()

    nuc_list = [nucname.zzaaam(nuc) for nuc in s.split()]
    nuc_list.sort()
    return nuc_list

def parse_for_all_isotopes(htmlfile):
    """Parses an elemental html file, returning a set of all occuring isotopes."""
    isos = set()
    with open(htmlfile, 'r') as f:
        for line in f:
            m = all_iso_regex.search(line)
            if m is not None:
                isos.add(nucname.zzaaam(m.group(1)))
    return isos

def sigma_t(nuc, T=300.0, E_g=None, E_n=None, phi_n=None):
    """Calculates the total neutron cross section for a nuclide. 

    .. math::
        \\sigma_{t, g} = \\sigma_{a, g} + \\sigma_{s, g}

    Parameters
    ----------
    nuc : int, str, Material, or dict-like 
        A nuclide or nuclide-atom fraction mapping for which to calculate the 
        total cross section.
    T : float, optional
        Tempurature of the target material [kelvin].
    E_g : array-like of floats, optional
        New, lower fidelity energy group structure [MeV] that is of length G+1. 
    E_n : array-like of floats, optional
        Higher resolution energy group structure [MeV] that is of length N+1. 
    phi_n : array-like of floats, optional
        The high-fidelity flux [n/cm^2/s] to collapse the fission cross-section over.  
        Length N.  

    Returns
    -------
    sig_t_g : ndarray 
        An array of the total cross section.

    """
    _prep_cache(E_g, E_n, phi_n)

    if isinstance(nuc, collections.Iterable) and not isinstance(nuc, basestring):
        return _atom_weight_channel(sigma_t, nuc, T)

    # Get the total XS
    nuc_zz = nucname.zzaaam(nuc)
    key_a = ('sigma_a_g', nuc_zz)
    key_s = ('sigma_t_g', nuc_zz, T)
    key_t = ('sigma_t_g', nuc_zz, T)

    # Don't recalculate anything if you don't have to
    if key_t in xs_cache:
        return xs_cache[key_t]

    # This calculation requires the abosorption cross-section
    if key_a not in xs_cache:
        xs_cache[key_a] = sigma_a(nuc, E_g, E_n, phi_n)

    # This calculation requires the scattering cross-section
    if key_s not in xs_cache:
        xs_cache[key_s] = sigma_s(nuc, T, E_g, E_n, phi_n)

    # Sum over all h indeces
    sig_t_g = xs_cache[key_a] + xs_cache[key_s]

    # Put this value back into the cache, with the appropriate label
    xs_cache[key_t] = sig_t_g
    return sig_t_g

def _to_zzaaam(nuc, m, s):
    nuc_zz = nucname.zzaaam(nuc.strip())
    if m == 'M':
        state = s.strip()
        if 0 < len(state):
            state = int(state)
        else:
            state = 1
        nuc_zz += state
    return nuc_zz

    def _load_reaction(self, nuc, rx, temp=300.0):
        """ 
        Note: EAF data does not use temperature information (temp)

        Parameters
        ----------
        nuc : int
            Nuclide in zzaaam form.
        rx : str 
            Reaction MT # in nnnm form.
            OR: (eventually)
            Reaction key: 'gamma', 'alpha', 'p', etc.

        See Also
        --------
        EAF_RX : list
            List of valid MT #s in the EAF data.
        EAF_RX_MAP : dict
            Dictionary for converting string reaction identifiers to MT #s.

        """
        nuc = nucname.zzaaam(nuc)

        # Munging the rx to an MT#
        try:
            int(rx)
        except ValueError:
            try:
                rx = EAF_RX_MAP[rx]
            except KeyError:
                pass

        # Check if usable rx #
        if rx is None:
            return None
        if str(rx) not in EAF_RX:
            msg = "the reaction '{rx}' is not valid.".format(rx=rx)
            raise IndexError(msg)

        # Grab data
        with tb.openFile(nuc_data, 'r') as f:
            cond = "(nuc_zz == {0}) & (rxnum == '{1}')".format(nuc, rx)
            node = f.root.neutron.eaf_xs.eaf_xs
            rows = [np.array(row['xs']) for row in node.where(cond)]

        if len(rows) == 0:
            rxdata = None
        elif 1 < len(rows):
            rows = np.array(rows)
            rxdata = rows.sum(axis=0)
        else:
            rxdata = rows[0]

        return rxdata

def sigma_f(nuc, E_g=None, E_n=None, phi_n=None):
    """Calculates the neutron fission cross-section for a nuclide for a new, 
    lower resolution group structure using a higher fidelity flux.  Note that 
    g indexes G, n indexes N, and G < N.  If any of these are None-valued, 
    values from the cache are used.  The energy groups and fluxes are normally 
    ordered from highest-to-lowest energy.

    Parameters
    ----------
    nuc : int, str, Material, or dict-like 
        A nuclide or nuclide-atom fraction mapping for which to calculate the 
        fission cross-section.
    E_g : array-like of floats, optional
        New, lower fidelity energy group structure [MeV] that is of length G+1. 
    E_n : array-like of floats, optional
        Higher resolution energy group structure [MeV] that is of length N+1. 
    phi_n : array-like of floats, optional
        The high-fidelity flux [n/cm^2/s] to collapse the fission cross-section over.  
        Length N.  

    Returns
    -------
    sigma_f_g : ndarray 
        An array of the collapsed fission cross-section.

    Notes
    -----
    This always pulls the fission cross-section out of nuc_data library.    

    """
    _prep_cache(E_g, E_n, phi_n)

    if isinstance(nuc, collections.Iterable) and not isinstance(nuc, basestring):
        return _atom_weight_channel(sigma_f, nuc)

    # Get the fission XS
    nuc_zz = nucname.zzaaam(nuc)
    sigma_f_n_nuc_zz = ('sigma_f_n', nuc_zz)
    sigma_f_g_nuc_zz = ('sigma_f_g', nuc_zz)

    # Don't recalculate anything if you don't have to
    if sigma_f_g_nuc_zz in xs_cache:
        return xs_cache[sigma_f_g_nuc_zz]
    else:
        sigma_f_n = xs_cache[sigma_f_n_nuc_zz]

    # Perform the group collapse, knowing that the right data is in the cache
    sigma_f_g = group_collapse(sigma_f_n, xs_cache['phi_n'], phi_g=xs_cache['phi_g'], 
                               partial_energies=xs_cache['partial_energy_matrix'])

    # Put this value back into the cache, with the appropriate label
    xs_cache[sigma_f_g_nuc_zz] = sigma_f_g

    return sigma_f_g

 def elemental_row(row):
     if re.match('[A-Z][a-z]?-?(\d{1,3})?$', row[0]):
         element = nucname.zzaaam(row[0])
         weight_frac = row[3]
         if nucname.name(element) in elemental_mats:
             composition.update(elemental_mats[row[0].upper()])
         else:
             composition[element] = float(weight_frac)
         nuc_zz.add(element)
     else:
         pass

def test_materials():
    keys = ["U235", "U238"]
    assert_equal(set(ff.materials.keys()), set(keys))

    for iso in keys:
        assert_equal(ff.materials[iso].mass, 1.0)
        assert_equal(ff.materials[iso].comp[nucname.zzaaam(iso)], 1.0)

    u235 = Material({922350: 1.0}, 1.0, name="U-235")
    u238 = Material({922380: 1.0}, 1.0, name="U-238")
    o16  = Material({80160:  1.0}, 1.0, name="O-16")
    mats = {"U235": u235, "U238": u238, "O16": o16}
    ff.materials = mats

    keys = ["U235", "U238", "O16"]
    assert_equal(set(ff.materials.keys()), set(keys))

    for iso in keys:
        assert_equal(ff.materials[iso].mass, 1.0)
        assert_equal(ff.materials[iso].comp[nucname.zzaaam(iso)], 1.0)

def test_sigma_s():
    if not hasattr(rx_h5.root, 'sigma_s'):
        raise nose.SkipTest

    for nuc in nucs:    
        nuc_zz = nucname.zzaaam(nuc)

        sig_t_arr, sig_t = read_array(rx_h5.root.sigma_t, nuc)
        sig_s_arr, sig_s = read_array(rx_h5.root.sigma_s, nuc)
        sig_s_gh_arr, sig_s_gh = read_array(rx_h5.root.sigma_s_gh, nuc)

        yield check_le, sig_s, sig_t, [sig_s_arr._v_pathname, sig_t_arr._v_pathname]
        yield check_array_almost_eq, sig_s, sig_s_gh.sum(axis=-2), [sig_s_arr._v_pathname, 'sum(' + sig_s_gh_arr._v_pathname + ')']

def sigma_s(nuc, T, E_g=None, E_n=None, phi_n=None):
    """Calculates the neutron scattering cross-section for a nuclide. 

    .. math::
        \\sigma_{s, g} = \\sum_{h} \\sigma_{s, g\\to h} 

    Parameters
    ----------
    nuc : int, str, Material, or dict-like 
        A nuclide or nuclide-atom fraction mapping for which to calculate the 
        scattering cross section.
    T : float
        Tempurature of the target material [kelvin].
    E_g : array-like of floats, optional
        New, lower fidelity energy group structure [MeV] that is of length G+1. 
    E_n : array-like of floats, optional
        Higher resolution energy group structure [MeV] that is of length N+1. 
    phi_n : array-like of floats, optional
        The high-fidelity flux [n/cm^2/s] to collapse the fission cross-section over.  
        Length N.  

    Returns
    -------
    sig_s_g : ndarray 
        An array of the scattering cross section.

    """
    _prep_cache(E_g, E_n, phi_n)

    if isinstance(nuc, collections.Iterable) and not isinstance(nuc, basestring):
        return _atom_weight_channel(sigma_s, nuc, T)

    nuc_zz = nucname.zzaaam(nuc)
    key_g = ('sigma_s_g', nuc_zz, T)
    key_gh = ('sigma_s_gh', nuc_zz, T)

    # Don't recalculate anything if you don't have to
    if key_g in xs_cache:
        return xs_cache[key_g]

    # This calculation requires the scattering kernel
    if key_gh not in xs_cache:
        xs_cache[key_gh] = sigma_s_gh(nuc, T, E_g, E_n, phi_n)

    # Sum over all h
    sig_s_g = xs_cache[key_gh].sum(axis=1)

    # Put this value back into the cache, with the appropriate label
    xs_cache[key_g] = sig_s_g

    return sig_s_g