Vectorization of `contains_lonlat`

[camposandro]

Would it be possible to have contains_lonlat handle multiple sky coordinates checks at a time, in a vectorized fashion, to speed up operations and prevent an excessive number of function calls?

What I am thinking of in terms of interface would be something similar to what we currently have, but the function call would allow for the specification of 1D arrays for ra and dec :

polygon.contains_lonlat([300, 300, 272, 272], [-50, -55, -55, -50])

Is there any particular reason why it is currently restricted to scalars?