Repository: atcollab/at (language: python)

Issue:
pyat detuning help message and lattice periodicity
Dear all,
the function `at.detuning` in pyat calculates the derivative wrt action $J=\frac{x^2\gamma}{2}$ where $\gamma=\frac{1+\alpha^2}{\beta}$.
The help message for `q1` should be `[[dQx/dJx, dQy/dJx], [dQx/dJy, dQy/dJy]]`.
Instead, the help message states it is the derivative wrt to $x$ and $y$.
```python
>>>help(at.detuning)
detuning(ring: at.lattice.lattice_object.Lattice, xm: Optional[float] = 3e-05, ym: Optional[float] = 3e-05, npoints:
...
        q1 (ndarray): Amplitude detuning coefficients with shape (2, 2)
                      ``[[dQx/dx, dQy/dx], [dQx/dy, dQy/dy]]``
```
it also seems to fail when periodicity is larger than one.
See https://github.com/atcollab/at/pull/821#issuecomment-2353087386


Discussion / hints:

Fix patch:
diff --git a/pyat/at/physics/linear.py b/pyat/at/physics/linear.py
index 9a81e7d41..ea27ca17b 100644
--- a/pyat/at/physics/linear.py
+++ b/pyat/at/physics/linear.py
@@ -1,13 +1,13 @@
 """
 Coupled or non-coupled 4x4 linear motion
 """
+
 from __future__ import annotations
 
 import warnings
 from collections.abc import Callable
 from math import sqrt, pi, sin, cos, atan2
 
-import numpy
 import numpy as np
 from scipy.linalg import solve
 
@@ -22,71 +22,86 @@
 from ..lattice import frequency_control
 from ..tracking import internal_lpass
 
-__all__ = ['linopt', 'linopt2', 'linopt4', 'linopt6', 'avlinopt',
-           'get_optics', 'get_tune', 'get_chrom']
+__all__ = [
+    "linopt",
+    "linopt2",
+    "linopt4",
+    "linopt6",
+    "avlinopt",
+    "get_optics",
+    "get_tune",
+    "get_chrom",
+]
 
 _S = jmat(1)
-_S2 = numpy.array([[0, 1], [-1, 0]], dtype=numpy.float64)
+_S2 = np.array([[0, 1], [-1, 0]], dtype=np.float64)
 
 # dtype for structured array containing linopt parameters
-_DATA2_DTYPE = [('alpha', numpy.float64, (2,)),
-                ('beta', numpy.float64, (2,)),
-                ('mu', numpy.float64, (2,))
-                ]
-
-_DATA4_DTYPE = [('alpha', numpy.float64, (2,)),
-                ('beta', numpy.float64, (2,)),
-                ('mu', numpy.float64, (2,)),
-                ('gamma', numpy.float64),
-                ('A', numpy.float64, (2, 2)),
-                ('B', numpy.float64, (2, 2)),
-                ('C', numpy.float64, (2, 2))
-                ]
-
-_DATA6_DTYPE = [('alpha', numpy.float64, (2,)),
-                ('beta', numpy.float64, (2,)),
-                ('mu', numpy.float64, (3,)),
-                ('R', numpy.float64, (3, 6, 6)),
-                ('A', numpy.float64, (6, 6)),
-                ('dispersion', numpy.float64, (4,)),
-                ]
-
-_DATAX_DTYPE = [('alpha', numpy.float64, (2,)),
-                ('beta', numpy.float64, (2,)),
-                ('mu', numpy.float64, (2,)),
-                ('R', numpy.float64, (2, 4, 4)),
-                ('A', numpy.float64, (4, 4)),
-                ]
-
-_W24_DTYPE = [('W', numpy.float64, (2,)),
-              ('Wp', numpy.float64, (2,)),
-              ('dalpha', numpy.float64, (2,)),
-              ('dbeta', numpy.float64, (2,)),
-              ('dmu', numpy.float64, (2,)),
-              ('ddispersion', numpy.float64, (4,)),
-              ]
-
-
-_W6_DTYPE = [('W', numpy.float64, (2,)),
-             ('Wp', numpy.float64, (2,)),
-             ('dalpha', numpy.float64, (2,)),
-             ('dbeta', numpy.float64, (2,)),
-             ('dmu', numpy.float64, (3,)),
-             ('dR', numpy.float64, (3, 6, 6)),
-             ('ddispersion', numpy.float64, (4,)),
-             ]
-
-_WX_DTYPE = [('W', numpy.float64, (2,)),
-             ('Wp', numpy.float64, (2,)),
-             ('dalpha', numpy.float64, (2,)),
-             ('dbeta', numpy.float64, (2,)),
-             ('dmu', numpy.float64, (2,)),
-             ('dR', numpy.float64, (2, 4, 4)),
-             ('ddispersion', numpy.float64, (4,)),
-             ]
-
-
-_IDX_DTYPE = [('idx', numpy.uint32)]
+_DATA2_DTYPE = [
+    ("alpha", np.float64, (2,)),
+    ("beta", np.float64, (2,)),
+    ("mu", np.float64, (2,)),
+]
+
+_DATA4_DTYPE = [
+    ("alpha", np.float64, (2,)),
+    ("beta", np.float64, (2,)),
+    ("mu", np.float64, (2,)),
+    ("gamma", np.float64),
+    ("A", np.float64, (2, 2)),
+    ("B", np.float64, (2, 2)),
+    ("C", np.float64, (2, 2)),
+]
+
+_DATA6_DTYPE = [
+    ("alpha", np.float64, (2,)),
+    ("beta", np.float64, (2,)),
+    ("mu", np.float64, (3,)),
+    ("R", np.float64, (3, 6, 6)),
+    ("A", np.float64, (6, 6)),
+    ("dispersion", np.float64, (4,)),
+]
+
+_DATAX_DTYPE = [
+    ("alpha", np.float64, (2,)),
+    ("beta", np.float64, (2,)),
+    ("mu", np.float64, (2,)),
+    ("R", np.float64, (2, 4, 4)),
+    ("A", np.float64, (4, 4)),
+]
+
+_W24_DTYPE = [
+    ("W", np.float64, (2,)),
+    ("Wp", np.float64, (2,)),
+    ("dalpha", np.float64, (2,)),
+    ("dbeta", np.float64, (2,)),
+    ("dmu", np.float64, (2,)),
+    ("ddispersion", np.float64, (4,)),
+]
+
+
+_W6_DTYPE = [
+    ("W", np.float64, (2,)),
+    ("Wp", np.float64, (2,)),
+    ("dalpha", np.float64, (2,)),
+    ("dbeta", np.float64, (2,)),
+    ("dmu", np.float64, (3,)),
+    ("dR", np.float64, (3, 6, 6)),
+    ("ddispersion", np.float64, (4,)),
+]
+
+_WX_DTYPE = [
+    ("W", np.float64, (2,)),
+    ("Wp", np.float64, (2,)),
+    ("dalpha", np.float64, (2,)),
+    ("dbeta", np.float64, (2,)),
+    ("dmu", np.float64, (2,)),
+    ("dR", np.float64, (2, 4, 4)),
+    ("ddispersion", np.float64, (4,)),
+]
+
+
+_IDX_DTYPE = [("idx", np.uint32)]
 
 
 def _twiss22(t12, alpha0, beta0):
@@ -94,26 +109,28 @@ def _twiss22(t12, alpha0, beta0):
     bbb = t12[:, 0, 1]
     aaa = t12[:, 0, 0] * beta0 - bbb * alpha0
     beta = (aaa * aaa + bbb * bbb) / beta0
-    alpha = -(aaa * (t12[:, 1, 0] * beta0 - t12[:, 1, 1] * alpha0) +
-              bbb * t12[:, 1, 1]) / beta0
-    mu = numpy.arctan2(bbb, aaa)
+    alpha = (
+        -(aaa * (t12[:, 1, 0] * beta0 - t12[:, 1, 1] * alpha0) + bbb * t12[:, 1, 1])
+        / beta0
+    )
+    mu = np.arctan2(bbb, aaa)
     # Unwrap negative jumps in betatron phase advance
