Skip to content

Commit 5925f0b

Browse files
jobovy-agentclaude
jobovy-agent
and
claude
authored
actionAngleStaeckel C: frequencies and angles near turning points (#1424)
Co-authored-by: Claude Fable 5 <noreply@anthropic.com>
1 parent 02d589b commit 5925f0b

2 files changed

Lines changed: 172 additions & 22 deletions

File tree

‎galpy/actionAngle/actionAngle_c_ext/actionAngleStaeckel.c‎

Lines changed: 106 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -295,7 +295,7 @@ void calcu0(int ndata,
295295
u_lo = gsl_min_fminimizer_x_lower (s);
296296
u_hi = gsl_min_fminimizer_x_upper (s);
297297
status = gsl_min_test_interval (u_lo, u_hi,
298-
9.9999999999999998e-13,
298+
4.4408920985006262e-16,
299299
4.4408920985006262e-16);
300300
}
301301
while (status == GSL_CONTINUE && iter < max_iter);
@@ -1333,8 +1333,8 @@ void calcAnglesStaeckel(int ndata,
13331333
midpoint= *(umin+ii)+ 0.5 * ( *(umax+ii) - *(umin+ii) );
13341334
if ( *(pux+ii) > 0. ) {
13351335
if ( *(ux+ii) > midpoint ) {
1336-
mid= 2. * asin( sqrt( ( *(umax+ii) - *(ux+ii) )
1337-
/ ( *(umax+ii) - *(umin+ii) ) ) );
1336+
mid= 2. * asin( sqrt( fmin ( fmax ( ( *(umax+ii) - *(ux+ii) )
1337+
/ ( *(umax+ii) - *(umin+ii) ), 0. ), 1. ) ) );
13381338
(AngleuInt+tid)->function = &dJRdEHighStaeckelIntegrand;
13391339
Or1= gsl_integration_glfixed (AngleuInt+tid,0.,mid,T);
13401340
(AngleuInt+tid)->function = &dJRdI3HighStaeckelIntegrand;
@@ -1347,8 +1347,8 @@ void calcAnglesStaeckel(int ndata,
13471347
I3r1= M_PI * *(dJRdI3+ii) - I3r1;
13481348
}
13491349
else {
1350-
mid= 2. * asin( sqrt( ( *(ux+ii) - *(umin+ii) )
1351-
/ ( *(umax+ii) - *(umin+ii) ) ) );
1350+
mid= 2. * asin( sqrt( fmin ( fmax ( ( *(ux+ii) - *(umin+ii) )
1351+
/ ( *(umax+ii) - *(umin+ii) ), 0. ), 1. ) ) );
13521352
(AngleuInt+tid)->function = &dJRdELowStaeckelIntegrand;
13531353
Or1= gsl_integration_glfixed (AngleuInt+tid,0.,mid,T);
13541354
(AngleuInt+tid)->function = &dJRdI3LowStaeckelIntegrand;
@@ -1361,8 +1361,8 @@ void calcAnglesStaeckel(int ndata,
13611361
}
13621362
else {
13631363
if ( *(ux+ii) > midpoint ) {
1364-
mid= 2. * asin( sqrt( ( *(umax+ii) - *(ux+ii) )
1365-
/ ( *(umax+ii) - *(umin+ii) ) ) );
1364+
mid= 2. * asin( sqrt( fmin ( fmax ( ( *(umax+ii) - *(ux+ii) )
1365+
/ ( *(umax+ii) - *(umin+ii) ), 0. ), 1. ) ) );
13661366
(AngleuInt+tid)->function = &dJRdEHighStaeckelIntegrand;
13671367
Or1= gsl_integration_glfixed (AngleuInt+tid,0.,mid,T);
13681368
Or1*= *(delta+ii*delta_stride) / sqrt(2.);
@@ -1375,8 +1375,8 @@ void calcAnglesStaeckel(int ndata,
13751375
*(Anglephi+ii)= M_PI * *(dJRdLz+ii) - *(Lz+ii) * gsl_integration_glfixed (AngleuInt+tid,0.,mid,T) / *(delta+ii*delta_stride) / sqrt(2.);
13761376
}
13771377
else {
1378-
mid= 2. * asin( sqrt( ( *(ux+ii) - *(umin+ii) )
1379-
/ ( *(umax+ii) - *(umin+ii) ) ) );
1378+
mid= 2. * asin( sqrt( fmin ( fmax ( ( *(ux+ii) - *(umin+ii) )
1379+
/ ( *(umax+ii) - *(umin+ii) ), 0. ), 1. ) ) );
13801380
(AngleuInt+tid)->function = &dJRdELowStaeckelIntegrand;
13811381
Or1= gsl_integration_glfixed (AngleuInt+tid,0.,mid,T);
13821382
Or1*= *(delta+ii*delta_stride) / sqrt(2.);
@@ -1405,10 +1405,10 @@ void calcAnglesStaeckel(int ndata,
14051405
if ( *(vx+ii) < midpoint || *(vx+ii) > (M_PI - midpoint) ) {
14061406
// chi of the current v, measured from the vmin turning point along
14071407
// the full loop [vmin, pi - vmin]; v beyond the midplane mirrors
1408-
mid = 2. * asin( sqrt( ( ( *(vx+ii) > 0.5 * M_PI )
1408+
mid = 2. * asin( sqrt( fmin ( fmax ( ( ( *(vx+ii) > 0.5 * M_PI )
14091409
? ( M_PI - *(vx+ii) - *(vmin+ii) )
14101410
: ( *(vx+ii) - *(vmin+ii) ) )
1411-
/ ( M_PI - 2. * *(vmin+ii) ) ) );
1411+
/ ( M_PI - 2. * *(vmin+ii) ), 0. ), 1. ) ) );
14121412
(AnglevInt+tid)->function = &dJzdELowStaeckelIntegrand;
14131413
Or2= gsl_integration_glfixed (AnglevInt+tid,0.,mid,T);
14141414
Or2*= *(delta+ii*delta_stride) / sqrt(2.);
@@ -1451,10 +1451,10 @@ void calcAnglesStaeckel(int ndata,
14511451
if ( *(vx+ii) < midpoint || *(vx+ii) > (M_PI - midpoint)) {
14521452
// chi of the current v, measured from the vmin turning point along
14531453
// the full loop [vmin, pi - vmin]; v beyond the midplane mirrors
1454-
mid = 2. * asin( sqrt( ( ( *(vx+ii) > 0.5 * M_PI )
1454+
mid = 2. * asin( sqrt( fmin ( fmax ( ( ( *(vx+ii) > 0.5 * M_PI )
14551455
? ( M_PI - *(vx+ii) - *(vmin+ii) )
14561456
: ( *(vx+ii) - *(vmin+ii) ) )
1457-
/ ( M_PI - 2. * *(vmin+ii) ) ) );
1457+
/ ( M_PI - 2. * *(vmin+ii) ), 0. ), 1. ) ) );
14581458
(AnglevInt+tid)->function = &dJzdELowStaeckelIntegrand;
14591459
Or2= gsl_integration_glfixed (AnglevInt+tid,0.,mid,T);
14601460
Or2*= *(delta+ii*delta_stride) / sqrt(2.);
@@ -1591,7 +1591,35 @@ void calcUminUmax(int ndata,
15911591
meps= GSL_FN_EVAL(JRRoot+tid,*(ux+ii)-0.000001);
15921592
if ( fabs(GSL_FN_EVAL(JRRoot+tid,*(ux+ii))) < 0.0000001 && peps*meps < 0. ){ //we are at umin or umax
15931593
if ( peps < 0. && meps > 0. ) {//umax
1594-
*(umax+ii)= *(ux+ii);
1594+
//The point lies within 1e-6 of umax. Adopting the point itself as
1595+
//the turning point is fine for the actions (an O(eps) endpoint
1596+
//error only enters them at O(eps^1.5)), but the frequency and
1597+
//angle integrands diverge as 1/sqrt(W) at the endpoint, so there
1598+
//it enters at O(sqrt(eps)); solve for the true root instead (the
1599+
//+-1e-6 bracket has the sign change by the branch condition),
1600+
//keeping the point itself only as the degenerate fallback
1601+
status = gsl_root_fsolver_set ((s+tid)->s, JRRoot+tid,
1602+
*(ux+ii) - 0.000001,
1603+
*(ux+ii) + 0.000001);
1604+
if (status == GSL_EINVAL)
1605+
// LCOV_EXCL_START
1606+
*(umax+ii)= *(ux+ii);//can't happen: the branch condition peps*meps<0 is the bracket
1607+
// LCOV_EXCL_STOP
1608+
else {
1609+
iter= 0;
1610+
do
1611+
{
1612+
iter++;
1613+
status = gsl_root_fsolver_iterate ((s+tid)->s);
1614+
u_lo = gsl_root_fsolver_x_lower ((s+tid)->s);
1615+
u_hi = gsl_root_fsolver_x_upper ((s+tid)->s);
1616+
status = gsl_root_test_interval (u_lo, u_hi,
1617+
4.4408920985006262e-16,
1618+
4.4408920985006262e-16);
1619+
}
1620+
while (status == GSL_CONTINUE && iter < max_iter);
1621+
*(umax+ii)= gsl_root_fsolver_root ((s+tid)->s);
1622+
}
15951623
u_lo= 0.9 * (*(ux+ii) - 0.000001);
15961624
u_hi= *(ux+ii) - 0.0000001;
15971625
while ( GSL_FN_EVAL(JRRoot+tid,u_lo) >= 0. && u_lo > 0.000000001){
@@ -1611,7 +1639,7 @@ void calcUminUmax(int ndata,
16111639
u_lo = gsl_root_fsolver_x_lower ((s+tid)->s);
16121640
u_hi = gsl_root_fsolver_x_upper ((s+tid)->s);
16131641
status = gsl_root_test_interval (u_lo, u_hi,
1614-
9.9999999999999998e-13,
1642+
4.4408920985006262e-16,
16151643
4.4408920985006262e-16);
16161644
}
16171645
while (status == GSL_CONTINUE && iter < max_iter);
@@ -1626,7 +1654,30 @@ void calcUminUmax(int ndata,
16261654
}
16271655
}
16281656
else {// JB: Should catch all: if ( peps > 0. && meps < 0. ){//umin
1629-
*(umin+ii)= *(ux+ii);
1657+
//Same as the umax case above: solve for the true root, with the
1658+
//point itself only as the degenerate fallback
1659+
status = gsl_root_fsolver_set ((s+tid)->s, JRRoot+tid,
1660+
*(ux+ii) - 0.000001,
1661+
*(ux+ii) + 0.000001);
1662+
if (status == GSL_EINVAL)
1663+
// LCOV_EXCL_START
1664+
*(umin+ii)= *(ux+ii);//can't happen: the branch condition peps*meps<0 is the bracket
1665+
// LCOV_EXCL_STOP
1666+
else {
1667+
iter= 0;
1668+
do
1669+
{
1670+
iter++;
1671+
status = gsl_root_fsolver_iterate ((s+tid)->s);
1672+
u_lo = gsl_root_fsolver_x_lower ((s+tid)->s);
1673+
u_hi = gsl_root_fsolver_x_upper ((s+tid)->s);
1674+
status = gsl_root_test_interval (u_lo, u_hi,
1675+
4.4408920985006262e-16,
1676+
4.4408920985006262e-16);
1677+
}
1678+
while (status == GSL_CONTINUE && iter < max_iter);
1679+
*(umin+ii)= gsl_root_fsolver_root ((s+tid)->s);
1680+
}
16301681
u_lo= *(ux+ii) + 0.000001;
16311682
u_hi= 1.1 * (*(ux+ii) + 0.000001);
16321683
while ( GSL_FN_EVAL(JRRoot+tid,u_hi) >= 0. && u_hi < asinh(37.5/ *(delta+ii*delta_stride))) {
@@ -1648,7 +1699,7 @@ void calcUminUmax(int ndata,
16481699
u_lo = gsl_root_fsolver_x_lower ((s+tid)->s);
16491700
u_hi = gsl_root_fsolver_x_upper ((s+tid)->s);
16501701
status = gsl_root_test_interval (u_lo, u_hi,
1651-
9.9999999999999998e-13,
1702+
4.4408920985006262e-16,
16521703
4.4408920985006262e-16);
16531704
}
16541705
while (status == GSL_CONTINUE && iter < max_iter);
@@ -1687,7 +1738,7 @@ void calcUminUmax(int ndata,
16871738
u_lo = gsl_root_fsolver_x_lower ((s+tid)->s);
16881739
u_hi = gsl_root_fsolver_x_upper ((s+tid)->s);
16891740
status = gsl_root_test_interval (u_lo, u_hi,
1690-
9.9999999999999998e-13,
1741+
4.4408920985006262e-16,
16911742
4.4408920985006262e-16);
16921743
}
16931744
while (status == GSL_CONTINUE && iter < max_iter);
@@ -1723,7 +1774,7 @@ void calcUminUmax(int ndata,
17231774
u_lo = gsl_root_fsolver_x_lower ((s+tid)->s);
17241775
u_hi = gsl_root_fsolver_x_upper ((s+tid)->s);
17251776
status = gsl_root_test_interval (u_lo, u_hi,
1726-
9.9999999999999998e-13,
1777+
4.4408920985006262e-16,
17271778
4.4408920985006262e-16);
17281779
}
17291780
while (status == GSL_CONTINUE && iter < max_iter);
@@ -1803,8 +1854,41 @@ void calcVmin(int ndata,
18031854
(JzRoot+tid)->function = &JzStaeckelIntegrandSquared;
18041855
(JzRoot+tid)->params = params+tid;
18051856
//Find starting points for minimum
1806-
if ( fabs(GSL_FN_EVAL(JzRoot+tid,*(vx+ii))) < 0.0000001) //we are at vmin
1807-
*(vmin+ii)= ( *(vx+ii) > 0.5 * M_PI ) ? M_PI - *(vx+ii): *(vx+ii);
1857+
if ( fabs(GSL_FN_EVAL(JzRoot+tid,*(vx+ii))) < 0.0000001) {//we are at vmin
1858+
//The point lies within 1e-6 of vmin. Adopting the point itself as
1859+
//the turning point is fine for the actions (an O(eps) endpoint
1860+
//error only enters them at O(eps^1.5)), but the frequency and
1861+
//angle integrands diverge as 1/sqrt(W) at the endpoint, so there
1862+
//it enters at O(sqrt(eps)); solve for the true root instead: W is
1863+
//negative below vmin and positive above, so the +-1e-6 bracket
1864+
//around the (mirrored) point holds unless the point is degenerate
1865+
//at roundoff level, where the point IS the turning point and the
1866+
//old assignment remains as the fallback
1867+
v_lo= ( ( *(vx+ii) > 0.5 * M_PI ) ? M_PI - *(vx+ii): *(vx+ii) )
1868+
- 0.000001;
1869+
if ( v_lo < 0.000000001 ) v_lo= 0.000000001;
1870+
v_hi= ( ( *(vx+ii) > 0.5 * M_PI ) ? M_PI - *(vx+ii): *(vx+ii) )
1871+
+ 0.000001;
1872+
if ( v_hi > 0.5 * M_PI ) v_hi= 0.5 * M_PI;
1873+
status = gsl_root_fsolver_set ((s+tid)->s, JzRoot+tid, v_lo, v_hi);
1874+
if (status == GSL_EINVAL)
1875+
*(vmin+ii)= ( *(vx+ii) > 0.5 * M_PI ) ? M_PI - *(vx+ii): *(vx+ii);
1876+
else {
1877+
iter= 0;
1878+
do
1879+
{
1880+
iter++;
1881+
status = gsl_root_fsolver_iterate ((s+tid)->s);
1882+
v_lo = gsl_root_fsolver_x_lower ((s+tid)->s);
1883+
v_hi = gsl_root_fsolver_x_upper ((s+tid)->s);
1884+
status = gsl_root_test_interval (v_lo, v_hi,
1885+
4.4408920985006262e-16,
1886+
4.4408920985006262e-16);
1887+
}
1888+
while (status == GSL_CONTINUE && iter < max_iter);
1889+
*(vmin+ii)= gsl_root_fsolver_root ((s+tid)->s);
1890+
}
1891+
}
18081892
else {
18091893
if ( *(vx+ii) > 0.5 * M_PI ){
18101894
v_lo= 0.9 * ( M_PI - *(vx+ii) );
@@ -1832,7 +1916,7 @@ void calcVmin(int ndata,
18321916
v_lo = gsl_root_fsolver_x_lower ((s+tid)->s);
18331917
v_hi = gsl_root_fsolver_x_upper ((s+tid)->s);
18341918
status = gsl_root_test_interval (v_lo, v_hi,
1835-
9.9999999999999998e-13,
1919+
4.4408920985006262e-16,
18361920
4.4408920985006262e-16);
18371921
}
18381922
while (status == GSL_CONTINUE && iter < max_iter);

‎tests/test_actionAngle.py‎

Lines changed: 66 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2521,6 +2521,72 @@ def test_actionAngleStaeckel_zerolz_actions_c():
25212521

25222522

25232523
# Check that precision increases with increasing Gauss-Legendre order
2524+
def test_actionAngleStaeckel_c_angles_freqs_near_turning_point():
2525+
# c=True frequencies and angles used to go wrong for points near a
2526+
# turning point: when |p^2| < 1e-7 there, the C code adopted the
2527+
# evaluation point ITSELF as the turning point instead of solving for
2528+
# it. That O(eps) endpoint error enters the actions only at
2529+
# O(eps^1.5) (invisible, ~1e-12) but the 1/sqrt(W)-divergent
2530+
# frequency and angle integrands at O(sqrt(eps)) -- theta_z errors up
2531+
# to ~1e-4 -- and, through ~1e-12-absolute root tolerances, also left
2532+
# an order-growing angle error on generic points. Frequencies are
2533+
# torus constants and the Python path is exact here, so both provide
2534+
# sharp regression checks.
2535+
import numpy
2536+
2537+
from galpy.actionAngle import actionAngleStaeckel
2538+
from galpy.orbit import Orbit
2539+
from galpy.potential import KuzminKutuzovStaeckelPotential
2540+
2541+
kkp = KuzminKutuzovStaeckelPotential(amp=4.0, ac=5.0, Delta=1.3)
2542+
o = Orbit([1.1, 0.35, 1.1, 0.3, 0.25, 0.0])
2543+
ts = numpy.linspace(0.0, 8.0, 17)
2544+
o.integrate(ts, kkp)
2545+
R, vR, vT, z, vz, phi = (
2546+
numpy.array([float(f(t)) for t in ts])
2547+
for f in (o.R, o.vR, o.vT, o.z, o.vz, o.phi)
2548+
)
2549+
# a point of this orbit within ~1e-8 of its upper vertical turning
2550+
# point (p_v^2 ~ 6e-8), given as a literal so the near-turning
2551+
# regime is hit deterministically on every platform
2552+
R[-1], vR[-1], vT[-1], z[-1], vz[-1], phi[-1] = (
2553+
1.0374878950211397,
2554+
0.26453976273320806,
2555+
1.166278667740467,
2556+
0.3280135549475248,
2557+
0.0331712140016195,
2558+
-2.550334704756344,
2559+
)
2560+
aAC = actionAngleStaeckel(pot=kkp, delta=1.3, c=True, order=100)
2561+
aAP = actionAngleStaeckel(pot=kkp, delta=1.3, c=False)
2562+
C = aAC.actionsFreqsAngles(R, vR, vT, z, vz, phi)
2563+
# frequencies are torus constants: they may not vary along the orbit
2564+
# (the near-turning point used to be off by ~2e-4 in Omega_z)
2565+
for k, name in ((3, "Omega_R"), (5, "Omega_z")):
2566+
spread = numpy.ptp(numpy.array(C[k])) / numpy.fabs(
2567+
numpy.median(numpy.array(C[k]))
2568+
)
2569+
assert spread < 1e-8, (
2570+
"c=True %s varies along an orbit by %g near a turning point"
2571+
% (name, spread)
2572+
)
2573+
# angles agree with the (exact) Python path pointwise, including at
2574+
# the near-turning sample (used to be off by ~4e-4 in theta_z)
2575+
P = aAP.actionsFreqsAngles(R, vR, vT, z, vz, phi)
2576+
for k, name in ((6, "theta_R"), (8, "theta_z")):
2577+
d = (
2578+
numpy.remainder(
2579+
numpy.array(C[k]) - numpy.array(P[k]) + numpy.pi, 2.0 * numpy.pi
2580+
)
2581+
- numpy.pi
2582+
)
2583+
assert numpy.max(numpy.fabs(d)) < 1e-6, (
2584+
"c=True %s disagrees with c=False by %g near a turning point"
2585+
% (name, numpy.max(numpy.fabs(d)))
2586+
)
2587+
return None
2588+
2589+
25242590
def test_actionAngleStaeckel_actions_order():
25252591
from galpy.actionAngle import actionAngleStaeckel
25262592
from galpy.orbit import Orbit

0 commit comments

Comments
 (0)