docs: def kip_cont(self, modstart, modstop, dcoeff_thresh=1.e12,
                   xres=1000, ylims=[0., 0.], xlims=[0., 0.], yres=2000,
                   ixaxis='log_time_left', outfile='',
                   landscape_plot=False, codev='KEP', kepswitch=12555,
                   outlines=False, write_preproc=False, hatching=False):
          """
          This function creates a Kippenhahn diagram as a contour plot of
          the se output using a convection flag or diffusion coefficient
          threshold.
  
          Parameters
          ----------
          modstart, modstop : integer
              First and last cycle numbers to plot.
          dcoeff_thresh : float, optional
              Diffusion coefficient threshold, above which assume
              convection is present (only for MESA, where
              convection_indicator is not there)).  The default is 1.e12.
          xres, yres : integer, optional
              x and y resolution.  The default for xres is 1000, the
              default for yres is 2000.
          xlims, ylims : list, optional
              x and y plot limits.  The default is [0., 0.].
          ixaxis : string, optional
              Choose one of 'age', 'model_number' or 'log_time_left'.  The
              default is 'log_time_left'.
          outfile : string, optional
              Name of output file including extension, which will be saved
              in 300 dpi.  The default is ''.
          landscape_plot : boolean, optional
              For a landscape plot.  The default is False.
          codev : string, optional
              Coose one of 'KEP', 'MES' or 'GNV'.  The default is 'KEP'.
          kepswitch : integer, optional
              The default is 12555.
          outlines : boolean, optional
              The default is False.
          write_preproc : boolean, optional
              The default is False.
          hatching : boolean, optional
              The default if False.
  
          Examples
          --------
          >>> pt.kip_cont(0, -1, xres=1000, yres=1000, ixaxis='log_time_left')
  
          Notes
          -----
          Unfortunately, although the implicit do loops are faster, I
          cannot implement a % bar, but this example call should not take
          more than 1 minute.
  
          """
  
          tmp_cycles=self.se.cycles
          original_cyclelist = [int(tmp_cycles[i]) for i in range(len(tmp_cycles))]
          nmodels=len(self.se.cycles[modstart:modstop])
          print('nmodels           = ', nmodels)
          sparse = int(max(1,old_div(nmodels,xres)))
          cyclelist = original_cyclelist[modstart:modstop:sparse]
  
          age_unit=self.get('age_unit')
          if codev == 'KEP':
              oneyear = 365.*24.*3600.
          elif codev == 'MES':
              oneyear=self.get('one_year')
          elif codev == 'GNV':
              oneyear=1.e0
  
  
          m_min=ylims[0]
          m_max=ylims[1]
          if m_min==0.:
              m_min=0.
              ylims[0] = m_min
          if m_max==0.:
              m_max=float(self.se.get('mini'))
              ylims[1] = m_max
  
  
          dy = old_div((m_max-m_min),float(yres))
          y = np.arange(m_min, m_max, dy)
  
          Z = np.zeros([len(y),len(cyclelist)],float)
  
          datatype='convection'
          if 'convection_indicator' not in self.se.dcols:
              datatype='dcoeff'
  
          def list_to_string(array):
              string=''
              for el in array:
                  string+='    '+str(el)+'    '
              string+='\n'
              return string
  
          if write_preproc==True:
              if os.path.exists(self.sedir+'/conv_data_preproc.txt'):
                  print('Preprocessor already exists (conv_data_preproc.txt).')
                  print('Please remove this file if you want to create a new preprocessor file')
                  sys.exit()
              print('write_preproc=True')
              print('getting complete dataset')
              print('getting mass coordinates')
              allcycs=[int(self.se.cycles[i]) for i in range(len(self.se.cycles))]
              mass=self.se.get(allcycs,'mass')
              print('getting conv')
              if datatype=='dcoeff':
                  conv=se.f.se.get(allcycs,'dcoeff')
              else:
                  conv=self.se.get(allcycs,'convection_indicator')
              print('getting cycles')
              models=self.se.cycles
              print('getting ages')
              if codev=='KEP':
                  deltat=self.se.get(allcycs,'deltat')
                  ages=np.cumsum(deltat)
              else:
                  ages=self.se.ages
              print('writing preprocessor file')
              f=open(self.sedir+'/conv_data_preproc.txt','w')
              #f=open('./conv_data_preproc.txt','w')
              f.write('FORMAT OF THIS FILE:\n')
              f.write('model_number\n')
              f.write('star_age/yr\n')
              f.write('mass/Mo\n')
              f.write('convection_indicator or dcoeff if unavailable \n')
              f.write('*************************************************************\n')
  
              for i in range(len(models)):
                  percent = int(i*100/len(self.se.cycles))
                  sys.stdout.flush()
                  sys.stdout.write("\rprogress " + "...%d%%" % percent)
                  f.write(str(models[i])+'\n')
                  f.write(str(ages[i])+'\n')
                  f.write(list_to_string(mass[i]))
                  f.write(list_to_string(conv[i]))
  
              f.close()
              print('preprocessor file written')
              sys.exit()
  
          def realarray(array):
              """
              Parameters
              ----------
              array : list
                  Takes array of strings
  
              Returns
              --------
              list
                  array of floats
  
              """
              floatarray=np.array([float(el) for el in array])
              return floatarray
  
          print('getting data...')
  
          if os.path.exists(self.sedir+'/conv_data_preproc.txt'):
              print('... from preprocessor file ...')
              lines=open(self.sedir+'/conv_data_preproc.txt','r').readlines()
              #lines=open('./conv_data_preproc.txt','r').readlines()
              lines=lines[6:]
              mod=np.array([float(line.split()[0]) for line in lines[::4]])
              age=old_div(np.array([float(line.split()[0]) for line in lines[1::4]]),oneyear)
              # masses and mix types will be read just as an array of strings (each string is a profile)
              # to be extracted later
              masses=lines[2::4]
              conv=lines[3::4]
  
              tmp=[realarray(massco.split()) for massco in masses]
              masses=tmp
  
              tmp=[realarray(mixvec.split()) for mixvec in conv]
              conv=tmp
  
              print(len(age),len(mod),len(masses),len(conv))
  
          else:
              print('... from HDF5 data ...')
              masses=self.se.get(cyclelist,'mass')
              mod=cyclelist
              if self.deltat == None:
                  print('getting deltat since not initialised')
                  self.deltat = self.se.get(self.se.cycles,'deltat')
              age=old_div(np.cumsum(self.deltat), (3600.*24.*365.))
  
              if datatype=='convection':
                  conv=self.se.get(cyclelist,'convection_indicator')
              else:
                  conv=self.se.get(cyclelist,'dcoeff')
  
          print('calculating convection matrix... ')
  
          print(len(masses))
          conv_i_vec_matrix=np.zeros([len(y),len(masses)],float)
  
          for i in range(len(masses)):
              percent = int(i*100/len(masses))
              sys.stdout.flush()
              sys.stdout.write("\rcreating color map mix " + "...%d%%" % percent)
  
              mass = masses[i]
              mix = conv[i]
  #            print i, len(mass),len(mix)
              # rearrange centre to surface:
              if mass[-1]<mass[0]:
                  mass=np.append(mass,[0.])[::-1]
                  mix=np.append(mix,[0.])[::-1]
              else:
                  mass=np.append(mass[::-1],[0.])[::-1]
                  mix=np.append(mix[::-1],[0.])[::-1]
              f = interpolate.interp1d(mass, mix, bounds_error=False, fill_value=0.)
              mixnew = f(y)
  
              for j in range(len(mixnew)):
                  conv_i_vec_matrix[j,i] = mixnew[j]
  
  
          print('getting stellar mass...')
          mtot = np.array([np.max(m) for m in masses])
  
          #########################################################################
          # PLOT:
          fsize=20
  
          params = {'axes.labelsize':  fsize,
            'text.fontsize':   fsize,
            'legend.fontsize': fsize,
            'xtick.labelsize': fsize*0.8,
            'ytick.labelsize': fsize*0.8,
            'xtick.major.pad': 8,          # distance to major tick label in points
            'ytick.major.pad': 6,          # distance to major tick label in points
            'text.usetex': False}
          pl.rcParams.update(params)
          fig = pl.figure(1)
          ax = pl.axes()
          if landscape_plot == True:
              fig.set_size_inches(9,4)
              fsize=20
              pl.gcf().subplots_adjust(bottom=0.15)
              pl.gcf().subplots_adjust(right=0.85)
  
          ax.set_ylabel('$\mathrm{Mass }(M_\odot)$',fontsize=fsize)
  
          xxx=[]
          age_array = age
  
          if ixaxis == 'log_time_left':
          # log of time left until core collapse
              gage= np.array(age_array)
              lage=np.zeros(len(gage))
              agemin = max(old_div(abs(gage[-1]-gage[-2]),5.),1.e-10)
              for i in np.arange(len(gage)):
                  if gage[-1]-gage[i]>agemin:
                      lage[i]=np.log10(gage[-1]-gage[i]+agemin)
                  else :
                      lage[i]=np.log10(agemin)
              xxx = np.array([lage[int(i)] for i in cyclelist]) # np.array([lage[i] for i in cyclelist])
              print(len(xxx))
              print('plot versus time left')
              ax.set_xlabel('$\mathrm{log}_{10}(t^*/\mathrm{yr})$',fontsize=fsize)
              if xlims == [0.,0.]:
                  xlims = [xxx[0],xxx[-1]]
          elif ixaxis =='model_number':
              cyclelist=mod[modstart:modstop:sparse] # np.array([int(cycle) for cycle in cyclelist])
              xxx= cyclelist
              print('plot versus model number')
              ax.set_xlabel('Model number',fontsize=fsize)
              if xlims == [0.,0.]:
                  xlims = [cyclelist[0],cyclelist[-1]]
          elif ixaxis =='age':
              xxx = old_div(np.array([age_array[int(i)] for i in cyclelist]),1.e6)
              print('plot versus age')
              ax.set_xlabel('Age [Myr]',fontsize=fsize)
              if xlims == [0.,0.]:
                  xlims = [xxx[0],xxx[-1]]
  
          #cmapMIX=mpl.colors.ListedColormap(['w','#8B8386']) # rose grey
          cmapMIX = mpl.colors.ListedColormap(['#C4C4C4']) # grey77
          #cmap = mpl.colors.ListedColormap(['w','b'])
  
          # some stuff for rasterizing only the contour part of the plot, for nice, but light, eps:
          class ListCollection(Collection):
              def __init__(self, collections, **kwargs):
                  Collection.__init__(self, **kwargs)
                  self.set_collections(collections)
              def set_collections(self, collections):
                  self._collections = collections
              def get_collections(self):
                  return self._collections
              @allow_rasterization
              def draw(self, renderer):
                  for _c in self._collections:
                      _c.draw(renderer)
  
          def insert_rasterized_contour_plot(c):
              collections = c.collections
              for _c in collections:
                  _c.remove()
              cc = ListCollection(collections, rasterized=True)
              ax = pl.gca()
              ax.add_artist(cc)
              return cc
  
  
          print('plotting contours')
          ax.autoscale(False)
          low_level = dcoeff_thresh - old_div(dcoeff_thresh,10)
          high_level= dcoeff_thresh + old_div(dcoeff_thresh,10)
          if datatype=='dcoeff':
              low_level = dcoeff_thresh - old_div(dcoeff_thresh,10)
              high_level= dcoeff_thresh + old_div(dcoeff_thresh,10)
              CMIX=ax.contourf(xxx,y,conv_i_vec_matrix, cmap=cmapMIX, alpha=1.0,levels=[low_level,1.e99])
          else:
              CMIX=ax.contourf(xxx,y,conv_i_vec_matrix, cmap=cmapMIX, alpha=1.0,levels=[0.5,1.5])
  
          if outlines == True:
              if datatype=='dcoeff':
                  CMIXoutlines=ax.contour(xxx,y,conv_i_vec_matrix, cmap=cmapMIX,levels=[low_level,high_level])
              else:
                  CMIX_outlines=ax.contour(xxx,y,conv_i_vec_matrix, cmap=cmapMIX,levels=[0.5,1.5])
          ax.plot(xxx,mtot,color='k')
          ax.axis([xlims[0],xlims[-1],ylims[0],ylims[1]])
  
          def contour_to_hatched_patches(cntrset, hatch_colors, hatch_patterns,
                                 remove_contour=True):
              from itertools import cycle
  
              patches_list = []
              for pathcollection in cntrset.collections:
                  patches_list.append([PathPatch(p1) for p1 in  pathcollection.get_paths()])
                  if remove_contour:
                      pathcollection.remove()
  
              for patches, hc, hp  in zip(patches_list,
                                          cycle(hatch_colors), cycle(hatch_patterns)):
                  for p in patches:
                      p.set_fc("none")
                      p.set_ec("k")
                      p.set_hatch(hp)
                      ax.add_patch(p)
  
          hatch_colors = "k"
          hatch_patterns = "/\|-+xoO.*"
          if hatching==True:
              contour_to_hatched_patches(CMIX, hatch_colors, hatch_patterns,
                      remove_contour=True)
  
          # uncomment or rasterizing the contour - for making nice eps figures that aren't huge files:
          #insert_rasterized_contour_plot(CMIX)
          #if outlines==True:
          #    insert_rasterized_contour_plot(CMIXoutlines)
  
          if outfile!='':
              if outfile[-4:] == '.png':
                  fig.savefig(outfile,dpi=300)
              elif outfile[-4:] == '.eps':
                  fig.savefig(outfile,format='eps')
              elif outfile[-3:] == '.ps':
                  fig.savefig(outfile,format='ps')
              elif outfile[-4:] == '.pdf':
                  fig.savefig(outfile,format='pdf')
              else:
                  print('file format not specified as .eps or .png, saving as eps:')
                  print(outfile+'.eps')
                  fig.savefig(outfile+'.eps',format='eps')
  
          pl.show()