docs: 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