code: def _prepare_ws(self, w0, mmap, n_steps):
          
          from ..dynamics import PhaseSpacePosition
          if not isinstance(w0, PhaseSpacePosition):
              w0 = PhaseSpacePosition.from_w(w0)
  
          arr_w0 = w0.w(self._func_units)
  
          self.ndim, self.norbits = arr_w0.shape
          self.ndim = self.ndim//2
  
          return_shape = (2*self.ndim, n_steps+1, self.norbits)
          if mmap is None:
              # create the return arrays
              ws = np.zeros(return_shape, dtype=float)
  
          else:
              if mmap.shape != return_shape:
                  raise ValueError("Shape of memory-mapped array doesn't match "
                                   "expected shape of return array ({} vs {})"
                                   .format(mmap.shape, return_shape))
  
              if not mmap.flags.writeable:
                  raise TypeError("Memory-mapped array must be a writable mode, "
                                  " not '{}'".format(mmap.mode))
  
              ws = mmap
  
          return w0, arr_w0, ws
docs: Decide how to make the return array. If mmap is False, this returns a
          full array of zeros, but with the correct shape as the output. If mmap
          is True, return a pointer to a memory-mapped array. The latter is
          particularly useful for integrating a large number of orbits or
          integrating a large number of time steps.