<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def zmax(self, return_times=False, func=np.mean, interp_kwargs=None, minimize_kwargs=None, approximate=False):
 """ Estimate the maximum ``z`` height of the orbit by identifying local maxima in the absolute value of the ``z`` position and interpolating between timesteps near the maxima. By default, this returns the mean of all local maxima. To get, e.g., the largest ``z`` excursion, pass in ``func=np.max``. To get all ``z`` maxima, pass in ``func=None``. Parameters func : func (optional) A function to evaluate on all of the identified z maximum times. return_times : bool (optional) Also return the times of maximum. interp_kwargs : dict (optional) Keyword arguments to be passed to :class:`scipy.interpolate.InterpolatedUnivariateSpline`. minimize_kwargs : dict (optional) Keyword arguments to be passed to :class:`scipy.optimize.minimize`. approximate : bool (optional) Compute approximate values by skipping interpolation. Returns ------- zs : float, :class:`~numpy.ndarray` Either a single number or an array of maximum z heights. times : :class:`~numpy.ndarray` (optional, see ``return_times``) If ``return_times=True``, also returns an array of the apocenter times. """


        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

        # time must increase
        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)