Similar fix to umin/umax root-finding to not skip edge case as was ne…

…cessary for actionAngleAdiabatic in 263afbd; also need to slightly change tolerance of one test, because now passes sometimes circular orbits now go through different code
jobovy · Nov 17, 2018 · 59eec24 · 59eec24
1 parent 263afbd
commit 59eec24
Show file tree

Hide file tree

Showing 3 changed files with 12 additions and 12 deletions.
diff --git a/galpy/actionAngle/actionAngle_c_ext/actionAngleAdiabatic.c b/galpy/actionAngle/actionAngle_c_ext/actionAngleAdiabatic.c
@@ -344,7 +344,7 @@ void calcRapRperi(int ndata,
     peps= GSL_FN_EVAL(JRRoot+tid,*(R+ii)+0.0000001);
     meps= GSL_FN_EVAL(JRRoot+tid,*(R+ii)-0.0000001);
     if ( fabs(GSL_FN_EVAL(JRRoot+tid,*(R+ii))) < 0.0000001 && peps*meps < 0 ){ //we are at rap or rperi
-      if ( peps < 0. && meps > 0. ) {//umax
+      if ( peps < 0. && meps > 0. ) {//rap
 	*(rap+ii)= *(R+ii);
 	R_lo= 0.9 * (*(R+ii) - 0.0000001);
 	R_hi= *(R+ii) - 0.00000001;
@@ -379,7 +379,7 @@ void calcRapRperi(int ndata,
 	// LCOV_EXCL_STOP
 	*(rperi+ii) = gsl_root_fsolver_root ((s+tid)->s);
       }
-      else if ( peps > 0. && meps < 0. ){//umin
+      else {// JB: Should catch all: if ( peps > 0. && meps < 0. ){//rperi
 	*(rperi+ii)= *(R+ii);
 	R_lo= *(R+ii) + 0.0000001;
 	R_hi= 1.1 * (*(R+ii) + 0.0000001);

diff --git a/galpy/actionAngle/actionAngle_c_ext/actionAngleStaeckel.c b/galpy/actionAngle/actionAngle_c_ext/actionAngleStaeckel.c
@@ -1553,14 +1553,10 @@ void calcUminUmax(int ndata,
     (JRRoot+tid)->function = &JRStaeckelIntegrandSquared;
     (JRRoot+tid)->params = params+tid;
     //Find starting points for minimum
-    if ( fabs(GSL_FN_EVAL(JRRoot+tid,*(ux+ii))) < 0.0000001){ //we are at umin or umax
-      peps= GSL_FN_EVAL(JRRoot+tid,*(ux+ii)+0.000001);
-      meps= GSL_FN_EVAL(JRRoot+tid,*(ux+ii)-0.000001);
-      if ( fabs(peps) < 0.00000001 && fabs(meps) < 0.00000001 ) {//circular
-	*(umin+ii) = *(ux+ii);
-	*(umax+ii) = *(ux+ii);
-      }
-      else if ( peps < 0. && meps > 0. ) {//umax
+    peps= GSL_FN_EVAL(JRRoot+tid,*(ux+ii)+0.000001);
+    meps= GSL_FN_EVAL(JRRoot+tid,*(ux+ii)-0.000001);
+    if ( fabs(GSL_FN_EVAL(JRRoot+tid,*(ux+ii))) < 0.0000001 && peps*meps < 0. ){ //we are at umin or umax
+      if ( peps < 0. && meps > 0. ) {//umax
 	*(umax+ii)= *(ux+ii);
 	u_lo= 0.9 * (*(ux+ii) - 0.000001);
 	u_hi= *(ux+ii) - 0.0000001;
@@ -1595,7 +1591,7 @@ void calcUminUmax(int ndata,
 	  *(umin+ii) = gsl_root_fsolver_root ((s+tid)->s);
 	}
       }
-      else if ( peps > 0. && meps < 0. ){//umin
+      else {// JB: Should catch all: if ( peps > 0. && meps < 0. ){//umin
 	*(umin+ii)= *(ux+ii);
 	u_lo= *(ux+ii) + 0.000001;
 	u_hi= 1.1 * (*(ux+ii) + 0.000001);
@@ -1632,6 +1628,10 @@ void calcUminUmax(int ndata,
 	*(umax+ii) = gsl_root_fsolver_root ((s+tid)->s);
       }
     }
+    else if ( fabs(peps) < 0.00000001 && fabs(meps) < 0.00000001 && peps <= 0 && meps <= 0 ) {//circular
+	*(umin+ii) = *(ux+ii);
+	*(umax+ii) = *(ux+ii);
+    }
     else {
       u_lo= 0.9 * *(ux+ii);
       u_hi= *(ux+ii);

diff --git a/tests/test_actionAngle.py b/tests/test_actionAngle.py
@@ -923,7 +923,7 @@ def test_actionAngleStaeckel_basic_actions_u0_interppot_c():
     #circular orbit
     R,vR,vT,z,vz= 1.,0.,1.,0.,0. 
     js= aAS(R,vR,vT,z,vz)
-    assert numpy.fabs(js[0][0]) < 10.**-16., 'Circular orbit in the MWPotential does not have Jr=0'
+    assert numpy.fabs(js[0][0]) < 10.**-12., 'Circular orbit in the MWPotential does not have Jr=0'
     assert numpy.fabs(js[2][0]) < 10.**-16., 'Circular orbit in the MWPotential does not have Jz=0'
     #Close-to-circular orbit
     R,vR,vT,z,vz= 1.01,0.01,1.,0.01,0.01