def kippenhahn(self, num_frame, xax, t0_model=0,
                   title='Kippenhahn diagram', tp_agb=0., t_eps=5.e2,
                   plot_star_mass=True, symbol_size=8, c12_bm=False,
                   print_legend=True):
        """Kippenhahn plot as a function of time or model.

        Parameters
        ----------
        num_frame : integer
            Number of frame to plot this plot into, if <0 open no new
            figure.
        xax : string
            Either 'model', 'time' or 'logtimerev' to indicate what is
            to be used on the x-axis.
        t0_model : integer, optional
            If xax = 'time' then model for the zero point in time, for
            AGB plots this would be usually the model of the 1st TP,
            which can be found with the Kippenhahn plot.  The default
            is 0.
        title : string, optional
            The figure title.  The default is "Kippenhahn diagram".
        tp_agb : float, optional
            If > 0. then,
            ylim=[h1_min*1.-tp_agb/100 : h1_max*1.+tp_agb/100] with
            h1_min, h1_max the min and max H-free core mass coordinate.
            The default is 0. .
        t_eps : float, optional
            Final time for logtimerev.  The default is '5.e2'.
        plot_star_mass : boolean, optional
            If True, then plot the stellar mass as a line as well.  The
            default is True.
        symbol_size : integer, optional
            Size of convection boundary marker.  The default is 8.
        c12_bm : boolean, optional
            If we plot c12_boundary_mass or not.  The default is False.
        print_legend : boolean, optionla
            Show or do not show legend.  The defalut is True.

        """

        if num_frame >= 0:
            pyl.figure(num_frame)

        t0_mod=[]

        if xax == 'time':
            xaxisarray = self.get('star_age')
            if t0_model > 0:
                ind=self.get('model_number')
                t0_model=where(ind>t0_model)[0][0]
                t0_mod=xaxisarray[t0_model]
            else:
                t0_mod = 0.
            print('zero time is '+str(t0_mod))
        elif xax == 'model':
            xaxisarray = self.get('model_number')
            #t0_mod=xaxisarray[t0_model]
            t0_mod = 0.
        elif xax == 'logtimerev':
            xaxi    = self.get('star_age')
            xaxisarray = np.log10(np.max(xaxi)+t_eps-xaxi)
            t0_mod = 0.
        else:
            print('kippenhahn_error: invalid string for x-axis selction.'+\
                  ' needs to be "time" or "model"')


        plot_bounds=True
        try:
            h1_boundary_mass  = self.get('h1_boundary_mass')
            he4_boundary_mass = self.get('he4_boundary_mass')
            if c12_bm:
                c12_boundary_mass = self.get('c12_boundary_mass')
        except:
            try:
                h1_boundary_mass  = self.get('he_core_mass')
                he4_boundary_mass = self.get('c_core_mass')
                if c12_bm:
                    c12_boundary_mass = self.get('o_core_mass')
            except:
                plot_bounds=False

        star_mass         = self.get('star_mass')
        mx1_bot           = self.get('mx1_bot')*star_mass
        mx1_top           = self.get('mx1_top')*star_mass
        mx2_bot           = self.get('mx2_bot')*star_mass
        mx2_top           = self.get('mx2_top')*star_mass


        if xax == 'time':
            if t0_model>0:
                pyl.xlabel('$t - t_0$ $\mathrm{[yr]}$')
            else:
                pyl.xlabel('t / yrs')
        elif xax == 'model':
            pyl.xlabel('model number')
        elif xax == 'logtimerev':
            pyl.xlabel('$\log(t_{final} - t)$  $\mathrm{[yr]}$')

        pyl.plot(xaxisarray[t0_model:]-t0_mod,mx1_bot[t0_model:],linestyle='None',color='blue',alpha=0.3,marker='o',markersize=symbol_size,label='convection zones')
        pyl.plot(xaxisarray[t0_model:]-t0_mod,mx1_top[t0_model:],linestyle='None',color='blue',alpha=0.3,marker='o',markersize=symbol_size)
        pyl.plot(xaxisarray[t0_model:]-t0_mod,mx2_bot[t0_model:],linestyle='None',color='blue',alpha=0.3,marker='o',markersize=symbol_size)
        pyl.plot(xaxisarray[t0_model:]-t0_mod,mx2_top[t0_model:],linestyle='None',color='blue',alpha=0.3,marker='o',markersize=symbol_size)

        if plot_bounds:
            pyl.plot(xaxisarray[t0_model:]-t0_mod,h1_boundary_mass[t0_model:],color='red',linewidth=2,label='H-free core')
            pyl.plot(xaxisarray[t0_model:]-t0_mod,he4_boundary_mass[t0_model:],color='green',linewidth=2,linestyle='dashed',label='He-free core')
            if c12_bm:
                pyl.plot(xaxisarray[t0_model:]-t0_mod,c12_boundary_mass[t0_model:],color='purple',linewidth=2,linestyle='dotted',label='C-free core')
        if plot_star_mass is True:
            pyl.plot(xaxisarray[t0_model:]-t0_mod,star_mass[t0_model:],label='$M_\star$')
        pyl.ylabel('$m_\mathrm{r}/\mathrm{M}_\odot$')
        if print_legend:
            pyl.legend(loc=2)
        if tp_agb > 0.:
            h1_min = min(h1_boundary_mass[t0_model:])
            h1_max = max(h1_boundary_mass[t0_model:])
            h1_min = h1_min*(1.-old_div(tp_agb,100.))
            h1_max = h1_max*(1.+old_div(tp_agb,100.))
            print('setting ylim to zoom in on H-burning:',h1_min,h1_max)
            pyl.ylim(h1_min,h1_max)