Based on MachineCurrents.m from MatDyn by Stijn Cole, developed at
        Katholieke Universiteit Leuven. See U{http://www.esat.kuleuven.be/
        electa/teaching/matdyn/} for more information.

        @param Xg: Generator state variables.
        @param U: Generator voltages.

        @rtype: tuple
        @return: Currents and electric power of generators.

def machineCurrents(self, Xg, U):
        """ Based on MachineCurrents.m from MatDyn by Stijn Cole, developed at
        Katholieke Universiteit Leuven. See U{http://www.esat.kuleuven.be/
        electa/teaching/matdyn/} for more information.

        @param Xg: Generator state variables.
        @param U: Generator voltages.

        @rtype: tuple
        @return: Currents and electric power of generators.
        """
        generators = self.dyn_generators

        # Initialise.
        ng = len(generators)
        Id = zeros(ng)
        Iq = zeros(ng)
        Pe = zeros(ng)

        typ1 = [g._i for g in generators if g.model == CLASSICAL]
        typ2 = [g._i for g in generators if g.model == FOURTH_ORDER]

        # Generator type 1: classical model
        delta = Xg[typ1, 0]
        Eq_tr = Xg[typ1, 2]

        xd = array([g.xd for g in generators])

        Pe[typ1] = \
            1 / xd * abs(U[typ1]) * abs(Eq_tr) * sin(delta - angle(U[typ1]))

        # Generator type 2: 4th order model
        delta = Xg[typ1, 0]
        Eq_tr = Xg[typ1, 2]
        Ed_tr = Xg[typ1, 3]

        xd_tr = array([g.xd_tr for g in generators])
        xq_tr = array([g.xq_tr for g in generators])

        theta = angle(U)

        # Transform U to rotor frame of reference.
        vd = -abs(U[typ2]) * sin(delta - theta[typ2])
        vq =  abs(U[typ2]) * cos(delta - theta[typ2])

        Id[typ2] =  (vq - Eq_tr) / xd_tr
        Iq[typ2] = -(vd - Ed_tr) / xq_tr

        Pe[typ2] = \
            Eq_tr * Iq[typ2] + Ed_tr * Id[typ2] + \
            (xd_tr - xq_tr) * Id[typ2] * Iq[typ2]

        return Id, Iq, Pe