  We consider the classical momentum- or velocity-dependent two-dimensional
Hamiltonian given by $$\mathcal H_N = p_1^2 + p_2^2 +\sum_{n=1}^N \gamma_n(q_1
p_1 + q_2 p_2)^n ,$$ where $q_i$ and $p_i$ are generic canonical variables,
$\gamma_n$ are arbitrary coefficients, and $N\in \mathbb N$. For $N=2$, being
both $\gamma_1,\gamma_2$ different from zero, this reduces to the classical
Zernike system. We prove that $\mathcal H_N$ always provides a superintegrable
system (for any value of $\gamma_n$ and $N$) by obtaining the corresponding
constants of the motion explicitly, which turn out to be of higher-order in the
momenta. Such generic results are not only applied to the Euclidean plane, but
also to the sphere and the hyperbolic plane. In the latter curved spaces,
$\mathcal H_N $ is expressed in geodesic polar coordinates showing that such a
new superintegrable Hamiltonian can be regarded as a superposition of the
isotropic 1:1 curved (Higgs) oscillator with even-order anharmonic curved
oscillators plus another superposition of higher-order momentum-dependent
potentials. Furthermore, the symmetry algebra determined by the constants of
the motion is also studied, giving rise to a $(2N-1)$th-order polynomial
algebra. As a byproduct, the Hamiltonian $\mathcal H_N $ is interpreted as a
family of superintegrable perturbations of the classical Zernike system.
Finally, it is shown that $\mathcal H_N$ (and so the Zernike system as well) is
endowed with a Poisson $\mathfrak{sl}(2,\mathbb R)$-coalgebra symmetry which
would allow for further possible generalizations that are also discussed.
