docs: 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)