I had that, it came from the fact that when we have a b-value=0 the corresponding bvec is [0,0,0] and the arccos is NaN. I change the NaN to zero because for the SHORE the \phi(0) SH are constant so the coordinate in which they are calculated is not important. 
It was the fastest solution, we can also check for b=0 and change the corresponding bvec into something else and then turn the radius in 0.  In this way we avoid the warning.
None of these solution are 100% clean of course...