<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def radiated_intensity(rho, i, j, epsilonp, rm, omega_level, xi, N, D, unfolding):
 r"""Return the radiated intensity in a given direction. 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