docs: def phase_transformation(Ne, Nl, rm, xi, return_equations=False):
      """Returns a phase transformation theta_i.
  
      The phase transformation is defined in a way such that
          theta1 + omega_level1 = 0.
  
      >>> xi = np.zeros((1, 2, 2))
      >>> xi[0, 1, 0] = 1.0
      >>> xi[0, 0, 1] = 1.0
      >>> rm = np.zeros((3, 2, 2))
      >>> rm[0, 1, 0] = 1.0
      >>> rm[1, 1, 0] = 1.0
      >>> rm[2, 1, 0] = 1.0
      >>> phase_transformation(2, 1, rm, xi)
      [-omega_1, -omega_1 - varpi_1]
  
      """
      # We first define the needed variables
      E0, omega_laser = define_laser_variables(Nl)
      theta = [Symbol('theta'+str(i+1)) for i in range(Ne)]
  
      # We check for the case of xi being a list of matrices.
      if type(xi) == list:
          xi = np.array([[[xi[l][i, j]
                         for j in range(Ne)] for i in range(Ne)]
                         for l in range(Nl)])
  
      # We find all the equations that the specified problem has to fulfil.
      eqs = []
      for i in range(Ne):
          for j in range(0, i):
              if (rm[0][i, j] != 0) or \
                 (rm[1][i, j] != 0) or \
                 (rm[2][i, j] != 0):
                  for l in range(Nl):
                      if xi[l, i, j] == 1:
                          eqs += [-omega_laser[l] + theta[j] - theta[i]]
  
      if return_equations:
          return eqs
  
      # We solve the system of equations.
      sol = sympy.solve(eqs, theta, dict=True)
      sol = sol[0]
      # We add any missing theta that may be left outside if the system is
      # under determined.
  
      extra_thetas = []
      for i in range(Ne):
          if theta[i] not in sol.keys():
              sol.update({theta[i]: theta[i]})
              extra_thetas += [theta[i]]
  
      # We make the solution such that theta1 + omega_level1 = 0.
      omega_level, omega, gamma = define_frequencies(Ne)
      eq_crit = sol[theta[0]] + omega_level[0]
      ss = sympy.solve(eq_crit, extra_thetas[0])[0]
      ss = {extra_thetas[0]: ss}
  
      sol_simple = [sol[theta[i]].subs(ss)
                    for i in range(Ne)]
  
      # sol = []
      # for i in range(Ne):
      #     soli = []
      #     for l in range(Nl):
      #         soli += [sympy.diff(sol_simple[theta[i]], omega_laser[l])]
      #     sol += [soli]
  
      return sol_simple