docs: def cores(self, incycle, **keyw):
          """
          This function uses the abundances as a function of time to
          return core masses.
  
          Parameters
          ----------
          incycle : integer
              incycle is the cycle to choose where to take the core
              masses.
          keyw : dict
              A dict of key word arguments.
  
          Returns
          -------
          list
              [cores, core_type, core_info]
              cores contains the core masses for the types found in
              core_type.  Core_info stores a variety of information about
              the abundances at certain mass coordinates.
  
          Notes
          -----
          These core masses are:
  
          +----------+-----------------------------------------+
          | M_alpha  | mass of the helium core                 |
          +----------+-----------------------------------------+
          | M_CO     | mass of the carbon-oxygen core          |
          +----------+-----------------------------------------+
          | M_si     | mass of the silicon core                |
          +----------+-----------------------------------------+
          | M_fe     | mass of the iron core                   |
          +----------+-----------------------------------------+
          | M_rem    | mass of the remnant star post-supernova |
          +----------+-----------------------------------------+
  
  
          The following keywords can also be used:
  
          +------------------+------------------+
          | Keyword Argument | Default Value    |
          +==================+==================+
          | abund            | "iso_massf"      |
          +------------------+------------------+
          | isoa             | "A"              |
          +------------------+------------------+
          | isom             | "isomeric_state" |
          +------------------+------------------+
          | isoz             | "Z"              |
          +------------------+------------------+
          | mass             | "mass"           |
          +------------------+------------------+
          | core_opt         | 1                |
          +------------------+------------------+
  
          All arguments change the name of fields used to read data from
          HDF5 files, other than core_opt.  core_opt controls which
          scheme to use for calculating the core masses.
  
              If core_opt is 0, then cores uses alpha-isotopes to
              calculate the silicon and iron cores.  Use this for stellar
              evolution output.
  
                  Si_core = Si28+S32+Ar36+Ca40+Ti44
  
                  Fe_core = Cr48+Fe52+Ni56
  
              If core_opt is 1, then cores uses all elements in a network
              with proton number Z=23 as a boundary.  Use this for MPPNP
              output.
  
                  Si_core = Sum of all isotopes with 14 <= Z <= 22
  
                  Fe_core = Sum of all isotopes with 23 <= Z <= 28
  
          """
  
          def infomod(core_opt, *inp):
              """
              This calls the correct infomod function depending on the
              value of core_opt.
  
              """
              if core_opt == 0:
                  return infomod1(inp[0], inp[1], inp[2])
              elif core_opt == 1:
                  return infomod2(inp[0], inp[1], inp[2], inp[3], inp[4])
  
          def infomod1(shell, yps, isoindex):
              """
              Function for defining data to print into a string.  This is
              used for the case of core_opt = 0.
  
              """
  
              xh = yps[shell][isoindex[0]]
              xhe = yps[shell][isoindex[1]]
              xc = yps[shell][isoindex[2]]
              xo = yps[shell][isoindex[3]]
              xsi = yps[shell][isoindex[4]]
              xs = yps[shell][isoindex[5]]
              xa = yps[shell][isoindex[6]]
              xca = yps[shell][isoindex[7]]
              xti = yps[shell][isoindex[8]]
              xcr = yps[shell][isoindex[9]]
              xfe = yps[shell][isoindex[10]]
              xni = yps[shell][isoindex[11]]
  
              xsicore = xsi + xs + xa + xca + xti
              xfecore = xcr + xfe + xni
  
              return ' h=' + "%12.4e"%(xh) + ' he=' + "%12.4e"%(xhe) + \
              ' c=' + "%12.4e"%(xc) + ' o=' + "%12.4e"%(xo) + \
              ' si=' + "%12.4e"%(xsicore) + ' fe=' + "%12.4e"%(xfecore)
  
          def infomod2(shell, yps, isoindex, xsicore, xfecore):
              """
              Function for defining data to print into a string.  This is
              used for the case of core_opt = 1.
  
              """
              xsicore = 0.
              xsicoren = 0.
              xfecore = 0.
              xfecoren = 0.
  
              xh = yps[shell][isoindex[0]]
              xhe = yps[shell][isoindex[1]]
              xc = yps[shell][isoindex[2]]
              xo = yps[shell][isoindex[3]]
  
              for j in iso_si:
                  xsicore += yps[i][j]
                  xsicoren += yps[i+iter][j]
              for j in iso_fe:
                  xfecore += yps[i][j]
                  xfecoren += yps[i+iter][j]
  
              return 'shellnb = ' + str(shell) + ' h=' + "%12.4e"%(xh) + \
              ' he=' + "%12.4e"%(xhe) + ' c=' + "%12.4e"%(xc) + \
              ' o=' + "%12.4e"%(xo) + ' si=' + "%12.4e"%(xsicore) + \
              ' fe=' + "%12.4e"%(xfecore)
  
          if ("isoa" in keyw) == False:
              keyw["isoa"] = "A"
          if ("isoz" in keyw) == False:
              keyw["isoz"] = "Z"
          if ("mass" in keyw) == False:
              keyw["mass"] = "mass"
          if ("abund" in keyw) == False:
              keyw["abund"] = "iso_massf"
          if ("coreopt" in keyw) == False:
              core_opt = 1
          else:
              core_opt = keyw["coreopt"]
  
          cores = [0,0,0,0,0,0]
          core_type = ["","","","","",""]
          core_info = []
          iso_si = []
          iso_fe = []
          first = True
  
          xm_list     = self.se.get(incycle, keyw["mass"])
          abund_list  = self.se.get(incycle, keyw["abund"])
          tables = self.se.Tables
          isoa_list   = tables[0]
          isoz_list   = tables[1]
  
          # Check the order of the input data
          xmass = xm_list
          yps = abund_list
          centre = 0
          surface = len(xmass)-1
          iter = -1
          if isinstance(xmass, list) == True:
              if xmass[0] > xmass[1]:
                  # mass is descending with shell number and the centre of the star
                  # is the last shell
                  centre = len(xmass)-1
                  surface = 0
                  iter = 1
  
          mco01 = 0.
          na    = 0
          nco   = 0
  
          h1index = -1
          he4index = -1
          c12index = -1
          o16index = -1
          si28index = -1
          s32index = -1
          a36index = -1
          ca40index = -1
          ti44index = -1
          cr48index = -1
          fe52index = -1
          ni56index = -1
  
          # Determine yps indices for certain abundances
          for i in range(len(isoa_list)):
              try:
                  A = int(isoa_list[i])
                  Z = int(isoz_list[i])
              except TypeError:
                  A = int(isoa_list[i][0])
                  Z = int(isoz_list[i][0])
  
              if A == 1 and Z == 1:
                  h1index = i
              elif A == 4 and Z == 2:
                  he4index = i
              elif A == 12 and Z == 6:
                  c12index = i
              elif A == 16 and Z == 8:
                  o16index = i
  
              if core_opt == 0:
                  if A == 28 and Z == 14:
                      si28index = i
                  elif A == 32 and Z == 16:
                      s32index = i
                  elif A == 36 and Z == 18:
                      a36index = i
                  elif A == 40 and Z == 20:
                      ca40index = i
                  elif A == 44 and Z == 22:
                      ti44index = i
                  elif A == 48 and Z == 24:
                      cr48index = i
                  elif A == 52 and Z == 26:
                      fe52index = i
                  elif A == 56 and Z == 38:
                      ni56index = i
  
          if h1index == -1 or he4index == -1 or c12index == -1 or o16index == -1:
              print("A key isotope(s) is not found in network!  Please check for the \
              presence of H1, He4, C12 and O16.")
              os.sys.exit()
  
          # Si_core = Si28+S32+Ar36+Ca40+Ti44
          # Fe_core = Cr48+Fe52+Ni56
          if core_opt == 0:
              if si28index == -1 or s32index == -1 or a36index == -1 or ca40index == -1 \
              or ti44index == -1 or cr48index == -1 or fe52index == -1 or ni56index == -1:
                  print("Key isotopes for measuring the core mass with core_opt = 0 \
                  are missing.  Setting core_opt = 1.")
                  core_opt = 1
  
          # Si_core = Sum of all isotopes with 14 <= Z <= 22
          # Fe_core = Sum of all isotopes with 23 <= Z <= 28
          if core_opt == 1:
              for i in range(len(isoa_list)):
                  try:
                      A = int(isoa_list[i])
                      Z = int(isoz_list[i])
                  except TypeError:
                      A = int(isoa_list[i][0])
                      Z = int(isoz_list[i][0])
                  if Z >= 14 and Z <= 22:
                      iso_si.append(i)
                  if Z >= 23 and Z <= 28:
                      iso_fe.append(i)
  
          isoindex = [h1index, he4index, c12index, o16index]
          if core_opt == 0:
              isoindex.extend([si28index, s32index, a36index, ca40index,
              ti44index, cr48index, fe52index, ni56index])
  
          if first == True:
              first = False
              cores[0] = xmass[surface]
              core_type[0] = "Final"
              core_info.append(infomod(core_opt, surface, yps, isoindex, iso_si, iso_fe))
  
          # Iterate over shells to determine the core masses
          for i in np.arange(surface, centre, iter):
  
              xsicore = 0.
              xfecore = 0.
              xsicoren = 0.
              xfecoren = 0.
  
              # Abundances of key isotopes at this shell (and the next one)
              xh = yps[i][h1index]
              xhe = yps[i][he4index]
              xc = yps[i][c12index]
              xo = yps[i][o16index]
  
              xhn = yps[i+iter][h1index]
              xhen = yps[i+iter][he4index]
              xcn = yps[i+iter][c12index]
              xon = yps[i+iter][o16index]
  
              xcocore = xc + xo
              xcocoren = xcn + xon
  
              if core_opt == 0:
                  xsi = yps[i][si28index]
                  xs = yps[i][s32index]
                  xa = yps[i][a36index]
                  xca = yps[i][ca40index]
                  xti = yps[i][ti44index]
                  xcr = yps[i][cr48index]
                  xfe = yps[i][fe52index]
                  xni = yps[i][ni56index]
  
                  xsin = yps[i+iter][si28index]
                  xsn = yps[i+iter][s32index]
                  xan = yps[i+iter][a36index]
                  xcan = yps[i+iter][ca40index]
                  xtin = yps[i+iter][ti44index]
                  xcrn = yps[i+iter][cr48index]
                  xfen = yps[i+iter][fe52index]
                  xnin = yps[i+iter][ni56index]
  
                  xsicore = xsi + xs + xa + xca + xti
                  xfecore = xcr + xfe + xni
  
                  xsicoren = xsin + xsn + xan + xcan + xtin
                  xfecoren = xcrn + xfen + xnin
  
              elif core_opt == 1:
                  for j in iso_si:
                      xsicore += yps[i][j]
                      xsicoren += yps[i+iter][j]
                  for j in iso_fe:
                      xfecore += yps[i][j]
                      xfecoren += yps[i+iter][j]
  
              # M_alpha:
              if xh >= 1.e-9:
                  if xhen > 0.75 and xhe <= 0.75:
                      if cores[1] == 0:
                          cores[1] = xmass[i]
                          core_type[1] = 'Malpha 75%'
                          core_info.append('Malpha 75% '+str(xmass[i-iter])+ \
                          infomod(core_opt, i-iter, yps, isoindex, iso_si, iso_fe))
                          core_info.append('Malpha 75% '+str(xmass[i])+ \
                          infomod(core_opt, i, yps, isoindex, iso_si, iso_fe))
                          core_info.append('Malpha 75% '+str(xmass[i+iter])+ \
                          infomod(core_opt, i+iter, yps, isoindex, iso_si, iso_fe))
                  if xhen > 0.5 and xhe <= 0.5:
                      core_info.append('Malpha 50% '+str(xmass[i])+ \
                      infomod(core_opt, i, yps, isoindex, iso_si, iso_fe))
                  if xhn < 1.e-2 and xh >= 1.e-2:
                      na=i
                      core_info.append('Malpha 1.e-2 '+str(xmass[i])+ \
                      infomod(core_opt, i, yps, isoindex, iso_si, iso_fe))
                  if xhn < 1.e-4 and xh >= 1.e-4:
                      core_info.append('Malpha 1.e-4 '+str(xmass[i])+ \
                      infomod(core_opt, i, yps, isoindex, iso_si, iso_fe))
                  if xhn < 1.e-5 and xh >= 1.e-5:
                      core_info.append('Malpha 1.e-5 '+str(xmass[i])+ \
                      infomod(core_opt, i, yps, isoindex, iso_si, iso_fe))
  
              # M_CO:
              if xhe >= 1.e-9 and xh < 1.e-9:
                  if xcocoren > 0.5 and xcocore <= 0.5:
                      core_info.append('Mco 50% '+str(xmass[i])+ \
                      infomod(core_opt, i, yps, isoindex, iso_si, iso_fe))
                  if xcocoren > 0.75 and xcocore <= 0.75:
                      core_info.append('Mco 75% '+str(xmass[i])+ \
                      infomod(core_opt, i, yps, isoindex, iso_si, iso_fe))
                  if xhen < 1.e-2 and xhe >= 1.e-2:
                      if cores[2] == 0:
                          mco01=xmass[i]
                          nco=i
                          cores[2] = xmass[i]
                          core_type[2] = 'Mco'
                          core_info.append('Mco 1.e-2 '+str(xmass[i-iter])+ \
                          infomod(core_opt, i-iter, yps, isoindex, iso_si, iso_fe))
                          core_info.append('Mco 1.e-2 '+str(xmass[i])+ \
                          infomod(core_opt, i, yps, isoindex, iso_si, iso_fe))
                          core_info.append('Mco 1.e-2 '+str(xmass[i+iter])+ \
                          infomod(core_opt, i+iter, yps, isoindex, iso_si, iso_fe))
                  if xhen < 1.e-4 and xhe >= 1.e-4:
                      core_info.append('Mco 1.e-4 '+str(xmass[i])+ \
                      infomod(core_opt, i, yps, isoindex, iso_si, iso_fe))
                  if xhen < 1.e-5 and xhe >= 1.e-5:
                      core_info.append('Mco 1.e-5 '+str(xmass[i])+ \
                      infomod(core_opt, i, yps, isoindex, iso_si, iso_fe))
  
              # M_Si:
              if xcocore >= 1.e-7:
                  if xsicoren > 0.5 and xsicore <= 0.5:
                      if cores[3] == 0:
                          core_info.append('Msi 50% '+str(xmass[i-iter])+ \
                          infomod(core_opt, i-iter, yps, isoindex, iso_si, iso_fe))
                          core_info.append('Msi 50% '+str(xmass[i])+ \
                          infomod(core_opt, i, yps, isoindex, iso_si, iso_fe))
                          core_info.append('Msi 50% '+str(xmass[i+iter])+ \
                          infomod(core_opt, i+iter, yps, isoindex, iso_si, iso_fe))
                          cores[3] = xmass[i]
                          core_type[3] = 'Msi 50%'
                  if xsicoren > 0.75 and xsicore <= 0.75:
                      core_info.append('Msi 75% '+str(xmass[i])+ \
                      infomod(core_opt, i, yps, isoindex, iso_si, iso_fe))
                  if xcocoren > 1.e-2 and xcocore <= 1.e-2:
                      core_info.append('Msi 1.e-2 '+str(xmass[i])+ \
                      infomod(core_opt, i, yps, isoindex, iso_si, iso_fe))
  
              # M_fe:
              if xfecoren > 0.5 and xfecore <= 0.5:
                  if cores[4] == 0:
                      core_info.append('Mfe 50% '+str(xmass[i-iter])+ \
                      infomod(core_opt, i-iter, yps, isoindex, iso_si, iso_fe))
                      core_info.append('Mfe 50% '+str(xmass[i])+ \
                      infomod(core_opt, i, yps, isoindex, iso_si, iso_fe))
                      core_info.append('Mfe 50% '+str(xmass[i+iter])+ \
                      infomod(core_opt, i+iter, yps, isoindex, iso_si, iso_fe))
                      cores[4] = xmass[i]
                      core_type[4] = 'Mfe 50%'
              if xfecoren > 0.75 and xfecore <= 0.75:
                  core_info.append('Mfe 75% '+str(xmass[i])+ \
                  infomod(core_opt, i, yps, isoindex, iso_si, iso_fe))
              if xsicoren > 1.e-2 and xsicore <= 1.e-2:
                  core_info.append('Mfe 1.e-2 '+str(xmass[i])+ \
                  infomod(core_opt, i, yps, isoindex, iso_si, iso_fe))
  
          mcof=mco01
          i=na
          while i < nco:
              i+=1
              xco = yps[i][c12index] + yps[i][o16index]
              xcop = yps[i-iter][c12index] + yps[i-iter][o16index]
              mcof+=(xco+xcop)/2.*(xmass[i-iter]-xmass[i])
          cores[2] = mcof
          core_type[2] = 'Mco int'
  
          # M_rem vs M_CO relationship
          xmco = np.array([0.72,1.01,1.38,2.10,2.97,4.94,7.20,13.88,24.75,38.12,58.73])
          xmre = np.array([0.72,1.01,1.23,1.42,1.65,2.15,2.72,4.3,7.6,11.4,17.1])
  
          if 'Mco' in core_type:
              for i in np.arange(len(xmco)):
                  if cores[2] < xmco[i]:
                      xco = yps[i][c12index] + yps[i][o16index]
                      xcop = yps[i-iter][c12index] + yps[i-iter][o16index]
                      cores[5] = xmre[i-1]+(xmre[i]-xmre[i-1])*(cores[j]-xmco[i-1])/(xmco[i]-xmco[i-1])
                      core_type[5] = 'Mrem'
                      break
  
          return [cores, core_type, core_info]