-    # dmu = numpy.diff(numpy.append([0], mu))
+    # dmu = np.diff(np.append([0], mu))
     # jumps = dmu < -1.0e-3
-    # mu += numpy.cumsum(jumps) * 2.0 * numpy.pi
+    # mu += np.cumsum(jumps) * 2.0 * np.pi
     return alpha, beta, mu
 
 
 def _closure(m22):
     diff = (m22[0, 0] - m22[1, 1]) / 2.0
     try:
-        sinmu = numpy.sign(m22[0, 1]) * sqrt(-m22[0, 1]*m22[1, 0] - diff*diff)
-        cosmu = 0.5 * numpy.trace(m22)
+        sinmu = np.sign(m22[0, 1]) * sqrt(-m22[0, 1] * m22[1, 0] - diff * diff)
+        cosmu = 0.5 * np.trace(m22)
         alpha = diff / sinmu
         beta = m22[0, 1] / sinmu
-        return alpha, beta, cosmu + sinmu*1j
-    except ValueError:          # Unstable motion
-        return numpy.nan, numpy.nan, numpy.nan
+        return alpha, beta, cosmu + sinmu * 1j
+    except ValueError:  # Unstable motion
+        return np.nan, np.nan, np.nan
 
 
 # noinspection PyShadowingNames,PyPep8Naming
@@ -127,11 +144,11 @@ def _tunes(ring, **kwargs):
         mt, _ = find_m44(ring, **kwargs)
     try:
         _, vps = a_matrix(mt)
-        tunes = numpy.mod(numpy.angle(vps) / 2.0 / pi, 1.0)
+        tunes = np.mod(np.angle(vps) / 2.0 / pi, 1.0)
     except AtError:
-        warnings.warn(AtWarning('Unstable ring'))
-        tunes = numpy.empty(nd)
-        tunes[:] = numpy.nan
+        warnings.warn(AtWarning("Unstable ring"))
+        tunes = np.empty(nd)
+        tunes[:] = np.nan
     return tunes
 
 
@@ -141,25 +158,28 @@ def _analyze2(mt, ms):
     B = mt[2:, 2:]
     alp0_a, bet0_a, vp_a = _closure(A)
     alp0_b, bet0_b, vp_b = _closure(B)
-    vps = numpy.array([vp_a, vp_b])
-    el0 = (numpy.array([alp0_a, alp0_b]), numpy.array([bet0_a, bet0_b]), 0.0)
+    vps = np.array([vp_a, vp_b])
+    el0 = (np.array([alp0_a, alp0_b]), np.array([bet0_a, bet0_b]), 0.0)
     alpha_a, beta_a, mu_a = _twiss22(ms[:, :2, :2], alp0_a, bet0_a)
     alpha_b, beta_b, mu_b = _twiss22(ms[:, 2:, 2:], alp0_b, bet0_b)
-    els = (numpy.stack((alpha_a, alpha_b), axis=1),
-           numpy.stack((beta_a, beta_b), axis=1),
-           numpy.stack((mu_a, mu_b), axis=1))
+    els = (
+        np.stack((alpha_a, alpha_b), axis=1),
+        np.stack((beta_a, beta_b), axis=1),
+        np.stack((mu_a, mu_b), axis=1),
+    )
     return vps, _DATA2_DTYPE, el0, els, _W24_DTYPE
 
 
 def _analyze4(mt, ms):
     """Analysis of a 4D 1-turn transfer matrix according to Sagan, Rubin"""
+
     def propagate(t12):
         M = t12[:2, :2]
         N = t12[2:, 2:]
         m = t12[:2, 2:]
         n = t12[2:, :2]
         ff = n @ C + g * N
-        gamma = sqrt(numpy.linalg.det(ff))
+        gamma = sqrt(np.linalg.det(ff))
         e12 = (g * M - m @ _S @ C.T @ _S.T) / gamma
         f12 = ff / gamma
         a12 = e12 @ A @ _S @ e12.T @ _S.T
@@ -172,82 +192,89 @@ def propagate(t12):
     m = mt[:2, 2:]
     n = mt[2:, :2]
     H = m + _S @ n.T @ _S.T
-    detH = numpy.linalg.det(H)
+    detH = np.linalg.det(H)
     if detH == 0.0:
         g = 1.0
         C = -H
         A = M
         B = N
     else:
-        t = numpy.trace(M - N)
+        t = np.trace(M - N)
         t2 = t * t
         t2h = t2 + 4.0 * detH
         g2 = (1.0 + sqrt(t2 / t2h)) / 2
         g = sqrt(g2)
-        C = -H * numpy.sign(t) / (g * sqrt(t2h))
-        A = g2 * M - g * (m @ _S @ C.T @ _S.T + C @ n) + \
-            C @ N @ _S @ C.T @ _S.T
-        B = g2 * N + g * (_S @ C.T @ _S.T @ m + n @ C) + \
-            _S @ C.T @ _S.T @ M @ C
+        C = -H * np.sign(t) / (g * sqrt(t2h))
+        A = g2 * M - g * (m @ _S @ C.T @ _S.T + C @ n) + C @ N @ _S @ C.T @ _S.T
+        B = g2 * N + g * (_S @ C.T @ _S.T @ m + n @ C) + _S @ C.T @ _S.T @ M @ C
     alp0_a, bet0_a, vp_a = _closure(A)
     alp0_b, bet0_b, vp_b = _closure(B)
-    vps = numpy.array([vp_a, vp_b])
-    el0 = (numpy.array([alp0_a, alp0_b]),
-           numpy.array([bet0_a, bet0_b]),
-           0.0, g, A, B, C)
+    vps = np.array([vp_a, vp_b])
+    el0 = (np.array([alp0_a, alp0_b]), np.array([bet0_a, bet0_b]), 0.0, g, A, B, C)
     if ms.shape[0] > 0:
         e, f, g, ai, bi, ci = zip(*[propagate(mi) for mi in ms])
-        alp_a, bet_a, mu_a = _twiss22(numpy.array(e), alp0_a, bet0_a)
-        alp_b, bet_b, mu_b = _twiss22(numpy.array(f), alp0_b, bet0_b)
-        els = (numpy.stack((alp_a, alp_b), axis=1),
-               numpy.stack((bet_a, bet_b), axis=1),
-               numpy.stack((mu_a, mu_b), axis=1), numpy.array(g),
-               numpy.stack(ai, axis=0), numpy.stack(bi, axis=0),
-               numpy.stack(ci, axis=0))
+        alp_a, bet_a, mu_a = _twiss22(np.array(e), alp0_a, bet0_a)
+        alp_b, bet_b, mu_b = _twiss22(np.array(f), alp0_b, bet0_b)
+        els = (
+            np.stack((alp_a, alp_b), axis=1),
+            np.stack((bet_a, bet_b), axis=1),
+            np.stack((mu_a, mu_b), axis=1),
+            np.array(g),
+            np.stack(ai, axis=0),
+            np.stack(bi, axis=0),
+            np.stack(ci, axis=0),
+        )
     else:
-        els = (numpy.empty((0, 2)), numpy.empty((0, 2)),
-               numpy.empty((0, 2)), numpy.empty((0,)),
-               numpy.empty((0, 2, 2)), numpy.empty((0, 2, 2)),
-               numpy.empty((0, 2, 2)))
+        els = (
+            np.empty((0, 2)),
+            np.empty((0, 2)),
+            np.empty((0, 2)),
+            np.empty((0,)),
+            np.empty((0, 2, 2)),
+            np.empty((0, 2, 2)),
+            np.empty((0, 2, 2)),
+        )
     return vps, _DATA4_DTYPE, el0, els, _W24_DTYPE
 
 
 def _analyze6(mt, ms):
     """Analysis of a 2D, 4D, 6D 1-turn transfer matrix
     according to Wolski"""
