fix division by zero warnings for values near sun - #10881

Use of erfa.pn() would take care of this also.

indeed, erfa.pn is a logical choice here - @StuartLittlefair - obviously not a big deal, but maybe an idea to change it? It will do the right thing, with xyz as the last dimension.

erfa.pn is a nice solution, but it doesn't handle shapes perfectly and I lack @mhvk's numpy wizardry to be sure of doing the right thing. e.g suppose the q array is (N, 3), the way I do it gives sundist a shape of (N, 1) so I can use it in np.where later.

If I use erfa.pn then sundist has shape (N,) - so I could use sundist[:, np.newaxis] but then that doesn't work for scalars? Does the solution below work?

sundist, q = erfa.pn(q)
sundist = np.atleast_1d(sundist)[:, np.newaxis]
q = np.where(sundist > 1.0e-10, q, pco)

@StuartLittlefair - it would be sundist[..., np.newaxis] to append at the end (works for scalars too):

erfa.pn([1.,2.,3.])[0][..., np.newaxis]
Out[15]: array([3.74165739])
erfa.pn(np.ones((2,2,3)))[0][..., np.newaxis]
Out[16]: 
array([[[1.73205081],
        [1.73205081]],

       [[1.73205081],
        [1.73205081]]])

God bless you @mhvk and your numpy-Fu. And thanks @mkbrewer for the erfa reminder...