def presnyields(self, *cycles, **keyw):
        """
        This function calculates the presupernova yields of a full
        structure profile from a remnant mass, mrem, to the surface.

        Parameters
        ----------
        cycles : variadic tuple
            cycle[0] is the cycle to perform the presupernova yields
            calculations on.  If cycle[1] is also specified, the yields
            are outputted using 'initial' abundances from cycle[1],
            otherwise the ejected masses are outputted.
        keyw : dict
            A dict of key word arguments.

        Notes
        -----
        The following keywords can be used:

        +------------------+---------------+
        | Keyword Argument | Default Value |
        +==================+===============+
        | abund            | "iso_massf"   |
        +------------------+---------------+
        | xm               | "mass"        |
        +------------------+---------------+
        | mrem             | 0             |
        +------------------+---------------+

        abund and xm are used when the variables within the input file
        differ in their names.  The default values are set to the
        output typically found in an MPPNP output file.  For example,
        if the table for the abundances is called "abundances" instead
        of the default value, use abund = "abundances" as a keyword
        argument.

        mrem is specified using a keyword argument and tells the program
        where to begin integrating.

        """
        abund_list  = []
        xm_list  = []

        if ("xm" in keyw) == False:
            keyw["xm"] = "mass"
        if ("abund" in keyw) == False:
            keyw["abund"] = "iso_massf"
        if ("mrem" in keyw) == False:
            mrem = 0.
        else:
            mrem = keyw["mrem"]

        # Only two cycles are required in this program.
        # Any more will be ignored
        cylen = len(cycles)
        if cylen > 2:
            cylen = 2
        for i in range(cylen):
            cycle = cycles[i]
            abund_list.append(self.se.get(cycle, keyw['abund']))
            xm_list.append(self.se.get(cycle, keyw['xm']))

        isoX = abund_list[0]
        xm = xm_list[0]
        if cylen == 2:
            isoXini = abund_list[1] # initial abundances

        niso = len(isoX[0,:])
        nshells = len(xm)

        X_i = np.zeros([niso], float)
        ybound = np.zeros([niso], float)
        xarray = np.zeros([nshells+1], float)
        yarray = np.zeros([nshells+1], float)

        # This part determines the index of the mass coordinate array which
        # is closest to the specified remnant mass.  This is used in the next
        # part, which is used to interpolate the abundances at the remnant mass
        for k in range(nshells):
            k1 = k
            if mrem<=xm[k]:
                break

        # This part is the interpolation of the isotopes found at the remnant
        # mass.
        for i in range(niso):
            if k1>=1:
                if isoX[k1-1,i]!=0.0:
                    m=old_div((isoX[k1,i]-isoX[k1-1,i]),(xm[k1]-xm[k1-1]))
                    ybound[i] = isoX[k1-1,i] +m*(mrem-xm[k1-1])
                else:
                    ybound[i]=1.e-99
            if k1==0:
                if isoX[k1,i]!=0.0:
                    ybound[i]=isoX[k1,i]
                else:
                    ybound[i]=1.e-99

        # This part merges the interpolated data and the existing arrays into
        # the arrays xarray and yarray.  Once this is done, the summation is
        # made.
        xarray[0] = mrem
        xarray[1:nshells-k1+1] = xm[k1:nshells]
        for i in range(niso):
            yarray[0] = ybound[i]
            for j in range(nshells-k1):
                if isoX[k1+j,i] != 0.0:
                    yarray[j+1] = isoX[k1+j,i]
                else:
                    yarray[j+1] = 1.e-99

            if cylen == 1:
                # Calculate the ejected masses
                for j in range(nshells-k1):
                    X_i[i] = X_i[i] + ((0.5*(yarray[j+1] + yarray[j])) * \
                    (xarray[j+1] - xarray[j]))
            elif cylen == 2:
                # Calculate the SN yield.
                for j in range(nshells-k1):
                    X_i[i] = X_i[i] + ((0.5*(yarray[j+1] + yarray[j]) - isoXini[-1,i]) * \
                    (xarray[j+1] - xarray[j]))


        return X_i