docs: def _tp_finder(self, dcycle):   # Private routine
          """
          Routine to find thermal pulses in given star and returns an
          index vector that gives the cycle number in which the thermal
          pulse occure.
  
          The routine looks for the C/O ratio jumping up and up, so only
          useful in TP-AGB star.  A vector is given back that indicates
          the position of the cycle that is at 95% of the thermal pulse
          (to make sure it's not in the next one and that most of the
          processing is done). The script also returns the co_ratio
          vector - the C/O ratio (number fraction) at the given thermal
          pulse.
  
          """
          # read in c and o isotopes for all cycles, regarding deltacycle
          last_cycle = int(self.se.cycles[len(self.se.cycles)-1])
          cyc_tp = list(range(1,last_cycle + dcycle, dcycle))
          all_data = array(self.get(cyc_tp,['C-12','C-13','O-16','O-17','O-18']))
          c_nf = np.zeros(len(all_data))
          o_nf = np.zeros(len(all_data))
          for i in range(len(all_data)):
              c_nf[i] = all_data[i][0] + all_data[i][1]
              o_nf[i] = all_data[i][2] + all_data[i][3] + all_data[i][4]
  
          # search for thermal pulses
          co_ratio = (old_div(c_nf, o_nf)) * 15.9994 / 12.0107
          tp_guess     = 200   # this should be an upper limit!
          tp_guess_max = 200   # to through an error
          # guess variables, i is the actual break criterion, n a max counter
          gi = 0
          gn = 0
  
          while gi != 1 and gn < 10000:
              tp_ind = list()
              i = 0
              while i < len(co_ratio)-2:
                  gcompar= old_div(1., (dcycle*tp_guess*100.))
                  slope1 = old_div((co_ratio[i+1]-co_ratio[i]),(dcycle))
                  slope2 = old_div((co_ratio[i+2]-co_ratio[i+1]),dcycle)
                  if slope1 > gcompar and slope2 < gcompar and co_ratio[i+1] > co_ratio[i]:
                      tp_ind.append(i+1)
                      i += 3   # jump three cycles to avoid defining a single cycle twice!
                  else:
                      i += 1
  
              if abs(len(tp_ind) - tp_guess) < old_div(tp_guess,2):   # gotta be within factor two of guess
                  gi = 1
              else:
                  gn += 1
                  tp_guess /= 2
          # check w/ maximum of thermal pulses allowed
          if len(tp_ind) > tp_guess_max:
              print('Problem detected with number of pulses')
          # create thermal pulse vector
          tp_startf = zeros(len(tp_ind))   # found start
          for i in range(len(tp_startf)):
              tp_startf[i] = cyc_tp[tp_ind[i]]
          # read out isotopic composition at 95% of the thermal pulse and the initial of the star
          # set up thermal pulse positions
          tp_limits = zeros(len(tp_startf)+1)
          for i in range(len(tp_startf)):
              tp_limits[i] = tp_startf[i]
          tp_limits[len(tp_limits)-1] = int(self.se.cycles[len(self.se.cycles)-1])
          # thermal pulse position (where to read the isotope ratio)
          tp_pos = list()
          for i in range(len(tp_startf)):
              tp_pos.append(int(tp_limits[i] + 0.95 * (tp_limits[i+1] - tp_limits[i])))
          # create co_ret vector to return c/o ratio vector
          co_return = zeros(len(tp_pos))
          for i in range(len(tp_pos)):
              co_return[i] = co_ratio[tp_ind[i]]
          # return the two vectors
          return tp_pos,co_return