+
     def get_phase(a22):
         """Return the phase for A standardization"""
         return atan2(a22[0, 1], a22[0, 0])
 
     def standardize(aa, slcs):
         """Apply rotation to put A in std form"""
+
         def rot2(slc):
             rot = -get_phase(aa[slc, slc])
             cs = cos(rot)
             sn = sin(rot)
-            return aa[:, slc] @ numpy.array([[cs, sn], [-sn, cs]])
+            return aa[:, slc] @ np.array([[cs, sn], [-sn, cs]])
 
-        return numpy.concatenate([rot2(slc) for slc in slcs], axis=1)
+        return np.concatenate([rot2(slc) for slc in slcs], axis=1)
 
     def r_matrices(ai):
         # Rk = A * S * Ik * inv(A) * S.T
         def mul2(slc):
             return ai[:, slc] @ tt[slc, slc]
 
-        ais = numpy.concatenate([mul2(slc) for slc in slices], axis=1)
+        ais = np.concatenate([mul2(slc) for slc in slices], axis=1)
         invai = solve(ai, ss.T)
-        ri = numpy.array([ais[:, sl] @ invai[sl, :] for sl in slices])
-        mui = numpy.array([get_phase(ai[sl, sl]) for sl in slices])
+        ri = np.array([ais[:, sl] @ invai[sl, :] for sl in slices])
+        mui = np.array([get_phase(ai[sl, sl]) for sl in slices])
         return mui, ri, ai
 
     def propagate4(ri, phi, ai):
-        betai = numpy.stack((ri[..., 0, 0, 0], ri[..., 1, 2, 2]), axis=-1)
-        alphai = -numpy.stack((ri[..., 0, 1, 0], ri[..., 1, 3, 2]), axis=-1)
+        betai = np.stack((ri[..., 0, 0, 0], ri[..., 1, 2, 2]), axis=-1)
+        alphai = -np.stack((ri[..., 0, 1, 0], ri[..., 1, 3, 2]), axis=-1)
         return alphai, betai, phi, ri, ai
 
     def propagate6(ri, phi, ai):
         alphai, betai, phi, ri, ai = propagate4(ri, phi, ai)
-        dispersion = ri[..., 2, :4, 4] / ri[..., 2, 4, 4, numpy.newaxis]
+        dispersion = ri[..., 2, :4, 4] / ri[..., 2, 4, 4, np.newaxis]
         return alphai, betai, phi, ri, ai, dispersion
 
     nv = mt.shape[0]
@@ -270,70 +297,95 @@ def propagate6(ri, phi, ai):
     el0 = propagate(r0, phi0, astd)
     if ms.shape[0] > 0:
         ps, rs, aas = zip(*[r_matrices(mi @ astd) for mi in ms])
-        els = propagate(numpy.array(rs), numpy.array(ps), numpy.array(aas))
+        els = propagate(np.array(rs), np.array(ps), np.array(aas))
     elif dms >= 3:
-        els = (numpy.empty((0, dms)), numpy.empty((0, dms)),
-               numpy.empty((0, dms)),
-               numpy.empty((0, dms, nv, nv)), numpy.empty((0, nv, nv)),
-               numpy.empty((0, 4)))
+        els = (
+            np.empty((0, dms)),
+            np.empty((0, dms)),
+            np.empty((0, dms)),
+            np.empty((0, dms, nv, nv)),
+            np.empty((0, nv, nv)),
+            np.empty((0, 4)),
+        )
     else:
-        els = (numpy.empty((0, dms)), numpy.empty((0, dms)),
-               numpy.empty((0, dms)),
-               numpy.empty((0, dms, nv, nv)), numpy.empty((0, nv, nv)))
+        els = (
+            np.empty((0, dms)),
+            np.empty((0, dms)),
+            np.empty((0, dms)),
+            np.empty((0, dms, nv, nv)),
+            np.empty((0, nv, nv)),
+        )
     return vps, dtype, el0, els, wtype
 
 
 # noinspection PyShadowingNames,PyPep8Naming
-def _linopt(ring: Lattice, analyze, refpts=None, dp=None, dct=None, df=None,
-            orbit=None, twiss_in=None, get_chrom=False, get_w=False,
-            keep_lattice=False, mname='M', add0=(), adds=(), cavpts=None,
-            **kwargs):
+def _linopt(
+    ring: Lattice,
+    analyze,
+    refpts=None,
+    dp=None,
+    dct=None,
+    df=None,
+    orbit=None,
+    twiss_in=None,
+    get_chrom=False,
+    get_w=False,
+    keep_lattice=False,
+    mname="M",
+    add0=(),
+    adds=(),
+    cavpts=None,
+    **kwargs,
+):
     """"""
+
     def build_sigma(orbit, dp=None):
         """Build the initial distribution at entrance of the transfer line"""
         try:
-            d0 = twiss_in['dispersion']
+            d0 = twiss_in["dispersion"]
         except (ValueError, KeyError):  # record arrays throw ValueError !
-            d0 = numpy.zeros((4,))
+            d0 = np.zeros((4,))
 
         try:
-            rmat = twiss_in['R']
+            rmat = twiss_in["R"]
         except (ValueError, KeyError):  # record arrays throw ValueError !
             rmat = None
 
         try:
-            alphas = twiss_in['alpha']
-            betas = twiss_in['beta']
+            alphas = twiss_in["alpha"]
+            betas = twiss_in["beta"]
         except (ValueError, KeyError):  # record arrays throw ValueError !
-            alphas = numpy.zeros((2,))
-            betas = numpy.ones((2,))
+            alphas = np.zeros((2,))
+            betas = np.ones((2,))
 
-        dorbit = numpy.hstack((d0, 1.0, 0.0))
+        dorbit = np.hstack((d0, 1.0, 0.0))
         if orbit is None:
             try:
-                orbit = twiss_in['closed_orbit']
+                orbit = twiss_in["closed_orbit"]
             except (ValueError, KeyError):  # record arrays throw ValueError !
-                orbit = numpy.zeros((6,))
+                orbit = np.zeros((6,))
             if dp is not None:
                 orbit = orbit + dorbit * dp
 
         if dp is not None:
             try:
-                drmat = twiss_in['dR']
+                drmat = twiss_in["dR"]
             except (ValueError, KeyError):  # record arrays throw ValueError !
                 drmat = None
 
             try:
-                dd0 = twiss_in['ddispersion']
-                dalpha = twiss_in['dalpha']
-                dbeta = twiss_in['dbeta']
+                dd0 = twiss_in["ddispersion"]
+                dalpha = twiss_in["dalpha"]
+                dbeta = twiss_in["dbeta"]
             except (ValueError, KeyError):  # record arrays throw ValueError !
-                msg = ("'get_w' option for a line requires 'twiss_in' calculated "
-                       "with 'get_w' activated")
+                msg = (
+                    "'get_w' option for a line requires 'twiss_in' calculated "
+                    "with 'get_w' activated"
+                )
                 raise AtError(msg)
 
-            orbit = orbit + numpy.hstack((dd0, 1.0, 0.0)) * dp * dp
-            dorbit = numpy.hstack((d0+dd0*dp, 1.0, 0.0))
+            orbit = orbit + np.hstack((dd0, 1.0, 0.0)) * dp * dp
+            dorbit = np.hstack((d0 + dd0 * dp, 1.0, 0.0))
 
             if (rmat is not None) and (drmat is not None):
                 rmat = rmat + drmat * dp
@@ -343,23 +395,23 @@ def build_sigma(orbit, dp=None):
 
         if rmat is not None:
             # For some reason, "emittances" must be different...
-            sigm = rmat[0, ...]+10.0*rmat[1, ...]
+            sigm = rmat[0, ...] + 10.0 * rmat[1, ...]
             if rmat.shape[0] >= 3:
-                sigm = sigm+0.1*rmat[2, ...]
-                dorbit = rmat[2, :, 4] / rmat[2, 4, 4, numpy.newaxis]
+                sigm = sigm + 0.1 * rmat[2, ...]
+                dorbit = rmat[2, :, 4] / rmat[2, 4, 4, np.newaxis]
         else:
             slices = [slice(2 * i, 2 * (i + 1)) for i in range(2)]
