docs: def paschen_back_energies(fine_state, Bz):
      r"""Return Paschen-Back regime energies for a given fine state and\
      magnetic field.
  
      >>> ground_state = State("Rb", 87, 5, 0, 1/Integer(2))
      >>> Bz = 200.0
      >>> Bz = Bz/10000
      >>> for f_group in paschen_back_energies(ground_state, Bz):
      ...     print(f_group)
      [1.51284728917866e-24 3.80485568127324e-25 -7.51876152924007e-25
       -1.88423787397534e-24]
      [-1.51229355210131e-24 -3.80300989101543e-25 7.51691573898227e-25
       1.88368413689800e-24]
  
      """
      element = fine_state.element
      isotope = fine_state.isotope
  
      N = fine_state.n
      L = fine_state.l
      J = fine_state.j
      II = Atom(element, isotope).nuclear_spin
      MJ = [-J+i for i in range(2*J+1)]
      MI = [-II+i for i in range(2*II+1)]
  
      Ahfs = fine_state.Ahfs
      Bhfs = fine_state.Bhfs
      gL, gS, gI, gJ = lande_g_factors(element, isotope, L, J)
  
      energiesPBack = []
      for mj in MJ:
          energiesMJ = []
          unperturbed_energy = hbar*State(element, isotope, N, L, J).omega
          for mi in MI:
              energyMI = unperturbed_energy
              energyMI += 2*pi*hbar*Ahfs*mi*mj
              if J != 1/Integer(2) and II != 1/Integer(2):
                  num = 9*(mi*mj)**2 - 3*J*(J+1)*mi**2
                  num += -3*II*(II+1)*mj**2 + II*(II+1)*J*(J+1)
                  den = 4*J*(2*J-1)*II*(2*II-1)
                  energyMI += 2*pi*hbar*Bhfs*num/den
  
              energyMI += muB*(gJ*mj+gI*mi)*Bz
  
              energiesMJ += [energyMI]
          energiesPBack += [energiesMJ]
      return array(energiesPBack)