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()
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'.