docs: def radiated_intensity(rho, i, j, epsilonp, rm, omega_level, xi,
                         N, D, unfolding):
      r"""Return the radiated intensity in a given direction.
  
      >>> from fast import State, Integer, split_hyperfine_to_magnetic
      >>> g = State("Rb", 87, 5, 1, 3/Integer(2), 0)
      >>> e = State("Rb", 87, 4, 2, 5/Integer(2), 1)
      >>> magnetic_states = split_hyperfine_to_magnetic([g, e])
      >>> omega0 = magnetic_states[0].omega
      >>> omega_level = [ei.omega - omega0 for ei in magnetic_states]
      >>> Ne = len(magnetic_states)
  
      >>> N = 4e6
      >>> D = 0.1
      >>> unfolding = Unfolding(Ne, True, True, True)
  
      >>> rho = np.zeros((Ne, Ne))
      >>> rho[0, 0] = 0.8
      >>> rho[3, 3] = 0.2
      >>> rho[3, 0] = 0.3
      >>> rho[0, 3] = 0.3
      >>> rho = unfolding(rho)
  
      >>> ep = np.array([1, 1j, 0])/np.sqrt(2.0)
      >>> ex = np.array([1, 0, 0])
      >>> r0 = 4.75278521538619e-11
      >>> rm = np.zeros((3, Ne, Ne), complex)
      >>> rm[0, 1, 0] = -r0
      >>> rm[0, 3, 0] = r0
      >>> rm[1, 1, 0] = -1j*r0
      >>> rm[1, 3, 0] = -1j*r0
      >>> rm[1, 2, 0] = -np.sqrt(2)*r0
  
      >>> xi = np.zeros((1, Ne, Ne))
      >>> xi[0, 1, 0] = 1
      >>> xi[0, 2, 0] = 1
      >>> xi[0, 3, 0] = 1
      >>> xi[0, :, :] += xi[0, :, :].transpose()
  
      >>> print(radiated_intensity(rho, 1, 0, ex, rm,
      ...                          omega_level, xi, N, D, unfolding))
      4.60125990174e-22
  
      """
      def inij(i, j, ilist, jlist):
          if (i in ilist) and (j in jlist):
              return 1
          else:
              return 0
  
      rm = np.array(rm)
      Nl = xi.shape[0]
      Ne = xi.shape[1]
      aux = define_simplification(omega_level, xi, Nl)
      u = aux[0]
      omega_levelu = aux[2]
  
      ui = u(i)
      uj = u(j)
      omegaij = omega_levelu[ui] - omega_levelu[uj]
  
      ilist = [ii for ii in range(Ne) if u(ii) == ui]
      jlist = [jj for jj in range(Ne) if u(jj) == uj]
  
      rp = np.array([rm[ii].conjugate().transpose() for ii in range(3)])
  
      rm = np.array([[[rm[p, ii, jj]*inij(ii, jj, ilist, jlist)
                     for jj in range(Ne)]
                     for ii in range(Ne)]
                     for p in range(3)])
      rp = np.array([[[rp[p, ii, jj]*inij(jj, ii, ilist, jlist)
                     for jj in range(Ne)]
                     for ii in range(Ne)]
                     for p in range(3)])
  
      epsilonm = epsilonp.conjugate()
  
      Adag = cartesian_dot_product(rm, epsilonp)
      A = cartesian_dot_product(rp, epsilonm)
  
      fact = alpha_num*N*hbar_num*omegaij**3/2/np.pi/c_num**2/D**2
  
      Iop = fact * np.dot(Adag, A)
      intensity = observable(Iop, rho, unfolding)
      intensity = float(np.real(intensity))
      return intensity