def kippenhahn_CO(self, num_frame, xax, t0_model=0,
                      title='Kippenhahn diagram', tp_agb=0.,
                      ylim_CO=[0,0]):
        """
        Kippenhahn plot as a function of time or model with CO ratio

        Parameters
        ----------
        num_frame : integer
            Number of frame to plot this plot into.
        xax : string
            Either model or time to indicate what is to be used on the
            x-axis.
        t0_model : integer, optional
            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
            Figure title.  The defalut 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 defalut is 0.
        ylim_CO : list
            if ylim_CO is [0,0], then it is automaticly set.  The
            default is [0,0].

        """

        pyl.figure(num_frame)

        if xax == 'time':
            xaxisarray = self.get('star_age')
        elif xax == 'model':
            xaxisarray = self.get('model_number')
        else:
            print('kippenhahn_error: invalid string for x-axis selction.'+\
                  ' needs to be "time" or "model"')

        t0_mod=xaxisarray[t0_model]

        plot_bounds=True
        try:
            h1_boundary_mass  = self.get('h1_boundary_mass')
            he4_boundary_mass = self.get('he4_boundary_mass')
        except:
            try:
                h1_boundary_mass  = self.get('he_core_mass')
                he4_boundary_mass = self.get('c_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
        surface_c12       = self.get('surface_c12')
        surface_o16       = self.get('surface_o16')

        COratio=old_div((surface_c12*4.),(surface_o16*3.))

        pyl.plot(xaxisarray[t0_model:]-t0_mod,COratio[t0_model:],'-k',label='CO ratio')
        pyl.ylabel('C/O ratio')
        pyl.legend(loc=4)
        if ylim_CO[0] is not 0 and  ylim_CO[1] is not 0:
            pyl.ylim(ylim_CO)
        if xax == 'time':
            pyl.xlabel('t / yrs')
        elif xax == 'model':
            pyl.xlabel('model number')

        pyl.twinx()
        if plot_bounds:
            pyl.plot(xaxisarray[t0_model:]-t0_mod,h1_boundary_mass[t0_model:],label='h1_boundary_mass')
            pyl.plot(xaxisarray[t0_model:]-t0_mod,he4_boundary_mass[t0_model:],label='he4_boundary_mass')
        pyl.plot(xaxisarray[t0_model:]-t0_mod,mx1_bot[t0_model:],',r',label='conv bound')
        pyl.plot(xaxisarray[t0_model:]-t0_mod,mx1_top[t0_model:],',r')
        pyl.plot(xaxisarray[t0_model:]-t0_mod,mx2_bot[t0_model:],',r')
        pyl.plot(xaxisarray[t0_model:]-t0_mod,mx2_top[t0_model:],',r')
        pyl.plot(xaxisarray[t0_model:]-t0_mod,star_mass[t0_model:],label='star_mass')
        pyl.ylabel('mass coordinate')
        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)