-            ab = numpy.stack((alphas, betas), axis=1)
-            sigm = numpy.zeros((4, 4))
+            ab = np.stack((alphas, betas), axis=1)
+            sigm = np.zeros((4, 4))
             for slc, (alpha, beta) in zip(slices, ab):
-                gamma = (1.0+alpha*alpha)/beta
-                sigm[slc, slc] = numpy.array([[beta, -alpha], [-alpha, gamma]])
+                gamma = (1.0 + alpha * alpha) / beta
+                sigm[slc, slc] = np.array([[beta, -alpha], [-alpha, gamma]])
 
         return orbit, sigm, dorbit
 
-    def chrom_w(ringup, ringdn, orbitup, orbitdn,
-                refpts=None, **kwargs):
+    def chrom_w(ringup, ringdn, orbitup, orbitdn, refpts=None, **kwargs):
         """Compute the chromaticity and W-functions"""
+
         # noinspection PyShadowingNames
         def off_momentum(rng, orb0, dp=None, **kwargs):
             if twiss_in is None:
@@ -379,34 +431,36 @@ def off_momentum(rng, orb0, dp=None, **kwargs):
                 dpdn = None
             vps, _, el0, els, wtype = analyze(mxx, ms)
             tunes = _tunes(rng, orbit=orb0)
-            o0up, oup = get_orbit(ring, refpts=refpts, guess=orb0, dp=dpup,
-                                  orbit=o0up, **kwargs)
-            o0dn, odn = get_orbit(ring, refpts=refpts, guess=orb0, dp=dpdn,
-                                  orbit=o0dn, **kwargs)
+            o0up, oup = get_orbit(
+                ring, refpts=refpts, guess=orb0, dp=dpup, orbit=o0up, **kwargs
+            )
+            o0dn, odn = get_orbit(
+                ring, refpts=refpts, guess=orb0, dp=dpdn, orbit=o0dn, **kwargs
+            )
             d0 = (o0up - o0dn)[:4] / dp_step
-            ds = numpy.array([(up - dn)[:4] / dp_step for up, dn in zip(oup, odn)])
+            ds = np.array([(up - dn)[:4] / dp_step for up, dn in zip(oup, odn)])
             return tunes, el0, els, d0, ds, wtype
 
         def wget(ddp, elup, eldn, has_r):
             """Compute the chromatic amplitude function"""
-            *data_up, = elup  # Extract alpha and beta
-            *data_dn, = eldn
+            (*data_up,) = elup  # Extract alpha and beta
+            (*data_dn,) = eldn
             alpha_up, beta_up, mu_up = data_up[:3]
             alpha_dn, beta_dn, mu_dn = data_dn[:3]
-            db = numpy.array(beta_up - beta_dn) / ddp
+            db = np.array(beta_up - beta_dn) / ddp
             mb = (beta_up + beta_dn) / 2
-            da = numpy.array(alpha_up - alpha_dn) / ddp
+            da = np.array(alpha_up - alpha_dn) / ddp
             ma = (alpha_up + alpha_dn) / 2
             wa = da - ma / mb * db
             wb = db / mb
-            ww = numpy.sqrt(wa ** 2 + wb ** 2)
-            wp = numpy.arctan2(wa, wb)
-            dmu = numpy.array(mu_up - mu_dn) / ddp
+            ww = np.sqrt(wa**2 + wb**2)
+            wp = np.arctan2(wa, wb)
+            dmu = np.array(mu_up - mu_dn) / ddp
             data_out = (ww, wp, da, db, dmu)
             if has_r:
                 r_up = data_up[3]
                 r_dn = data_dn[3]
-                data_out += (numpy.array(r_up - r_dn) / ddp, )
+                data_out += (np.array(r_up - r_dn) / ddp,)
             return data_out
 
         deltap = orbitup[4] - orbitdn[4]
@@ -416,9 +470,9 @@ def wget(ddp, elup, eldn, has_r):
         tunesdn, el0dn, elsdn, d0dn, dsdn, _ = data_dn
         has_r = len(wtype) == 7
         # in 6D, dp comes out of find_orbit6
-        chrom = (tunesup-tunesdn) / deltap
-        dd0 = numpy.array(d0up - d0dn) / deltap
-        dds = numpy.array(dsup - dsdn) / deltap
+        chrom = (tunesup - tunesdn) / deltap
+        dd0 = np.array(d0up - d0dn) / deltap
+        dds = np.array(dsup - dsdn) / deltap
         data0 = wget(deltap, el0up, el0dn, has_r)
         datas = wget(deltap, elsup, elsdn, has_r)
         data0 = data0 + (dd0,)
@@ -427,13 +481,12 @@ def wget(ddp, elup, eldn, has_r):
 
     def unwrap(mu):
         """Remove the phase jumps"""
-        dmu = numpy.diff(numpy.concatenate((numpy.zeros((1, dms)),
-                                            mu)), axis=0)
-        jumps = dmu < -1.e-3
-        mu += numpy.cumsum(jumps, axis=0) * 2.0 * numpy.pi
+        dmu = np.diff(np.concatenate((np.zeros((1, dms)), mu)), axis=0)
+        jumps = dmu < -1.0e-3
+        mu += np.cumsum(jumps, axis=0) * 2.0 * np.pi
 
-    dp_step = kwargs.get('DPStep', DConstant.DPStep)
-    addtype = kwargs.pop('addtype', [])
+    dp_step = kwargs.get("DPStep", DConstant.DPStep)
+    addtype = kwargs.pop("addtype", [])
 
     if ring.is_6d:
         get_matrix = find_m66
@@ -444,15 +497,16 @@ def unwrap(mu):
 
     o0up = None
     o0dn = None
-    if twiss_in is None:   # Ring
+    if twiss_in is None:  # Ring
         if orbit is None:
-            orbit, _ = get_orbit(ring, dp=dp, dct=dct, df=df,
-                                 keep_lattice=keep_lattice, **kwargs)
+            orbit, _ = get_orbit(
+                ring, dp=dp, dct=dct, df=df, keep_lattice=keep_lattice, **kwargs
+            )
             keep_lattice = True
         # Get 1-turn transfer matrix
         mt, ms = get_matrix(ring, refpts=refpts, orbit=orbit, **kwargs)
         mxx = mt
-    else:                       # Transfer line
+    else:  # Transfer line
         orbit, sigma, dorbit = build_sigma(orbit)
         # Get 1-turn transfer matrix
         mt, ms = get_matrix(ring, refpts=refpts, orbit=orbit, **kwargs)
@@ -482,32 +536,38 @@ def unwrap(mu):
 
     nrefs = orbs.shape[0]
     dms = vps.size
-    if dms >= 3:            # 6D processing
-        dtype = dtype + [('closed_orbit', numpy.float64, (6,)),
-                         ('M', numpy.float64, (2*dms, 2*dms)),
-                         ('s_pos', numpy.float64)]
-        data0 = (orb0, numpy.identity(2*dms), 0.0)
+    if dms >= 3:  # 6D processing
+        dtype = dtype + [
+            ("closed_orbit", np.float64, (6,)),
+            ("M", np.float64, (2 * dms, 2 * dms)),
+            ("s_pos", np.float64),
+        ]
+        data0 = (orb0, np.identity(2 * dms), 0.0)
         datas = (orbs, ms, spos)
         length = ring.get_s_pos(len(ring))[0]
-        damping_rates = -numpy.log(numpy.absolute(vps))
+        damping_rates = -np.log(np.absolute(vps))
         damping_times = length / clight / damping_rates
-    else:               # 4D processing
-        kwargs['keep_lattice'] = True
+    else:  # 4D processing
+        kwargs["keep_lattice"] = True
         dpup = orb0[4] + 0.5 * dp_step
         dpdn = orb0[4] - 0.5 * dp_step
