docs: def find_TPs_and_DUPs(self, percent=5., makefig=False):
          """
          Function which finds TPs and uses the calc_DUP_parameter
          function.  To calculate DUP parameter evolution dependent of
          the star or core mass.
  
          Parameters
          ----------
          fig : integer
              Figure number to plot.
          t0_model : integer
              First he-shell lum peak.
          percent : float
              dredge-up is defined as when the mass dredged up is a certain
              percent of the total mass dredged up during that event, which
              is set by the user in this variable.
              The default is 5.
          makefig :
              do you want a figure to be made?
  
          Returns
          -------
          TPmods : array
              model numbers at the peak of each thermal pulse
          DUPmods : array
              model numbers at the dredge-up, where dredge-up is defined as
              when the mass dredged up is a certain percent of the total mass
              dredged up during that event, which is set by the user
          TPend : array
              model numbers at the end of the PDCZ for each TP
          lambda : array
              DUP efficiency for each pulse
          """
  
          t0_model=self.find_first_TP()
          t0_idx=(t0_model-self.get("model_number")[0])
          first_TP_he_lum=10**(self.get("log_LHe")[t0_idx])
          he_lum=10**(self.get("log_LHe")[t0_idx:])
          h_lum=10**(self.get("log_LH")[t0_idx:])
          model=self.get("model_number")[t0_idx:]
  
          try:
             h1_bndry=self.get("h1_boundary_mass")[t0_idx:]
          except:
             try:
                  h1_bndry=self.get('he_core_mass')[t0_idx:]
             except:
                  pass
          # SJ find TPs by finding local maxima in He-burning luminosity and
          # checking that the he_lum is greater than the h_lum:
          maxima=[0]
          for i in range(2,len(model)-1):
              if he_lum[i] > he_lum[i-1] and he_lum[i] > he_lum[i+1]:
                  if he_lum[i-1] > he_lum[i-2] and he_lum[i+1] > he_lum[i+2]:
                      if he_lum[i] > h_lum[i]:
                          maxima.append(i)
  
          # find DUPs when h-boundary first decreases by more than XX% of the total DUP
          # depth:
          DUPs=[]
          TPend=[]
          maxDUPs=[]
          for i in range(len(maxima)):
              idx1=maxima[i]
              try:
                  idx2=maxima[i+1]
              except IndexError:
                  idx2=-1
              bound=h1_bndry[idx1:idx2]
              bound0=bound[0]
              if bound0==min(bound) or bound0 < min(bound): # then no DUP
                  DUP=idx1
                  DUPs.append(DUP)
                  maxDUPs.append(DUP)
              else:
                  maxDUPs.append(idx1+bound.argmin()) # model number of deepest extend of 3DUP
                  maxDUP=bound0-min(bound) # total mass dredged up in DUP
                  db=bound - bound[0]
                  db_maxDUP = old_div(db, maxDUP)
                  DUP=np.where(db_maxDUP <= old_div(-float(percent),100.))[0][0]
                  DUPs.append(DUP+idx1)
  #                # Alternative definition, where envelope reaches mass coordinate
  #                # where top of PDCZ had resided during the TP:
  #                top=self.get('mx2_top')[idx1]
  #                DUP=np.abs(bound-top).argmin()
  #                DUPs.append(DUP+idx1)
  
          # find end of PDCZ by seeking from TP peak and checking mx2_bot:
              mx2b=self.get('mx2_bot')[t0_idx:][idx1:idx2]
              for i in range(len(mx2b)):
                  if mx2b[i]==0.:
                      endTP=i+idx1
                      TPend.append(endTP)
                      break
  
          # 3DUP efficiency:
          lambd=[0.]
          for i in range(1,len(maxima)):
              dmenv = h1_bndry[maxima[i]] - h1_bndry[maxDUPs[i-1]]
              dmdredge = h1_bndry[maxima[i]] - h1_bndry[maxDUPs[i]]
              lambd.append(old_div(dmdredge,dmenv))
  
          TPmods = maxima + t0_idx
          DUPmods = DUPs + t0_idx
          TPend = TPend + t0_idx
  
          return TPmods, DUPmods, TPend, lambd