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)