code: def align_circulation_with_z(self, circulation=None):
          if circulation is None:
              circulation = self.circulation()
          circulation = atleast_2d(circulation, insert_axis=1)
          cart = self.cartesian
          pos = cart.xyz
          vel = np.vstack((cart.v_x.value[None],
                           cart.v_y.value[None],
                           cart.v_z.value[None])) * cart.v_x.unit
          if pos.ndim < 3:
              pos = pos[...,np.newaxis]
              vel = vel[...,np.newaxis]
          if (circulation.shape[0] != self.ndim or
                  circulation.shape[1] != pos.shape[2]):
              raise ValueError("Shape of 'circulation' array should match the "
                               "shape of the position/velocity (minus the time "
                               "axis).")
          new_pos = pos.copy()
          new_vel = vel.copy()
          for n in range(pos.shape[2]):
              if circulation[2,n] == 1 or np.all(circulation[:,n] == 0):
                  continue
              if sum(circulation[:,n]) > 1:
                  logger.warning("Circulation about multiple axes - are you sure "
                                 "the orbit has been integrated for long enough?")
              if circulation[0,n] == 1:
                  circ = 0
              elif circulation[1,n] == 1:
                  circ = 1
              else:
                  raise RuntimeError("Should never get here...")
              new_pos[circ,:,n] = pos[2,:,n]
              new_pos[2,:,n] = pos[circ,:,n]
              new_vel[circ,:,n] = vel[2,:,n]
              new_vel[2,:,n] = vel[circ,:,n]
          return self.__class__(pos=new_pos.reshape(cart.xyz.shape),
                                vel=new_vel.reshape(cart.xyz.shape),
                                t=self.t,
                                hamiltonian=self.hamiltonian)