vectorizing functions with multiple return values?

My goal is to provide a vectorized Python binding to C++ functions that return more than a single value, like

void ang2vec (double theta, double phi, double &x, double &y, double &z)
  {
  x=sin(theta)*cos(phi);
  y=sin(theta)*sin(phi);
  z=cos(theta);
  }

py::vectorize(ang2vec) does not work, producing error messages like
"pybind11/numpy.h:732:5: error: forming pointer to reference type ‘double&’".

I tried to translate this into a more Python-like interface:

tuple<double,double,double> ang2vec2 (double theta, double phi)
  {
  return tuple<double,double,double>(sin(theta)*cos(phi),sin(theta)*sin(phi),cos(theta));
  }

but unfortunately vectorizing this fails as well.

Are there fundamental obstacles to supporting this scenario, or am I just the first one trying to do this? :)

Start by tracing py::vectorize and the failure reported at pybind11/numpy.h:732 for both the reference-argument and tuple-return examples. Determine whether multiple outputs can be supported and identify the tests needed to establish the expected Python-facing behavior; done means a documented, tested resolution of this use case.