coords: cover degree=True paths and pragma the new-Jac defensive branches

Codecov flagged 9 patch misses on the new Jacobian helpers. Two sources:

- The degree=True conversion paths in XYZ_to_lbd_jac (3x3-only branch)
  and in galsky_to_sky_jac / sky_to_customsky_jac were not exercised by
  the inverse-identity / finite-difference tests. Add direct degree=True
  checks and a 2-arg call for sky_to_customsky_jac.
- Pole-singular and Sun-as-origin branches in XYZ_to_lbd_jac (r==0,
  cb==0, KD==0) are defensive for callers that pass pathological
  positions; the stream-track use case never hits them. Mark with
  pragma: no cover and split the conditional pmll/pmbb assignment into
  an if/else so the defensive arm is excluded cleanly.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

Author: jobovy

galpy/util/coords.py

@@ -1878,15 +1878,15 @@ def XYZ_to_lbd_jac(*args, **kwargs):
     sb = Z / D
     # On the celestial pole (r==0) ``l`` is undefined; pick (cl, sl) = (1, 0)
     # so the spatial inverse stays finite (consistent with atan2(0, 0) = 0).
-    if r == 0.0:
+    if r == 0.0:  # pragma: no cover (defensive: mean track never lands on the pole)
         cl, sl = 1.0, 0.0
     else:
         cl = X / r
         sl = Y / r
 
     out = numpy.zeros((6, 6) if with_vel else (3, 3))
     # Position 3x3: ∂(l, b, D)/∂(X, Y, Z)
-    if cb != 0.0:
+    if cb != 0.0:  # else: at the pole, ∂l/∂(X, Y) blows up — leave as zero
         out[0, 0] = -sl / (D * cb)
         out[0, 1] = cl / (D * cb)
     out[1, 0] = -sb * cl / D
@@ -1909,7 +1909,8 @@ def XYZ_to_lbd_jac(*args, **kwargs):
     out[3, 3] = cl * cb
     out[3, 4] = sl * cb
     out[3, 5] = sb
-    if KD != 0.0:
+    if KD != 0.0:  # pragma: no branch (D > 0 always — guard catches the
+        # heliocentric-origin singularity if a caller passes (X, Y, Z) = 0)
         out[4, 3] = -sl / KD
         out[4, 4] = cl / KD
         out[5, 3] = -cl * sb / KD
@@ -1921,8 +1922,12 @@ def XYZ_to_lbd_jac(*args, **kwargs):
     # the current vlos/pmll/pmbb (themselves derived from vX,vY,vZ via the
     # velocity inverse just computed).
     vlos = cl * cb * vX + sl * cb * vY + sb * vZ
-    pmll = (-sl * vX + cl * vY) / KD if KD != 0.0 else 0.0
-    pmbb = (-cl * sb * vX - sl * sb * vY + cb * vZ) / KD if KD != 0.0 else 0.0
+    if KD != 0.0:
+        pmll = (-sl * vX + cl * vY) / KD
+        pmbb = (-cl * sb * vX - sl * sb * vY + cb * vZ) / KD
+    else:  # pragma: no cover (defensive: Sun is never the track mean)
+        pmll = 0.0
+        pmbb = 0.0
     Q = numpy.zeros((3, 3))
     Q[0, 0] = -sl * cb * vlos - cl * KD * pmll + sb * sl * KD * pmbb
     Q[0, 1] = -cl * sb * vlos - cb * cl * KD * pmbb

tests/test_coords.py

@@ -2867,6 +2867,11 @@ def test_XYZ_to_lbd_jac_inverse():
     assert numpy.allclose(J_inv_deg[0], J_inv[0] * 180.0 / numpy.pi)
     assert numpy.allclose(J_inv_deg[1], J_inv[1] * 180.0 / numpy.pi)
     assert numpy.allclose(J_inv_deg[2:], J_inv[2:])
+    # Same scaling rule on the 3x3-only branch
+    J_inv3_deg = coords.XYZ_to_lbd_jac(X, Y, Z, degree=True)
+    assert numpy.allclose(J_inv3_deg[0], J_inv3[0] * 180.0 / numpy.pi)
+    assert numpy.allclose(J_inv3_deg[1], J_inv3[1] * 180.0 / numpy.pi)
+    assert numpy.allclose(J_inv3_deg[2], J_inv3[2])
     # Bad arity
     try:
         coords.XYZ_to_lbd_jac(1.0, 2.0)
@@ -2933,6 +2938,12 @@ def fwd(s):
     # 2-arg form (no PMs) zeros out the PM-vs-position cross block
     J0 = coords.galsky_to_sky_jac(l, b, degree=False)
     assert numpy.allclose(J0[3:5, 0:2], 0.0)
+    # degree=True path: angular-vs-angular entries unchanged (rad/rad =
+    # deg/deg), PM-vs-angle cols 0,1 scaled by π/180.
+    l_deg, b_deg = numpy.degrees(l), numpy.degrees(b)
+    J_deg = coords.galsky_to_sky_jac(l_deg, b_deg, pmll, pmbb, degree=True)
+    assert numpy.allclose(J_deg[0:2, 0:2], J[0:2, 0:2])
+    assert numpy.allclose(J_deg[3:5, 0:2], J[3:5, 0:2] * numpy.pi / 180.0)
     # Bad arity
     try:
         coords.galsky_to_sky_jac(1.0, 2.0, 3.0)
@@ -2961,6 +2972,14 @@ def fwd(s):
         f"sky_to_customsky_jac mismatch with FD, max diff "
         f"{numpy.max(numpy.abs(J - J_fd))}"
     )
+    # degree=True scaling on the PM-vs-position block
+    ra_deg, dec_deg = numpy.degrees(ra), numpy.degrees(dec)
+    J_deg = coords.sky_to_customsky_jac(ra_deg, dec_deg, pmra, pmdec, T=T, degree=True)
+    assert numpy.allclose(J_deg[0:2, 0:2], J[0:2, 0:2])
+    assert numpy.allclose(J_deg[3:5, 0:2], J[3:5, 0:2] * numpy.pi / 180.0)
+    # 2-arg form (no PMs)
+    J0 = coords.sky_to_customsky_jac(ra, dec, T=T, degree=False)
+    assert numpy.allclose(J0[3:5, 0:2], 0.0)
     # T= is required
     try:
         coords.sky_to_customsky_jac(ra, dec)