docs: def kip_cont(self, ifig=110, modstart=0, modstop=-1,t0_model=0,
                   outfile='out.png', xlims=None, ylims=None,
                   xres=1000, yres=1000, ixaxis='model_number',
                   mix_zones=20, burn_zones=20, plot_radius=False,
                   engenPlus=True, engenMinus=False,
                   landscape_plot=False, rad_lines=False, profiles=[],
                   showfig=True, outlines=True, boundaries=True,
                   c12_boundary=False, rasterise=False, yscale='1.',
                   engenlevels=None,CBM=False,fsize=14):
          """
          This function creates a Kippenhahn plot with energy flux using
          contours.
  
          This plot uses mixing_regions and burning_regions written to
          your history.data or star.log. Set both variables in the
          log_columns.list file to 20 as a start.
  
          The output log file should then contain columns called
          "mix_type_n", "mix_qtop_n", "burn_type_n" and "burn_qtop_n".
          The number of columns (i.e. the biggest value of n) is what
          goes in the arguments as mix_zones and burn_zones.  DO NOT
          WORRY! if you do not have these columns, just leave the default
          values alone and the script should recognise that you do not
          have these columns and make the most detailed plot that is
          available to you.
  
          Defaults are set to get some plot, that may not look great if
          you zoom in interactively. Play with xres and yres as well as
          setting the xlims to ylims to the region you are interested in.
  
          Parameters
          ----------
          ifig : integer, optional
              Figure frame number.  The default is 110.
          modstart : integer, optional
              Model from which you want to plot (be careful if your
              history.data or star.log output is sparse...).  If it is
              0 then it starts from the beginning, works even if
              log_cnt > 1.
              The default is 0.
          modstop : integer, optional
              Model to which you wish to plot, -1 corresponds to end
              [if log_cnt>1, devide modstart and modstop by log_cnt,
              this needs to be improved! SJ: this should be ficed now].
              The defalut is -1.
          t0_model : integer, optional
              Model number from which to reset the time to 0. Typically,
              if modstart!=0, t0_model=modstart is a good choice, but we
              leave the choice to the user in case the time is wished to
              start from 0 at a different key point of the evolution.
              The default value is 0.
          outfile : sting, optional
              'filename + extension' where you want to save the figure.
              The defalut is "out.png".
          xlims, ylims : list, optional
              Plot limits, however these are somewhat obsolete now that
              we have modstart and modstop.  Leaving them as 0. is
              probably no slower, and you can always zoom in afterwards
              in mpl.  ylims is important for well resolved thermal
              pulse etc plots; it's best to get the upper and lower limits
              of he-intershell using s.kippenhahn_CO(1,'model') first.
              The default is [0., 0.].
          xres, yres : integer, optional
              plot resolution. Needless to say that increasing these
              values will yield a nicer plot with some slow-down in
              plotting time.  You will most commonly change xres.  For a
              prelim plot, try xres~200, then bump it up to anywhere from
              1000-10000 for real nicely resolved, publication quality
              plots.  The default is 1000.
          ixaxis : string, optional
              Either 'log_time_left', 'age', or 'model_number'.  The
              default is "model_number".
          mix_zones, burn_zones : integer, optional
              As described above, if you have more detailed output about
              your convection and energy generation boundaries in columns
              mix_type_n, mix_qtop_n, burn_type_n and burn_qtop_n, you
              need to specify the total number of columns for mixing zones
              and burning zones that you have.  Can't work this out from
              your history.data or star.log file?  Check the
              history_columns.list that you used, it'll be the number
              after "mixing regions" and "burning regions".  Can't see
              these columns?  leave it and 2 conv zones and 2 burn zones
              will be drawn using other data that you certainly should
              have in your history.data or star.log file.  The defalut for
              both is 20.
          plot_radius : boolean, optional
              Whether on a second y-axis you want to plot the radius of
              the surface and the he-free core.  The default is False.
          engenPlus : boolean
              Plot energy generation contours for eps_nuc>0.  The default
              is True.
          endgenMinus : boolean, optional
              Plot energy generation contours for eos_nuc<0.  The default
              is True.
          landscape_plot : boolean, optionla
              The default is False.
          rad_lines : boolean, optional
              The deafault is False.
          profiles : list, optional
              The default is [].
          showfig : boolean, optional
              The default is True.
          outlines : boolean, optional
              Whether or not to plot outlines of conv zones in darker
              colour.
          boundaries : boolean, optional
              Whether or not to plot H-, He- and C-free boundaries.
          c12_boundary : boolean, optional
              The default is False.
          rasterise : boolean, optional
              Whether or not to rasterise the contour regions to make
              smaller vector graphics figures.  The default is False.
          yscale : string, optional
              Re-scale the y-axis by this amount
          engenlevels : list
              Give cusstom levels to the engenPlus contour. If None,
              the levels are chosen automatically.
              The default is None.
          CBM : boolean, optional
              plot contours for where CBM is active?
  
          fsize : integer
              font size for labels
  
          Notes
          -----
          The parameter xlims is depricated.
  
          """
          if ylims is None:
              ylims=[0.,0.]
          if xlims is None:
              xlims=[0.,0.]
          
          # Find correct modstart and modstop:
          mod=np.array([int(i) for i in self.get('model_number')])
          mod1=np.abs(mod-modstart).argmin()
          mod2=np.abs(mod-modstop).argmin()
          if modstart != 0 : modstart=mod1
          if modstop != -1 : modstop=mod2
  
          xxyy=[self.get('star_age')[modstart:modstop],self.get('star_age')[modstart:modstop]]
          mup = max(float(self.get('star_mass')[0])*1.02,1.0)
          nmodels=len(self.get('model_number')[modstart:modstop])
  
          if ylims == [0.,0.]:
              mup   = max(float(self.get('star_mass')[0])*1.02,1.0)
              mDOWN = 0.
          else:
              mup = ylims[1]
              mDOWN = ylims[0]
  
          # y-axis resolution
          ny=yres
          #dy=mup/float(ny)
          dy = old_div((mup-mDOWN),float(ny))
  
          # x-axis resolution
          maxpoints=xres
          dx=int(max(1,old_div(nmodels,maxpoints)))
  
          #y = np.arange(0., mup, dy)
          y = np.arange(mDOWN, mup, dy)
          x = np.arange(0, nmodels, dx)
          Msol=1.98892E+33
  
          engenstyle = 'full'
  
          B1=np.zeros([len(y),len(x)],float)
          B2=np.zeros([len(y),len(x)],float)
          try:
              self.get('burn_qtop_1')
          except:
              engenstyle = 'twozone'
          if engenstyle == 'full' and (engenPlus == True or engenMinus == True):
              ulimit_array = np.array([self.get('burn_qtop_'+str(j))[modstart:modstop:dx]*\
                             self.get('star_mass')[modstart:modstop:dx] for j in range(1,burn_zones+1)])
              #ulimit_array = np.around(ulimit_array,decimals=len(str(dy))-2)
              llimit_array = np.delete(ulimit_array,-1,0)
              llimit_array = np.insert(ulimit_array,0,0.,0)
              #llimit_array = np.around(llimit_array,decimals=len(str(dy))-2)
              btype_array = np.array([self.get('burn_type_'+str(j))[modstart:modstop:dx] for j in range(1,burn_zones+1)])
              old_percent = 0
              for i in range(len(x)):
                  # writing status
                  percent = int(i*100/(len(x) - 1))
                  if percent >= old_percent + 5:
                      sys.stdout.flush()
                      sys.stdout.write("\r creating color map burn " + "...%d%%" % percent)
                      old_percent = percent
  
                  for j in range(burn_zones):
                      if btype_array[j,i] > 0. and abs(btype_array[j,i]) < 99.:
                          B1[(np.abs(y-llimit_array[j][i])).argmin():(np.abs(y-ulimit_array[j][i])).argmin()+1,i] = 10.0**(btype_array[j,i])
                      elif btype_array[j,i] < 0. and abs(btype_array[j,i]) < 99.:
                          B2[(np.abs(y-llimit_array[j][i])).argmin():(np.abs(y-ulimit_array[j][i])).argmin()+1,i] = 10.0**(abs(btype_array[j,i]))
              print(' \n')
  
          if engenstyle == 'twozone' and (engenPlus == True or engenMinus == True):
              V=np.zeros([len(y),len(x)],float)
              old_percent = 0
              for i in range(len(x)):
                  # writing status
                  percent = int(i*100/(len(x) - 1))
                  if percent >= old_percent + 5:
                      sys.stdout.flush()
                      sys.stdout.write("\r creating color map1 " + "...%d%%" % percent)
                      old_percent = percent
                  llimitl1=old_div(self.get('epsnuc_M_1')[modstart:modstop][i*dx],Msol)
                  ulimitl1=old_div(self.get('epsnuc_M_4')[modstart:modstop][i*dx],Msol)
                  llimitl2=old_div(self.get('epsnuc_M_5')[modstart:modstop][i*dx],Msol)
                  ulimitl2=old_div(self.get('epsnuc_M_8')[modstart:modstop][i*dx],Msol)
                  llimith1=old_div(self.get('epsnuc_M_2')[modstart:modstop][i*dx],Msol)
                  ulimith1=old_div(self.get('epsnuc_M_3')[modstart:modstop][i*dx],Msol)
                  llimith2=old_div(self.get('epsnuc_M_6')[modstart:modstop][i*dx],Msol)
                  ulimith2=old_div(self.get('epsnuc_M_7')[modstart:modstop][i*dx],Msol)
                  # lower thresh first, then upper thresh:
                  if llimitl1!=ulimitl1:
                      for k in range(ny):
                          if llimitl1<=y[k] and ulimitl1>y[k]:
                              V[k,i]=10.
                  if llimitl2!=ulimitl2:
                      for k in range(ny):
                          if llimitl2<=y[k] and ulimitl2>y[k]:
                              V[k,i]=10.
                  if llimith1!=ulimith1:
                      for k in range(ny):
                          if llimith1<=y[k] and ulimith1>y[k]:
                              V[k,i]=30.
                  if llimith2!=ulimith2:
                      for k in range(ny):
                          if llimith2<=y[k] and ulimith2>y[k]:
                              V[k,i]=30.
              print(' \n')
  
          mixstyle = 'full'
          try:
              self.get('mix_qtop_1')
          except:
              mixstyle = 'twozone'
          if mixstyle == 'full':
              old_percent = 0
              Z=np.zeros([len(y),len(x)],float)
              if CBM:
                  Zcbm=np.zeros([len(y),len(x)],float)
              ulimit_array = np.array([self.get('mix_qtop_'+str(j))[modstart:modstop:dx]*self.get('star_mass')[modstart:modstop:dx] for j in range(1,mix_zones+1)])
              llimit_array = np.delete(ulimit_array,-1,0)
              llimit_array = np.insert(ulimit_array,0,0.,0)
              mtype_array = np.array([self.get('mix_type_'+str(j))[modstart:modstop:dx] for j in range(1,mix_zones+1)])
              for i in range(len(x)):
                  # writing status
                  percent = int(i*100/(len(x) - 1))
                  if percent >= old_percent + 5:
                      sys.stdout.flush()
                      sys.stdout.write("\r creating color map mix " + "...%d%%" % percent)
                      old_percent = percent
                  for j in range(mix_zones):
                      if mtype_array[j,i] == 1.:
                          Z[(np.abs(y-llimit_array[j][i])).argmin():(np.abs(y-ulimit_array[j][i])).argmin()+1,i] = 1.
                      if CBM:
                          if mtype_array[j,i] == 2.:
                              Zcbm[(np.abs(y-llimit_array[j][i])).argmin():(np.abs(y-ulimit_array[j][i])).argmin()+1,i] = 1.
              print(' \n')
  
          if mixstyle == 'twozone':
              Z=np.zeros([len(y),len(x)],float)
              old_percent = 0
              for i in range(len(x)):
              # writing reading status
                  # writing status
                  percent = int(i*100/(len(x) - 1))
                  if percent >= old_percent + 5:
                      sys.stdout.flush()
                      sys.stdout.write("\r creating color map mix " + "...%d%%" % percent)
                      old_percent = percent
  
                  ulimit=self.get('conv_mx1_top')[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                  llimit=self.get('conv_mx1_bot')[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                  if llimit!=ulimit:
                      for k in range(ny):
                          if llimit<=y[k] and ulimit>y[k]:
                              Z[k,i]=1.
                  ulimit=self.get('conv_mx2_top')[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                  llimit=self.get('conv_mx2_bot')[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                  if llimit!=ulimit:
                      for k in range(ny):
                          if llimit<=y[k] and ulimit>y[k]:
                              Z[k,i]=1.
              print(' \n')
  
          if rad_lines == True:
              masses = np.arange(0.1,1.5,0.1)
              rads=[[],[],[],[],[],[],[],[],[],[],[],[],[],[]]
              modno=[]
              for i in range(len(profiles)):
                  p=mesa_profile('./LOGS',profiles[i])
                  modno.append(p.header_attr['model_number'])
                  for j in range(len(masses)):
                      idx=np.abs(p.get('mass')-masses[j]).argmin()
                      rads[j].append(p.get('radius')[idx])
  
          print('engenstyle was ', engenstyle)
          print('mixstyle was ', mixstyle)
          print('\n finished preparing color map')
  
          ########################################################################
          #----------------------------------plot--------------------------------#
          fig = pyl.figure(ifig)
          #fsize=20
          if landscape_plot == True:
              fig.set_size_inches(9,4)
              pl.gcf().subplots_adjust(bottom=0.2)
              pl.gcf().subplots_adjust(right=0.85)
  
          params = {'axes.labelsize':  fsize,
            'axes.labelsize': fsize,
            'font.size': fsize,
            'legend.fontsize': fsize,
            'xtick.labelsize': fsize,
            'ytick.labelsize': fsize,
                    'text.usetex': False}
          pyl.rcParams.update(params)
  
          #ax=pl.axes([0.1,0.1,0.9,0.8])
  
          #fig=pl.figure()
          ax=pl.axes()
  
          if ixaxis == 'log_time_left':
          # log of time left until core collapse
              gage= self.get('star_age')
              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 = lage[modstart:modstop]
              print('plot versus time left')
              ax.set_xlabel('$ \\log_{10}(t-t_\mathrm{end})\ /\  \mathrm{[yr]}$',fontsize=fsize)
              if xlims[1] == 0.:
                  xlims = [xxx[0],xxx[-1]]
          elif ixaxis =='model_number':
              xxx= self.get('model_number')[modstart:modstop]
              print('plot versus model number')
              ax.set_xlabel('Model number',fontsize=fsize)
              if xlims[1] == 0.:
                  xlims = [self.get('model_number')[modstart],self.get('model_number')[modstop]]
          elif ixaxis =='age':
              if t0_model != 0:
                  t0_mod=np.abs(mod-t0_model).argmin()
                  xxx= self.get('star_age')[modstart:modstop] - self.get('star_age')[t0_mod]
                  print('plot versus age')
                  ax.set_xlabel('t - %.5e / [yr]' %self.get('star_age')[modstart],fontsize=fsize)
              else:
                  xxx= old_div(self.get('star_age')[modstart:modstop],1.e6)
                  ax.set_xlabel('t [Myr]',fontsize=fsize)
              if xlims[1] == 0.:
                  xlims = [xxx[0],xxx[-1]]
  
          ax.set_ylabel('$\mathrm{enclosed\ mass\ /\ [M_\odot]}$',fontsize=fsize)
  
          # 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
  
          cmapMIX = matplotlib.colors.ListedColormap(['w','#8B8386']) # rose grey
          if CBM:
              cmapCBM = matplotlib.colors.ListedColormap(['w','g']) # green
          cmapB1  = pyl.cm.get_cmap('Blues')
          cmapB2  = pl.cm.get_cmap('Reds')
  
          ylims1=[0.,0.]
          ylims1[0]=ylims[0]
          ylims1[1]=ylims[1]
          if ylims == [0.,0.]:
              ylims[0] = 0.
              ylims[1] = mup
              #print("Setting ylims[1] to mup="+str(mup))
          if ylims[0] != 0.:
              ylab='$(\mathrm{Mass }$ - '+str(ylims[0])
              if yscale!='1.':
                  ylab+=') / '+yscale+' $M_\odot$'
              else:
                  ylab+=') / $M_\odot$'
              ax.set_ylabel(ylab)
              y = y - ylims[0]
              y = y*float(yscale) # SJONES tweak
              ylims[0] = y[0]
              ylims[1] = y[-1]
  
          print('plotting contours')
          CMIX    = ax.contourf(xxx[::dx],y,Z, cmap=cmapMIX,alpha=0.6,levels=[0.5,1.5])
          #CMIX    = ax.pcolor(xxx[::dx],y,Z, cmap=cmapMIX,alpha=0.6,vmin=0.5,vmax=1.5)
          if rasterise==True:
              insert_rasterized_contour_plot(CMIX)
          if outlines == True:
              CMIX_outlines    = ax.contour(xxx[::dx],y,Z, cmap=cmapMIX)
              if rasterise==True:
                  insert_rasterized_contour_plot(CMIX_outlines)
  
          if CBM:
              CCBM    = ax.contourf(xxx[::dx],y,Zcbm, cmap=cmapCBM,alpha=0.6,levels=[0.5,1.5])
              if rasterise==True:
                  insert_rasterized_contour_plot(CCBM)
              if outlines == True:
                  CCBM_outlines    = ax.contour(xxx[::dx],y,Zcbm, cmap=cmapCBM)
                  if rasterise==True:
                      insert_rasterized_contour_plot(CCBM_outlines)
  
          if engenstyle == 'full' and engenPlus == True:
              if engenlevels!= None:
                  CBURN1  = ax.contourf(xxx[::dx],y,B1, cmap=cmapB1, alpha=0.5,\
                          locator=matplotlib.ticker.LogLocator(),levels=engenlevels)
                  if outlines:
                      CB1_outlines  = ax.contour(xxx[::dx],y,B1, cmap=cmapB1, alpha=0.7, \
                          locator=matplotlib.ticker.LogLocator(),levels=engenlevels)
              else:
                  CBURN1  = ax.contourf(xxx[::dx],y,B1, cmap=cmapB1, alpha=0.5, \
                                        locator=matplotlib.ticker.LogLocator())
                  if outlines:
                      CB1_outlines  = ax.contour(xxx[::dx],y,B1, cmap=cmapB1, alpha=0.7, \
                                                 locator=matplotlib.ticker.LogLocator())
              CBARBURN1 = pyl.colorbar(CBURN1)
              CBARBURN1.set_label('$|\epsilon_\mathrm{nuc}-\epsilon_{\\nu}| \; (\mathrm{erg\,g}^{-1}\mathrm{\,s}^{-1})$',fontsize=fsize)
              if rasterise==True:
                  insert_rasterized_contour_plot(CBURN1)
                  if outlines:
                      insert_rasterized_contour_plot(CB1_outlines)
  
          if engenstyle == 'full' and engenMinus == True:
              CBURN2  = ax.contourf(xxx[::dx],y,B2, cmap=cmapB2, alpha=0.5, locator=matplotlib.ticker.LogLocator())
              if outlines:
                  CBURN2_outlines  = ax.contour(xxx[::dx],y,B2, cmap=cmapB2, alpha=0.7, locator=matplotlib.ticker.LogLocator())
              CBARBURN2 = pl.colorbar(CBURN2)
              if engenPlus == False:
                  CBARBURN2.set_label('$|\epsilon_\mathrm{nuc}-\epsilon_{\\nu}| \; (\mathrm{erg\,g}^{-1}\mathrm{\,s}^{-1})$',fontsize=fsize)
              if rasterise==True:
                  insert_rasterized_contour_plot(CBURN2)
                  if outlines:
                      insert_rasterized_contour_plot(CB2_outlines)
  
          if engenstyle == 'twozone' and (engenPlus == True or engenMinus == True):
              ax.contourf(xxx[::dx],y,V, cmap=cmapB1, alpha=0.5)
  
          print('plotting patches')
          mtot=self.get('star_mass')[modstart:modstop][::dx]
          mtot1=(mtot-ylims1[0])*float(yscale)
          ax.plot(xxx[::dx],mtot1,'k-')
  
          if boundaries == True:
              print('plotting abund boundaries')
              try:
                  bound=self.get('h1_boundary_mass')[modstart:modstop]
                  bound1=(bound-ylims1[0])*float(yscale)
                  ax.plot(xxx,bound1,label='H boundary',linestyle='-')
  
                  bound=self.get('he4_boundary_mass')[modstart:modstop]
                  bound1=(bound-ylims1[0])*float(yscale)
                  ax.plot(xxx,bound1,label='He boundary',linestyle='--')
  
                  bound=self.get('c12_boundary_mass')[modstart:modstop]
                  bound1=(bound-ylims1[0])*float(yscale)
                  ax.plot(xxx,bound1,label='C boundary',linestyle='-.')
  
              except:
                  try:
                      bound=self.get('he_core_mass')[modstart:modstop]
                      bound1=(bound-ylims1[0])*float(yscale)
                      ax.plot(xxx,bound1,label='H boundary',linestyle='-')
  
                      bound=self.get('c_core_mass')[modstart:modstop]-ylims[0]
                      bound1=(bound-ylims1[0])*float(yscale)
                      ax.plot(xxx,bound1,label='He boundary',linestyle='--')
  
                      bound=self.get('o_core_mass')[modstart:modstop]-ylims[0]
                      bound1=(bound-ylims1[0])*float(yscale)
                      ax.plot(xxx,bound1,label='C boundary',linestyle='-.')
  
                      bound=self.get('si_core_mass')[modstart:modstop]-ylims[0]
                      bound1=(bound-ylims1[0])*float(yscale)
                      ax.plot(xxx,bound1,label='C boundary',linestyle='-.')
  
                      bound=self.get('fe_core_mass')[modstart:modstop]-ylims[0]
                      bound1=(bound-ylims1[0])*float(yscale)
                      ax.plot(xxx,bound1,label='C boundary',linestyle='-.')
  
                  except:
  #                    print 'problem to plot boundaries for this plot'
                      pass
  
          ax.axis([xlims[0],xlims[1],ylims[0],ylims[1]])
  
          if plot_radius == True:
              ax2=pyl.twinx()
              ax2.plot(xxx,np.log10(self.get('he4_boundary_radius')[modstart:modstop]),label='He boundary radius',color='k',linewidth=1.,linestyle='-.')
              ax2.plot(xxx,self.get('log_R')[modstart:modstop],label='radius',color='k',linewidth=1.,linestyle='-.')
              ax2.set_ylabel('log(radius)')
          if rad_lines == True:
              ax2=pyl.twinx()
              for i in range(len(masses)):
                  ax2.plot(modno,np.log10(rads[i]),color='k')
  
          if outfile[-3:]=='png':
              fig.savefig(outfile,dpi=300)
          elif outfile[-3:]=='eps':
              fig.savefig(outfile,format='eps')
          elif outfile[-3:]=='pdf':
              fig.savefig(outfile,format='pdf')
          if showfig == True:
              pyl.show()