Fix to prevent error in gc tracing when vpar = 0.

hiddenSymmetries · Jan 26, 2024 · f8a5668 · f8a5668
1 parent db10d50
commit f8a5668
Showing 1 changed file with 6 additions and 2 deletions.
diff --git a/src/simsoptpp/tracing.cpp b/src/simsoptpp/tracing.cpp
@@ -418,6 +418,8 @@ class GuidingCenterNoKBoozerRHS {
             double v_perp2 = 2*mu*modB;
             double fak1 = m*v_par*v_par/modB + m*mu;
             double D = ((q + m*v_par*dIdpsi/modB)*G - (-q*iota + m*v_par*dGdpsi/modB)*I)/iota;
+            double F = (q + m*v_par*dIdpsi/modB);
+            double C = (-q*iota + m*v_par*dGdpsi/modB);
 
             double sdot = (I*dmodBdzeta - G*dmodBdtheta)*fak1/(D*iota*psi0);
             double tdot = (G*dmodBdpsi*fak1 - (-q*iota + m*v_par*dGdpsi/modB)*v_par*modB)/(D*iota);
@@ -432,7 +434,8 @@ class GuidingCenterNoKBoozerRHS {
                 dydt[1] = tdot;
             }
             dydt[2] = ((q + m*v_par*dIdpsi/modB)*v_par*modB - dmodBdpsi*fak1*I)/(D*iota);
-            dydt[3] = - (mu / v_par) * (dmodBdpsi * sdot * psi0 + dmodBdtheta * tdot + dmodBdzeta * dydt[2]);
+            dydt[3] = modB*mu*(dmodBdtheta*C - dmodBdzeta*F)/(F*G-C*I);
+            // dydt[3] = - (mu / v_par) * (dmodBdpsi * sdot * psi0 + dmodBdtheta * tdot + dmodBdzeta * dydt[2]);
         }
 };
 
@@ -517,7 +520,8 @@ class GuidingCenterBoozerRHS {
                 dydt[1] = tdot;
             }
             dydt[2] = (F*v_par*modB - dmodBdpsi*fak1*I + K*fak1*dmodBdtheta)/(D*iota);
-            dydt[3] = - (mu / v_par) * (dmodBdpsi * sdot * psi0 + dmodBdtheta * tdot + dmodBdzeta * dydt[2]);
+            // dydt[3] = - (mu / v_par) * (dmodBdpsi * sdot * psi0 + dmodBdtheta * tdot + dmodBdzeta * dydt[2]);
+            dydt[3] = modB*mu*(dmodBdtheta*C - dmodBdzeta*F)/(F*G-C*I);
         }
 };