diff --git a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp
index 830033b59..e77ebb5c7 100644
--- a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp
+++ b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp
@@ -43,17 +43,17 @@ namespace GridKit
   {
     assert(refframe_ || node_ref_->size() == 1);
     assert(node_bus_->size() == 2);
-    // internals [\delta_i, Pi, Qi, phi_di, phi_qi, gamma_di, gamma_qi, il_di, il_qi, vo_di, vo_qi, io_di, io_qi]
+    // internals [Pi, Qi, phi_di, phi_qi, gamma_di, gamma_qi, il_di, il_qi, vo_di, vo_qi, io_di, io_qi, \delta_i]
     // externals [\omega_ref, vba_out, vbb_out]
     size_     = 16;
     n_intern_ = refframe_ ? 12 : 13;
     n_extern_ = refframe_ ? 4 : 3;
     idc_      = id;
-    nnz_      = refframe_ ? 80 : 78;
+    nnz_      = refframe_ ? 77 : 78;
 
     if (refframe_)
     {
-      extern_indices_ = {0, 1, 2, 3};
+      extern_indices_ = {0, 1, 2, 15};
     }
     else
     {
@@ -91,54 +91,56 @@ namespace GridKit
   template <class ScalarT, typename IdxT>
   int DistributedGenerator<ScalarT, IdxT>::evaluateInternalResidual()
   {
-    ScalarT omega = wb_ - mp_ * y_[4];
+    ScalarT omega = wb_ - mp_ * y_[3];
 
     // Take incoming voltages to current rotator reference frame
-    ScalarT vbd_in = std::cos(y_[3]) * y_[1] + std::sin(y_[3]) * y_[2];
-    ScalarT vbq_in = -std::sin(y_[3]) * y_[1] + std::cos(y_[3]) * y_[2];
+    ScalarT vbd_in = std::cos(y_[15]) * y_[1] + std::sin(y_[15]) * y_[2];
+    ScalarT vbq_in = -std::sin(y_[15]) * y_[1] + std::cos(y_[15]) * y_[2];
 
     // ### Internal Componenets ##
-    // Rotor difference angle
-    f_[3] = -yp_[3] + omega - y_[0];
-
     // P and Q equations
-    f_[4] = -yp_[4] + wc_ * (y_[12] * y_[14] + y_[13] * y_[15] - y_[4]);
-    f_[5] = -yp_[5] + wc_ * (-y_[12] * y_[15] + y_[13] * y_[14] - y_[5]);
+    f_[3] = -yp_[3] + wc_ * (y_[11] * y_[13] + y_[12] * y_[14] - y_[3]);
+    f_[4] = -yp_[4] + wc_ * (-y_[11] * y_[14] + y_[12] * y_[13] - y_[4]);
 
     // Voltage control
-    ScalarT vod_star = Vn_ - nq_ * y_[5];
+    ScalarT vod_star = Vn_ - nq_ * y_[4];
     ScalarT voq_star = static_cast<ScalarT>(0.0);
 
-    f_[6] = -yp_[6] + vod_star - y_[12];
-    f_[7] = -yp_[7] + voq_star - y_[13];
+    f_[5] = -yp_[5] + vod_star - y_[11];
+    f_[6] = -yp_[6] + voq_star - y_[12];
 
-    ScalarT ild_star = F_ * y_[14] - wb_ * Cf_ * y_[13] + Kpv_ * (vod_star - y_[12]) + Kiv_ * y_[6];
-    ScalarT ilq_star = F_ * y_[15] + wb_ * Cf_ * y_[12] + Kpv_ * (voq_star - y_[13]) + Kiv_ * y_[7];
+    ScalarT ild_star = F_ * y_[13] - wb_ * Cf_ * y_[12] + Kpv_ * (vod_star - y_[11]) + Kiv_ * y_[5];
+    ScalarT ilq_star = F_ * y_[14] + wb_ * Cf_ * y_[11] + Kpv_ * (voq_star - y_[12]) + Kiv_ * y_[6];
 
     // Current control
-    f_[8] = -yp_[8] + ild_star - y_[10];
-    f_[9] = -yp_[9] + ilq_star - y_[11];
+    f_[7] = -yp_[7] + ild_star - y_[9];
+    f_[8] = -yp_[8] + ilq_star - y_[10];
 
-    ScalarT vid_star = -wb_ * Lf_ * y_[11] + Kpc_ * (ild_star - y_[10]) + Kic_ * y_[8];
-    ScalarT viq_star = wb_ * Lf_ * y_[10] + Kpc_ * (ilq_star - y_[11]) + Kic_ * y_[9];
+    ScalarT vid_star = -wb_ * Lf_ * y_[10] + Kpc_ * (ild_star - y_[9]) + Kic_ * y_[7];
+    ScalarT viq_star = wb_ * Lf_ * y_[9] + Kpc_ * (ilq_star - y_[10]) + Kic_ * y_[8];
 
     // Output LC Filter
-    f_[10] = -yp_[10] - (rLf_ / Lf_) * y_[10] + omega * y_[11] + (vid_star - y_[12]) / Lf_;
-    f_[11] = -yp_[11] - (rLf_ / Lf_) * y_[11] - omega * y_[10] + (viq_star - y_[13]) / Lf_;
+    f_[9]  = -yp_[9] - (rLf_ / Lf_) * y_[9] + omega * y_[10] + (vid_star - y_[11]) / Lf_;
+    f_[10] = -yp_[10] - (rLf_ / Lf_) * y_[10] - omega * y_[9] + (viq_star - y_[12]) / Lf_;
 
-    f_[12] = -yp_[12] + omega * y_[13] + (y_[10] - y_[14]) / Cf_;
-    f_[13] = -yp_[13] - omega * y_[12] + (y_[11] - y_[15]) / Cf_;
+    f_[11] = -yp_[11] + omega * y_[12] + (y_[9] - y_[13]) / Cf_;
+    f_[12] = -yp_[12] - omega * y_[11] + (y_[10] - y_[14]) / Cf_;
 
     // Output Connector
-    f_[14] = -yp_[14] - (rLc_ / Lc_) * y_[14] + omega * y_[15] + (y_[12] - vbd_in) / Lc_;
-    f_[15] = -yp_[15] - (rLc_ / Lc_) * y_[15] - omega * y_[14] + (y_[13] - vbq_in) / Lc_;
+    f_[13] = -yp_[13] - (rLc_ / Lc_) * y_[13] + omega * y_[14] + (y_[11] - vbd_in) / Lc_;
+    f_[14] = -yp_[14] - (rLc_ / Lc_) * y_[14] - omega * y_[13] + (y_[12] - vbq_in) / Lc_;
+
+    // Rotor difference angle
+    if (!refframe_)
+      f_[15] = -yp_[15] + omega - y_[0];
+
     return 0;
   }
 
   template <class ScalarT, typename IdxT>
   int DistributedGenerator<ScalarT, IdxT>::evaluateExternalResidual()
   {
-    ScalarT omega = wb_ - mp_ * y_[4];
+    ScalarT omega = wb_ - mp_ * y_[3];
     // ref common ref motor angle
     if (refframe_)
     {
@@ -151,8 +153,8 @@ namespace GridKit
 
     // output
     // current transformed to common frame
-    f_[1] = std::cos(y_[3]) * y_[14] - std::sin(y_[3]) * y_[15];
-    f_[2] = std::sin(y_[3]) * y_[14] + std::cos(y_[3]) * y_[15];
+    f_[1] = std::cos(y_[15]) * y_[13] - std::sin(y_[15]) * y_[14];
+    f_[2] = std::sin(y_[15]) * y_[13] + std::cos(y_[15]) * y_[14];
     return 0;
   }
 
@@ -195,39 +197,34 @@ namespace GridKit
     // Create dF/dy
     // r = 1
 
-    ctemp = {3, 14, 15};
+    ctemp = {13, 14, 15};
     rtemp.clear();
     for (size_t i = 0; i < ctemp.size(); i++)
       rtemp.push_back(1);
-    valtemp = {-sin(static_cast<RealT>(y_[3])) * static_cast<RealT>(y_[14]) - cos(static_cast<RealT>(y_[3])) * static_cast<RealT>(y_[15]),
-               cos(static_cast<RealT>(y_[3])),
-               -sin(static_cast<RealT>(y_[3]))};
+    valtemp = {
+        cos(static_cast<RealT>(y_[15])),
+        -sin(static_cast<RealT>(y_[15])),
+        -sin(static_cast<RealT>(y_[15])) * static_cast<RealT>(y_[13]) - cos(static_cast<RealT>(y_[15])) * static_cast<RealT>(y_[14]),
+    };
     this->setJacValues(rtemp, ctemp, valtemp);
 
     // r = 2
 
-    ctemp = {3, 14, 15};
+    ctemp = {13, 14, 15};
     rtemp.clear();
     for (size_t i = 0; i < ctemp.size(); i++)
       rtemp.push_back(2);
-    valtemp = {cos(static_cast<RealT>(y_[3])) * static_cast<RealT>(y_[14]) - sin(static_cast<RealT>(y_[3])) * static_cast<RealT>(y_[15]),
-               sin(static_cast<RealT>(y_[3])),
-               cos(static_cast<RealT>(y_[3]))};
-    this->setJacValues(rtemp, ctemp, valtemp);
-
-    // r = 3
-
-    ctemp = {0, 3, 4};
-    rtemp.clear();
-    for (size_t i = 0; i < ctemp.size(); i++)
-      rtemp.push_back(3);
-    valtemp = {-1.0, -alpha_, -mp_};
+    valtemp = {
+        sin(static_cast<RealT>(y_[15])),
+        cos(static_cast<RealT>(y_[15])),
+        cos(static_cast<RealT>(y_[15])) * static_cast<RealT>(y_[13]) - sin(static_cast<RealT>(y_[15])) * static_cast<RealT>(y_[14]),
+    };
     this->setJacValues(rtemp, ctemp, valtemp);
 
     // r = 0
     if (refframe_)
     {
-      ctemp = {0, 4};
+      ctemp = {0, 3};
       rtemp.clear();
       for (size_t i = 0; i < ctemp.size(); i++)
         rtemp.push_back(0);
@@ -235,140 +232,155 @@ namespace GridKit
       this->setJacValues(rtemp, ctemp, valtemp);
     }
 
+    // r = 3
+    ctemp = {3, 11, 12, 13, 14};
+    rtemp.clear();
+    for (size_t i = 0; i < ctemp.size(); i++)
+      rtemp.push_back(3);
+    valtemp = {-wc_ - alpha_,
+               wc_ * static_cast<RealT>(y_[13]),
+               wc_ * static_cast<RealT>(y_[14]),
+               wc_ * static_cast<RealT>(y_[11]),
+               wc_ * static_cast<RealT>(y_[12])};
+    this->setJacValues(rtemp, ctemp, valtemp);
+
     // r = 4
-    ctemp = {4, 12, 13, 14, 15};
+    ctemp = {4, 11, 12, 13, 14};
     rtemp.clear();
     for (size_t i = 0; i < ctemp.size(); i++)
       rtemp.push_back(4);
     valtemp = {-wc_ - alpha_,
-               wc_ * static_cast<RealT>(y_[14]),
-               wc_ * static_cast<RealT>(y_[15]),
+               -wc_ * static_cast<RealT>(y_[14]),
+               wc_ * static_cast<RealT>(y_[13]),
                wc_ * static_cast<RealT>(y_[12]),
-               wc_ * static_cast<RealT>(y_[13])};
+               -wc_ * static_cast<RealT>(y_[11])};
     this->setJacValues(rtemp, ctemp, valtemp);
 
     // r = 5
-    ctemp = {5, 12, 13, 14, 15};
+    ctemp = {4, 5, 11};
     rtemp.clear();
     for (size_t i = 0; i < ctemp.size(); i++)
       rtemp.push_back(5);
-    valtemp = {-wc_ - alpha_,
-               -wc_ * static_cast<RealT>(y_[15]),
-               wc_ * static_cast<RealT>(y_[14]),
-               wc_ * static_cast<RealT>(y_[13]),
-               -wc_ * static_cast<RealT>(y_[12])};
+    valtemp = {-nq_, -alpha_, -1.0};
     this->setJacValues(rtemp, ctemp, valtemp);
 
     // r = 6
-    ctemp = {5, 6, 12};
+    ctemp = {6, 12};
     rtemp.clear();
     for (size_t i = 0; i < ctemp.size(); i++)
       rtemp.push_back(6);
-    valtemp = {-nq_, -alpha_, -1.0};
+    valtemp = {-alpha_, -1.0};
     this->setJacValues(rtemp, ctemp, valtemp);
 
     // r = 7
-    ctemp = {7, 13};
+    ctemp = {4, 5, 7, 9, 11, 12, 13};
     rtemp.clear();
     for (size_t i = 0; i < ctemp.size(); i++)
       rtemp.push_back(7);
-    valtemp = {-alpha_, -1.0};
+    valtemp = {-Kpv_ * nq_, Kiv_, -alpha_, -1.0, -Kpv_, -Cf_ * wb_, F_};
     this->setJacValues(rtemp, ctemp, valtemp);
 
     // r = 8
-    ctemp = {5, 6, 8, 10, 12, 13, 14};
+    ctemp = {6, 8, 10, 11, 12, 14};
     rtemp.clear();
     for (size_t i = 0; i < ctemp.size(); i++)
       rtemp.push_back(8);
-    valtemp = {-Kpv_ * nq_, Kiv_, -alpha_, -1.0, -Kpv_, -Cf_ * wb_, F_};
+    valtemp = {Kiv_, -alpha_, -1.0, Cf_ * wb_, -Kpv_, F_};
     this->setJacValues(rtemp, ctemp, valtemp);
 
     // r = 9
-    ctemp = {7, 9, 11, 12, 13, 15};
+    ctemp = {3, 4, 5, 7, 9, 10, 11, 12, 13};
     rtemp.clear();
     for (size_t i = 0; i < ctemp.size(); i++)
       rtemp.push_back(9);
-    valtemp = {Kiv_, -alpha_, -1.0, Cf_ * wb_, -Kpv_, F_};
-    this->setJacValues(rtemp, ctemp, valtemp);
-
-    // r = 10
-    ctemp = {4, 5, 6, 8, 10, 11, 12, 13, 14};
-    rtemp.clear();
-    for (size_t i = 0; i < ctemp.size(); i++)
-      rtemp.push_back(10);
-    valtemp = {-mp_ * static_cast<RealT>(y_[11]),
+    valtemp = {-mp_ * static_cast<RealT>(y_[10]),
                -(Kpc_ * Kpv_ * nq_) / Lf_,
                (Kpc_ * Kiv_) / Lf_,
                Kic_ / Lf_,
                -(Kpc_ + rLf_) / Lf_ - alpha_,
-               -mp_ * static_cast<RealT>(y_[4]),
+               -mp_ * static_cast<RealT>(y_[3]),
                -(Kpc_ * Kpv_ + 1.0) / Lf_,
                -(Cf_ * Kpc_ * wb_) / Lf_,
                (F_ * Kpc_) / Lf_};
     this->setJacValues(rtemp, ctemp, valtemp);
 
-    // r = 11
-    ctemp = {4, 7, 9, 10, 11, 12, 13, 15};
+    // r = 10
+    ctemp = {3, 6, 8, 9, 10, 11, 12, 14};
     rtemp.clear();
     for (size_t i = 0; i < ctemp.size(); i++)
-      rtemp.push_back(11);
-    valtemp = {mp_ * static_cast<RealT>(y_[10]),
+      rtemp.push_back(10);
+    valtemp = {mp_ * static_cast<RealT>(y_[9]),
                (Kiv_ * Kpc_) / Lf_,
                Kic_ / Lf_,
-               mp_ * static_cast<RealT>(y_[4]),
+               mp_ * static_cast<RealT>(y_[3]),
                -(Kpc_ + rLf_) / Lf_ - alpha_,
                (Cf_ * Kpc_ * wb_) / Lf_,
                -(Kpc_ * Kpv_ + 1.0) / Lf_,
                (F_ * Kpc_) / Lf_};
     this->setJacValues(rtemp, ctemp, valtemp);
 
-    // r = 12
-    ctemp = {4, 10, 12, 13, 14};
+    // r = 11
+    ctemp = {3, 9, 11, 12, 13};
     rtemp.clear();
     for (size_t i = 0; i < ctemp.size(); i++)
-      rtemp.push_back(12);
-    valtemp = {-mp_ * static_cast<RealT>(y_[13]),
+      rtemp.push_back(11);
+    valtemp = {-mp_ * static_cast<RealT>(y_[12]),
                1.0 / Cf_,
                -alpha_,
-               wb_ - mp_ * static_cast<RealT>(y_[4]),
+               wb_ - mp_ * static_cast<RealT>(y_[3]),
                -1.0 / Cf_};
     this->setJacValues(rtemp, ctemp, valtemp);
 
+    // r = 12
+    ctemp = {3, 10, 11, 12, 14};
+    rtemp.clear();
+    for (size_t i = 0; i < ctemp.size(); i++)
+      rtemp.push_back(12);
+    valtemp = {mp_ * static_cast<RealT>(y_[11]), 1.0 / Cf_, -wb_ + mp_ * static_cast<RealT>(y_[3]), -alpha_, -1.0 / Cf_};
+    this->setJacValues(rtemp, ctemp, valtemp);
+
     // r = 13
-    ctemp = {4, 11, 12, 13, 15};
+    ctemp = {1, 2, 3, 11, 13, 14, 15};
     rtemp.clear();
     for (size_t i = 0; i < ctemp.size(); i++)
       rtemp.push_back(13);
-    valtemp = {mp_ * static_cast<RealT>(y_[12]), 1.0 / Cf_, -wb_ + mp_ * static_cast<RealT>(y_[4]), -alpha_, -1.0 / Cf_};
+    valtemp = {
+        (1.0 / Lc_) * -cos(static_cast<RealT>(y_[15])),
+        (1.0 / Lc_) * -sin(static_cast<RealT>(y_[15])),
+        -mp_ * static_cast<RealT>(y_[14]),
+        1.0 / Lc_,
+        -rLc_ / Lc_ - alpha_,
+        wb_ - mp_ * static_cast<RealT>(y_[3]),
+        (1.0 / Lc_) * (sin(static_cast<RealT>(y_[15])) * static_cast<RealT>(y_[1]) - cos(static_cast<RealT>(y_[15])) * static_cast<RealT>(y_[2])),
+    };
     this->setJacValues(rtemp, ctemp, valtemp);
 
     // r = 14
-    ctemp = {1, 2, 3, 4, 12, 14, 15};
+    ctemp = {1, 2, 3, 12, 13, 14, 15};
     rtemp.clear();
     for (size_t i = 0; i < ctemp.size(); i++)
       rtemp.push_back(14);
-    valtemp = {(1.0 / Lc_) * -cos(static_cast<RealT>(y_[3])),
-               (1.0 / Lc_) * -sin(static_cast<RealT>(y_[3])),
-               (1.0 / Lc_) * (sin(static_cast<RealT>(y_[3])) * static_cast<RealT>(y_[1]) - cos(static_cast<RealT>(y_[3])) * static_cast<RealT>(y_[2])),
-               -mp_ * static_cast<RealT>(y_[15]),
-               1.0 / Lc_,
-               -rLc_ / Lc_ - alpha_,
-               wb_ - mp_ * static_cast<RealT>(y_[4])};
+    valtemp = {
+        (1.0 / Lc_) * sin(static_cast<RealT>(y_[15])),
+        (1.0 / Lc_) * -cos(static_cast<RealT>(y_[15])),
+        mp_ * static_cast<RealT>(y_[13]),
+        1.0 / Lc_,
+        -wb_ + mp_ * static_cast<RealT>(y_[3]),
+        -rLc_ / Lc_ - alpha_,
+        (1.0 / Lc_) * (cos(static_cast<RealT>(y_[15])) * static_cast<RealT>(y_[1]) + sin(static_cast<RealT>(y_[15])) * static_cast<RealT>(y_[2])),
+    };
     this->setJacValues(rtemp, ctemp, valtemp);
 
-    // r = 15
-    ctemp = {1, 2, 3, 4, 13, 14, 15};
-    rtemp.clear();
-    for (size_t i = 0; i < ctemp.size(); i++)
-      rtemp.push_back(15);
-    valtemp = {(1.0 / Lc_) * sin(static_cast<RealT>(y_[3])),
-               (1.0 / Lc_) * -cos(static_cast<RealT>(y_[3])),
-               (1.0 / Lc_) * (cos(static_cast<RealT>(y_[3])) * static_cast<RealT>(y_[1]) + sin(static_cast<RealT>(y_[3])) * static_cast<RealT>(y_[2])),
-               mp_ * static_cast<RealT>(y_[14]),
-               1.0 / Lc_,
-               -wb_ + mp_ * static_cast<RealT>(y_[4]),
-               -rLc_ / Lc_ - alpha_};
-    this->setJacValues(rtemp, ctemp, valtemp);
+    if (!refframe_)
+    {
+      // r = 15
+      ctemp = {0, 3, 15};
+      rtemp.clear();
+      for (size_t i = 0; i < ctemp.size(); i++)
+        rtemp.push_back(15);
+      valtemp = {-1.0, -mp_, -alpha_};
+      this->setJacValues(rtemp, ctemp, valtemp);
+    }
 
     return 0;
   }
@@ -384,7 +396,7 @@ namespace GridKit
 
     if (refframe_)
     {
-      this->setExternalConnectionNodes(3, INVALID_INDEX<IdxT>);
+      this->setExternalConnectionNodes(15, INVALID_INDEX<IdxT>);
     }
 
     return 0;
diff --git a/examples/Enzyme/PowerElectronics/main.cpp b/examples/Enzyme/PowerElectronics/main.cpp
index 22f447c0f..bb5f2d245 100644
--- a/examples/Enzyme/PowerElectronics/main.cpp
+++ b/examples/Enzyme/PowerElectronics/main.cpp
@@ -48,7 +48,13 @@ void evaluateResidual(std::vector<double> y_, std::vector<double> yp_, std::vect
 
   constexpr bool ref_frame_ = true;
 
-  double omega = wb_ - mp_ * y_[4];
+  using ScalarT = double;
+
+  ScalarT omega = wb_ - mp_ * y_[3];
+
+  ScalarT vbd_in = std::cos(y_[15]) * y_[1] + std::sin(y_[15]) * y_[2];
+  ScalarT vbq_in = -std::sin(y_[15]) * y_[1] + std::cos(y_[15]) * y_[2];
+
   if (ref_frame_)
   {
     f_[0] = omega - y_[0];
@@ -58,39 +64,47 @@ void evaluateResidual(std::vector<double> y_, std::vector<double> yp_, std::vect
     f_[0] = 0.0;
   }
 
-  f_[1] = cos(y_[3]) * y_[14] - sin(y_[3]) * y_[15];
-  f_[2] = sin(y_[3]) * y_[14] + cos(y_[3]) * y_[15];
+  // output
+  // current transformed to common frame
+  f_[1] = std::cos(y_[15]) * y_[13] - std::sin(y_[15]) * y_[14];
+  f_[2] = std::sin(y_[15]) * y_[13] + std::cos(y_[15]) * y_[14];
 
-  double vbd_in = cos(y_[3]) * y_[1] + sin(y_[3]) * y_[2];
-  double vbq_in = -sin(y_[3]) * y_[1] + cos(y_[3]) * y_[2];
+  // ### Internal Componenets ##
+  // P and Q equations
+  f_[3] = -yp_[3] + wc_ * (y_[11] * y_[13] + y_[12] * y_[14] - y_[3]);
+  f_[4] = -yp_[4] + wc_ * (-y_[11] * y_[14] + y_[12] * y_[13] - y_[4]);
 
-  f_[3] = -yp_[3] + omega - y_[0];
-  f_[4] = -yp_[4] + wc_ * (y_[12] * y_[14] + y_[13] * y_[15] - y_[4]);
-  f_[5] = -yp_[5] + wc_ * (-y_[12] * y_[15] + y_[13] * y_[14] - y_[5]);
+  // Voltage control
+  ScalarT vod_star = Vn_ - nq_ * y_[4];
+  ScalarT voq_star = static_cast<ScalarT>(0.0);
 
-  double vod_star = Vn_ - nq_ * y_[5];
-  double voq_star = 0.0;
+  f_[5] = -yp_[5] + vod_star - y_[11];
+  f_[6] = -yp_[6] + voq_star - y_[12];
 
-  f_[6] = -yp_[6] + vod_star - y_[12];
-  f_[7] = -yp_[7] + voq_star - y_[13];
+  ScalarT ild_star = F_ * y_[13] - wb_ * Cf_ * y_[12] + Kpv_ * (vod_star - y_[11]) + Kiv_ * y_[5];
+  ScalarT ilq_star = F_ * y_[14] + wb_ * Cf_ * y_[11] + Kpv_ * (voq_star - y_[12]) + Kiv_ * y_[6];
 
-  double ild_star = F_ * y_[14] - wb_ * Cf_ * y_[13] + Kpv_ * (vod_star - y_[12]) + Kiv_ * y_[6];
-  double ilq_star = F_ * y_[15] + wb_ * Cf_ * y_[12] + Kpv_ * (voq_star - y_[13]) + Kiv_ * y_[7];
+  // Current control
+  f_[7] = -yp_[7] + ild_star - y_[9];
+  f_[8] = -yp_[8] + ilq_star - y_[10];
 
-  f_[8] = -yp_[8] + ild_star - y_[10];
-  f_[9] = -yp_[9] + ilq_star - y_[11];
+  ScalarT vid_star = -wb_ * Lf_ * y_[10] + Kpc_ * (ild_star - y_[9]) + Kic_ * y_[7];
+  ScalarT viq_star = wb_ * Lf_ * y_[9] + Kpc_ * (ilq_star - y_[10]) + Kic_ * y_[8];
 
-  double vid_star = -wb_ * Lf_ * y_[11] + Kpc_ * (ild_star - y_[10]) + Kic_ * y_[8];
-  double viq_star = wb_ * Lf_ * y_[10] + Kpc_ * (ilq_star - y_[11]) + Kic_ * y_[9];
+  // Output LC Filter
+  f_[9]  = -yp_[9] - (rLf_ / Lf_) * y_[9] + omega * y_[10] + (vid_star - y_[11]) / Lf_;
+  f_[10] = -yp_[10] - (rLf_ / Lf_) * y_[10] - omega * y_[9] + (viq_star - y_[12]) / Lf_;
 
-  f_[10] = -yp_[10] - (rLf_ / Lf_) * y_[10] + omega * y_[11] + (vid_star - y_[12]) / Lf_;
-  f_[11] = -yp_[11] - (rLf_ / Lf_) * y_[11] - omega * y_[10] + (viq_star - y_[13]) / Lf_;
+  f_[11] = -yp_[11] + omega * y_[12] + (y_[9] - y_[13]) / Cf_;
+  f_[12] = -yp_[12] - omega * y_[11] + (y_[10] - y_[14]) / Cf_;
 
-  f_[12] = -yp_[12] + omega * y_[13] + (y_[10] - y_[14]) / Cf_;
-  f_[13] = -yp_[13] - omega * y_[12] + (y_[11] - y_[15]) / Cf_;
+  // Output Connector
+  f_[13] = -yp_[13] - (rLc_ / Lc_) * y_[13] + omega * y_[14] + (y_[11] - vbd_in) / Lc_;
+  f_[14] = -yp_[14] - (rLc_ / Lc_) * y_[14] - omega * y_[13] + (y_[12] - vbq_in) / Lc_;
 
-  f_[14] = -yp_[14] - (rLc_ / Lc_) * y_[14] + omega * y_[15] + (y_[12] - vbd_in) / Lc_;
-  f_[15] = -yp_[15] - (rLc_ / Lc_) * y_[15] - omega * y_[14] + (y_[13] - vbq_in) / Lc_;
+  // Rotor difference angle
+  if (!ref_frame_)
+    f_[15] = -yp_[15] + omega - y_[0];
 }
 
 // Function that computes the Jacobian via automatic differentiation
diff --git a/examples/PowerElectronics/DistributedGeneratorTest/DGTest.cpp b/examples/PowerElectronics/DistributedGeneratorTest/DGTest.cpp
index b2d067985..fb927c0be 100644
--- a/examples/PowerElectronics/DistributedGeneratorTest/DGTest.cpp
+++ b/examples/PowerElectronics/DistributedGeneratorTest/DGTest.cpp
@@ -49,7 +49,24 @@ int main(int /* argc */, char const** /* argv */)
   GridKit::DistributedGenerator<double, size_t> dg(0, parms, true, &dg_signal, &bus);
 
   std::vector<double> t1(16, 0.0);
-  std::vector<double> t2{0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5};
+  std::vector<double> t2{
+      0.0,
+      0.1,
+      0.2,
+      0.4,
+      0.5,
+      0.6,
+      0.7,
+      0.8,
+      0.9,
+      1.0,
+      1.1,
+      1.2,
+      1.3,
+      1.4,
+      1.5,
+      0.3,
+  };
 
   dg_signal.allocate();
   bus.allocate();
@@ -64,7 +81,6 @@ int main(int /* argc */, char const** /* argv */)
   std::vector<double> true_vec{3.141592277589793e+02,
                                8.941907747838389e-01,
                                1.846733023014284e+00,
-                               3.141592277589793e+02,
                                1.014543000000000e+02,
                                -1.507680000000000e+01,
                                3.787993500000000e+02,
diff --git a/examples/PowerElectronics/Microgrid/Microgrid.cpp b/examples/PowerElectronics/Microgrid/Microgrid.cpp
index c6c0690be..3b20c14e0 100644
--- a/examples/PowerElectronics/Microgrid/Microgrid.cpp
+++ b/examples/PowerElectronics/Microgrid/Microgrid.cpp
@@ -177,14 +177,14 @@ int main(int /* argc */, char const** /* argv */)
   sysmodel->yp()[2]      = parms1.Vn_;
   sysmodel->yp()[4]      = parms1.Kpv_ * parms1.Vn_;
   sysmodel->yp()[6]      = (parms1.Kpc_ * parms1.Kpv_ * parms1.Vn_) / parms1.Lf_;
-  sysmodel->yp()[12 + 3] = parms1.Vn_;
-  sysmodel->yp()[12 + 5] = parms1.Kpv_ * parms1.Vn_;
-  sysmodel->yp()[12 + 7] = (parms1.Kpc_ * parms1.Kpv_ * parms1.Vn_) / parms1.Lf_;
+  sysmodel->yp()[12 + 2] = parms1.Vn_;
+  sysmodel->yp()[12 + 4] = parms1.Kpv_ * parms1.Vn_;
+  sysmodel->yp()[12 + 6] = (parms1.Kpc_ * parms1.Kpv_ * parms1.Vn_) / parms1.Lf_;
   for (size_t i = 2; i < 4; i++)
   {
-    sysmodel->yp()[13 * i - 1 + 3] = parms2.Vn_;
-    sysmodel->yp()[13 * i - 1 + 5] = parms2.Kpv_ * parms2.Vn_;
-    sysmodel->yp()[13 * i - 1 + 7] = (parms2.Kpc_ * parms2.Kpv_ * parms2.Vn_) / parms2.Lf_;
+    sysmodel->yp()[13 * i - 1 + 2] = parms2.Vn_;
+    sysmodel->yp()[13 * i - 1 + 4] = parms2.Kpv_ * parms2.Vn_;
+    sysmodel->yp()[13 * i - 1 + 6] = (parms2.Kpc_ * parms2.Kpv_ * parms2.Vn_) / parms2.Lf_;
   }
 
   // since the intial P_com = 0
@@ -254,7 +254,6 @@ int main(int /* argc */, char const** /* argv */)
       -2.668928293656362e-06,
       6.321941919221522e+01,
       -3.509200178595996e+01,
-      -7.555954467454730e-03,
       2.297580486511343e+04,
       8.742028429066131e+03,
       3.710079564796484e-02,
@@ -267,7 +266,7 @@ int main(int /* argc */, char const** /* argv */)
       3.465673854181523e-05,
       6.232933406188410e+01,
       -2.371564475187742e+01,
-      -8.273939686941580e-02,
+      -7.555954467454730e-03,
       1.727775042678524e+04,
       1.649365247247288e+04,
       3.116555157570849e-02,
@@ -280,7 +279,7 @@ int main(int /* argc */, char const** /* argv */)
       -1.496407194199739e-04,
       4.861823504694532e+01,
       -4.642797132602495e+01,
-      -8.445727984408551e-02,
+      -8.273939686941580e-02,
       1.727723725566433e+04,
       9.182386962936238e+03,
       3.024959333190777e-02,
@@ -293,6 +292,7 @@ int main(int /* argc */, char const** /* argv */)
       1.076423957830039e-04,
       4.718938116520511e+01,
       -2.507094256286497e+01,
+      -8.445727984408551e-02,
       -1.881248349415025e+01,
       2.114714832305742e+01,
       4.329946674909793e+01,
diff --git a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp
index 12e685c3d..07246e83a 100644
--- a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp
+++ b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp
@@ -251,9 +251,9 @@ int test(index_type Nsize, real_type error_tol, bool debug_output)
   // Create initial derivatives specifics generated in MATLAB
   for (index_type i = 0; i < 2 * Nsize; i++)
   {
-    sys_model->yp()[13 * i - 1 + 3] = DGParams_list[i].Vn_;
-    sys_model->yp()[13 * i - 1 + 5] = DGParams_list[i].Kpv_ * DGParams_list[i].Vn_;
-    sys_model->yp()[13 * i - 1 + 7] = (DGParams_list[i].Kpc_ * DGParams_list[i].Kpv_ * DGParams_list[i].Vn_) / DGParams_list[i].Lf_;
+    sys_model->yp()[13 * i - 1 + 2] = DGParams_list[i].Vn_;
+    sys_model->yp()[13 * i - 1 + 4] = DGParams_list[i].Kpv_ * DGParams_list[i].Vn_;
+    sys_model->yp()[13 * i - 1 + 6] = (DGParams_list[i].Kpc_ * DGParams_list[i].Kpv_ * DGParams_list[i].Vn_) / DGParams_list[i].Lf_;
   }
 
   // since the intial P_com = 0, the set the intial vector to the reference frame
diff --git a/examples/PowerElectronics/ScaleMicrogrid/SolutionKeys.hpp b/examples/PowerElectronics/ScaleMicrogrid/SolutionKeys.hpp
index f604b5e9e..bc3b94bac 100644
--- a/examples/PowerElectronics/ScaleMicrogrid/SolutionKeys.hpp
+++ b/examples/PowerElectronics/ScaleMicrogrid/SolutionKeys.hpp
@@ -43,7 +43,6 @@ const std::vector<real_type> answer_key_N2 = {
     -3.01001804456668e-06,
     63.2193618667302,
     -35.0920143829141,
-    -0.00755583999797958,
     22975.8457840701,
     8742.0172166684,
     0.0371009878165673,
@@ -56,7 +55,7 @@ const std::vector<real_type> answer_key_N2 = {
     3.79041261513032e-05,
     62.3296560841534,
     -23.7157298601546,
-    -0.0827401584095097,
+    -0.00755583999797958,
     17277.712521392,
     16493.7578328327,
     0.0311645357391787,
@@ -69,7 +68,7 @@ const std::vector<real_type> answer_key_N2 = {
     -0.000167008045114112,
     48.6166481095298,
     -46.4283022775928,
-    -0.0844566209033113,
+    -0.0827401584095097,
     17277.2493364959,
     9182.29479881977,
     0.0302503981389185,
@@ -82,6 +81,7 @@ const std::vector<real_type> answer_key_N2 = {
     0.000121054462396047,
     47.19063792117,
     -25.0705499054599,
+    -0.0844566209033113,
     -18.8125403122072,
     21.1471334522508,
     43.2997340692497,
@@ -115,7 +115,6 @@ const std::vector<real_type> answer_key_N4 = {
     -0.000167782960864797,
     76.0120710728117,
     -28.9413214479885,
-    -0.0154989336357971,
     27832.0772154743,
     7838.80352597054,
     0.0448129806183038,
@@ -128,7 +127,7 @@ const std::vector<real_type> answer_key_N4 = {
     -5.17440472732025e-05,
     75.285919083362,
     -21.1922719962368,
-    -0.142596248886844,
+    -0.0154989336357971,
     20950.1519748442,
     16334.3924580166,
     0.0377787488551238,
@@ -141,7 +140,7 @@ const std::vector<real_type> answer_key_N4 = {
     -8.66167274613643e-05,
     58.9348242715342,
     -45.9457664944992,
-    -0.164628885385433,
+    -0.142596248886844,
     20956.2846481359,
     12345.7855779785,
     0.0371673110634281,
@@ -154,7 +153,7 @@ const std::vector<real_type> answer_key_N4 = {
     8.24481547667206e-05,
     57.9810719040216,
     -34.1573355283559,
-    -0.245126576498858,
+    -0.164628885385433,
     20978.9206932923,
     17875.6825958711,
     0.0380785081131898,
@@ -167,7 +166,7 @@ const std::vector<real_type> answer_key_N4 = {
     5.52721744209165e-05,
     59.4025238741703,
     -50.6117204294148,
-    -0.253672269778405,
+    -0.245126576498858,
     20990.6142554844,
     11779.5685373245,
     0.0371243965775962,
@@ -180,7 +179,7 @@ const std::vector<real_type> answer_key_N4 = {
     0.000209247259473183,
     57.9140782926947,
     -32.5257339504705,
-    -0.288747285297145,
+    -0.253672269778405,
     21004.5418902462,
     16686.981801999,
     0.0379131316203435,
@@ -193,7 +192,7 @@ const std::vector<real_type> answer_key_N4 = {
     8.14114753778542e-05,
     59.1444846638236,
     -47.0195494892644,
-    -0.285876486857154,
+    -0.288747285297145,
     21005.895228757,
     8274.13913126834,
     0.0366267978764156,
@@ -206,6 +205,7 @@ const std::vector<real_type> answer_key_N4 = {
     0.000263351556352378,
     57.1378679159653,
     -22.518140603782,
+    -0.285876486857154,
     -6.57009481597077,
     28.1053117219546,
     68.341829672125,
@@ -259,7 +259,6 @@ const std::vector<real_type> answer_key_N8 = {
     -0.00135963018489767,
     79.8611351716573,
     -26.7318789917162,
-    -0.018074876389402,
     29296.9645459408,
     7415.51099954903,
     0.0471789367393934,
@@ -272,7 +271,7 @@ const std::vector<real_type> answer_key_N8 = {
     -0.00065261697325357,
     79.261388452934,
     -19.9860446882469,
-    -0.161178046356866,
+    -0.018074876389402,
     22178.7337755581,
     15750.48869542,
     0.0399575496599213,
@@ -285,7 +284,7 @@ const std::vector<real_type> answer_key_N8 = {
     -0.000888421043595241,
     62.3337589017493,
     -44.1546192236143,
-    -0.191206982234175,
+    -0.161178046356866,
     22245.3587939821,
     12476.0312983393,
     0.0395296431152022,
@@ -298,7 +297,7 @@ const std::vector<real_type> answer_key_N8 = {
     0.000606091049765824,
     61.6667589580497,
     -34.5611135316627,
-    -0.300475535038271,
+    -0.191206982234175,
     22471.354929649,
     17410.9400696732,
     0.0407475331761201,
@@ -311,7 +310,7 @@ const std::vector<real_type> answer_key_N8 = {
     -0.00050508442766251,
     63.5661867600733,
     -49.1745383652835,
-    -0.321406735230742,
+    -0.300475535038271,
     22562.1063953332,
     12532.040739496,
     0.0400730650460856,
@@ -324,7 +323,7 @@ const std::vector<real_type> answer_key_N8 = {
     0.00123813270414561,
     62.5146115608051,
     -34.7565572443792,
-    -0.398765881347133,
+    -0.321406735230742,
     22822.9625175864,
     17526.9123376124,
     0.0413769691613727,
@@ -337,7 +336,7 @@ const std::vector<real_type> answer_key_N8 = {
     -0.000255137918253834,
     64.548036917491,
     -49.5382706457328,
-    -0.409294621792759,
+    -0.398765881347133,
     22917.4640781459,
     11700.1980576741,
     0.0405370737861658,
@@ -350,7 +349,7 @@ const std::vector<real_type> answer_key_N8 = {
     0.00143605100750772,
     63.2384503029344,
     -32.3525635247183,
-    -0.460038982115987,
+    -0.409294621792759,
     23175.0870418607,
     17429.0372806116,
     0.0419717095772374,
@@ -363,7 +362,7 @@ const std::vector<real_type> answer_key_N8 = {
     -0.000116590766441278,
     65.4757901489451,
     -49.2486679790579,
-    -0.461848591689221,
+    -0.460038982115987,
     23264.6732405303,
     10919.0360806903,
     0.0409947737541995,
@@ -376,7 +375,7 @@ const std::vector<real_type> answer_key_N8 = {
     0.00154556746783549,
     63.9523565934322,
     -30.1044436917781,
-    -0.491529580804623,
+    -0.461848591689221,
     23488.6137650364,
     17313.8654584652,
     0.0424952351974662,
@@ -389,7 +388,7 @@ const std::vector<real_type> answer_key_N8 = {
     0.000169934086992562,
     66.2924545883603,
     -48.9139721676174,
-    -0.486548093474613,
+    -0.491529580804623,
     23557.6136117665,
     10241.0790360359,
     0.0413774884672506,
@@ -402,7 +401,7 @@ const std::vector<real_type> answer_key_N8 = {
     0.00142420279185394,
     64.5492984953198,
     -28.1556597465189,
-    -0.499726325044563,
+    -0.486548093474613,
     23722.4168131373,
     16958.7765974266,
     0.0428487548646152,
@@ -415,7 +414,7 @@ const std::vector<real_type> answer_key_N8 = {
     0.000573219240561573,
     66.844044569212,
     -47.8519918014414,
-    -0.490219448053285,
+    -0.499726325044563,
     23769.8747006891,
     9267.5206785042,
     0.0415639300026818,
@@ -428,7 +427,7 @@ const std::vector<real_type> answer_key_N8 = {
     0.00128332197988872,
     64.8400012051758,
     -25.3838027935624,
-    -0.491586583798841,
+    -0.490219448053285,
     23854.0240861864,
     15160.7060140903,
     0.0427385845553072,
@@ -441,7 +440,7 @@ const std::vector<real_type> answer_key_N8 = {
     0.000734666517031364,
     66.6722602990682,
     -42.4779079553274,
-    -0.482907936930927,
+    -0.491586583798841,
     23870.5548740189,
     6166.60742794455,
     0.0412253604361869,
@@ -454,6 +453,7 @@ const std::vector<real_type> answer_key_N8 = {
     0.00130490570318174,
     64.3117765573399,
     -16.6904284244068,
+    -0.482907936930927,
     -2.9267653635869,
     30.5878043390112,
     75.9238910016987,