Background
The BunchTuneAnalysis class "normalizes" the phase space coordinates using a set of Twiss parameters (alpha, beta, dispersion) in each plane. The phase space coordinates $\mathbf{x} = (x, x')$ are linearly related to the normalized coordinates $\tilde{\mathbf{x}} = (\tilde{x}, \tilde{x}')$:
$$
\mathbf{x} = \mathbf{V} \tilde{\mathbf{x}}.
$$
Ignoring dispersion, the symplectic $2 \times 2$ (un)normalization matrix $\mathbf{V}$ is written
$$
\mathbf{V} = \mathbf{V}(\alpha, \beta) =
\frac{1}{\beta}
\begin{pmatrix}
+\beta & +0 \\ -\alpha & +1
\end{pmatrix}.
$$
We choose ${ \alpha, \beta }$ so that the one-turn transfer matrix $\mathbf{M}$ decomposes as
$$
\mathbf{M} = \mathbf{V} \mathbf{R} \mathbf{V}^{-1},
$$
where $\mathbf{R} = \mathbf{R}(\phi)$ is a rotation matrix by phase advance $\phi$. We define the tune $\nu$ as the number of oscillations per turn:
$$
\nu = \frac{\phi}{2 \pi}.
$$
This class calculates the tunes from the bunch coordinates on subsequent turns. After each set of coordinates is normalized, the change in phase angle is computed for each particle. The change in phase angle is computed from $\tan(\phi) = \tilde{x}' / \tilde{x}$.
Beam-based normalization
The method above depends crucially on the input lattice parameters ${ \alpha, \beta }$. The tune estimates will be incorrect if the beam is not matched to the lattice optics. This could be the case for strong space charge intensities, where the linear component of the space charge forces effectively modifies the lattice functions. In this case, if we assume the beam covariance matrix $\mathbf{\Sigma} = \langle \mathbf{x} \mathbf{x}^T \rangle$ is periodic, we can compute the correct normalization matrix directly from $\mathbf{\Sigma}$. See here and the paper by Lebedev and Bogacz. The method is to compute the eigenvectors of $\mathbf{\Sigma} \mathbf{U}$, where $\mathbf{U}$ is the unit symplectic matrix.
This method will only work if there enough particles in the bunch to compute the covariance matrix, which requires a few thousand particles. One option for implementation would be to keep the same class, but add an option to use the covariance matrix instead of predefined twiss parameters.
Coupling
If there is $x-y$ coupling in the accelerator lattice, the $4 \times 4$ normalization matrix $\mathbf{V}$ will have off-block-diagonal elements. Although the method of estimating the tunes would remain the same, but the parameterization of $\mathbf{V}$ would change. There are multiple parameterizations available. A new tune calculator class could be developed that adopts one of these parameterizations, or multiple classes to provide options for different parameterizations.
Beam-based normalization is straightforward using the eigenvector method described above. Note that linear space charge can generate coupling even without coupled focusing elements.
Other tune calculation methods
Might it make sense to also include frequency-based tune estimation methods in PyORBIT? This can always be done offline by collecting turn-by-turn coordinates, but since this is a common operation, it might make sense to include. These methods can also utilize the normalization methods described above.
