The flow around different prolate (needle-like) and oblate (disc-like) spheroids is studied using a multi-relaxation-time lattice Boltzmann method. We compute the mean drag coefficient $C_{D,\unicode[STIX]{x1D719}}$ at different incident angles $\unicode[STIX]{x1D719}$ for a wide range of Reynolds numbers ($\mathit{Re}$). We show that the sine-squared drag law $C_{D,\unicode[STIX]{x1D719}}=C_{D,\unicode[STIX]{x1D719}=0^{\circ }}+(C_{D,\unicode[STIX]{x1D719}=90^{\circ }}-C_{D,\unicode[STIX]{x1D719}=0^{\circ }})\sin ^{2}\unicode[STIX]{x1D719}$ holds up to large Reynolds numbers, $\mathit{Re}=2000$. Further, we explore the physical origin behind the sine-squared law, and reveal that, surprisingly, this does not occur due to linearity of flow fields. Instead, it occurs due to an interesting pattern of pressure distribution contributing to the drag at higher $\mathit{Re}$ for different incident angles. The present results demonstrate that it is possible to perform just two simulations at $\unicode[STIX]{x1D719}=0^{\circ }$ and $\unicode[STIX]{x1D719}=90^{\circ }$ for a given $\mathit{Re}$ and obtain particle-shape-specific $C_{D}$ at arbitrary incident angles. However, the model has limited applicability to flatter oblate spheroids, which do not exhibit the sine-squared interpolation, even for $\mathit{Re}=100$, due to stronger wake-induced drag. Regarding lift coefficients, we find that the equivalent theoretical equation can provide a reasonable approximation, even at high $\mathit{Re}$, for prolate spheroids.