-        o0up, oup = get_orbit(ring, refpts=refpts, guess=orb0, dp=dpup,
-                              orbit=o0up, **kwargs)
-        o0dn, odn = get_orbit(ring, refpts=refpts, guess=orb0, dp=dpdn,
-                              orbit=o0dn, **kwargs)
+        o0up, oup = get_orbit(
+            ring, refpts=refpts, guess=orb0, dp=dpup, orbit=o0up, **kwargs
+        )
+        o0dn, odn = get_orbit(
+            ring, refpts=refpts, guess=orb0, dp=dpdn, orbit=o0dn, **kwargs
+        )
         d0 = (o0up - o0dn)[:4] / dp_step
-        ds = numpy.array([(up - dn)[:4] / dp_step for up, dn in zip(oup, odn)])
-        dtype = dtype + [('dispersion', numpy.float64, (4,)),
-                         ('closed_orbit', numpy.float64, (6,)),
-                         (mname, numpy.float64, (2*dms, 2*dms)),
-                         ('s_pos', numpy.float64)]
+        ds = np.array([(up - dn)[:4] / dp_step for up, dn in zip(oup, odn)])
+        dtype = dtype + [
+            ("dispersion", np.float64, (4,)),
+            ("closed_orbit", np.float64, (6,)),
+            (mname, np.float64, (2 * dms, 2 * dms)),
+            ("s_pos", np.float64),
+        ]
         data0 = (d0, orb0, mt, get_s_pos(ring, len(ring))[0])
         datas = (ds, orbs, ms, spos)
-        damping_times = numpy.nan
+        damping_times = np.nan
 
     if get_w:
         dtype = dtype + wtype
@@ -520,19 +580,22 @@ def unwrap(mu):
         deltap = o0up[4] - o0dn[4]
         chrom = (tunesup - tunesdn) / deltap
     else:
-        chrom = numpy.nan
+        chrom = np.nan
 
-    beamdata = numpy.array((tunes, chrom, damping_times),
-                           dtype=[('tune', numpy.float64, (dms,)),
-                                  ('chromaticity', numpy.float64, (dms,)),
-                                  ('damping_time', numpy.float64, (dms,))
-                                  ]).view(numpy.recarray)
+    beamdata = np.array(
+        (tunes, chrom, damping_times),
+        dtype=[
+            ("tune", np.float64, (dms,)),
+            ("chromaticity", np.float64, (dms,)),
+            ("damping_time", np.float64, (dms,)),
+        ],
+    ).view(np.recarray)
 
     dtype = dtype + addtype
-    elemdata0 = numpy.array(el0+data0+add0, dtype=dtype).view(numpy.recarray)
-    elemdata = numpy.recarray((nrefs,), dtype=dtype)
+    elemdata0 = np.array(el0 + data0 + add0, dtype=dtype).view(np.recarray)
+    elemdata = np.recarray((nrefs,), dtype=dtype)
     if nrefs > 0:
-        for name, value in zip(numpy.dtype(dtype).names, els+datas+adds):
+        for name, value in zip(np.dtype(dtype).names, els + datas + adds):
             elemdata[name] = value
         unwrap(elemdata.mu)
     return elemdata0, beamdata, elemdata
@@ -542,6 +605,9 @@ def unwrap(mu):
 def linopt2(ring: Lattice, *args, **kwargs):
     r"""Linear analysis of an uncoupled lattice
 
+    :py:func:`linopt2` computes the linear optics parameters on a single cell
+    (*periodicity* is not taken into account).
+
     Parameters:
         ring:   Lattice description.
 
@@ -653,6 +719,9 @@ def linopt4(ring: Lattice, *args, **kwargs):
 
     4D-analysis of coupled motion following Sagan/Rubin [7]_
 
+    :py:func:`linopt4` computes the linear optics parameters on a single cell
+    (*periodicity* is not taken into account).
+
     Parameters:
         ring:   Lattice description.
 
@@ -774,6 +843,9 @@ def linopt6(ring: Lattice, *args, **kwargs):
     *  a field **R**, as provided by ATLINOPT6, or
     * the fields **beta** and **alpha**, as provided by linopt and linopt6
 
+    :py:func:`linopt6` computes the linear optics parameters on a single cell
+    (*periodicity* is not taken into account).
+
     Parameters:
         ring:   Lattice description.
 
@@ -907,18 +979,24 @@ def linopt_auto(ring: Lattice, *args, **kwargs):
         The output varies depending whether :py:func:`.linopt2` or
         :py:func:`.linopt6` is called. To be used with care!
     """
-    if not (kwargs.pop('coupled', True) or ring.is_6d):
+    if not (kwargs.pop("coupled", True) or ring.is_6d):
         return linopt2(ring, *args, **kwargs)
     else:
         return linopt6(ring, *args, **kwargs)
 
 
-def get_optics(ring: Lattice, refpts: Refpts = None,
-               dp: float = None,
-               method: Callable = linopt6,
-               **kwargs):
+def get_optics(
+    ring: Lattice,
+    refpts: Refpts = None,
+    dp: float = None,
+    method: Callable = linopt6,
+    **kwargs,
+):
     """Linear analysis of a fully coupled lattice
 
+    :py:func:`get_optics` computes the linear optics parameters on a single cell
+    (*periodicity* is not taken into account).
+
     Parameters:
         ring:                   Lattice description.
         refpts:                 Elements at which data is returned.
@@ -1003,8 +1081,13 @@ def get_optics(ring: Lattice, refpts: Refpts = None,
 
 # noinspection PyPep8Naming
 @check_6d(False)
-def linopt(ring: Lattice, dp: float = 0.0, refpts: Refpts = None,
-           get_chrom: bool = False, **kwargs):
+def linopt(
+    ring: Lattice,
+    dp: float = 0.0,
+    refpts: Refpts = None,
+    get_chrom: bool = False,
+    **kwargs,
+):
     """Linear analysis of a H/V coupled lattice (deprecated)
 
     Parameters:
