docs: def find_first_TP(self):
  
          """
          Find first TP of the TPAGB phase and returns the model
          number at its LHe maximum.
  
          Parameters
          ----------
  
          """
  
          star_mass         = self.get('star_mass')
          he_lumi           = self.get('log_LHe')
          h_lumi            = self.get('log_LH')
          mx2_bot           = self.get('mx2_bot')*star_mass
          try:
             h1_boundary_mass  = self.get('h1_boundary_mass')
             he4_boundary_mass = self.get('he4_boundary_mass')
          except:
             try:
                  h1_boundary_mass  = self.get('he_core_mass')
                  he4_boundary_mass = self.get('c_core_mass')
             except:
                  pass
  
  
          TP_bot=np.array(self.get('conv_mx2_bot'))*np.array(self.get('star_mass'))
          TP_top=np.array(self.get('conv_mx2_top'))*np.array(self.get('star_mass'))
          lum_array=[]
          activate=False
          models=[]
          pdcz_size=[]
          for i in range(len(h1_boundary_mass)):
              if (h1_boundary_mass[i]-he4_boundary_mass[i] <0.2) and (he4_boundary_mass[i]>0.2):
                  if (mx2_bot[i]>he4_boundary_mass[i]) and (he_lumi[i]>h_lumi[i]):
                          if TP_top[i]>he4_boundary_mass[i]:
                                  pdcz_size.append(TP_top[i]-TP_bot[i])
                                  activate=True
                                  lum_array.append(he_lumi[i])
                                  models.append(i)
                                  #print(TP_bot[i],TP_top[i])
                  if (activate == True) and (he_lumi[i]<h_lumi[i]):
                          #if fake tp
                          if max(pdcz_size)<1e-5:
                                  active=False
                                  lum_array=[]
                                  models=[]
                                  print('fake tp')
                          else:
                                  break
          t0_model = models[np.argmax(lum_array)]
          return t0_model