code: def find_actions(results, t, N_matrix=8, use_box=False, ifloop=False, ifprint = True):
      
  
      # Determine orbit class
      loop = assess_angmom(results)
      arethereloops = np.any(loop>0)
      if(arethereloops and not use_box):
          L = loop_actions(flip_coords(results,loop),t,N_matrix, ifprint)
          if(L==None):
              if(ifprint):
                  print("Failed to find actions for this orbit")
              return
          # Used for switching J_2 and J_3 for long-axis loop orbits
          # This is so the orbit classes form a continuous plane in action space
          # if(loop[0]):
          #     L[0][1],L[0][2]=L[0][2],L[0][1]
          #     L[1][1],L[1][2]=L[1][2],L[1][1]
          #     L[1][4],L[1][5]=L[1][5],L[1][4]
          #     L[3].T[1],L[3].T[2]=L[3].T[2],L[3].T[1]
      else:
          L = box_actions(results,t,N_matrix, ifprint)
          if(L==None):
              if(ifprint):
                  print("Failed to find actions for this orbit")
              return
      if(ifloop):
          return L,loop
      else:
          return L
docs: Main routine:
          Takes a series of phase-space points from an orbit integration at times t and returns
          L = (act,ang,n_vec,toy_aa, pars) where act is the actions, ang the initial angles and
          frequencies, n_vec the n vectors of the Fourier modes, toy_aa the toy action-angle
          coords, and pars are the toy potential parameters
          N_matrix sets the maximum |n| of the Fourier modes used,
          use_box forces the routine to use the triaxial harmonic oscillator as the toy potential,
          ifloop=True returns orbit classification,
          ifprint=True prints progress messages.