@@ -1078,11 +1161,19 @@ def linopt(ring: Lattice, dp: float = 0.0, refpts: Refpts = None,
 
     :meta private:
     """
-    analyze = _analyze4 if kwargs.pop('coupled', True) else _analyze2
-    eld0, bd, eld = _linopt(ring, analyze, refpts, dp=dp, get_chrom=get_chrom,
-                            add0=(0,), adds=(get_uint32_index(ring, refpts),),
-                            addtype=[('idx', numpy.uint32)],
-                            mname='m44', **kwargs)
+    analyze = _analyze4 if kwargs.pop("coupled", True) else _analyze2
+    eld0, bd, eld = _linopt(
+        ring,
+        analyze,
+        refpts,
+        dp=dp,
+        get_chrom=get_chrom,
+        add0=(0,),
+        adds=(get_uint32_index(ring, refpts),),
+        addtype=[("idx", np.uint32)],
+        mname="m44",
+        **kwargs,
+    )
     return eld0, bd.tune, bd.chromaticity, eld
 
 
@@ -1094,6 +1185,10 @@ def avlinopt(ring: Lattice, dp: float = 0.0, refpts: Refpts = None, **kwargs):
     :py:func:`avlinopt` returns average beta, mu, dispersion over the lattice
     elements.
 
+    :py:func:`avlinopt` returns average beta, mu, dispersion over the lattice
+    elements. :py:func:`avlinopt` acts on a single cell (*periodicity* is not taken
+    into account).
+
     Parameters:
         ring:       Lattice description.
         dp:         Momentum deviation.
@@ -1323,11 +1418,20 @@ def displong(disp0, ir, k2, lg):
 
 
 @frequency_control
-def get_tune(ring: Lattice, *, method: str = 'linopt',
-             dp: float = None, dct: float = None, df: float = None,
-             orbit: Orbit = None, **kwargs):
+def get_tune(
+    ring: Lattice,
+    *,
+    method: str = "linopt",
+    dp: float = None,
+    dct: float = None,
+    df: float = None,
+    orbit: Orbit = None,
+    **kwargs,
+):
     r"""Computes the tunes using several available methods
 
+    :py:func:`get_tune` may use several methods depending on a *method* keyword.
+
     Parameters:
         ring:                   Lattice description
         method:                 ``'linopt'`` returns the tunes from the
@@ -1362,6 +1466,7 @@ def get_tune(ring: Lattice, *, method: str = 'linopt',
     Returns:
         tunes (ndarray):        array([:math:`\nu_x,\nu_y`])
     """
+
     # noinspection PyShadowingNames
     def gen_centroid(ring, ampl, nturns, remove_dc, ld):
         nv = ld.A.shape[0]
@@ -1370,38 +1475,50 @@ def gen_centroid(ring, ampl, nturns, remove_dc, ld):
         p0[2] += ampl
         if nv >= 6:
             p0[4] += ampl
-        p1 = numpy.squeeze(internal_lpass(ring, p0, nturns, len(ring)))
+        p1 = np.squeeze(internal_lpass(ring, p0, nturns, len(ring)))
         if remove_dc:
-            p1 -= numpy.mean(p1, axis=1, keepdims=True)
+            p1 -= np.mean(p1, axis=1, keepdims=True)
         p2 = solve(ld.A, p1[:nv, :])
-        return numpy.conjugate(p2.T.view(dtype=complex).T)
+        return np.conjugate(p2.T.view(dtype=complex).T)
 
     get_integer = kwargs.pop("get_integer", False)
     if get_integer:
-        assert method == 'linopt', 'Integer tune only accessible with method=linopt'
-    if method == 'linopt':
+        assert method == "linopt", "Integer tune only accessible with method=linopt"
+    if method == "linopt":
         if get_integer:
-            _, _, c = get_optics(ring, refpts=range(len(ring)+1),
-                                 dp=dp, dct=dct, df=df, orbit=orbit)
-            tunes = c.mu[-1]/(2*numpy.pi)
+            _, _, c = get_optics(
+                ring, refpts=range(len(ring) + 1), dp=dp, dct=dct, df=df, orbit=orbit
+            )
+            tunes = c.mu[-1] / (2 * np.pi) * ring.periodicity
         else:
             tunes = _tunes(ring, dp=dp, dct=dct, df=df, orbit=orbit)
+            tunes, _ = np.modf(tunes * ring.periodicity)
     else:
-        nturns = kwargs.pop('nturns', 512)
-        ampl = kwargs.pop('ampl', 1.0e-6)
-        remove_dc = kwargs.pop('remove_dc', True)
+        nturns = kwargs.pop("nturns", 512)
+        ampl = kwargs.pop("ampl", 1.0e-6)
+        remove_dc = kwargs.pop("remove_dc", True)
         ld, _, _ = linopt6(ring, dp=dp, dct=dct, df=df, orbit=orbit)
         cents = gen_centroid(ring, ampl, nturns, remove_dc, ld)
         tunes = get_tunes_harmonic(cents, method=method, **kwargs)
+        tunes, _ = np.modf(tunes * ring.periodicity)
     return tunes
 
 
 @frequency_control
-def get_chrom(ring: Lattice, *, method: str = 'linopt',
-              dp: float = None, dct: float = None, df: float = None,
-              cavpts: Refpts = None, **kwargs):
+def get_chrom(
+    ring: Lattice,
+    *,
+    method: str = "linopt",
+    dp: float = None,
+    dct: float = None,
+    df: float = None,
+    cavpts: Refpts = None,
+    **kwargs,
+):
     r"""Computes the chromaticities using several available methods
 
+    :py:func:`get_tune` may use several methods depending on a *method* keyword.
+
     Parameters:
         ring:               Lattice description.
         method:             ``'linopt'`` returns the tunes from the
@@ -1441,9 +1558,9 @@ def get_chrom(ring: Lattice, *, method: str = 'linopt',
         chromaticities (ndarray):   array([:math:`\xi_x,\xi_y`])
     """
 
-    dp_step = kwargs.pop('DPStep', DConstant.DPStep)
-    if method == 'fft':
-        print('Warning fft method not accurate to get the chromaticity')
+    dp_step = kwargs.pop("DPStep", DConstant.DPStep)
+    if method == "fft":
+        print("Warning fft method not accurate to get the chromaticity")
 
     if ring.is_6d:
         f0 = ring.get_rf_frequency(cavpts=cavpts)
diff --git a/pyat/at/physics/nonlinear.py b/pyat/at/physics/nonlinear.py
index 7599fb237..c7036a266 100644
--- a/pyat/at/physics/nonlinear.py
+++ b/pyat/at/physics/nonlinear.py
@@ -1,7 +1,8 @@
 """
 Non-linear optics
 """
-import numpy
+
+import numpy as np
 from typing import Optional, Sequence
 from scipy.special import factorial
 from ..lattice import Element, Lattice
@@ -10,12 +11,16 @@
 from .linear import get_tune, get_chrom, linopt6
 from .harmonic_analysis import get_tunes_harmonic
 
-__all__ = ['detuning', 'chromaticity', 'gen_detuning_elem', 'tunes_vs_amp']
+__all__ = ["detuning", "chromaticity", "gen_detuning_elem", "tunes_vs_amp"]
 
 
-def tunes_vs_amp(ring: Lattice, amp: Optional[Sequence[float]] = None,
-                 dim: Optional[int] = 0, nturns: Optional[int] = 512,
-                 **kwargs) -> numpy.ndarray:
+def tunes_vs_amp(
+    ring: Lattice,
+    amp: Optional[Sequence[float]] = None,
+    dim: Optional[int] = 0,
+    nturns: Optional[int] = 512,
+    **kwargs,
+) -> np.ndarray:
     r"""Generates a range of tunes for varying x, y, or z amplitudes
 
     Parameters:
@@ -37,49 +42,57 @@ def tunes_vs_amp(ring: Lattice, amp: Optional[Sequence[float]] = None,
     """
 
     def _gen_part(ring, amp, dim, orbit, ld, nturns):
-        invx = numpy.array([[1/numpy.sqrt(ld['beta'][0]), 0],
-                            [ld['alpha'][0]/numpy.sqrt(ld['beta'][0]),
-                            numpy.sqrt(ld['beta'][0])]])
-
-        invy = numpy.array([[1/numpy.sqrt(ld['beta'][1]), 0],
-                            [ld['alpha'][1]/numpy.sqrt(ld['beta'][1]),
-                            numpy.sqrt(ld['beta'][1])]])
-        part = numpy.array([orbit, ] * len(amp)).T + 1.0e-6
+        invx = np.array(
+            [
+                [1 / np.sqrt(ld["beta"][0]), 0],
+                [ld["alpha"][0] / np.sqrt(ld["beta"][0]), np.sqrt(ld["beta"][0])],
+            ]
+        )
+
+        invy = np.array(
+            [
+                [1 / np.sqrt(ld["beta"][1]), 0],
+                [ld["alpha"][1] / np.sqrt(ld["beta"][1]), np.sqrt(ld["beta"][1])],
+            ]
+        )
+        part = np.array([orbit] * len(amp)).T + 1.0e-6
         part[dim, :] += amp
         part = internal_lpass(ring, part, nturns=nturns)
         sh = part.shape
-        partx = numpy.reshape(part[0, :], (sh[1], sh[3]))
-        partxp = numpy.reshape(part[1, :], (sh[1], sh[3]))
-        party = numpy.reshape(part[2, :], (sh[1], sh[3]))
-        partyp = numpy.reshape(part[3, :], (sh[1], sh[3]))
-        px = numpy.array([numpy.matmul(invx, [partx[i], partxp[i]])
-                          for i in range(len(amp))])
-        py = numpy.array([numpy.matmul(invy, [party[i], partyp[i]])
-                          for i in range(len(amp))])
-        return px[:, 0, :] - 1j*px[:, 1, :], py[:, 0, :] - 1j*py[:, 1, :]
+        partx = np.reshape(part[0, :], (sh[1], sh[3]))
+        partxp = np.reshape(part[1, :], (sh[1], sh[3]))
+        party = np.reshape(part[2, :], (sh[1], sh[3]))
+        partyp = np.reshape(part[3, :], (sh[1], sh[3]))
+        px = np.array([np.matmul(invx, [partx[i], partxp[i]]) for i in range(len(amp))])
+        py = np.array([np.matmul(invy, [party[i], partyp[i]]) for i in range(len(amp))])
+        return px[:, 0, :] - 1j * px[:, 1, :], py[:, 0, :] - 1j * py[:, 1, :]
 
     l0, bd, _ = linopt6(ring)
-    orbit = l0['closed_orbit']
-    tunes = bd['tune']
+    orbit = l0["closed_orbit"]
+    tunes = bd["tune"]
 
     if amp is not None:
         partx, party = _gen_part(ring, amp, dim, orbit, l0, nturns)
         qx = get_tunes_harmonic(partx, **kwargs)
         qy = get_tunes_harmonic(party, **kwargs)
-        tunes = numpy.vstack((qx, qy)).T
+        tunes = np.vstack((qx, qy)).T
 
-    return numpy.array(tunes)
+    return np.array(tunes)
 
 
-def detuning(ring: Lattice,
-             xm: Optional[float] = 0.3e-4, ym: Optional[float] = 0.3e-4,
-             npoints: Optional[int] = 3,
-             nturns: Optional[int] = 512, **kwargs):
+def detuning(
+    ring: Lattice,
+    xm: Optional[float] = 0.3e-4,
+    ym: Optional[float] = 0.3e-4,
+    npoints: Optional[int] = 3,
+    nturns: Optional[int] = 512,
+    **kwargs,
+):
     """Computes the tunes for a sequence of amplitudes
 
     This function uses :py:func:`tunes_vs_amp` to compute the tunes for
     the specified amplitudes. Then it fits this data and returns
-    the detuning coefficiant dQx/dx, dQy/dx, dQx/dy, dQy/dy and the
+    the detuning coefficiant dQx/Jx, dQy/Jx, dQx/Jy, dQy/Jy and the
     qx, qy arrays versus x, y arrays
 
     Parameters:
@@ -97,44 +110,56 @@ def detuning(ring: Lattice,
 
     Returns:
         q0 (ndarray): qx, qy from horizontal and vertical amplitude scans.
-                      Fit results with shape (2, 2)
+                      Fit results with shape (2, 2). Only the fractional part
+                      is returned
         q1 (ndarray): Amplitude detuning coefficients with shape (2, 2)
-                      ``[[dQx/dx, dQy/dx], [dQx/dy, dQy/dy]]``
+                      ``[[dQx/dJx, dQy/dJx], [dQx/dJy, dQy/dJy]]``
         x (ndarray): x amplitudes (npoints, )
         q_dx (ndarray): qx, qy tunes as a function of x amplitude (npoints, 2)
+                        Only the fractional parts are returned
         y (ndarray): y amplitudes (npoints, )
         q_dy (ndarray): qx, qy tunes as a function of y amplitude (npoints, 2)
+                        Only the fractional parts are returned
     """
     lindata0, _, _ = linopt6(ring)
     gamma = (1 + lindata0.alpha * lindata0.alpha) / lindata0.beta
 
-    x = numpy.linspace(-xm, xm, npoints)
-    y = numpy.linspace(-ym, ym, npoints)
+    x = np.linspace(-xm, xm, npoints)
+    y = np.linspace(-ym, ym, npoints)
     x2 = x * x
     y2 = y * y
 
     q_dx = tunes_vs_amp(ring, amp=x, dim=0, nturns=nturns, **kwargs)
     q_dy = tunes_vs_amp(ring, amp=y, dim=2, nturns=nturns, **kwargs)
-
-    idx = numpy.isfinite(q_dx[:, 0]) & numpy.isfinite(q_dx[:, 1])
-    idy = numpy.isfinite(q_dy[:, 0]) & numpy.isfinite(q_dy[:, 1])
-
-    fx = numpy.polyfit(x2[idx], q_dx[idx], 1)
-    fy = numpy.polyfit(y2[idy], q_dy[idy], 1)
-
-    q0 = [[fx[1, 0], fx[1, 1]],
-          [fy[1, 0], fy[1, 1]]]
-    q1 = [[2 * fx[0, 0] / gamma[0], 2 * fx[0, 1] / gamma[0]],
-          [2 * fy[0, 0] / gamma[1], 2 * fy[0, 1] / gamma[1]]]
-
-    return numpy.array(q0), numpy.array(q1), x, q_dx, y, q_dy
-
-
-def chromaticity(ring: Lattice, method: Optional[str] = 'linopt',
-                 dpm: Optional[float] = 0.02, npoints: Optional[int] = 11,
-                 order: Optional[int] = 3,
-                 dp: Optional[float] = 0,
-                 **kwargs):
+    q_dx, _ = np.modf(q_dx * ring.periodicity)
+    q_dy, _ = np.modf(q_dy * ring.periodicity)
+
+    idx = np.isfinite(q_dx[:, 0]) & np.isfinite(q_dx[:, 1])
+    idy = np.isfinite(q_dy[:, 0]) & np.isfinite(q_dy[:, 1])
+
+    fx = np.polyfit(x2[idx], q_dx[idx], 1)
+    fy = np.polyfit(y2[idy], q_dy[idy], 1)
+
+    q0 = np.array([[fx[1, 0], fx[1, 1]], [fy[1, 0], fy[1, 1]]])
+    q1 = np.array(
+        [
+            [2 * fx[0, 0] / gamma[0], 2 * fx[0, 1] / gamma[0]],
+            [2 * fy[0, 0] / gamma[1], 2 * fy[0, 1] / gamma[1]],
+        ]
+    )
+
+    return q0, q1, x, q_dx, y, q_dy
+
+
+def chromaticity(
+    ring: Lattice,
+    method: Optional[str] = "linopt",
+    dpm: Optional[float] = 0.02,
+    npoints: Optional[int] = 11,
+    order: Optional[int] = 3,
+    dp: Optional[float] = 0,
+    **kwargs,
+):
     r"""Computes the non-linear chromaticities
 
     This function computes the tunes for the specified momentum offsets.
@@ -157,24 +182,25 @@ def chromaticity(ring: Lattice, method: Optional[str] = 'linopt',
 
     Returns:
         fit (ndarray): Horizontal Vertical polynomial coefficients
-                       from ``numpy.polyfit`` of shape (2, order+1)
+                       from ``np.polyfit`` of shape (2, order+1)
         dpa: dp array of shape (npoints,)
         qz: tune array of shape (npoints, 2)
     """
     if order == 0:
         return get_chrom(ring, dp=dp, **kwargs)
     elif order > npoints - 1:
-        raise ValueError('order should be smaller than npoints-1')
+        raise ValueError("order should be smaller than npoints-1")
     else:
-        dpa = numpy.linspace(-dpm, dpm, npoints)
+        dpa = np.linspace(-dpm, dpm, npoints)
         qz = []
         for dpi in dpa:
-            qz.append(get_tune(ring, method=method, dp=dp+dpi, remove_dc=True,
-                      **kwargs))
-        fit = numpy.polyfit(dpa, qz, order)[::-1]
-        fitx = fit[:, 0]/factorial(numpy.arange(order + 1))
-        fity = fit[:, 1]/factorial(numpy.arange(order + 1))
-        return numpy.array([fitx, fity]), dpa, numpy.array(qz)
+            qz.append(
+                get_tune(ring, method=method, dp=dp + dpi, remove_dc=True, **kwargs)
+            )
+        fit = np.polyfit(dpa, qz, order)[::-1]
+        fitx = fit[:, 0] / factorial(np.arange(order + 1))
+        fity = fit[:, 1] / factorial(np.arange(order + 1))
+        return np.array([fitx, fity]), dpa, np.array(qz)
 
 
 def gen_detuning_elem(ring: Lattice, orbit: Optional[Orbit] = None) -> Element:
@@ -193,12 +219,22 @@ def gen_detuning_elem(ring: Lattice, orbit: Optional[Orbit] = None) -> Element:
         orbit, _ = find_orbit(ring)
     lindata0, _, _ = ring.linopt6(get_chrom=False, orbit=orbit)
     xsi = get_chrom(ring.radiation_off(copy=True))
-    r0, r1, x, q_dx, y, q_dy = detuning(ring.radiation_off(copy=True),
-                                        xm=1.0e-4, ym=1.0e-4, npoints=3)
-    nonlin_elem = Element('NonLinear', PassMethod='DeltaQPass',
-                          Betax=lindata0.beta[0], Betay=lindata0.beta[1],
-                          Alphax=lindata0.alpha[0], Alphay=lindata0.alpha[1],
-                          Qpx=xsi[0], Qpy=xsi[1],
-                          A1=r1[0][0], A2=r1[0][1], A3=r1[1][1],
-                          T1=-orbit, T2=orbit)
+    r0, r1, x, q_dx, y, q_dy = detuning(
+        ring.radiation_off(copy=True), xm=1.0e-4, ym=1.0e-4, npoints=3
+    )
+    nonlin_elem = Element(
+        "NonLinear",
+        PassMethod="DeltaQPass",
+        Betax=lindata0.beta[0],
+        Betay=lindata0.beta[1],
+        Alphax=lindata0.alpha[0],
+        Alphay=lindata0.alpha[1],
+        Qpx=xsi[0],
+        Qpy=xsi[1],
+        A1=r1[0][0],
+        A2=r1[0][1],
+        A3=r1[1][1],
+        T1=-orbit,
+        T2=orbit,
+    )
     return nonlin_elem


Test patch:
diff --git a/atmat/attests/pytests.m b/atmat/attests/pytests.m
index 4b6531ce1..d087c64e4 100644
--- a/atmat/attests/pytests.m
+++ b/atmat/attests/pytests.m
@@ -164,23 +164,25 @@ function offmomdf(testCase, dp)
         function tunechrom4(testCase,lat,dp)
             % test on and off-momentum tunes of 4D lattices
             lattice=testCase.ring4.(lat);
+            periodicity = atGetRingProperties(lattice.m,'Periodicity');
             [mtune,mchrom]=tunechrom(lattice.m,'get_chrom',dp=dp);
             ptune=double(lattice.p.get_tune(pyargs(dp=dp)));
             pchrom=double(lattice.p.get_chrom(pyargs(dp=dp)));
-            testCase.verifyEqual(mtune,ptune,AbsTol=2.e-9);
-            testCase.verifyEqual(mchrom,pchrom,AbsTol=3.e-4);
+            testCase.verifyEqual(mod(mtune*periodicity,1),ptune,AbsTol=2.e-9);
+            testCase.verifyEqual(mchrom*periodicity,pchrom,AbsTol=3.e-4);
         end
 
         function tunechrom6(testCase,lat2,dp)
             % test on and off-momentum tunes of 6D lattices
             lattice=testCase.ring6.(lat2);
+            periodicity = atGetRingProperties(lattice.m,'Periodicity');
             mlat=atsetcavity(lattice.m,frequency='nominal',dp=dp);
             plat=lattice.p.set_rf_frequency(pyargs(dp=dp,copy=true));
             [mtune,mchrom]=tunechrom(mlat,'get_chrom');
             ptune=double(plat.get_tune());
             pchrom=double(plat.get_chrom());
-            testCase.verifyEqual(mtune,ptune,AbsTol=1.e-9);
-            testCase.verifyEqual(mchrom,pchrom,AbsTol=2.e-4);
+            testCase.verifyEqual(mod(mtune*periodicity,1),ptune,AbsTol=1.e-9);
+            testCase.verifyEqual(mchrom*periodicity,pchrom,AbsTol=2.e-4);
         end
 
         function linopt1(testCase,dp)
diff --git a/pyat/test/test_physics.py b/pyat/test/test_physics.py
index 8d98710d3..cd2519acf 100644
--- a/pyat/test/test_physics.py
+++ b/pyat/test/test_physics.py
@@ -264,10 +264,11 @@ def test_get_tune_chrom(hmba_lattice):
     qplin = hmba_lattice.get_chrom()
     qharm = hmba_lattice.get_tune(method='interp_fft')
     qpharm = hmba_lattice.get_chrom(method='interp_fft')
-    assert_close(qlin, [0.38156245, 0.85437541], rtol=1e-8)
-    assert_close(qharm, [0.38156245, 0.85437541], rtol=1e-8)
-    assert_close(qplin, [0.1791909, 0.12242558], rtol=1e-5)
-    assert_close(qpharm, [0.17919145, 0.12242622], rtol=1e-5)
+    print(qlin, qharm)
+    assert_close(qlin, [0.2099983, 0.34001317],atol=1e-8)
+    assert_close(qharm, [0.20999833, 0.34001324], atol=1e-8)
+    assert_close(qplin, [5.734099, 3.917612], atol=1e-8)
+    assert_close(qpharm, [5.734123, 3.917639], atol=1e-8)
 
 
 def test_nl_detuning_chromaticity(hmba_lattice):
@@ -276,19 +277,36 @@ def test_nl_detuning_chromaticity(hmba_lattice):
                                                method='interp_fft', npoints=11)
     q0, q1, _, _, _, _ = at.nonlinear.detuning(hmba_lattice,
                                                npoints=11)
-    assert_close(nlqplin, [[0.38156741, 0.17908231, 1.18656034, -16.47368694],
-                           [0.85437409, 0.1224062, 2.01744075, -3.06407764]],
+    assert_close(nlqplin, [[0.2101570,  5.730634,  37.96993, -527.1578],
+                           [0.3399707,  3.916998,  64.55810, -98.05048]],
                  atol=1e-12, rtol=1e-5)
-    assert_close(nlqpharm, [[0.38156741, 0.17908228, 1.18656178, -16.47370342],
-                            [0.85437409, 0.12240619, 2.01744051, -3.06407046]],
+    assert_close(nlqpharm, [[0.2101570,  5.730630,  37.96992, -527.1587],
+                            [0.3399708,  3.916997,  64.55809, -98.05020]],
                  atol=1e-12, rtol=1e-5)
-    assert_close(q0, [[0.38156263, 0.85437553], [0.38156263, 0.85437553]],
+    assert_close(q0, [[0.210004, 0.340017], [0.210004, 0.340017]],
                  atol=1e-12, rtol=1e-5)
-    assert_close(q1, [[3005.74776344, -3256.81838517],
-                      [-3258.24669916,  1615.13729938]],
+    assert_close(q1, [[96183.925683, -104218.18371],
+                      [-104263.908197,   51684.400417]],
                  atol=1e-12, rtol=1e-5)
 
 
+def test_periodicity(hmba_lattice):
+    q32 = at.linear.get_tune(hmba_lattice)
+    qp32 = at.linear.get_chrom(hmba_lattice)
+    dq32, *_ = at.nonlinear.detuning(hmba_lattice)
+    dqp32, *_ = at.nonlinear.chromaticity(hmba_lattice, npoints=11)
+    hmba_lattice.periodicity = 1
+    hmba_lattice = hmba_lattice.repeat(32)
+    q1 = at.linear.get_tune(hmba_lattice)
+    qp1 = at.linear.get_chrom(hmba_lattice)
+    dq1, *_ = at.nonlinear.detuning(hmba_lattice)
+    dqp1, *_ = at.nonlinear.chromaticity(hmba_lattice, npoints=11)
+    assert_close(q32, q1, atol=1e-12)
+    assert_close(qp32, qp1, atol=1e-12)
+    assert_close(dq32, dq1, atol=1e-7)
+    assert_close(dqp32, dqp1, atol=1e-7)
+
+
 def test_quantdiff(hmba_lattice):
     hmba_lattice = hmba_lattice.radiation_on(copy=True)
     dmat = physics.radiation.quantdiffmat(hmba_lattice)