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