The following families of integrals involving a product of two Bessel functions \(J_n(\nu)\) are considered: \begin{align*} I(m,n,k,\alpha) &= \int_0^\infty J_{m+1/2}(\nu) J_{n+1/2}(\nu) \frac{1}{\nu^k\sqrt{\alpha^2-\nu^2}} \,\mathrm{d}\nu, \\ J(m,n,k,\alpha) &= \int_0^\infty J_{m+1/2}(\nu) J_{n+1/2}(\nu) \frac{\sqrt{\alpha^2-\nu^2}}{\nu^k} \,\mathrm{d}\nu. \end{align*} Due to the half-integer in the index, these Bessel functions can be written more conveniently in terms of spherical Bessel functions. The authors are interested in evaluating these integrals numerically for particular choices of the parameters \(m,n,k,\alpha\), which is important in studying the acoustic and electromagnetic scattering from circular disks and apertures. Their solution is a rapidly converging series that is based on the Mellin-Barnes integral representation of the product of Bessel functions. They implemented their method in Matlab and considered several examples to illustrate their results.