Because the interval of integration is finite we can tackle this integral with the method of steepest descent. We'll write $\DeclareMathOperator{rre}{Re} \DeclareMathOperator{iim}{Im}$
$$
\begin{align}
S(p) &= \int_0^1 \frac{y \sqrt{1-y^2}}{(\varepsilon^2-1)y^2+1} \sin(py) dy \\
&= \frac{1}{2i} \left( \int_0^1 \frac{y \sqrt{1-y^2}}{(\varepsilon^2-1)y^2+1} e^{ipy} dy - \int_0^1 \frac{y \sqrt{1-y^2}}{(\varepsilon^2-1)y^2+1} e^{-ipy} dy \right) \\
&= \frac{1}{2i} \Bigl( I_1(p) - I_2(p) \Bigr) \tag{1}
\end{align}
$$
and study $I_1$ and $I_2$ independently. (If $\varepsilon^2$ were real we would only need to consider $\iim I_1(p)$, but let's study both to be safe.)
These integrals exist as long as $\varepsilon^2$ does not lie in the real interval $(-\infty,0)$. Let $\pm y_\varepsilon$ be the roots of the equation
$$
(\varepsilon^2-1)y^2+1 = 0
$$
(assuming $\varepsilon \neq \pm 1$) and define the real number
$$
a = a(\varepsilon) = \begin{cases}
1 & \text{if } -1 \leq \varepsilon \leq 1, \\
\tfrac{1}{2} |\iim y_\varepsilon| & \text{otherwise.}
\end{cases}
$$
By definition we have $a > 0$.
Let's put $I_1$ into more standard terms by writing
$$
I_1(p) = \int_0^1 f(y) e^{pg(y)}\,dy.
$$
The current integration contour runs along a path of constant altitude on the surface described by
$$
\rre g(y) = \rre(iy) = -\iim y,
$$
and at the endpoints of the contour (namely $y=0$ and $y=1$) the paths of steepest descent are rays parallel to the positive imaginary axis. We therefore deform the contour according to the following picture, where the original contour (the segment from $y=0$ to $y=1$) is shown in blue:
By our choice of $a$ we avoid the poles of the integrand and may thus conclude that
$$
I_1(p) = \int_0^1 = \int_{\gamma_1} + \int_{\gamma_2} - \int_{\gamma_3}. \tag{2}
$$
The contours $\gamma_1$ and $\gamma_3$ follow the paths of steepest descent, and the integral over the contour $\gamma_2$ is exponentially small:
$$
\left|\int_{\gamma_2} \frac{y \sqrt{1-y^2}}{(\varepsilon^2-1)y^2+1} e^{ipy} dy \right| \leq e^{-ap} \int_{\gamma_2} \left|\frac{y \sqrt{1-y^2}}{(\varepsilon^2-1)y^2+1}\right| |dy|. \tag{3}
$$
We can parameterize the contour $\gamma_1$ by $y = it$, $dy = idt$, to get
$$
\int_{\gamma_1} = -\int_0^a \frac{t \sqrt{1+t^2}}{1-(\varepsilon^2-1)t^2} e^{-pt}\,dt,
$$
which can be handled using Watson's lemma. Indeed,
$$
\int_{\gamma_1} = -\frac{1}{p^2} + O\left(\frac{1}{p^4}\right). \tag{4}
$$
Similarly, $\gamma_3$ can be parameterized by $y = 1+it$, $dy = idt$, to get
$$
\int_{\gamma_3} = e^{i(p+\pi/2)} \int_0^a \frac{(1+it)\sqrt{1-(1+it)^2}}{(\varepsilon^2-1)(1+it)^2 + 1} e^{-pt}\,dt,
$$
which can again be handled using Watson's lemma. We calculate
$$
\int_{\gamma_3} = e^{i(p+\pi/4)} \frac{\sqrt{\pi/2}}{\varepsilon^2 p^{3/2}} + O\left(\frac{1}{p^{5/2}}\right). \tag{5}
$$
Combining $(3)$, $(4)$, and $(5)$ in $(2)$ we find that
$$
I_1(p) = e^{i(p-3\pi/4)} \frac{\sqrt{\pi/2}}{\varepsilon^2 p^{3/2}} - \frac{1}{p^2} + O\left(\frac{1}{p^{5/2}}\right). \tag{6}
$$
An analogous method can be applied to $I_2$ using the contour in the lower half-plane joining the points $0, -ia, 1-ia, 1$. Doing this, we obtain the asymptotic
$$
I_2(p) = e^{-i(p-3\pi/4)} \frac{\sqrt{\pi/2}}{\varepsilon^2 p^{3/2}} - \frac{1}{p^2} + O\left(\frac{1}{p^{5/2}}\right). \tag{7}
$$
Substituting $(6)$ and $(7)$ into $(1)$ we conclude that
$$
S(p) = \sin\left(p - \frac{3\pi}{4}\right) \frac{\sqrt{\pi/2}}{\varepsilon^2 p^{3/2}} + O\left(\frac{1}{p^{5/2}}\right). \tag{8}
$$
Here is a plot of the numerical evaluation of $S(p)$ in blue and the asymptotic in $(8)$ in red. In this we've taken $\varepsilon = 3$.
The asymptotic $(8)$ is slightly out of phase with $S(p)$. This is corrected in the next term of the asymptotic expansion, given by
$$
S(p) = \sin\left(p - \frac{3\pi}{4}\right) \frac{\sqrt{\pi/2}}{\varepsilon^2 p^{3/2}} + \cos\left(p - \frac{3\pi}{4}\right) \frac{3(8-3\varepsilon^2)\sqrt{\pi/2}}{8 \varepsilon^4 p^{5/2}} + O\left(\frac{1}{p^{7/2}}\right).
$$
$\tag{9}$
