code: def get_xyz(self, xyz_axis=0):
          
          # Add new axis in x, y, z so one can concatenate them around it.
          # NOTE: just use np.stack once our minimum numpy version is 1.10.
          result_ndim = self.ndim + 1
          if not -result_ndim <= xyz_axis < result_ndim:
              raise IndexError('xyz_axis {0} out of bounds [-{1}, {1})'
                               .format(xyz_axis, result_ndim))
          if xyz_axis < 0:
              xyz_axis += result_ndim
  
          # Get components to the same units (very fast for identical units)
          # since np.concatenate cannot deal with quantity.
          unit = self._x1.unit
  
          sh = self.shape
          sh = sh[:xyz_axis] + (1,) + sh[xyz_axis:]
          components = [getattr(self, '_'+name).reshape(sh).to(unit).value
                        for name in self.attr_classes]
          xs_value = np.concatenate(components, axis=xyz_axis)
          return u.Quantity(xs_value, unit=unit, copy=False)
docs: Return a vector array of the x, y, and z coordinates.
  
          Parameters
          ----------
          xyz_axis : int, optional
              The axis in the final array along which the x, y, z components
              should be stored (default: 0).
  
          Returns
          -------
          xs : `~astropy.units.Quantity`
              With dimension 3 along ``xyz_axis``.