def zmax(self, return_times=False, func=np.mean,
             interp_kwargs=None, minimize_kwargs=None,
             approximate=False):
        if return_times and func is not None:
            raise ValueError("Cannot return times if reducing "
                             "using an input function. Pass `func=None` if "
                             "you want to return all individual values "
                             "and times.")
        if func is None:
            reduce = False
            func = lambda x: x
        else:
            reduce = True
        if self.t[-1] < self.t[0]:
            self = self[::-1]
        vals = []
        times = []
        for orbit in self.orbit_gen():
            v, t = orbit._max_helper(np.abs(orbit.cylindrical.z),
                                     interp_kwargs=interp_kwargs,
                                     minimize_kwargs=minimize_kwargs,
                                     approximate=approximate)
            vals.append(func(v))
            times.append(t)
        return self._max_return_helper(vals, times, return_times, reduce)