docs: def write_deltatable(self, filename='default', decayed=True,
                           dcycle=500, iniabufile='../../frames/mppnp/USEEPP/iniab2.0E-02GN93.ppn'):
          """
          This subroutine is to write out tables with delta values for
          cosmochemists to use in comparison with their data.
  
          No options are necessarily needed to load this routine, however,
          it might be useful to specify a filename.  This file furthermore
          searches for thermal pulses, hence, it's only useful for
          non-explosive TP-AGB stars!
  
          Parameters
          ----------
          filename : string, optional
              Choose the filename you want.  The default is "default".
          decayed : boolean, optional
              Value if decayed massfractions of isotope should be taken
              or not.  The default is False.
          dcycle : integer, optional
              Difference between cycles to search for thermal pulses.  The
              default is 500.
          iniabufile : string, optional
              File with initial abundances.  As a standard value, the
              GN93 file is used (in USEEPP folder).  Important, input
              file has to be USEEPP conforming, see NuGrid book! The
              default is '../../frames/mppnp/USEEPP/iniab2.0E-02GN93.ppn'.
  
          Notes
          -----
          Regardless of what value is used as filename it will be replaced
          with 'delta_outfile.txt'.
  
          """
          print('This is a preliminary version - contact Reto for more information!')
          # make filename
          filename = 'delta_outfile.txt'
          # which iso_massmf?
          if decayed:
              isomassmf = 'iso_massf_decay'
          else:
              isomassmf = 'iso_massf'
          # read in thermal pulse position and co_ratio (from private routine)
          tp_pos, co_return = self._tp_finder(dcycle)
          # define isotopes to read out and calculate delta values
          noele = 1   # there is at least one element available
          elelist = list()
          isolisttmp = list()
          isolist = list()
          isoabulist = list()
          isoabulisttmp = list()
          # load file w/ initial abundances
          inut = iniabu(iniabufile)
          # list of elements:
          itmp = inut.names[0][0:2].replace(' ','')
          elelist.append([itmp, inut.z[0]])
          isolisttmp.append(inut.a[0])
          isoabulisttmp.append(inut.abu[0])
  
          for i in range(1,len(inut.names)):
              if inut.names[i][0:2].replace(' ','') != itmp:
                  # add to elelist
                  itmp = inut.names[i][0:2].replace(' ','')
                  elelist.append([itmp, inut.z[i]])
                  # append to isolist
                  isolist.append(isolisttmp)
                  isolisttmp = list()
                  isolisttmp.append(inut.a[i])
                  # append to isoabulist
                  isoabulist.append(isoabulisttmp)
                  isoabulisttmp = list()
                  isoabulisttmp.append(inut.abu[i])
                  # increase number of elements
                  noele += 1
              else:
                  # add to isolisttmp
                  isolisttmp.append(inut.a[i])
                  isoabulisttmp.append(inut.abu[i])
          isolist.append(isolisttmp)
          isoabulist.append(isoabulisttmp)
  
          ### CHOP TEMPORARILY ###
          elelist = elelist[5:len(elelist)-2]
          isolist = isolist[5:len(elelist)-2]
          isoabulist = isoabulist[5:len(elelist)-2]
          # make an index vector for at which position for each line the most abundand isotope is
          isoabuindex = zeros(len(isoabulist),dtype=int)   # as integer!
          for i in range(len(isoabulist)):
              tmp = 0
              for j in range(len(isoabulist[i])):
                  if isoabulist[i][j] > tmp:
                      isoabuindex[i] = j
                      tmp = isoabulist[i][j]
  
          # and now - finally - make a list with all delta values that need to be calculated
          # 0: isotope name, 1: nominator isotope, 2: denominator isotope
          # also make solar system ratios list
          deltalist = list()
          solsysratio = list()
          for i in range(len(isolist)):
              if len(isolist[i]) > 1:
                  for j in range(len(isolist[i])):
                      if j != isoabuindex[i]:
                          deltalist.append([elelist[i][0], isolist[i][j], isolist[i][isoabuindex[i]]])
                          solsysratio.append(old_div(isoabulist[i][j], isoabulist[i][isoabuindex[i]]))
          solsysratio = array(solsysratio)
          # Make an array to write all data into. first column, cycle number, second column, co ratio, then all delta values
          write_out = zeros((len(tp_pos), len(deltalist) + 2))
  
          # now add first three columns to write_out
          for i in range(len(write_out)):
              write_out[i][0] = tp_pos[i]
              write_out[i][1] = co_return[i]
  
          # header for delta values
          deltaheader = list()
  
          ### BIG DELTA LOOP ###
          for deli in range(len(deltalist)):
              ## read in the two isotopes of interest, iso1 and iso2 (iso1 as nominator)
              iso1 = zeros(len(tp_pos))
              iso2 = zeros(len(tp_pos))
              # isotope names:
              if len(deltalist[deli][0]) == 1:
                  iso1name = deltalist[deli][0].upper() + '-' + str(int(deltalist[deli][1]))
                  iso2name = deltalist[deli][0].upper() + '-' + str(int(deltalist[deli][2]))
              else:
                  iso1name = deltalist[deli][0][0].upper() + deltalist[deli][0][1] + '-' + str(int(deltalist[deli][1]))
                  iso2name = deltalist[deli][0][0].upper() + deltalist[deli][0][1] + '-' + str(int(deltalist[deli][2]))
  
              # read in
              iso_data = self.get(tp_pos,isomassmf,[iso1name,iso2name])
              iso1 = np.zeros(len(iso_data))
              iso2 = np.zeros(len(iso_data))
              for i in range(len(iso_data)):
                  iso1[i] = iso_data[i][0]
                  iso2[i] = iso_data[i][1]
              # readin loop
              # for i in range(len(tp_pos)):
              #     iso1[i] = self.get(int(tp_pos[i]),isomassmf,iso1name)
              #     iso2[i] = self.get(int(tp_pos[i]),isomassmf,iso2name)
  
              ## calculate deltavalue array
              deltavalues = (iso1/iso2 / solsysratio[deli] - 1.) * 1000.
  
              ## append values to write_out
              for i in range(len(deltavalues)):
                  write_out[i][deli+2] = deltavalues[i]
  
              ## header
              deltaheader.append('del_' + deltalist[deli][0] + str(int(deltalist[deli][1])) + '_' + deltalist[deli][0] + str(int(deltalist[deli][2])))
  
  
          ### WRITE OUT DATAFILE ###
          outf = open(filename,'w')
  
          # write header
          outf.writelines('cycle_no,'.ljust(18) + 'c_o_ratio,'.ljust(18))
          for i in range(len(deltaheader)):
              if i != len(deltaheader)-1:
                  outfstr = deltaheader[i] + ','
              else:
                  outfstr = deltaheader[i]
              outf.writelines(outfstr.ljust(18))
          outf.writelines('\n')
  
          # write data
          for i in range(len(write_out)):
              for j in range(len(write_out[i])):
                  if j != len(write_out[i])-1:
                      outfstr = format(write_out[i][j],'.7E') + ','
                  else:
                      outfstr = format(write_out[i][j],'.7E')
                  outf.writelines(outfstr.ljust(18))
              outf.writelines('\n')
  
          # close file
          outf.close()