Age-depth plot

def plot(self, agebins=50, p=(2.5, 97.5), ax=None):
        """Age-depth plot"""
        if ax is None:
            ax = plt.gca()
        ax.hist2d(np.repeat(self.depth, self.age_ensemble.shape[1]), self.age_ensemble.flatten(),
                   (len(self.depth), agebins), cmin=1)
        ax.step(self.depth, self.age_median(), where='mid', color='red')
        ax.step(self.depth, self.age_percentile(p[0]), where='mid', color='red', linestyle=':')
        ax.step(self.depth, self.age_percentile(p[1]), where='mid', color='red', linestyle=':')
        ax.set_ylabel('Age (cal yr BP)')
        ax.set_xlabel('Depth (cm)')
        ax.grid(True)
        return ax