code: def generator(self, Xgen, Xexc, Xgov, Vgen):
          generators = self.dyn_generators
          omegas = 2 * pi * self.freq
          F = zeros(Xgen.shape)
          typ1 = [g._i for g in generators if g.model == CLASSICAL]
          typ2 = [g._i for g in generators if g.model == FOURTH_ORDER]
          omega = Xgen[typ1, 1]
          Pm0 = Xgov[typ1, 0]
          H = array([g.h for g in generators])[typ1]
          D = array([g.d for g in generators])[typ1]
          Pe = Vgen[typ1, 2]
          ddelta = omega = omegas
          domega = pi * self.freq / H * (-D * (omega - omegas) + Pm0 - Pe)
          dEq = zeros(len(typ1))
          F[typ1, :] = c_[ddelta, domega, dEq]
          omega = Xgen[typ2, 1]
          Eq_tr = Xgen[typ2, 2]
          Ed_tr = Xgen[typ2, 3]
          H = array([g.h for g in generators])
          D = array([g.d for g in generators])
          xd = array([g.xd for g in generators])
          xq = array([g.xq for g in generators])
          xd_tr = array([g.xd_tr for g in generators])
          xq_tr = array([g.xq_tr for g in generators])
          Td0_tr = array([g.td for g in generators])
          Tq0_tr = array([g.tq for g in generators])
          Id = Vgen[typ2, 0]
          Iq = Vgen[typ2, 1]
          Pe = Vgen[typ2, 2]
          Efd = Xexc[typ2, 0]
          Pm = Xgov[typ2, 0]
          ddelta = omega - omegas
          domega = pi * self.freq / H * (-D * (omega - omegas) + Pm - Pe)
          dEq = 1 / Td0_tr * (Efd - Eq_tr + (xd - xd_tr) * Id)
          dEd = 1 / Tq0_tr * (-Ed_tr - (xq - xq_tr) * Iq)
          F[typ2, :] = c_[ddelta, domega, dEq, dEd]
          return F