For very narrow spherical Gaussians, numerical issues can arise in finite precision floating calculations. We found the following special cases useful in our implementation. To evaluate the $M$ function when the variance is low, $v < 0.1$, we use $$ M\left(v,\theta _{\text{cone}},\theta _o\right) \approx \exp \left(\log \left(I_0\left(a\right)\right)+b-\frac{1}{v}+0.6931+\log \left(\frac{1}{2 v}\right)\right) $$ where $a = \cos \theta _{\text{cone}} \cos \theta_o / v$ and $b = \sin \theta _{\text{cone}} \sin \theta_o / v$. The log of the Bessel function can also be problematic for large arguments for we use the following special case for $\log( I_0(x))$ when $x > 12$: $$ \log( I_0(x)) \approx x + 0.5 \left( -\log ( 2 \pi ) + \log \frac{1}{x} + \frac{1}{8x} \right). $$